MODULE DDislocation ! ! Computes benchmark offsets/velocities due to shear dislocations ! from earthquakes/frictional slip in the brittle part of each fault, ! at the rate predicted by a finite-element model. ! ! Originally, exactly parallel to MODULE Dislocation.f90, ! except that all REAL variables were ! replaced with DOUBLE PRECISION (REAL*8) variables, ! and subprogram names had "D" added as a ! sign that they expect DOUBLE PRECISION ! arguments and return DOUBLE PRECISION answers. ! By Peter Bird, UCLA, gradually developed 1980-2015; ! major logic revision on 2002.08.13; ! REAL*8 added (for all REALs) 2015.02.12; ! minor numerical improvement on 2015.11.03. ! THEN, in 2020.09 there was a major upgrade to routines DChange, Halo, Aura, ! so that fault-displacment (Burgers) vectors could have an opening component, ! whose remote surface displacments are computed from the solution of ! Okada [1985, BSSA]. ! Note that this upgrade was NOT retroactively added to Dislocation.f90. ! ! Usage from a F-E program like Shells/OrbScore, with fault elements: ! make a single call for the entire F-E grid: ! ! CALL DCoseis (farg,fdip,geothe,geophi,mxnode, & ! & mxfel,numgeo,nfl,nodef,radius, & ! & v,wedge,xnode,ynode,zmnode,ztranf,& ! inputs ! & geouth,geouph) ! outputs ! ! For a typical year with no earthquakes, ! correct the predicted benchmark velocities ! by subtracting the coseismic part: ! ! DO i = 1, numgeo ! predic(1, i) = predic(1, i) - geouth(i) ! predic(2, i) = predic(2, i) - geouph(i) ! END DO ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Usage from a program (like NeoKinema) where the faults are not ! represented by elements: make many calls, one for each combination ! of fault segment and benchmark. Call DChange directly: ! ! benchmark_unlocked_vw = benchmark_vw ! whole array; initialize before sum ! DO ..... fault segments ........ ! ....translate fault data into Change variables... ! DO i = 1, internal_benchmarks ! btheta = ... ! bphi = ... ! CALL DChange (argume, & ! & btheta, bphi, & ! & dipf, lf, & ! & ftheta, fphi, & ! & radius, & ! & slip, & ! & wedge, & ! & ztop, zbot, & ! inputs ! & duthet, duphi) ! output ! For a typical year with no earthquakes, ! correct the observed benchmark velocities ! by adding the (estimated) missing coseismic part: ! benchmark_unlocked_vw(2 * i - 1) = benchmark_unlocked_vw(2 * i - 1) + duthet ! benchmark_unlocked_vw(2 * i ) = benchmark_unlocked_vw(2 * i ) + duphi ! END DO ! i = 1, internal_benchmarks ! END DO ! fault segments !NOTE that this typical usage will NOT make use of the new [2020] mode-1 Okada [1985] sections. !================================================================== CONTAINS SUBROUTINE DCoseis (farg, fdip, geothe, geophi, & & mxnode, mxfel, numgeo, nfl, nodef, radius, & & v, wedge, xnode, ynode, zmnode, ztranf, & ! inputs & uTheta, uPhi) ! outputs ! Computes the horizontal-plane velocity components ! (uTheta = South, uPhi = East) for each of NUMGEO benchmark sites, ! due to earthquake dislocations on the brittle parts of faults. ! (The usual reason for doing this is so the coseismic part ! of the benchmark velocities can be subtracted from the ! velocities predicted by the F-E model before comparing to ! geodetic data in an aseismic year.) ! GEOTHE and GEOPHI should be the theta and phi coordinates of ! the benchmarks, in radians. Theta is colatitude, measured ! Southward from the North pole. Phi is latitude, measured ! Eastward from the prime meridian. The same convention ! should be used for node positions in XNODE (= theta) and ! YNODE (= phi); the use of x for theta and y for phi is ! a memnonic convention adopted when flat-planet Plates ! was rewritten as spherical-planet Shells ! Similarly, the 1st (x, or theta) component of V should be ! the Southward velocity, while the 2nd (y, or phi) component ! of V should be the Eastward velocity at the nodes. ! Most of the actual work is done by the CALL DChange ! statement; this routine is only reponsible for cutting ! each fault into several segments (for greater accuracy) ! and summing the corrections due to each segment of ! each fault. ! The following parameter NSEGME determines how many segments ! each fault element is divided into; higher values may ! give more accuracy for benchmarks close to faults with ! spatially-varying slip rates. IMPLICIT NONE INTEGER, PARAMETER :: nsegme = 3 ! Arguments (see text above): INTEGER, INTENT(IN) :: mxfel, mxnode, nfl, numgeo INTEGER, DIMENSION(4, nfl), INTENT(IN) :: nodef DOUBLE PRECISION, DIMENSION(2, mxnode), INTENT(IN) :: v REAL*8, INTENT(IN) :: radius, wedge REAL*8, DIMENSION(2,mxfel), INTENT(IN) :: farg, fdip, ztranf REAL*8, DIMENSION(mxnode), INTENT(IN) :: xnode, ynode, zmnode REAL*8, DIMENSION(numgeo), INTENT(IN) :: geothe, geophi REAL*8, DIMENSION(numgeo), INTENT(OUT) :: uTheta, uPhi DOUBLE PRECISION :: lf INTEGER :: i, ilayer, j, m, n1, n2, n3, n4 REAL*8 :: al, argume, & & crossx, crossy, & & delvx, delvy, dip, duthet, duphi, & & fhival, height, moho, & & s, smid, unitx, unity, & & vx1, vx2, vx3, vx4, vy1, vy2, vy3, vy4, & & x1, x2, x3, x4, xbegin, xend, xmid, & & y1, y2, y3, y4, ybegin, yend, ymid, & & zbot, ztop REAL*8, DIMENSION(3) :: slip !--------------------------------------------------------------------- ! STATEMENT FUNCTIONS: ! Interpolation within one fault element: fhival(s, x1, x2) = x1 + s * (x2 - x1) !--------------------------------------------------------------------- ! Initialize to zero before beginning sums: DO j = 1, numgeo uTheta(j) = 0.0D0 uPhi(j) = 0.0D0 END DO ! Sum contributions from fault elements: DO 100 i = 1, nfl n1 = nodef(1, i) n2 = nodef(2, i) n3 = nodef(3, i) n4 = nodef(4, i) x1 = xnode(n1) x2 = xnode(n2) x3 = xnode(n3) x4 = xnode(n4) y1 = ynode(n1) y2 = ynode(n2) y3 = ynode(n3) y4 = ynode(n4) vx1 = v(1, n1) vx2 = v(1, n2) vx3 = v(1, n3) vx4 = v(1, n4) vy1 = v(2, n1) vy2 = v(2, n2) vy3 = v(2, n3) vy4 = v(2, n4) xbegin = x1 ybegin = y1 ! Sum over nseg segments within each fault: DO 90 m = 1, nsegme smid = ((m * 1.0D0) - 0.5D0) / (nsegme * 1.0D0) CALL DOnarc(smid, x1, y1, x2, y2, & ! input & xmid, ymid) ! output s = (m * 1.0D0) / (nsegme * 1.0D0) IF(m == nsegme) THEN xend = x2 yend = y2 ELSE CALL DOnarc(s, x1, y1, x2, y2, & ! input & xend, yend) ! output END IF al = DFltLen(ybegin, yend, radius, xbegin, xend) lf = 0.5D0 * al ! FARG is the direction from node 1 to node 2, ! in radians counterclockwise from the ! +theta (South) direction, and so is ARGUME: argume = DChord(farg(1, i), smid * 1.0D0, farg(2, i)) dip = fhival(smid, fdip(1, i), fdip(2, i)) ! We consider that the node-1,2 side of the fault ! moves while the node-3,4 side is fixed: delvx = fhival(smid, vx1, vx2) - fhival(smid, vx4, vx3) delvy = fhival(smid, vy1, vy2) - fhival(smid, vy4, vy3) ! UNIT is a horizontal unit vector pointing ! along the fault from the node-1 end to the node-2 end; ! this will be the same as the Z1 direction of HALO if ! the fault is vertical, and equal or opposite to the ! X1 direction of AURA if the fault is dipping: unitx = DCOS(argume) unity = DSIN(argume) ! CROSS is a horizontal unit vector pointing ! across the fault, 90 degrees clockwise from UNIT ! (when viewed from outside the planet); ! this will be the same as the Z2 direction of HALO ! if the fault is exactly vertical or ! equal or opposite to the X2 direction of AURA ! if the fault is dipping. crossx = + unity crossy = -unitx slip(1) = delvx * unitx + delvy * unity ! Note: positive for sinistral slip slip(2) = delvx * crossx + delvy * crossy ! Note: positive for opening across trace ! Vertical component is positive when down: IF (ABS(1.570796D0 - dip) > wedge) THEN slip(3) = slip(2) * DTAN(dip) ELSE slip(3) = 0.0D0 END IF ! Consider both crustal and mantle patches: moho = fhival(smid, zmnode(n1), zmnode(n2)) DO 80 ilayer = 1, 2 height = ztranf(ilayer, i) IF (height > 0.0D0) THEN IF (ilayer == 1) THEN ztop = 0.0D0 ELSE ztop = moho END IF zbot = ztop + height ! Compute effects of patch at all benchmarks: DO 70 j = 1, numgeo ! ------------------------------------------- CALL DChange(argume, & & geothe(j), geophi(j), & & dip, lf, & & xmid, ymid, & & radius, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output ! ------------------------------------------- uTheta(j) = uTheta(j) + duthet uPhi(j) = uPhi(j) + duphi 70 CONTINUE END IF 80 CONTINUE ! Prepare to loop to next segment: xbegin = xend ybegin = yend 90 CONTINUE 100 CONTINUE END SUBROUTINE DCoseis !----------------------------------------------------------------------- SUBROUTINE DChange (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & radius, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi, dur) ! output ! This routine is a "driver" or "wrapper" for the Cartesian (flat-Earth) ! *Mansinha & Smylie routine Halo (vertical shear-dislocation patch) ! or Aura (dipping shear-dislocation patch) ! *and tensile (mode-1; opening) crack routine of Okada [1985, BSSA], ! which permits them to be called from spherical-planet ! programs because it converts coordinates from ! (theta, phi) = (colatitude, longitude) in radians ! to the (z1,z2,z3) Cartesian units of a locally flat planet. ! Also, output (horizontal) vector (DUTHET, DUPHI) is expressed ! in the spherical-planet (theta, phi) system, so components are ! (South, East). Units are the same as input vector SLIP, ! which may be either a relative displacement or ! a relative velocity of the two sides of the fault. ! The third component of the output displacement[-rate] vector, ! DUR (radial, or up) is OPTIONAL. ! FTHETA and FPHI should give (theta, phi) coordinates of the ! midpoint of the trace of the plane containing the ! rectangular dislocation patch. ! BTHETA and BPHI should give (theta, phi) coordinates of the ! benchmark or test point at which the displacement (-rate) ! is desired. ! LF is the half-length (from center to end) of the dislocation ! patch, measured along a horizontal strike line. ! Units are the same as RADIUS, ZTOP, ZBOT (below). ! ARGUME is the argument of the trace of the dislocation patch, ! measured in radians counterclockwise from +theta (from South). ! DIPF is fault dip in radians measured clockwise ! (initially, down) from horizontal on the right side of the ! fault (when looking in direction ARGUME). ! The fault/crack dislocation vector (Burgers vector) ! SLIP must be already rotated into fault-trace-centered coordinates: ! SLIP(1) is the component parallel to the fault, ! and it is positive for sinistral sense of slip. ! SLIP(2) is the horizontal component in the direction ! perpendicular to the trace of the fault, ! and it is positive for divergence across the trace. ! SLIP(3) is the relative vertical component, and it ! is positive when the right side of the fault ! (when looking along direction ARGUME) is down. ! Note that SLIP(2) > 0 with SLIP(3) > 0 is associated with normal faulting; ! while SLIP(2) < 0 with SLIP(3) < 0 is associated with thrust faulting; ! and SLIP(2) > 0 with SLIP(3) < 0 is associated with tensile (mode-1) opening. ! WEDGE is a tolerance (in radians) for dip; if DIPF is ! within WEDGE of Pi/2, then fault is considered vertical ! and routine Halo is called; otherwise, Aura is called. ! ZTOP and ZBOT are (positive) depths to top and bottom of ! the dislocation patch, respectively. Units are the same ! as RADIUS, which gives the size of the planet. ! Note that a 0.0 input value will be accepted for ZTOP, ! but in the actual calculation a small positive depth ! will be substituted for this zero. ! NOTE: If distance from "F" point to "B" point exceeds ! one planetary radius, result is approximated as zero and ! Halo/Aura are never called. IMPLICIT NONE ! ------------------------------------------------------- ! Arguments (see text above): REAL*8, INTENT(IN) :: argume, btheta, bphi, dipf, & & ftheta, fphi, radius, wedge, zbot, ztop REAL*8, DIMENSION(3), INTENT(IN) :: slip REAL*8, INTENT(OUT) :: duthet, duphi ! S and E components REAL*8, INTENT(OUT), OPTIONAL :: dur ! radial or up component ! ------------------------------------------------------- ! Following internal variables must be DP to agree with Halo, Aura: DOUBLE PRECISION :: dbot, dtop, lf, theta DOUBLE PRECISION, DIMENSION(3) :: u, ucap, x, z ! Following are DP for precision in determination of ! variables FAR and ALPHA. (At short distances, ! unit vectors must be DP or location precision is ! lost and benchmark can end up on wrong side of fault!) DOUBLE PRECISION :: crossp, dot1, dot2, dotpro DOUBLE PRECISION, DIMENSION(3) :: buvec, dvec, fuvec, tvec, uPhi, uTheta INTEGER :: i LOGICAL :: vert REAL*8 :: alpha, cosaz, far, sinaz, tazimf ! Compute position relative to fault trace; ! this is where we convert from a spherical-Earth ! to a flat-Earth model, implicitly using a gnomonic ! projection (which preserves distance and azimuth ! from the center of the trace of the dislocation plane). fuvec(1) = DSIN(ftheta) * DCOS(fphi) fuvec(2) = DSIN(ftheta) * DSIN(fphi) fuvec(3) = DCOS(ftheta) buvec(1) = DSIN(btheta) * DCOS(bphi) buvec(2) = DSIN(btheta) * DSIN(bphi) buvec(3) = DCOS(btheta) dotpro = fuvec(1) * buvec(1) + fuvec(2) * buvec(2) + fuvec(3) * buvec(3) tvec(1) = fuvec(2) * buvec(3) - fuvec(3) * buvec(2) tvec(2) = fuvec(3) * buvec(1) - fuvec(1) * buvec(3) tvec(3) = fuvec(1) * buvec(2) - fuvec(2) * buvec(1) crossp = DSQRT(tvec(1)**2 + tvec(2)**2 + tvec(3)**2) far = radius * DATAN2(crossp, dotpro) ! !Notice of important change: !Prior to October 2015 the test condition was: !IF (far > (40.0D0 * zbot)) THEN !However, while plotting some very high-resolution (0.01 degree x 0.01 degree) !models of coseismic strains, created with StrainRates2015, I noticed that !this gives an objectionable "cutoff" artifact along a circular arc. !Thus, the outer radius for application of these solutions needs to be much bigger. !On the other hand, these flat-Earth solutions should certainly NOT be trusted !in the opposite hemisphere! As a compromise, I now cut them off at an arc !distance of one radian from the center of the trace of the dislocation patch: !The result will NOT be accurate in the far-field, but at least it will be !smoothly going to zero, so that it does not contaminate computed strain-rates. IF (far > radius) THEN duthet = 0.0D0 duphi = 0.0D0 ELSE dvec(1) = buvec(1) - fuvec(1) dvec(2) = buvec(2) - fuvec(2) dvec(3) = buvec(3) - fuvec(3) uTheta(1) = DCOS(ftheta) * DCOS(fphi) uTheta(2) = DCOS(ftheta) * DSIN(fphi) uTheta(3) = -DSIN(ftheta) dot1 = dvec(1) * uTheta(1) + dvec(2) * uTheta(2) + dvec(3) * uTheta(3) uPhi(1) = -DSIN(fphi) uPhi(2) = DCOS(fphi) uPhi(3) = 0.0D0 dot2 = dvec(1) * uPhi(1) + dvec(2) * uPhi(2) + dvec(3) * uPhi(3) alpha = DATAN2(dot2, dot1) z(1) = far * DCOS(alpha - argume) z(2) = -far * DSIN(alpha - argume) ! Note: Z(3) = 0 because we believe it is more important ! to convey that the benchmark is on the free surface ! than it is to convey the precise angular relationship ! to the dislocation. However, some day one might test ! a positive Z(3), in proportion to square of "far". z(3) = 0.0D0 vert = ABS(dipf - 1.5708) <= wedge IF (vert) THEN ! vertical fault or crack, so use Halo: tazimf = argume theta = dipf ! Z values computed above are still good ucap(1:3) = slip(1:3) dtop = ztop dbot = zbot ! Kludge to insure that input parameter DTOP is not zero, ! because, in that case, HALO may give numerically-unstable ! ("bad") results for near-field test points: dtop = MAX(dtop, 0.0001D0 * dbot) ! - - - - - - - - - - - - - - - - - - - - - - - CALL Halo (ucap, z, lf, dbot, dtop, & ! inputs & u) ! output ! - - - - - - - - - - - - - - - - - - - - - - - ELSE ! non-vertical, dipping fault, so use Aura: IF (dipf <= 1.5708D0) THEN tazimf = argume theta = dipf DO 10 i = 1, 3 x(i) = z(i) ucap(i) = slip(i) 10 CONTINUE ELSE ! Look at this fault from the other end/side, ! so that dip THETA will be less than Pi/2: tazimf = argume + 3.14159265358979D0 theta = 3.14159265358979D0 - dipf x(1) = -z(1) x(2) = -z(2) x(3) = z(3) ucap(1) = slip(1) ucap(2) = slip(2) ucap(3) = -slip(3) ! Note: Reversal of Z1, Z2 axes (but not Z3) ! during name-change to X1,X2,X3 axes (for AURA) ! combined with change of "moving" block ! leaves "sinistral" and "opening" components ! unchanged, while "relative vertical" component ! reverses. END IF dtop = ztop / DSIN(theta) dbot = zbot / DSIN(theta) ! Kludge to insure that input parameter DTOP is not zero, ! because, in that case, AURA may give numerically-unstable ! ("bad") results for near-field test points: dtop = MAX(dtop, 0.0001D0 * dbot) ! - - - - - - - - - - - - - - - - - - - - - - - CALL Aura (ucap, theta, x, lf, dbot, dtop, & ! inputs & u) ! output ! - - - - - - - - - - - - - - - - - - - - - - - END IF sinaz = DSIN(tazimf) cosaz = DCOS(tazimf) duthet = + cosaz * u(1) + sinaz * u(2) duphi = + sinaz * u(1) - cosaz * u(2) IF (PRESENT(dur)) THEN dur = -u(3) ! changing sign because both Halo and Aura use ! the Mansinha & Smylie positive-downwards convention, ! but this routine uses a positive-upwards convention. END IF END IF END SUBROUTINE DChange !----------------------------------------------------------------- SUBROUTINE Halo (ucap, z, l, dbot, dtop, & ! inputs & u) ! output ! Computes displacements in a uniform elastic half-space ! of Poisson's ratio 0.25 due to uniform slip (and/or opening) on a ! rectangular dislocation patch in a vertical plane, ! whose horizontal and vertical edges form a rectangle. ! Equations for fault-shear components based on Mansinha and Smylie (1967), ! J. Geophys. Res., v. 72, p. 4731-4743. ! Equations for opening/closing (mode-1) components based on Okada [1985], ! Bull. Seismol. Soc. Am., 74(4), 1135-1154. USE DSphere ! provided by Peter Bird as file DSphere.f90 IMPLICIT NONE ! Arguments: DOUBLE PRECISION, INTENT(IN) :: l, dbot, dtop DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: ucap, z DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: u INTEGER :: i DOUBLE PRECISION :: d_Okada, del_Okada, dTilde_Okada, epsilon_Okada, factor, I1_Okada, I3_Okada, I5_Okada, l_Okada, & & mu_over_lambda_plus_mu, nu_Okada, p_Okada, q, q_Okada, r, R_Okada, R2_Okada, & & term1, term2, term3, term4, term5, term6, u3_Okada, & & ux_Okada, uy_Okada, uz_Okada, & & w_Okada, x_Okada, XCap_Okada, XCap2_Okada, y_Okada, yTilde_Okada, z_Okada, & & zeta1, zeta3 ! The coordinate system is Cartesian and right-handed, ! with z3 down, z1 parallel to the fault/crack trace, and z2 ! perpendicular. The dislocation surface is at z2 = 0, ! DTOP <= z3 <= DBOT (note that DTOP > 0.), and ! -L <= z1 <= L. ! UCAP contains the three components of the displacement of the ! positive-z2 side relative to the other side; ! U contains the three components of the displacement ! at the observation point Z. u = 0.0D0 ! all 3 components initialized DO 100 i = 1, 4 factor = 1.0D0 IF ((i == 2).OR.(i == 3)) factor = -1.0D0 zeta1 = l IF ((i == 3).OR.(i == 4)) zeta1 = -l zeta3 = dbot IF ((i == 2).OR.(i == 4)) zeta3 = dtop r = DSQRT((z(1) - zeta1)**2 + z(2)**2 + (z(3) - zeta3)**2) q = DSQRT((z(1) - zeta1)**2 + z(2)**2 + (z(3) + zeta3)**2) term1 = z(2) * (z(1) - zeta1) * (2.0D0 / (r * (r + z(3) - zeta3)) & & - (5.0D0 * q + 8.0D0 * zeta3) / (2.0D0 * q * (q + z(3) + zeta3)**2) & & + 4.0D0 * z(3) * zeta3 * (2.0D0 * q + z(3) + zeta3) / & & (q**3 * (q + z(3) + zeta3)**2)) & & + 3.0D0 * DAtanPV((z(1) - zeta1) * (z(3) - zeta3), (r * z(2))) & & - 3.0D0 * DAtanPV((z(1) - zeta1) * (z(3) + zeta3), (q * z(2))) term2 = -DLOG(r + z(3) - zeta3) & & + 0.5D0 * DLOG(q + z(3) + zeta3) & & - (5.0D0 * z(3) - 3.0D0 * zeta3) / (2.0D0 * (q + z(3) + zeta3)) & & - 4.0D0 * z(3) * zeta3 / (q * (q + z(3) + zeta3)) & & + z(2)**2 * (2.0D0 / (r * (r + z(3) - zeta3)) & & - (5.0D0 * q + 8.0D0 * zeta3) / (2.0D0 * q * (q + z(3) + zeta3)**2) & & + 4.0D0 * z(3) * zeta3 * (2.0D0 * q + z(3) + zeta3) / & & (q**3 * (q + z(3) + zeta3)**2)) TERM3=Z(2)*(2.0D0/R-2.0D0/Q+4.0D0*Z(3)*ZETA3/Q**3 & & +3.0D0/(Q+Z(3)+ZETA3)+2.0D0*(Z(03)+3.0D0*ZETA3)/ & & (Q*(Q+Z(3)+ZETA3))) TERM4=Z(2)*(2.0D0/R+4.0D0/Q-4.0D0*Z(3)*ZETA3/Q**3 & & -6.0D0*Z(3)/(Q*(Q+Z(3)+ZETA3))) TERM5= -DLOG(R+Z(1)-ZETA1)+DLOG(Q+Z(1)-ZETA1) & & +(6.0D0*Z(3)**2+10.0D0*Z(3)*ZETA3)/(Q*(Q+Z(1)-ZETA1)) & & +(6.0D0*Z(3)*(Z(1)-ZETA1))/(Q*(Q+Z(3)+ZETA3)) & & +Z(2)**2*(2.0D0/(R*(R+Z(1)-ZETA1)) & & +4.0D0/(Q*(Q+Z(1)-ZETA1)) & & -4.0D0*Z(3)*ZETA3*(2.0D0*Q+Z(1)-ZETA1)/(Q**3*(Q+Z(1)-ZETA1)**2)) TERM6=Z(2)*(2.0D0*(Z(3)-ZETA3)/(R*(R+Z(1)-ZETA1)) & & -2.0D0*(Z(3)+2.0D0*ZETA3)/(Q*(Q+Z(1)-ZETA1)) & & -4.0D0*Z(3)*ZETA3*(Z(3)+ZETA3)*(2.0D0*Q+Z(1)-ZETA1)/ & & (Q**3*(Q+Z(1)-ZETA1)**2)) & & +3.0D0*DATANPV((Z(1)-ZETA1)*(Z(3)-ZETA3),(R*Z(2))) & & -3.0D0*DATANPV((Z(1)-ZETA1)*(Z(3)+ZETA3),(Q*Z(2))) !Begin with strike-slip component (most common application): u(1) = u(1) + factor * term1 * ucap(1) / 37.69911184D0 u(2) = u(2) + factor * term2 * ucap(1) / 37.69911184D0 U(3)=U(3)+FACTOR*TERM3*UCAP(1)/37.69911184D0 !Is there any vertical slip on this dislocation? !NOTE that the vertical component U(3) is positive-downwards. IF (ucap(3) /= 0.0D0) THEN U(1)=U(1)+FACTOR*TERM4*UCAP(3)/37.69911184D0 U(2)=U(2)+FACTOR*TERM5*UCAP(3)/37.69911184D0 U(3)=U(3)+FACTOR*TERM6*UCAP(3)/37.69911184D0 END IF 100 CONTINUE !------end of Mansinha & Smylie [1967] section--------------------------- !------beginning of Okada [1985, BSSA] section -------------------------- !Is there any opening/closing (mode-1 component) on this dislocation? IF (ucap(2) /= 0.0D0) THEN ! tensile / mode-1 / opening component: !Note that Okada uses a right-handed Cartesian coordinate system (like M&S) !but that his has z increasing UP, and therefore if we make his x axis (along-strike) !equal to the z1 axis of M&S, then Okada's y axis is opposite to the z2 axis of M&S. !Another distinction is that M&S have their (z1, z2) origin at a central point on !the fault trace (projected to the free surface), but Okada's (x, y) origin !is located over the deepest part of the opening patch. !Dip of crack/fault: del_Okada = Pi_over_2 ! defined in DSphere.f90 (== 90 degrees; only here in Halo) !Observer (test) location [which must be on the surface, using formulas of Okada, 1985]: x_Okada = z(1) y_Okada = -z(2) !N.B. Fortunately, for a vertical patch, the formulas above do NOT require any offset terms ! between M&S's (z1, z2) origin point and Okada's (x, y) origin point. ! This would not be true for other dips of the fault patch! z_Okada = 0.0D0 ! (On free surface. To change this, use formulas of Okada, 1992 instead of Okada, 1985.) !Dimensions of dislocation patch: d_Okada = dbot w_Okada = dbot - dtop ! should be positive. Note formula would be more complex for non-vertical dip. L_Okada = l ! half-length, along strike (same as M&S) !Okada dimensional parameters dependent on the above: p_Okada = d_Okada ! special simple form for vertical dip, when del_Okada = Pi_over_2 q_Okada = y_Okada ! only for vertical dip !Amount of opening {rate?}: u3_Okada = ucap(2) !Poisson's ratio term: mu_over_lambda_plus_mu = 0.50D0 ! (assumed here in Halo; equivalent to Poisson's ratio of 1/4) !Initialize before sum of 4 terms: ux_Okada = 0.0D0 uy_Okada = 0.0D0 uz_Okada = 0.0D0 DO 200 i = 1, 4 ! sum over 4 terms (from rectangular integration limits in 2-D): !Integration limits and associated sign-changing factor: factor = 1.0D0 ! unless... IF ((i == 2).OR.(i == 3)) factor = -1.0D0 ! replacing the choice made in the line above epsilon_Okada = x_Okada + L_Okada ! unless... IF (i >= 3) epsilon_Okada = x_Okada - L_Okada ! replacing the choice made in the line above nu_Okada = p_Okada ! unless... IF ((i == 2).OR.(i == 4)) nu_Okada = p_Okada - w_Okada ! replacing the choice made in the line above !Geometric factors dependent on the above: yTilde_Okada = q_Okada ! only for vertical dip dTilde_Okada = nu_Okada R2_Okada = epsilon_Okada**2 + yTilde_Okada**2 + dTilde_Okada**2 R_Okada = SQRT(R2_Okada) XCap2_Okada = epsilon_Okada**2 + q_Okada**2 XCap_Okada = SQRT(XCap2_Okada) I1_Okada = -0.5D0 * mu_over_lambda_plus_mu * epsilon_Okada*q_Okada/(R_Okada+dTilde_Okada)**2 I3_Okada = +0.5D0 * mu_over_lambda_plus_mu * ( (nu_Okada/(R_Okada+dTilde_Okada)) + & & (yTilde_Okada*q_Okada/(R_Okada+dTilde_Okada)**2) - & & LOG(R_Okada+nu_Okada) ) I5_Okada = -mu_over_lambda_plus_mu * epsilon_Okada/(R_Okada+dTilde_Okada) ux_Okada = ux_Okada + factor * (u3_Okada/Two_Pi) * ( (q_Okada**2/(R_Okada*(R_Okada+nu_Okada))) - & & I3_Okada ) uy_Okada = uy_Okada + factor * (u3_Okada/Two_Pi) * ( (-dTilde_Okada*q_Okada/(R_Okada*(R_Okada+epsilon_Okada))) - & & ((epsilon_Okada*q_Okada/(R_Okada*(R_Okada+nu_Okada))) - (ATAN((epsilon_Okada*nu_Okada) / (q_Okada*R_Okada)))) - & & I1_Okada ) uz_Okada = uz_Okada + factor * (u3_Okada/Two_Pi) * ( (yTilde_Okada*q_Okada/(R_Okada*(R_Okada+epsilon_Okada))) - & & I5_Okada ) 200 CONTINUE ! sum over 4 terms !Final contribution (reversing signs of two terms, per reversal of axis-directions between M&S versus O): u(1) = u(1) + ux_Okada u(2) = u(2) - uy_Okada u(3) = u(3) - uz_Okada END IF ! there is a tensile / mode-1 / opening component of source motion !------end of Okada [1985, BSSA] section -------------------------------- END SUBROUTINE Halo !----------------------------------------------------------------- SUBROUTINE Aura (u4, theta4, x4, l4, dbot4, dtop4, & ! inputs & v4) ! output ! Computes displacements in an infinite uniform elastic halfspace ! of Poisson's ratio 0.25 caused by uniform slip over the surface ! of a buried rectangular dislocation with horizontal and plunging ! sides. To calculate displacements caused by more general slip, ! use this routine many times in the kernel of an integral, ! assigning the local average slip over each very small rectangle. ! Equations for shear-displactions based on Mansinha and Smylie (1971), ! Bull. Seism. Soc. Amer., v. 61, p. 1433-1440. ! This code has been extended [2020.09] to also model tensile/mode-1/opening ! components of crack dislocations, according to the solution of ! Okada [1985], Bull. Seismol. Soc. Amer., v. 75, #4, p. 1135-1154. USE DSphere ! provided by Peter Bird as file DSphere.f90 IMPLICIT NONE ! Arguments: DOUBLE PRECISION, INTENT(IN) :: dbot4, dtop4, l4, theta4 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: u4, x4 DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: v4 ! Internal variables: INTEGER :: i DOUBLE PRECISION :: costh, d_Okada, dbot, del_Okada, deld, dTilde_Okada, dtop, & & epsilon_Okada, factor, h, & & I1_Okada, I3_Okada, I4_Okada, I5_Okada, & ! N.B. I2_Okada (for strike-slip) is not needed. & k, l, L_Okada, & & mu_over_lambda_plus_mu, & & nu_Okada, p_Okada, q, q_Okada, q2, q3, qsquar, & & r, R_Okada, r2, R2_Okada, r3, rsquar, & & secth, sinth, tanth, & & term1, term2, term3, term4, term5, term6, theta, & & ucap, ucap1, u3_Okada, ux_Okada, uy_Okada, uz_Okada, w_Okada, & & x_Okada, xCap_Okada, xCap2_Okada, xi, xi1, xi2, xi3, y_Okada, yTilde_Okada, z_Okada DOUBLE PRECISION, DIMENSION(3) :: u(3), v(3), x(3) !Copy all input values to new variables, to avoid inadvertent changes! DO i = 1, 3 u(i) = u4(i) x(i) = x4(i) END DO theta = theta4 l = l4 dbot = dbot4 dtop = dtop4 ! U(3) is the dislocation vector (i.e. the total offset vector, ! sometimes called the Burgers vector), ! in a Cartesian coordinate system where x1 is parallel to the ! fault strike, x2 is also horizontal and pointing 90 degrees to the ! right of +x1 (that is, pointing to the hanging wall, ! assuming that dip is less than Pi/2 radians!), ! and x3 is down. (Take care to make these right-handed.) ! THETA is the fault dip angle (in radians) from +x2 axis. ! X(3) gives the observation point in these coordinates. ! L is the half-length of the rectangle, along strike. ! DBOT and DTOP are oblique distances in the fault plane (and ! down the dip) from the surface to the bottom and top of the ! rectangular dislocation, respectively. ! Note: Fault does not break the surface, so DTOP > 0., ! although it may be small. If 0.0 is supplied for DTOP4, ! a small positive depth will be substituted in the calculation. ! V(3) is the answer: the 3-component displacement at the ! observation point in same units as U and in x1,x2,x3 coordinates. !Fault/crack dip section: sinth = SIN(theta) costh = COS(theta) tanth = sinth / costh secth = 1.0D0 / costh ! NOTE: FAULT MAY NOT BE VERTICAL IN THIS ROUTINE. ! (To handle vertically-dipping faults, CALL Halo instead.) !Dislocation (Burgers) vector section: ! ucap1 is the strike-slip component, with sinistral (left-lateral) considered positive ucap1 = u(1) !ucap is the dip-slip (shear) component, with normal-faulting considered positive: ucap = costh * u(2) + sinth * u(3) !u3_Okada (not used until the 2nd half of this code, in the Okada [1985] section) ! is the tensile / mode-1 / crack-opening component. u3_Okada = sinth * u(2) - costh * u(3) ! (Of course, crack-closing can be represented by negative crack-opening.) !NOTE that Aura would typically execute faster if some of these components ! were exactly zero, so that code-sections could be skipped. ! However, that rarely happens in REAL*8 floating-point computations! ! Extra LOGICAL arguments could be added to control this. ! Initialize remote displacement {rate} to zero, before adding terms. DO i = 1, 3 v(i) = 0.0D0 END DO ! Note that v(1:3) follows Mansinha & Smylie [1971] sign-conventions, ! with downwards displacement considered positive. !---------- begin Mansinha & Smylie [1971] component (shearing only) ----------------- r2 = x(2) * sinth - x(3) * costh r3 = x(2) * costh + x(3) * sinth q2 = x(2) * sinth + x(3) * costh q3 = -x(2) * costh + x(3) * sinth DO 100 i = 1, 4 factor = 1.0D0 IF ((i == 2).OR.(i == 3)) factor = -1.0D0 xi1 = l IF (i >= 3) xi1 = -l xi = dbot IF ((i == 2).OR.(i == 4)) xi = dtop xi2 = xi * costh xi3 = xi * sinth rsquar = (x(1) - xi1)**2 + (x(2) - xi2)**2 + (x(3) - xi3)**2 qsquar = (x(1) - xi1)**2 + (x(2) - xi2)**2 + (x(3) + xi3)**2 k = DSQRT(qsquar - (q3 + xi)**2) h = DSQRT(qsquar - (x(1) - xi1)**2) q = DSQRT(qsquar) r = DSQRT(rsquar) !Impose arbitrary limit on how close test point can be to the (extended) fault plane; !without this limit the equations below may give 0./0. indeterminacy and an ABEND. !The following version of the limit was historically used in Dislocation.f90 through 2015, !q = MAX(q, 1.00001D0 * ABS(x(1) - xi1)) !r = MAX(r, 1.00001D0 * ABS(x(1) - xi1)) !but on 2015.11.03 I found that the second version is preferable, for both !Dislocation.f90 (which uses DOUBLE PRECISION in this routine) and for DDislocation.f90: q = MAX(q, 1.0000000001D0 * ABS(x(1) - xi1)) r = MAX(r, 1.0000000001D0 * ABS(x(1) - xi1)) term1 = (x(1) - xi1) * (2.0D0 * r2 / (r * (r + r3 - xi)) & & - (4.0D0 * q2 - 2.0D0 * x(3) * costh) / (q * (q + q3 + xi)) & & - 3.0D0 * tanth / (q + x(3) + xi3) & & + 4.0D0 * q2 * x(3) * sinth / q**3 & & - 4.0D0 * q2 * q3 * x(3) * sinth * (2.0D0 * q + q3 + xi) / & & (q**3 * (q + q3 + xi)**2)) & & - 6.0D0 * tanth**2 * DAtanPV(((k - q2 * costh) * (q - k) + & & (q3 + xi) * k * sinth), ((x(1) - xi1) * (q3 + xi) * costh)) & & + 3.0D0 * DAtanPV((x(1) - xi1) * (r3 - xi), (r2 * r)) & & - 3.0D0 * DAtanPV((x(1) - xi1) * (q3 + xi), (q2 * q)) term2 = sinth * (3.0D0 * tanth * secth * Dlog(q + x(3) + xi3) - & & Dlog(r + r3 - xi) - (1.0D0 + 3.0D0 * tanth**2) * Dlog & & (q + q3 + xi)) + 2.0D0 * r2**2 * sinth / (r * (r + r3 - xi)) & & + 2.0D0 * r2 * costh / r - 2.0D0 * sinth * & & (2.0D0 * x(3) * (q2 * costh - q3 * sinth) + q2 * (q2 + x(2) * sinth)) & & / (q * (q + q3 + xi)) - 3.0D0 * tanth * (x(2) - xi2) / & & (q + x(3) + xi3) + 2.0D0 * (q2 * costh - q3 * sinth - & & x(3) * sinth**2) / q + 4.0D0 * q2 * x(3) * sinth * & & ((x(2) - xi2) + q3 * costh) / q**3 - 4.0D0 * q2**2 * q3 * x(3) & & * sinth**2 * (2.0D0 * q + q3 + xi) / (q**3 * (q + q3 + xi)**2) TERM3=COSTH*(DLOG(R+R3-XI)+(1.0D0+3.0D0*TANTH**2)*DLOG( & & Q+Q3+XI)-3.0D0*TANTH*SECTH*DLOG(Q+X(3)+XI3)) & & +2.0D0*R2*SINTH/R+2.0D0*SINTH*(Q2+X(2)*SINTH)/Q & & -2.0D0*R2**2*COSTH/(R*(R+R3-XI))+ & & (4.0D0*Q2*X(3)*SINTH**2-2.0D0*(Q2+X(2)*SINTH)*(X(3)+Q3*SINTH))/ & & (Q*(Q+Q3+XI))+ & & 4.0D0*Q2*X(3)*SINTH*(X(3)+XI3-Q3*SINTH)/Q**3 & & -4.0D0*Q2**2*Q3*X(3)*COSTH*SINTH*(2.0D0*Q+Q3+XI)/ & & (Q**3*(Q+Q3+XI)**2) term4 = (x(2) - xi2) * sinth * (2.0D0 / r + 4.0D0 / q - 4.0D0 * xi3 * x(3) / q**3 & & - 3.0D0 / (q + x(3) + xi3)) - costh * (3.0D0 * Dlog(q + x(3) + xi3) & & + 2.0D0 * (x(3) - xi3) / r + 4.0D0 * (x(3) - xi3) / q + 4.0D0 * xi3 * x(3) * (x(3) + xi3) & & / q**3) + 3.0D0 * (Dlog(q + x(3) + xi3) - sinth * Dlog(q + & & q3 + xi)) / costh + 6.0D0 * x(3) * (costh / q - q2 * sinth / & & (q * (q + q3 + xi))) term5 = sinth * (-Dlog(r + x(1) - xi1) + Dlog(q + x(1) - xi1) & & + 4.0D0 * xi3 * x(3) / (q * (q + x(1) - xi1)) + 3.0D0 * (x(1) - xi1) / & & (q + x(3) + xi3) + (x(2) - xi2)**2 * (2.0D0 / (r * (r + x(1) - xi1)) & & + 4.0D0 / (q * (q + x(1) - xi1)) - 4.0D0 * xi3 * x(3) * (2.0D0 * q + x(1) - xi1) / & & (q**3 * (q + x(1) - xi1)**2))) & & - costh * ((x(2) - xi2) * (2.0D0 * (x(3) - xi3) / (r * (r + x(1) - xi1)) & & + 4.0D0 * (x(3) - xi3) / (q * (q + x(1) - xi1)) + 4.0D0 * xi3 * x(3) * & & (x(3) + xi3) * (2.0D0 * q + x(1) - xi1) / (q**3 * (q + x(1) - xi1)**2)) & & + 6.0D0 * DAtanPV((x(1) - xi1) * (x(2) - xi2), ((h + x(3) + xi3) * (q + h))) & & - 3.0D0 * DAtanPV((x(1) - xi1) * (r3 - xi), (r2 * r)) + & & 6.0D0 * DAtanPV((x(1) - xi1) * (q3 + xi), (q2 * q))) + & & 6.0D0 * (secth * DAtanPV(((k - q2 * costh) * (q - k) + (q3 + xi) * k * sinth) & & , ((x(1) - xi1) * (q3 + xi) * costh)) & & + x(3) * (((sinth**2 - costh**2) * (q3 + xi) + 2.0D0 * q2 * & & costh * sinth) / (q * (q + x(1) - xi1)) + (x(1) - xi1) * sinth**2 & & / (q * (q + q3 + xi)))) TERM6=SINTH*((X(2)-XI2)*(2.0D0*(X(3)-XI3)/(R*(R+X(1)-XI1)) & & +4.0D0*(X(3)-XI3)/(Q*(Q+X(1)-XI1))-4.0D0*XI3*X(3)* & & (X(3)+XI3)*(2.0D0*Q+X(1)-XI1)/(Q**3*(Q+X(1)-XI1)**2)) & & -6.0D0*DATANPV((X(1)-XI1)*(X(2)-XI2),((H+X(3)+XI3)*(Q+H))) & & +3.0D0*DATANPV((X(1)-XI1)*(R3-XI),(R2*R))- & & 6.0D0*DATANPV((X(1)-XI1)*(Q3+XI),(Q2*Q))) & & +COSTH*(DLOG(R+X(1)-XI1)-DLOG(Q+X(1)-XI1)- & & 2.0D0*(X(3)-XI3)**2/(R*(R+X(1)-XI1)) & & -4.0D0*((X(3)+XI3)**2-XI3*X(3))/(Q*(Q+X(1)-XI1)) & & -4.0D0*XI3*X(3)*(X(3)+XI3)**2*(2.0D0*Q+X(1)-XI1)/ & & (Q**3*(Q+X(1)-XI1)**2)) & & +6.0D0*X(3)*(COSTH*SINTH*(2.0D0*(Q3+XI)/(Q*(Q+X(1)-XI1)) & & +(X(1)-XI1)/(Q*(Q+Q3+XI))) & & -Q2*(SINTH**2-COSTH**2)/(Q*(Q+X(1)-XI1))) deld = dbot - dtop v(1) = v(1) + factor * (term1 * ucap1 + term4 * ucap) / 37.69911184D0 v(2) = v(2) + factor * (term2 * ucap1 + term5 * ucap) / 37.69911184D0 V(3) = V(3) + FACTOR * (TERM3 * UCAP1 + TERM6 * UCAP) / 37.69911184D0 100 CONTINUE !---------- end Mansinha & Smylie [1971] component (shearing only) ------------------- !---------- begin Okada [1985] component (tensile/mode-1/opening cracks) ------------- IF (u3_Okada /= 0.0D0) THEN ! Note this this test is ALMOST ALWAYS passed, and therefore this Okada section almost always runs. ! However, its contribution may be infinitesimal if the dislocation is basically of shear/fault type. ! Aura would (typically) run faster if the test were for (ABS(u3_Okada) > n.nnDnn), ! or if extra LOGICAL arguments were added to control which code sections are run! !Note that Okada uses a right-handed Cartesian coordinate system (like M&S) !but that his has z increasing UP, and therefore if we make his x axis (along-strike) !equal to the z1 axis of M&S, then Okada's y axis is opposite to the z2 axis of M&S. !Another distinction is that M&S have their (z1, z2) origin at a central point on !the fault trace (projected to the free surface), but Okada's (x, y) origin !is located over the deepest part of the opening patch. !Poisson's ratio term: mu_over_lambda_plus_mu = 0.50D0 ! (assumed here in Aura; equivalent to Poisson's ratio of 1/4) !Dip of crack/fault: del_Okada = theta !Dimensions of dislocation patch: d_Okada = dbot * sinth ! measured vertically, not down-dip w_Okada = dbot - dtop ! measured down-dip; should be positive. L_Okada = l ! half-length, along strike (same as M&S) ! NOTE that Figure 1 of Okada [1985] shows half of the fault-patch with dashed-lines; ! we are adopting the convention that this "half" is part of the source; ! Okada explains in his text how to adjust his formulas [see * below] if this is so. !Observer (test) location [which must be on the surface, using formulas of Okada, 1985]: x_Okada = +x(1) ! reflecting parallelism of Okada's x-axis with Mansinha & Smylie's x1-axis y_Okada = -x(2) + & ! reflecting anti-parallelism of Okada's y-axis with Mansinha & Smylie's x2-axis; & d_Okada/tanth ! also note OFFSET between Okada and M&S origins along the dip-direction! z_Okada = 0.0D0 ! (On free surface. To change this, use formulas of Okada, 1992 instead of Okada, 1985.) !Distances used in Okada formulas to describe position of observer in a vertical cross-section perpendicular to fault/crack strike: p_Okada = y_Okada * costh + d_Okada * sinth ! distance from lower edge of dislocation patch; component measured parallel to the fault/crack q_Okada = y_Okada * sinth - d_Okada * costh ! component of distance from fault/crack plane (extended beyond dislocation patch, if necessary) !Initialize remote displacement {rate} before sum of 4 terms: ux_Okada = 0.0D0 uy_Okada = 0.0D0 uz_Okada = 0.0D0 DO 200 i = 1, 4 ! sum over 4 terms (from rectangular integration limits in 2-D): !Integration limit sign-changing factors, and integration limits: factor = 1.0D0 ! unless... IF ((i == 2).OR.(i == 3)) factor = -1.0D0 ! replacing the choice made in the previous line epsilon_Okada = x_Okada + L_Okada ! *This is where we choose the double-length, symmetrical source option. IF (i >= 3) epsilon_Okada = x_Okada - L_Okada ! replacing the choice made in the previous line nu_Okada = p_Okada ! unless... IF ((i == 2).OR.(i == 4)) nu_Okada = p_Okada - w_Okada ! replacing the choice made in the previous line !Geometric factors dependent on the above: yTilde_Okada = nu_Okada * costh + q_Okada * sinth dTilde_Okada = nu_Okada * sinth - q_Okada * costh R2_Okada = epsilon_Okada**2 + nu_Okada**2 + q_Okada**2 !R2_Okada = epsilon_Okada**2 + yTilde_Okada**2 + dTilde_Okada**2 ! (alternate formula given by Okada) R_Okada = SQRT(R2_Okada) XCap2_Okada = epsilon_Okada**2 + q_Okada**2 XCap_Okada = SQRT(XCap2_Okada) !Cryptic Okada factors/terms: Note that I2_Okada is not used, because it is associated with strike-slip component of source dislocation. ! Also note that factors must be computed in reverse order (5, 4, 3, 1), due to sequential self-referencing. I5_Okada = mu_over_lambda_plus_mu * (2.0D0/costh) * ATAN( (nu_Okada*(XCap_Okada+q_Okada*costh) + XCap_Okada*(R_Okada+XCap_Okada)*sinth) / & & (epsilon_Okada*(R_Okada+XCap_Okada)*costh) ) I4_Okada = mu_over_lambda_plus_mu * (1.0D0/costh) * ( LOG(R_Okada+dTilde_Okada) - sinth*LOG(R_Okada+nu_Okada) ) I3_Okada = mu_over_lambda_plus_mu * ( (1.0D0/costh)*(yTilde_Okada/(R_Okada+dTilde_Okada)) - LOG(R_Okada+nu_Okada) ) + tanth*I4_Okada I1_Okada = mu_over_lambda_plus_mu * ( (-1.0D0/costh) * (epsilon_Okada/(R_Okada+dTilde_Okada)) ) - tanth*I5_Okada !Final Okada formula for remote displacement due to a finite rectangular opening (mode-1) source dislocation: ux_Okada = ux_Okada + factor * (u3_Okada/Two_Pi) * ( (q_Okada**2/(R_Okada*(R_Okada+nu_Okada))) - & & I3_Okada*(sinth**2) ) uy_Okada = uy_Okada + factor * (u3_Okada/Two_Pi) * ( (-dTilde_Okada*q_Okada/(R_Okada*(R_Okada+epsilon_Okada))) - & & sinth*( (epsilon_Okada*q_Okada/(R_Okada*(R_Okada+nu_Okada))) - ATAN((epsilon_Okada*nu_Okada) / (q_Okada*R_Okada)) ) - & & I1_Okada*(sinth**2) ) uz_Okada = uz_Okada + factor * (u3_Okada/Two_Pi) * ( (yTilde_Okada*q_Okada/(R_Okada*(R_Okada+epsilon_Okada))) + & & costh*( (epsilon_Okada*q_Okada/(R_Okada*(R_Okada+nu_Okada))) - ATAN((epsilon_Okada*nu_Okada) / (q_Okada*R_Okada)) ) - & & I5_Okada*(sinth**2) ) 200 CONTINUE ! sum over 4 terms !Final contribution to remote displacment {rate}, reversing signs !of two terms, per reversal of two axis-directions between M&S versus O: v(1) = v(1) + ux_Okada v(2) = v(2) - uy_Okada v(3) = v(3) - uz_Okada END IF ! non-zero opening component !---------- end Okada [1985] component (tensile/mode-1/opening cracks) --------------- !Copy results to output variables: v4(1) = v(1) v4(2) = v(2) v4(3) = v(3) ! NOTE that this vertical component of remote displacment is positive-downwards, ! according to the Mansinha & Smylie [1971] convention. END SUBROUTINE Aura !------------------------------------------------------------------------------ SUBROUTINE DOnArc (s, theta1, phi1, theta2, phi2, & ! input & theta, phi) ! output ! Computes coordinates (THETA, PHI) = ! (colatitude, longitude) in radians ! for a point which lies on the great circle arc from ! (THETA1, PHI1) to (THETA2, PHI2). ! Input parameter S expresses the fractional distance, ! so input S = 0.0D0 returns (THETA1,PHI1) and ! input S = 1.0D0 returns (THETA2,PHI2). IMPLICIT NONE ! Arguments (see text above): REAL*8, INTENT(IN) :: s, theta1, phi1, theta2, phi2 REAL*8, INTENT(OUT) :: theta, phi REAL*8 :: comple, equat, size REAL*8, DIMENSION(3) :: tvec, uvec, uvec1, uvec2 ! These are all unit vectors (in the unit sphere) ! in a Cartesian coordinate system with ! x = (0E, 0N), y = (90E, 0N), z = North pole. uvec1(1) = DSIN(theta1) * DCOS(phi1) uvec1(2) = DSIN(theta1) * DSIN(phi1) uvec1(3) = DCOS(theta1) uvec2(1) = DSIN(theta2) * DCOS(phi2) uvec2(2) = DSIN(theta2) * DSIN(phi2) uvec2(3) = DCOS(theta2) comple = 1.0D0 - s tvec(1) = comple * uvec1(1) + s * uvec2(1) tvec(2) = comple * uvec1(2) + s * uvec2(2) tvec(3) = comple * uvec1(3) + s * uvec2(3) size = DSQRT(tvec(1)**2 + tvec(2)**2 + tvec(3)**2) uvec(1) = tvec(1) / size uvec(2) = tvec(2) / size uvec(3) = tvec(3) / size equat = DSQRT(uvec(1)**2 + uvec(2)**2) theta = DATAN2(equat, uvec(3)) phi = DATan2F(uvec(2), uvec(1)) END SUBROUTINE DOnArc !------------------------------------------------------------ REAL*8 FUNCTION DATan2F(y, x) ! Works like ATAN2 but corrects for case of (0.0D0, 0.0D0). ! Returns inverse tangent in radians. IMPLICIT NONE REAL*8, INTENT(IN) :: y, x IF (y == 0.0D0) THEN IF (x == 0.0D0) THEN DATAN2F = 0.0D0 ELSE DATAN2F = DATAN2(y, x) END IF ELSE DATAN2F = DATAN2(y, x) END IF END FUNCTION DATan2F !--------------------------------------------------------------- REAL*8 FUNCTION DFltLen (phi1, phi2, radius, theta1, theta2) ! Calculates length of great circle segment between ! point (THETA1,PHI1) and point (THETA2,PHI2), ! in physical length units (radians*planet RADIUS). IMPLICIT NONE REAL*8, INTENT(IN) :: phi1, phi2, radius, theta1, theta2 DOUBLE PRECISION :: ab ab = DSIN(theta1) * DSIN(theta2) * DCOS(phi1) * DCOS(phi2) + & & DSIN(theta1) * DSIN(theta2) * DSIN(phi1) * DSIN(phi2) + & & DCOS(theta1) * DCOS(theta2) ab = DACOS(ab) DFltLen = ab * radius END FUNCTION DFltLen !------------------------------------------------------------------- REAL*8 FUNCTION DChord (angle1, s, angle2) ! Returns an angle obtained by interpolation between ANGLE1 ! and ANGLE2. The interpolation method is not sensitive to ! possible cycle shifts (of 2*n*PI) between ANGLE1 and ANGLE2. ! Unit vectors are constructed for ANGLE1 and ANGLE2, and a ! linear chord is drawn between their tips. ! Double-precision S is the internal coordinate along the chord; ! it is dimensionless, with value 0.0 at ANGLE1 and 1.0 at ! ANGLE2. (The user may input S values outside this range ! to get results outside the (smaller) angle between ANGLE1 and ! ANGLE2, if desired.) The angle returned is that from the ! origin to this chord point. ! This algorithm should work equally well for angles measured ! either clockwise or counterclockwise from any reference, as ! long as the usage is consistent. ! Both the input angles and the result "Chord" are in radians. DOUBLE PRECISION, INTENT(IN) :: s REAL*8, INTENT(IN) :: angle1, angle2 REAL*8, DIMENSION(2) :: uvec1, uvec2, vecs uvec1(1) = DCOS(angle1) uvec1(2) = DSIN(angle1) uvec2(1) = DCOS(angle2) uvec2(2) = DSIN(angle2) vecs(1) = (1.0D0 - s) * uvec1(1) + s * uvec2(1) vecs(2) = (1.0D0 - s) * uvec1(2) + s * uvec2(2) DChord = DATAN2F(vecs(2), vecs(1)) END FUNCTION DChord DOUBLE PRECISION FUNCTION DAtanPV (y, x) ! Returns principal value of inverse tangent of (y,x) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: x, y DOUBLE PRECISION :: xx, yy IF (x >= 0.0D0) THEN xx = x yy = y ELSE xx = -x yy = -y END IF IF((yy == 0.0D0).OR.(xx == 0.0D0)) THEN IF(yy == 0.0D0) DAtanPV = 0.0D0 IF(xx == 0.0D0) DAtanPV = 1.57079632679490D0 ELSE DAtanPV = DATAN2(yy, xx) END IF END FUNCTION DAtanpv END MODULE DDislocation