MODULE Dislocation ! by Peter Bird, UCLA, gradually developed 1980-2002; ! last revision on 2002.08.13. ! Compute benchmark velocities due to elastic dislocations ! from earthquakes in the brittle part of each fault, at the rate ! predicted by a finite-element model: ! Usage from a F-E program like Shells/OrbScore, with fault elements: ! make a single call for the entire F-E grid: ! ! CALL Coseis (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 Change 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 Change (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 !================================================================== CONTAINS SUBROUTINE Coseis (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 Change ! 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, INTENT(IN) :: radius, wedge REAL, DIMENSION(2,mxfel), INTENT(IN) :: farg, fdip, ztranf REAL, DIMENSION(mxnode), INTENT(IN) :: xnode, ynode, zmnode REAL, DIMENSION(numgeo), INTENT(IN) :: geothe, geophi REAL, DIMENSION(numgeo), INTENT(OUT) :: utheta, uphi DOUBLE PRECISION :: lf INTEGER :: i, ilayer, j, m, n1, n2, n3, n4 REAL :: 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, 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.0 uphi(j) = 0.0 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 = (REAL(m) - 0.5) / REAL(nsegme) CALL Onarc(smid, x1, y1, x2, y2, & ! input & xmid, ymid) ! output s = REAL(m) / REAL(nsegme) IF(m == nsegme) THEN xend = x2 yend = y2 ELSE CALL Onarc(s, x1, y1, x2, y2, & ! input & xend, yend) ! output END IF al = fltlen(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 = chord(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 = COS(argume) unity = SIN(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.570796 - dip) > wedge) THEN slip(3) = slip(2) * TAN(dip) ELSE slip(3) = 0.0 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.0) THEN IF (ilayer == 1) THEN ztop = 0.0 ELSE ztop = moho END IF zbot = ztop + height ! Compute effects of patch at all benchmarks: DO 70 j = 1, numgeo ! ------------------------------------------- CALL Change(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 Coseis !----------------------------------------------------------------------- SUBROUTINE Change (argume, & & btheta, bphi, & & dipf, lf, & & ftheta, fphi, & & radius, & & slip, & & wedge, & & ztop, zbot, & ! inputs & duthet, duphi) ! output ! This routine is a "driver" or "wrapper" for the ! Mansinha & Smylie routine Halo (vertical dislocation patch) ! or Aura (dipping dislocation patch) ! 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, ! 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). ! SLIP must be already rotated into fault 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. ! 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: If distance from "F" point to "B" point exceeds ! (50 * ZBOT), result is approximated as zero and ! Halo/Aura are never called. IMPLICIT NONE ! ------------------------------------------------------- ! Arguments (see text above): REAL, INTENT(IN) :: argume, btheta, bphi, dipf, & & ftheta, fphi, radius, wedge, zbot, ztop REAL, DIMENSION(3), INTENT(IN) :: slip REAL, INTENT(OUT) :: duthet, duphi ! ------------------------------------------------------- ! 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 :: 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) = SIN(dble(ftheta)) * COS(dble(fphi)) fuvec(2) = SIN(dble(ftheta)) * SIN(dble(fphi)) fuvec(3) = COS(dble(ftheta)) buvec(1) = SIN(dble(btheta)) * COS(dble(bphi)) buvec(2) = SIN(dble(btheta)) * SIN(dble(bphi)) buvec(3) = COS(dble(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 = SQRT(tvec(1)**2 + tvec(2)**2 + tvec(3)**2) far = radius * ATAN2(crossp, dotpro) IF (far > (50.D0 * zbot)) THEN duthet = 0.0 duphi = 0.0 ELSE dvec(1) = buvec(1) - fuvec(1) dvec(2) = buvec(2) - fuvec(2) dvec(3) = buvec(3) - fuvec(3) utheta(1) = COS(dble(ftheta)) * COS(dble(fphi)) utheta(2) = COS(dble(ftheta)) * SIN(dble(fphi)) utheta(3) = -SIN(dble(ftheta)) dot1 = dvec(1) * utheta(1) + dvec(2) * utheta(2) + dvec(3) * utheta(3) uphi(1) = -SIN(dble(fphi)) uphi(2) = COS(dble(fphi)) uphi(3) = 0.0D0 dot2 = dvec(1) * uphi(1) + dvec(2) * uphi(2) + dvec(3) * uphi(3) alpha = ATAN2(dot2, dot1) z(1) = far * COS(alpha - argume) z(2) = -far * SIN(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.0 vert = ABS(dipf - 1.5708) <= wedge IF (vert) THEN tazimf = argume theta = dipf ! Z values computed above are still good ucap(1) = slip(1) ucap(2) = 0.0 ucap(3) = 0.0 dtop = ztop dbot = zbot ! - - - - - - - - - - - - - - - - - - - - - - - CALL Halo (ucap, z, lf, dbot, dtop, & ! inputs & u) ! output ! - - - - - - - - - - - - - - - - - - - - - - - ELSE ! non-vertical, dipping fault: IF (dipf <= 1.5708) 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.141593 theta = 3.141593 - 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 / SIN(theta) dbot = zbot / SIN(theta) ! Kludge to insure that DTOP is not zero ! because AURA gives erroneous results in that case: dtop = MAX(dtop, 0.001D0 * dbot) ! - - - - - - - - - - - - - - - - - - - - - - - CALL Aura (ucap, theta, x, lf, dbot, dtop, & ! inputs & u) ! output ! - - - - - - - - - - - - - - - - - - - - - - - END IF sinaz = SIN(tazimf) cosaz = COS(tazimf) duthet = + cosaz * u(1) + sinaz * u(2) duphi = + sinaz * u(1) - cosaz * u(2) END IF END SUBROUTINE Change !----------------------------------------------------------------- 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 on a ! rectangular dislocation in a vertical plane. ! Equations based on Mansinha and Smylie (1967), ! J. Geophys. Res., v. 72, p. 4731-4743. 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 :: factor, q, r, term1, term2, zeta1, zeta3 ! The coordinate system is Cartesian and right-handed, ! with z3 down, z1 parallel to the fault 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 slip of the ! positive-z2 side relative to the other side; UCAP(2) ! should therefore be zerp (but it is not used here). ! U contains the three components of the displacement ! at the observation point Z. DO 10 i = 1, 3 10 u(i) = 0. DO 100 i = 1, 4 factor = 1. IF(i == 2.OR.i == 3) factor = -1. zeta1 = l IF(i == 3.OR.i == 4) zeta1 = -l zeta3 = dbot IF(i == 2.OR.i == 4) zeta3 = dtop r = SQRT((z(1) - zeta1)**2 + z(2)**2 + (z(3) - zeta3)**2) q = SQRT((z(1) - zeta1)**2 + z(2)**2 + (z(3) + zeta3)**2) term1 = z(2) * (z(1) - zeta1) * (2. / (r * (r + z(3) - zeta3)) & & - (5. * q + 8. * zeta3) / (2. * q * (q + z(3) + zeta3)**2) & & + 4. * z(3) * zeta3 * (2. * q + z(3) + zeta3) / & & (q**3 * (q + z(3) + zeta3)**2)) & & + 3. * atanpv((z(1) - zeta1) * (z(3) - zeta3), (r * z(2))) & & - 3. * atanpv((z(1) - zeta1) * (z(3) + zeta3), (q * z(2))) term2 = -log(r + z(3) - zeta3) & & + 0.5 * log(q + z(3) + zeta3) & & - (5. * z(3) - 3. * zeta3) / (2. * (q + z(3) + zeta3)) & & - 4. * z(3) * zeta3 / (q * (q + z(3) + zeta3)) & & + z(2)**2 * (2. / (r * (r + z(3) - zeta3)) & & - (5. * q + 8. * zeta3) / (2. * q * (q + z(3) + zeta3)**2) & & + 4. * z(3) * zeta3 * (2. * q + z(3) + zeta3) / & & (q**3 * (q + z(3) + zeta3)**2)) ! TERM3=Z(2)*(2./R-2./Q+4.*Z(3)*ZETA3/Q**3 ! 1 +3./(Q+Z(3)+ZETA3)+2.*(Z(3)+3.*ZETA3)/ ! 2 (Q*(Q+Z(3)+ZETA3))) ! TERM4=Z(2)*(2./R+4./Q-4.*Z(3)*ZETA3/Q**3 ! 1 -6.*Z(3)/(Q*(Q+Z(3)+ZETA3))) ! TERM5= -LOG(R+Z(1)-ZETA1)+LOG(Q+Z(1)-ZETA1) ! 1 +(6.*Z(3)**2+10.*Z(3)*ZETA3)/(Q*(Q+Z(1)-ZETA1)) ! 2 +(6.*Z(3)*(Z(1)-ZETA1))/(Q*(Q+Z(3)+ZETA3)) ! 3 +Z(2)**2*(2./(R*(R+Z(1)-ZETA1)) ! 4 +4./(Q*(Q+Z(1)-ZETA1)) ! 5 -4.*Z(3)*ZETA3*(2.*Q+Z(1)-ZETA1)/(Q**3*(Q+Z(1)-ZETA1)**2)) ! TERM6=Z(2)*(2.*(Z(3)-ZETA3)/(R*(R+Z(1)-ZETA1)) ! 1 -2.*(Z(3)+2.*ZETA3)/(Q*(Q+Z(1)-ZETA1)) ! 2 -4.*Z(3)*ZETA3*(Z(3)+ZETA3)*(2.*Q+Z(1)-ZETA1)/ ! 3 (Q**3*(Q+Z(1)-ZETA1)**2)) ! 4 +3.*ATANPV((Z(1)-ZETA1)*(Z(3)-ZETA3),(R*Z(2))) ! 5 -3.*ATANPV((Z(1)-ZETA1)*(Z(3)+ZETA3),(Q*Z(2))) u(1) = u(1) + factor * term1 * ucap(1) / 37.7 u(2) = u(2) + factor * term2 * ucap(1) / 37.7 ! U(1)=U(1)+FACTOR*TERM4*UCAP(3)/37.7 ! U(2)=U(2)+FACTOR*TERM5*UCAP(3)/37.7 ! U(3)=U(3)+FACTOR*(TERM3*UCAP(1)+TERM6*UCAP(3))/37.7 100 CONTINUE 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 based on Mansinha and Smylie (1971), ! Bull.Seism.Soc.Amer., v. 61, p. 1433-1440. 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 INTEGER :: i DOUBLE PRECISION :: costh, dbot, deld, dtop, factor, & & h, k, l, q, q2, q3, qsquar, & & r, r2, r3, rsquar, & & secth, sinth, tanth, & & term1, term2, term4, term5, theta, & & ucap, ucap1, & & xi, xi1, xi2, xi3 DOUBLE PRECISION :: term5_mute ! added 2015.03.26 DOUBLE PRECISION, DIMENSION(3) :: u(3), v(3), x(3) DO 1 i = 1, 3 u(i) = u4(i) x(i) = x4(i) 1 CONTINUE theta = theta4 l = l4 dbot = dbot4 dtop = dtop4 ! U(3) is the dislocation vector (i.e. the total slip vector) ! of the hanging wall relative to the footwall, ! in a Cartesian coordinate system where x1 is parallel to the ! fault strike, x2 is also horizontal and pointing to hanging ! wall, 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. ! V(3) is the answer: the 3-component displacement at the ! observation point in same units as U and in x1,x2,x3 coordinates. sinth = SIN(theta) costh = COS(theta) tanth = sinth / costh ! NOTE: FAULT MAY NOT BE VERTICAL IN THIS ROUTINE. secth = 1. / costh 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 10 i = 1, 3 10 v(i) = 0.D0 DO 100 i = 1, 4 factor = 1.00D0 IF (i == 2.OR.i == 3) factor = -1.00D0 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 = SQRT(qsquar - (q3 + xi)**2) h = SQRT(qsquar - (x(1) - xi1)**2) ucap1 = u(1) ucap = SQRT(u(2)**2 + u(3)**2) IF(u(2) < 0.)ucap = -ucap q = SQRT(qsquar) r = SQRT(rsquar) q = MAX(q, 1.00001D0 * ABS(x(1) - xi1)) r = MAX(r, 1.00001D0 * ABS(x(1) - xi1)) term1 = (x(1) - xi1) * (2. * r2 / (r * (r + r3 - xi)) & & - (4. * q2 - 2. * x(3) * costh) / (q * (q + q3 + xi)) & & - 3. * tanth / (q + x(3) + xi3) & & + 4. * q2 * x(3) * sinth / q**3 & & - 4. * q2 * q3 * x(3) * sinth * (2. * q + q3 + xi) / & & (q**3 * (q + q3 + xi)**2)) & & - 6. * tanth**2 * atanpv(((k - q2 * costh) * (q - k) + & & (q3 + xi) * k * sinth), ((x(1) - xi1) * (q3 + xi) * costh)) & & + 3. * atanpv((x(1) - xi1) * (r3 - xi), (r2 * r)) & & - 3. * atanpv((x(1) - xi1) * (q3 + xi), (q2 * q)) term2 = sinth * (3. * tanth * secth * log(q + x(3) + xi3) - & & log(r + r3 - xi) - (1. + 3. * tanth**2) * log & & (q + q3 + xi)) + 2. * r2**2 * sinth / (r * (r + r3 - xi)) & & + 2. * r2 * costh / r - 2. * sinth * & & (2. * x(3) * (q2 * costh - q3 * sinth) + q2 * (q2 + x(2) * sinth)) & & / (q * (q + q3 + xi)) - 3. * tanth * (x(2) - xi2) / & & (q + x(3) + xi3) + 2. * (q2 * costh - q3 * sinth - & & x(3) * sinth**2) / q + 4. * q2 * x(3) * sinth * & & ((x(2) - xi2) + q3 * costh) / q**3 - 4. * q2**2 * q3 * x(3) & & * sinth**2 * (2. * q + q3 + xi) / (q**3 * (q + q3 + xi)**2) ! TERM3=COSTH*(LOG(R+R3-XI)+(1.+3.*TANTH**2)*LOG( ! 1 Q+Q3+XI)-3.*TANTH*SECTH*LOG(Q+X(3)+XI3)) ! 2 +2.*R2*SINTH/R+2.*SINTH*(Q2+X(2)*SINTH)/Q ! 3 -2.*R2**2*COSTH/(R*(R+R3-XI))+ ! 4 (4.*Q2*X(3)*SINTH**2-2.*(Q2+X(2)*SINTH)*(X(3)+Q3*SINTH))/ ! 5 (Q*(Q+Q3+XI))+ ! 6 4.*Q2*X(3)*SINTH*(X(3)+XI3-Q3*SINTH)/Q**3 ! 7 -4.*Q2**2*Q3*X(3)*COSTH*SINTH*(2.*Q+Q3+XI)/ ! 8 (Q**3*(Q+Q3+XI)**2) term4 = (x(2) - xi2) * sinth * (2. / r + 4. / q - 4. * xi3 * x(3) / q**3 & & - 3. / (q + x(3) + xi3)) - costh * (3. * log(q + x(3) + xi3) & & + 2. * (x(3) - xi3) / r + 4. * (x(3) - xi3) / q + 4. * xi3 * x(3) * (x(3) + xi3) & & / q**3) + 3. * (log(q + x(3) + xi3) - sinth * log(q + & & q3 + xi)) / costh + 6. * x(3) * (costh / q - q2 * sinth / & & (q * (q + q3 + xi))) term5 = sinth * (-log(r + x(1) - xi1) + log(q + x(1) - xi1) & & + 4. * xi3 * x(3) / (q * (q + x(1) - xi1)) + 3. * (x(1) - xi1) / & & (q + x(3) + xi3) + (x(2) - xi2)**2 * (2. / (r * (r + x(1) - xi1)) & & + 4. / (q * (q + x(1) - xi1)) - 4. * xi3 * x(3) * (2. * q + x(1) - xi1) / & & (q**3 * (q + x(1) - xi1)**2))) & & - costh * ((x(2) - xi2) * (2. * (x(3) - xi3) / (r * (r + x(1) - xi1)) & & + 4. * (x(3) - xi3) / (q * (q + x(1) - xi1)) + 4. * xi3 * x(3) * & & (x(3) + xi3) * (2. * q + x(1) - xi1) / (q**3 * (q + x(1) - xi1)**2)) & & + 6. * atanpv((x(1) - xi1) * (x(2) - xi2), ((h + x(3) + xi3) * (q + h))) & & - 3. * atanpv((x(1) - xi1) * (r3 - xi), (r2 * r)) + & & 6. * atanpv((x(1) - xi1) * (q3 + xi), (q2 * q))) + & & 6. * (secth * atanpv(((k - q2 * costh) * (q - k) + (q3 + xi) * k * sinth) & & , ((x(1) - xi1) * (q3 + xi) * costh)) & & + x(3) * (((sinth**2 - costh**2) * (q3 + xi) + 2. * q2 * & & costh * sinth) / (q * (q + x(1) - xi1)) + (x(1) - xi1) * sinth**2 & & / (q * (q + q3 + xi)))) term5_mute = 1.D0 / ( 1.D0 + 20.D0 * (ABS(x(1))/l) * EXP(-(10.D0 * x(2)/MAX(ABS(x(1)), dtop))**2) ) term5 = term5 * term5_mute ! added 2015.03.26 ! TERM6=SINTH*((X(2)-XI2)*(2.*(X(3)-XI3)/(R*(R+X(1)-XI1)) ! 1 +4.*(X(3)-XI3)/(Q*(Q+X(1)-XI1))-4.*XI3*X(3)* ! 2 (X(3)+XI3)*(2.*Q+X(1)-XI1)/(Q**3*(Q+X(1)-XI1)**2)) ! 3 -6.*ATANPV((X(1)-XI1)*(X(2)-XI2),((H+X(3)+XI3)*(Q+H))) ! 4 +3.*ATANPV((X(1)-XI1)*(R3-XI),(R2*R))- ! 5 6.*ATANPV((X(1)-XI1)*(Q3+XI),(Q2*Q))) ! 6 +COSTH*(LOG(R+X(1)-XI1)-LOG(Q+X(1)-XI1)- ! 7 2.*(X(3)-XI3)**2/(R*(R+X(1)-XI1)) ! 8 -4.*((X(3)+XI3)**2-XI3*X(3))/(Q*(Q+X(1)-XI1)) ! 9 -4.*XI3*X(3)*(X(3)+XI3)**2*(2.*Q+X(1)-XI1)/ ! A (Q**3*(Q+X(1)-XI1)**2)) ! B +6.*X(3)*(COSTH*SINTH*(2.*(Q3+XI)/(Q*(Q+X(1)-XI1)) ! C +(X(1)-XI1)/(Q*(Q+Q3+XI))) ! D -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 v4(1) = v(1) v4(2) = v(2) ! V4(3)=V(3) END SUBROUTINE Aura !------------------------------------------------------------------------------ SUBROUTINE OnArc (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.0 returns (THETA1,PHI1) and ! input S = 1.0 returns (THETA2,PHI2). IMPLICIT NONE ! Arguments (see text above): REAL, INTENT(IN) :: s, theta1, phi1, theta2, phi2 REAL, INTENT(OUT) :: theta, phi REAL :: comple, equat, size REAL, 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) = SIN(theta1) * COS(phi1) uvec1(2) = SIN(theta1) * SIN(phi1) uvec1(3) = COS(theta1) uvec2(1) = SIN(theta2) * COS(phi2) uvec2(2) = SIN(theta2) * SIN(phi2) uvec2(3) = COS(theta2) comple = 1.0 - 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 = SQRT(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 = SQRT(uvec(1)**2 + uvec(2)**2) theta = ATAN2(equat, uvec(3)) phi = atan2f(uvec(2), uvec(1)) END SUBROUTINE OnArc !------------------------------------------------------------ REAL FUNCTION ATAN2F(y, x) ! Works like ATAN2 but corrects for case of (0, 0). ! Returns inverse tangent in radians. IMPLICIT NONE REAL, INTENT(IN) :: y, x IF (y == 0.) THEN IF (x == 0.) THEN ATAN2F = 0. ELSE ATAN2F = ATAN2(y, x) END IF ELSE ATAN2F = ATAN2(y, x) END IF END FUNCTION ATAN2F !--------------------------------------------------------------- REAL FUNCTION FltLen (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, INTENT(IN) :: phi1, phi2, radius, theta1, theta2 DOUBLE PRECISION :: ab ab = SIN(theta1) * SIN(theta2) * COS(phi1) * COS(phi2) + & & SIN(theta1) * SIN(theta2) * SIN(phi1) * SIN(phi2) + & & COS(theta1) * COS(theta2) ab = ACOS(ab) FltLen = ab * radius END FUNCTION FltLen !------------------------------------------------------------------- REAL FUNCTION Chord (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, INTENT(IN) :: angle1, angle2 REAL, DIMENSION(2) :: uvec1, uvec2, vecs uvec1(1) = COS(angle1) uvec1(2) = SIN(angle1) uvec2(1) = COS(angle2) uvec2(2) = SIN(angle2) vecs(1) = (1.0D0 - s) * uvec1(1) + s * uvec2(1) vecs(2) = (1.0D0 - s) * uvec1(2) + s * uvec2(2) Chord = ATAN2F(vecs(2), vecs(1)) END FUNCTION Chord DOUBLE PRECISION FUNCTION Atanpv (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.D0) 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) Atanpv = 0.0D0 IF(xx == 0.0D0) Atanpv = 1.570796327D0 ELSE Atanpv = ATAN2(yy, xx) END IF END FUNCTION Atanpv END MODULE Dislocation