anchor.diff:  C----------------------------------------------------------------------C
   C======================================================================C
   
  -      SUBROUTINE dipterms(dipin)
  -C----
  -C Pre-calculate dip terms (sin(dip), cos(dip), etc.)
  -C
  -C Note: double precision math can tolerate fault dips close to 90, but
  -C not equal to 90 (ends up with divide by zero problems):
  -C     
  -C     if (dip.gt.89.9999d0) dip = 89.9999d0
  -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
  -
  +      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
   
  
positive.code: 
        SUBROUTINE o92rect(ux,uy,uz,x,y,stdp,evdp,dipin,rakin,wid,len,
       1                   slip,vp,vs,dens)
  C----
  C Static elastic displacement at a location (internal or surface) due
  C to a finite rectangular source shear dislocation in a uniform, 
  
negative.code:       mu = dens*vs*vs
        lamda = dens*vp*vp - 2.0d0*mu
        a   = (lamda+mu)/(lamda+2.0d0*mu)
        CA1 = (1.0d0-a)*0.5d0
        CA2 = a*0.5d0
        CB  = (1.0d0-a)/a
        CC  = 1.0d0-a
  
  