PROGRAM Spherical_Lurch ! Computes the coseismic displacement ("lurch") at the surface ! due to a set of rectangular seismic-dislocation patches at depth. ! This program is "spherical" because all input and output files ! use (longitude, latitude) coordinates to express locations. ! An input file with dislocation patches and their slips and rakes ! is required. The user will specify the desired grid for output. ! This program will produce a set of 6 GRD files with 3 components ! and 3 other descriptors (magnitude, trend, & plunge) ! of the surface displacement. Note that topography is not ! modeled; output values are on a flat surface, and input depths ! of dislocation patches should be relative to this surface. ! By Peter Bird, UCLA, 2017.03.17 USE DSphere ! provided by P. Bird as file DSphere.f90. ! DSphere is a library of geometric operations on a dimensionless unit sphere. ! The "D" in "DSphere" indicates DOUBLE PRECISION (or REAL*8) arguments and accuracy. USE DDislocation ! provided by P. Bird as file DDislocation.f90. ! DDislocation contains the rectangular-dislocation-patch solutions of Mansinha & Smylie [1967 JGR, 1971 BSSA]. ! Be sure to get the latest version (2017.03.17 or later), as earlier versions do not provide ! the vertical component of displacement, and therefore will not link properly with this program. USE DTriangular ! provided by P. Bird as file DTriangular.f90. ! DTriangular contains the triangular-dislocation-patch solutions of Nikkhoo & Walter [2015, GJI], ! translated from MatLab to Fortran 90 by P. Bird. IMPLICIT NONE CHARACTER*72 :: input_filename, short_line INTEGER :: columns, dislocation_type, i, input_style, ios, j, patches, rows, uAzimuth INTEGER, DIMENSION(:, :), ALLOCATABLE :: uPlunge, uTrend REAL*8, PARAMETER :: R_Earth = 6371.0D3 ! meters REAL*8, PARAMETER :: Poisson = 0.25D0 REAL*8 :: argume, & & bPhi, bTheta, & & cellsize_degrees, & & depth_km, depth_m, dip, dipF, downdip_depth_m, Ds, dUPhi, dURadial, dUTheta, dx, dy, & & Elon, Elon_at_SE_corner, Elon_at_SW_corner, Elon_surface_center, & & halfdip_surface_arc_radians, halfstrike_surface_arc_radians, & & lF, & & Nlat, Nlat_at_SW_corner, Nlat_at_NW_corner, Nlat_surface_center, & & phi_surface_center, & & rake, & & slip_m, strike, Ss, & & test_Elon, test_Nlat, test_depth_m, theta_surface_center, Ts, & & uEos, uHorizontal, updip_depth_m, uRadial, uSouth, & & wedge, wholedip_surface_arc_radians, & & x_km, & & y_km REAL*8, DIMENSION(3) :: slipVector ! input vector for DChange REAL*8, DIMENSION(3) :: center_uvec, & & downdip_center_uvec, downdip_start_uvec, downdip_end_uvec, & & omega_uvec, test_uvec, & & surface_center_uvec, & & updip_center_uvec, updip_start_uvec, updip_end_uvec ! where "start" and "end" refer to moving along the strike of the fault, in the direction given by azimuth "strike". REAL*8, DIMENSION(:, :), ALLOCATABLE :: uEast, uNorth, uUp, uMagnitude !Introduce this program: WRITE (*, "(' PROGRAM Spherical_Lurch:')") WRITE (*, *) WRITE (*, "(' Computes the coseismic displacement (""lurch"") at the surface')") WRITE (*, "(' due to a set of rectangular seismic-dislocation patches at depth.')") WRITE (*, "(' This program is ""spherical"" because all input and output files')") WRITE (*, "(' use (longitude, latitude) coordinates to express locations.')") WRITE (*, "(' An input file with dislocation patches and their slips and rakes')") WRITE (*, "(' is required. The user will specify the desired grid for output.')") WRITE (*, "(' This program will produce a set of 6 GRD files with 3 components')") WRITE (*, "(' and 3 other descriptors (magnitude, trend, & plunge)')") WRITE (*, "(' of the surface displacement. Note that topography is not')") WRITE (*, "(' modeled; output values are on a flat surface, and input depths')") WRITE (*, "(' of dislocation patches should be relative to this surface.')") WRITE (*, "(' By Peter Bird, UCLA, 2017.03.17')") CALL Pause() !Select input format & filename: 10 WRITE (*, *) WRITE (*, "(' Choice of input formats for the file with seismic dislocation patches:')") WRITE (*, "(' 1: Tinti et al. [2016, GRL] Amatrice earthquake of 2016.08.24')") WRITE (*, "(' ----------------------------------------------------------------------')") WRITE (*, "(' Enter an integer from the list above: '\)") READ (*, *) input_style IF ((input_style < 1).OR.(input_style > 1)) THEN WRITE (*, "(' ERROR: Please try again...')") CALL Pause() GO TO 10 END IF 20 WRITE (*, *) WRITE (*, "(' Enter filename for (existing) input file in this format:')") READ (*, "(A)") input_filename OPEN (UNIT = 1, FILE = TRIM(input_filename), STATUS = "OLD", PAD = "YES", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: Named file could not be read (from current folder).')") WRITE (*, "(' Please move the input file, and try again...')") CALL Pause() GO TO 20 END IF !Display headers in the dislocation input file: IF (input_style == 1) THEN WRITE (*, *) WRITE (*, "(' Here are the first 72 bytes of each header line from this input file:')") WRITE (*, "(' ------------------------------------------------------------------------')") DO i = 1, 13 READ (1, "(A)") short_line WRITE (*, "(' ',A)") TRIM(short_line) END DO WRITE (*, "(' ------------------------------------------------------------------------')") END IF CALL Pause() !Prompt user to define the grid of output points: WRITE (*, *) WRITE (*, "(' Next, you must define the grid for the output values:')") WRITE (*, "(' (HINT: You might want to copy these values from a matching ALOS.asc file.)')") 30 WRITE (*, "(' Enter number of grid columns: '\)") READ (*, *) columns IF (columns < 1) THEN WRITE (*, "(' ERROR: Please enter 1 or more.')") CALL Pause() GO TO 30 END IF 40 WRITE (*, "(' Enter number of grid rows: '\)") READ (*, *) rows IF (rows < 1) THEN WRITE (*, "(' ERROR: Please enter 1 or more.')") CALL Pause() GO TO 40 END IF 50 WRITE (*, "(' East longitude (use - for W) at SW corner of grid, degrees: '\)") READ (*, *) Elon_at_SW_corner 60 WRITE (*, "(' North latitude (use - for S) at SW corner of grid, degrees: '\)") READ (*, *) Nlat_at_SW_corner IF (ABS(Nlat_at_SW_corner) > 90.0D0) THEN WRITE (*, "(' ERROR: Please enter -90. <= latitude < 90.')") CALL Pause() GO TO 60 END IF 70 WRITE (*, "(' Enter grid cell-size, in degrees: '\)") READ (*, *) cellsize_degrees IF (cellsize_degrees <= 0.0D0) THEN WRITE (*, "(' ERROR: Enter a positive cell-size.')") CALL Pause() GO TO 70 END IF Elon_at_SE_corner = Elon_at_SW_corner + (columns - 1) * cellsize_degrees Nlat_at_NW_corner = Nlat_at_SW_corner + (rows - 1) * cellsize_degrees WRITE (*, *) WRITE (*, "(' Grid longitudes will range from ',F10.4,' to ',F10.4,' degrees East.')") Elon_at_SW_corner, Elon_at_SE_corner WRITE (*, "(' Grid latitudes will range from ',F10.4,' to ',F10.4,' degrees North.')") Nlat_at_SW_corner, Nlat_at_NW_corner !Ask user to decide how the problem will be solved? 80 WRITE (*, *) WRITE (*, "(' Available solution methods:')") WRITE (*, "(' ----------------------------------------------------------------------------')") WRITE (*, "(' 1: rectangular dislocation patches [Manshinha & Smylie, 1967 JGR, 1971 BSSA]')") WRITE (*, "(' 2: triangular dislocation patches [Nikkhoo & Walter, 2015 GJI]')") WRITE (*, "(' ----------------------------------------------------------------------------')") WRITE (*, "(' Enter an integer from the list above: '\)") READ (*, *) dislocation_type IF ((dislocation_type < 1).OR.(dislocation_type > 2)) THEN WRITE (*, "(' ERROR: Please try again...')") CALL Pause() GO TO 80 END IF !Allocate arrays for output GRDs: ALLOCATE ( uEast (rows, columns) ) ALLOCATE ( uNorth (rows, columns) ) ALLOCATE ( uUp (rows, columns) ) ALLOCATE ( uMagnitude(rows, columns) ) ALLOCATE ( uTrend (rows, columns) ) ALLOCATE ( uPlunge (rows, columns) ) !Initialize arrays that will be computed as a sum over dislocation-patch effects: uEast = 0.0D0 uNorth = 0.0D0 uUp = 0.0D0 !Begin indefinite loop over seismic dislocation patches (from input file): patches = 0 patching: DO READ (1, *, IOSTAT = ios) Elon, Nlat, depth_km, x_km, y_km, dx, dy, slip_m, rake, dip, strike IF (ios /= 0) EXIT patching patches = patches + 1 ! N.B. One of the header lines in the Tinti et al. [2016, GRL] file notes "(values in the center)" {of each sub-fault patch}. ! For some ODD reason they have negative signs on all "depths", so what they really mean is "relative elevation" ! with respect to the ficticious observation plane at the mean elevation of all used seismic stations. depth_km = -depth_km ! correcting this oddity (see above). depth_m = depth_km * 1000.0D0 !Characterize depth range of this patch (using non-negative numbers): updip_depth_m = depth_m - 0.5D0 * ABS(dy) * SIN(dip * radians_per_degree) downdip_depth_m = depth_m + 0.5D0 * ABS(dy) * SIN(dip * radians_per_degree) !========================================================================================================= IF (dislocation_type == 1) THEN ! rectangular dislocations of Mansinha & Smylie [1967 JGR, 1971 BSSA]: !Characterize center and surface_center with uvecs: CALL DLonLat_2_Uvec(Elon, Nlat, center_uvec) IF (ABS(dip) < 5.0D0) THEN ! flat or very gently-dipping dislocation patch; wholedip_surface_arc_radians would be very large or infinite! WRITE (*, "(' ERROR: Your input data contains dislocation patch(es) dipping less than 5 degrees.')") WRITE (*, "(' The rectangular-dislocation-method of Mansinha & Smylie [1971] cannot solve cases')") WRITE (*, "(' of zero dip, and is not very accurate when the dip is very small.')") WRITE (*, "(' Therefore, please re-start this program Spherical_Lurch,')") WRITE (*, "(' and select dislocation type 2, for triangular dislocation patches.')") CALL Pause() STOP ELSE IF (ABS(90.0D0 - dip) < 1.0D0) THEN ! vertical fault; avoid infinite TAN(dip)! surface_center_uvec = center_uvec ELSE ! most common case (5 < dip < 89): wholedip_surface_arc_radians = ABS(depth_m) / (TAN(dip * radians_per_degree) * R_Earth) ! defined as non-negative CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) - Pi_over_2), base_uvec = center_uvec, far_radians = wholedip_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = surface_center_uvec) ! outputs END IF !Find (lon, lat) and (theta, phi) coordinates of surface_center point (midpoint of virtual fault trace, when dislocation-patch plane is extended): CALL DUvec_2_LonLat(surface_center_uvec, Elon_surface_center, Nlat_surface_center) theta_surface_center = (90.0D0 - Nlat_surface_center) * radians_per_degree phi_surface_center = Elon_surface_center * radians_per_degree !Translate other properties of the dislocation patch into terms understood by DChange: argume = (180.0D0 - strike) * radians_per_degree dipF = dip * radians_per_degree lF = 0.5D0 * dx wedge = 15.0D0 * radians_per_degree ! (as in NeoKinema) slipVector(1) = slip_m * COS(rake * radians_per_degree) ! strike-slip, positive when sinistral slipVector(2) = slip_m * (-SIN(rake * radians_per_degree)) * COS(dipF) ! heave, positive for divergence slipVector(3) = slip_m * (-SIN(rake * radians_per_degree)) * SIN(dipF) ! throw, positive when the right side of the fault (when looking along direction ARGUME) is down. !Survey all test-points (virtual benchmarks) in the GRD array: DO i = 1, rows test_Nlat = Nlat_at_NW_corner - (i - 1) * cellsize_degrees bTheta = (90.0D0 - test_Nlat) * radians_per_degree DO j = 1, columns test_Elon = Elon_at_SW_corner + (j - 1) * cellsize_degrees bPhi = test_Elon * radians_per_degree CALL DChange (argume = argume, & ! inputs btheta = bTheta, bphi = bPhi, & ! inputs dipf = dipF, lf = lF, & ! inputs ftheta = theta_surface_center, fphi = phi_surface_center, & ! inputs radius = R_Earth, & ! inputs slip = slipVector, & ! inputs wedge = wedge, & ! inputs ztop = updip_depth_m, zbot = downdip_depth_m, & ! inputs duthet = dUTheta, duphi = dUPhi, dur = dURadial) ! output uEos = dUPhi uSouth = dUTheta uRadial = dURadial uEast(i, j) = uEast(i, j) + uEos uNorth(i, j) = uNorth(i, j) - uSouth uUp(i, j) = uUp(i, j) + uRadial END DO END DO !========================================================================================================= ELSE IF (dislocation_type == 2) THEN ! triangular dislocation patches of Nikkhoo & Walter [2015 GJI]: !Characterize center and corners of this dislocation path with uvecs (which identify surface points directly above those deep points): CALL DLonLat_2_Uvec(Elon, Nlat, center_uvec) halfdip_surface_arc_radians = 0.5D0 * ABS(dy) * COS(dip * radians_per_degree) / R_Earth ! defined as non-negative CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) + Pi_over_2), base_uvec = center_uvec, far_radians = halfdip_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = downdip_center_uvec) ! outputs CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) - Pi_over_2), base_uvec = center_uvec, far_radians = halfdip_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = updip_center_uvec) ! outputs halfstrike_surface_arc_radians = 0.5D0 * ABS(dx) / R_Earth ! defined as non-negative CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) + Pi), base_uvec = updip_center_uvec, far_radians = halfstrike_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = updip_start_uvec) ! outputs CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) ), base_uvec = updip_center_uvec, far_radians = halfstrike_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = updip_end_uvec) ! outputs CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) + Pi), base_uvec = downdip_center_uvec, far_radians = halfstrike_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = downdip_start_uvec) ! outputs CALL DTurn_To (azimuth_radians = ((strike * radians_per_degree) ), base_uvec = downdip_center_uvec, far_radians = halfstrike_surface_arc_radians, & ! inputs & omega_uvec = omega_uvec, result_uvec = downdip_end_uvec) ! outputs !Break slip into components as in TDdispHS: Ds = dip-slip (thrust-faulting positive); ! Ss = strike-slip (sinistral or left-lateral positive); ! Ts = tensile "slip" (opening; dike or sill intrusion positive): Ds = slip_m * SIN(rake * radians_per_degree) Ss = slip_m * COS(rake * radians_per_degree) Ts = 0.0D0 ! which will only increase if there is dike intrusion during the earthquake. !Scan all GRD observation points, and add displacements due to this one small seismic dislocation patch: DO i = 1, rows test_Nlat = Nlat_at_NW_corner - (i - 1) * cellsize_degrees DO j = 1, columns test_Elon = Elon_at_SW_corner + (j - 1) * cellsize_degrees CALL DLonLat_2_Uvec(test_Elon, test_Nlat, test_uvec) test_depth_m = 0.0D0 CALL Gore(R_m = R_Earth, P1_uvec = updip_start_uvec, P2_uvec = downdip_start_uvec, P3_uvec = downdip_end_uvec, & ! inputs & P1_depth_m = updip_depth_m, P2_depth_m = downdip_depth_m, P3_depth_m = downdip_depth_m, & ! inputs & Ds = Ds, Ss = Ss, Ts = Ts, Poisson = Poisson, & ! inputs & test_uvec = test_uvec, test_depth_m = test_depth_m, & ! inputs & uSouth = uSouth, uEast = uEos, uRadial = uRadial) ! outputs uEast(i, j) = uEast(i, j) + uEos uNorth(i, j) = uNorth(i, j) - uSouth uUp(i, j) = uUp(i, j) + uRadial CALL Gore(R_m = R_Earth, P1_uvec = downdip_end_uvec, P2_uvec = updip_end_uvec, P3_uvec = updip_start_uvec, & ! inputs & P1_depth_m = downdip_depth_m, P2_depth_m = updip_depth_m, P3_depth_m = updip_depth_m, & ! inputs & Ds = Ds, Ss = Ss, Ts = Ts, Poisson = Poisson, & ! inputs & test_uvec = test_uvec, test_depth_m = test_depth_m, & ! inputs & uSouth = uSouth, uEast = uEos, uRadial = uRadial) ! outputs uEast(i, j) = uEast(i, j) + uEos uNorth(i, j) = uNorth(i, j) - uSouth uUp(i, j) = uUp(i, j) + uRadial END DO END DO !========================================================================================================= END IF ! dislocation_type == 1, or ==2 ? WRITE (*, "(' Finished computing surface displacements due to patch #',I10)") patches END DO patching WRITE (*, *) WRITE (*, "(' ',I10,' seismic dislocation patches used in total solution.')") patches !Compute derived fields: DO i = 1, rows DO j = 1, columns uMagnitude(i, j) = SQRT(uEast(i, j)**2 + uNorth(i, j)**2 + uUp(i, j)**2) uHorizontal = SQRT(uEast(i, j)**2 + uNorth(i, j)**2) ! temporary scalar uPlunge(i, j) = NINT(degrees_per_radian * ATan2F(-uUp(i, j), uHorizontal)) uAzimuth = NINT(degrees_per_radian * ATan2F(uEast(i, j), uNorth(i, j))) IF (uAzimuth < 0) uAzimuth = uAzimuth + 360 uTrend(i, j) = uAzimuth END DO END DO !Create output GRD files: WRITE (*, *) WRITE (*, "(' Writing uEast.grd...')") OPEN (UNIT = 2, FILE = "uEast.grd") ! unconditional OPEN; overwrites any existing file of that name WRITE (2, "(3F16.10)") Elon_at_SW_corner, cellsize_degrees, Elon_at_SE_corner WRITE (2, "(3F16.10)") Nlat_at_SW_corner, cellsize_degrees, Nlat_at_NW_corner WRITE (2, "(10ES10.2)") ((uEast(i, j), j = 1, columns), i = 1, rows) CLOSE(2) WRITE (*, "(' Writing uNorth.grd...')") OPEN (UNIT = 3, FILE = "uNorth.grd") ! unconditional OPEN; overwrites any existing file of that name WRITE (3, "(3F16.10)") Elon_at_SW_corner, cellsize_degrees, Elon_at_SE_corner WRITE (3, "(3F16.10)") Nlat_at_SW_corner, cellsize_degrees, Nlat_at_NW_corner WRITE (3, "(10ES10.2)") ((uNorth(i, j), j = 1, columns), i = 1, rows) CLOSE(3) WRITE (*, "(' Writing uUp.grd...')") OPEN (UNIT = 4, FILE = "uUp.grd") ! unconditional OPEN; overwrites any existing file of that name WRITE (4, "(3F16.10)") Elon_at_SW_corner, cellsize_degrees, Elon_at_SE_corner WRITE (4, "(3F16.10)") Nlat_at_SW_corner, cellsize_degrees, Nlat_at_NW_corner WRITE (4, "(10ES10.2)") ((uUp(i, j), j = 1, columns), i = 1, rows) CLOSE(4) WRITE (*, "(' Writing uMagnitude.grd...')") OPEN (UNIT = 11, FILE = "uMagnitude.grd") ! unconditional OPEN; overwrites any existing file of that name WRITE (11, "(3F16.10)") Elon_at_SW_corner, cellsize_degrees, Elon_at_SE_corner WRITE (11, "(3F16.10)") Nlat_at_SW_corner, cellsize_degrees, Nlat_at_NW_corner WRITE (11, "(110ES9.2)") ((uMagnitude(i, j), j = 1, columns), i = 1, rows) CLOSE(11) WRITE (*, "(' Writing uTrend.grd...')") OPEN (UNIT = 12, FILE = "uTrend.grd") ! unconditional OPEN; overwrites any existing file of that name WRITE (12, "(3F16.10)") Elon_at_SW_corner, cellsize_degrees, Elon_at_SE_corner WRITE (12, "(3F16.10)") Nlat_at_SW_corner, cellsize_degrees, Nlat_at_NW_corner WRITE (12, "(25I4)") ((uTrend(i, j), j = 1, columns), i = 1, rows) CLOSE(12) WRITE (*, "(' Writing uPlunge.grd...')") OPEN (UNIT = 13, FILE = "uPlunge.grd") ! unconditional OPEN; overwrites any existing file of that name WRITE (13, "(3F16.10)") Elon_at_SW_corner, cellsize_degrees, Elon_at_SE_corner WRITE (13, "(3F16.10)") Nlat_at_SW_corner, cellsize_degrees, Nlat_at_NW_corner WRITE (13, "(25I4)") ((uPlunge(i, j), j = 1, columns), i = 1, rows) CLOSE(13) WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS REAL*8 FUNCTION ATan2F(y, x) !Fixes problem of undefined result for arguments of (0.0D0, 0.0D0) IMPLICIT NONE REAL*8, INTENT(IN) :: y, x IF ((x /= 0.0D0).OR.(y /= 0.0D0)) THEN ATan2F = ATAN2(y, x) ELSE ATan2F = 0.0D0 END IF END FUNCTION ATan2F SUBROUTINE Pause () IMPLICIT NONE WRITE (*, "(/' Press [Enter] to continue...')") READ (*, *) RETURN END SUBROUTINE Pause END PROGRAM Spherical_Lurch