anchor.diff:            logRe = dlog(R+eta)
         endif
         Y0  = Y11 - ksi*ksi*Y32
  -      
  +
         Z32 = sd/R3 - h*Y32
         Z53 = 3.0d0*sd/R5 - h*Y53
         Z0  = Z32 - ksi*ksi*Z53
  
positive.code: 
  C-------------------
  
        SUBROUTINE disp1stn(f,ksi,eta,q,z)
  C----
  C Components of displacement due to finite source for calculating
  
negative.code: C======================================================================C
  C Subroutines to compute deformation from point and finite sources in
  C an elastic half-space.
  C
  C Reference:
  C Okada, Y., 1992, Internal deformation due to shear and tensile
  C   faults in a half-space. Bulletin of the Seismological Society of
  C   America 82, pp. 1018-1040.
  C
  C Notes:
  C - All variables are assumed to be in SI units
  C - Angles input into subroutines are in degrees
  C - Coordinate system is relative to the fault orientation:
  C   - X points in the along-strike direction
  C   - Y points in the horizontal up-dip direction
  C   - Depths are defined positive DOWN
  C======================================================================C
  
  C      SUBROUTINE PTDISP(disp,x,y,stdp,evdp,dipin,rakin,area,slip,
  C     1                 vp,vs,dens)
        SUBROUTINE o92pt(ux,uy,uz,x,y,stdp,evdp,dipin,rakin,area,slip,
       1                 vp,vs,dens)
  C----
  C Compute displacement vector at a point in anywhere in a half-space
  C due to a point source
  C----
  
        IMPLICIT none
        REAL*8 x,y,stdp,evdp,z,c,d             ! Coordinates
        REAL*8 dipin,rakin,area,slip,Mss,Mds   ! Fault parameters
        REAL*8 vp,vs,dens                      ! Half-space parameters
        INTEGER i,j
        REAL*8 u(6,3),f(6,3),ux,uy,uz
  C      REAL*8 disp(3)
        REAL*8 zro
        DATA zro/0.0d0/
  
  