anchor.diff:  C----
   C Source and halfspace constants
   C----
  -      call dipterms(dipin)
  +      call DIPVAR(dipin)
         call hafspc(vp,vs,dens)
         call moment0(Mss,Mds,slip,rakin,area)
   C----
  
positive.code: C----------------------------------------------------------------------C
  C Okada, Y., 1992, Internal deformation due to shear and tensile       C
  C   faults in a half-space. Bulletin of the Seismological Society of   C
  C   America 82, pp. 1018-1040.                                         C
  C                                                                      C
  C MODIFICATIONS:                                                       C
  C   Fall 2012: Original file created                                   C
  C   Summer 2013: Instead of 4 large subroutines, compartmentalize       C
  C                calculations into several smaller subroutines.        C
  C                Pass variables to and from subroutines via common     C
  C                blocks to make variable list shorter.                 C
  C   2013-10-16: Began commenting code thoroughly, renaming variables   C
  C               and subroutines to be more consistent throughout code. C
  C----------------------------------------------------------------------C
  
        SUBROUTINE o92pt(ux,uy,uz,x,y,stdp,evdp,dipin,rakin,area,slip,
       1                 vp,vs,dens)
  C----
  C Static elastic displacement at a location (internal or surface) due to
  C a point source shear dislocation in an isotropic halfspace.
  C
  C INPUTS (ALL IN SI UNITS):
  C   - Geometric Parameters (relative to source location):
  C         x (along-strike), y (horizontal up-dip), station depth,
  C         event depth
  C   - Source parameters:
  C         fault dip (deg), fault rake (deg), fault area, slip
  C   - Physical parameters of half-space:
  C         vp, vs, density
  C OUTPUT: x, y, z displacements at station (in x-y coordinates, not N-E)
  C----
        IMPLICIT none
        REAL*8 pi,d2r,eps
        PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2,eps=1.0d-4)
        REAL*8 x,y,stdp,evdp,z,c,d
        REAL*8 dipin,rakin,area,slip,Mss,Mds
        REAL*8 vp,vs,dens
  
negative.code: C----------------------------------------------------------------------C
  C======================================================================C
  
        SUBROUTINE DIPVAR(DIPIN)
  C----
  C PRE-CALCULATE DIP VARIABLES.
  C IF (DIP.GT.89.9999D0) DIP = 89.9999D0 (TO AVOID DIVIDE BY ZERO ERROR)
  C----
        IMPLICIT NONE
        REAL*8 PI,D2R
        PARAMETER (PI=4.0D0*DATAN(1.0D0),D2R=PI/1.8D2)
        REAL*8 DIP,DIPIN
        REAL*8 SD,CD,S2D,C2D,CDCD,SDSD,CDSD
        COMMON /SOURCE/ SD,CD,S2D,C2D,CDCD,SDSD,CDSD
        IF (DIPIN.GT.89.9999D0) THEN
            DIP = 89.9999D0*D2R
        ELSE
            DIP = DIPIN*D2R
        ENDIF
        SD = DSIN(DIP)
        CD = DCOS(DIP)
        S2D = DSIN(2.0D0*DIP)
        C2D = DCOS(2.0D0*DIP)
        CDCD = CD*CD
        SDSD = SD*SD
        CDSD = SD*CD
        RETURN
        END
  
  