initial commit

5 years ago · 320dc65ae3
56 changed files with 9268 additions and 0 deletions
--- a/Makefile.am
+++ b/Makefile.am
@ -0,0 +1,2 @@
+AUTOMAKE_OPTIONS = foreign
+SUBDIRS = src/psgrn src/pscmp
--- a/README.md
+++ b/README.md
@ -0,0 +1,29 @@
+# PSGRN and PSCMP (packaged as fomosto backend)
+
+Code to calculate synthetic stress/strain/tilt/gravitational fields on a
+layered viscoelastic halfspace.
+
+PSGRN and PSCMP have been written by Rongjiang Wang.
+
+Packaging has been done by Hannes Vasyura-Bathke.
+
+## References
+
+- Wang, R., F. Lorenzo-Martín and F. Roth (2003), Computation of deformation
+  induced by earthquakes in a multi-layered elastic crust - FORTRAN programs
+  EDGRN/EDCMP, Computer and Geosciences, 29(2), 195-207.
+- Wang, R., F. Lorenzo-Martin and F. Roth (2006), PSGRN/PSCMP - a new code for
+  calculating co- and post-seismic deformation, geoid and gravity changes
+  based on the viscoelastic-gravitational dislocation theory, Computers and
+  Geosciences, 32, 527-541. DOI:10.1016/j.cageo.2005.08.006.
+- Wang, R. (2005), The dislocation theory: a consistent way for including the
+  gravity effect in (visco)elastic plane-earth models, Geophysical Journal
+  International, 161, 191-196.
+
+# Compile and install PSGRN and PSCMP
+```
+autoreconf -i   # only if 'configure' script is missing
+F77=gfortran FFLAGS=-mcmodel=medium ./configure
+make
+sudo make install
+```
--- a/configure.ac
+++ b/configure.ac
@ -0,0 +1,13 @@
+#                                               -*- Autoconf -*-
+# Process this file with autoconf to produce a configure script.
+
+AC_PREREQ([2.68])
+AC_INIT([fomosto_psgrn], [2008a])
+AM_INIT_AUTOMAKE([-Wall -Werror foreign])
+AC_PROG_F77
+AC_CONFIG_FILES([
+                 Makefile
+                 src/psgrn/Makefile
+                 src/pscmp/Makefile
+                 ])
+AC_OUTPUT
--- a/examples/pscmp08-wenchuan.inp
+++ b/examples/pscmp08-wenchuan.inp
@ -0,0 +1,424 @@
+#===============================================================================
+# This is input file of FORTRAN77 program "pscmp08" for modeling post-seismic
+# deformation induced by earthquakes in multi-layered viscoelastic media using
+# the Green's function approach. The earthquke source is represented by an
+# arbitrary number of rectangular dislocation planes. For more details, please
+# read the accompanying READ.ME file.
+#
+# written by Rongjiang Wang
+# GeoForschungsZentrum Potsdam
+# e-mail: wang@gfz-potsdam.de
+# phone +49 331 2881209
+# fax +49 331 2881204
+#
+# Last modified: Potsdam, July, 2008
+#
+#################################################################
+##                                                             ##
+## Green's functions should have been prepared with the        ##
+## program "psgrn08" before the program "pscmp08" is started.  ##
+##                                                             ##
+## For local Cartesian coordinate system, the Aki's convention ##
+## is used, that is, x is northward, y is eastward, and z is   ##
+## downward.                                                   ##
+##                                                             ##
+## If not specified otherwise, SI Unit System is used overall! ##
+##                                                             ##
+#################################################################
+#===============================================================================
+# OBSERVATION ARRAY
+# =================
+# 1. selection for irregular observation positions (= 0) or a 1D observation
+#    profile (= 1) or a rectangular 2D observation array (= 2): iposrec
+#
+#    IF (iposrec = 0 for irregular observation positions) THEN
+#
+# 2. number of positions: nrec
+#
+# 3. coordinates of the observations: (lat(i),lon(i)), i=1,nrec
+#
+#    ELSE IF (iposrec = 1 for regular 1D observation array) THEN
+#
+# 2. number of position samples of the profile: nrec
+#
+# 3. the start and end positions: (lat1,lon1), (lat2,lon2)
+#
+#    ELSE IF (iposrec = 2 for rectanglular 2D observation array) THEN
+#
+# 2. number of x samples, start and end values: nxrec, xrec1, xrec2
+#
+# 3. number of y samples, start and end values: nyrec, yrec1, yrec2
+#
+#    sequence of the positions in output data: lat(1),lon(1); ...; lat(nx),lon(1);
+#    lat(1),lon(2); ...; lat(nx),lon(2); ...; lat(1),lon(ny); ...; lat(nx),lon(ny).
+#
+#    Note that the total number of observation positions (nrec or nxrec*nyrec)
+#    should be <= NRECMAX (see pecglob.h)!
+#===============================================================================
+  0
+  180
+( 31.8010, 104.4430) ( 32.1820, 104.8720) ( 31.0600, 103.6910) ( 31.4860, 104.2250) ( 32.5700, 105.2200)
+( 31.3530, 104.1860) ( 31.7050, 104.4430) ( 31.0080, 103.1450) ( 32.3600, 104.8100) ( 30.9680, 103.7400)
+( 30.8800, 103.6200) ( 32.4050, 104.5710) ( 31.1560, 104.4400) ( 31.4860, 104.7600) ( 31.4860, 104.7810)
+( 32.0600, 103.5800) ( 31.8500, 102.6700) ( 32.0750, 103.1650) ( 30.7320, 104.0770) ( 32.0200, 105.8300)
+( 32.0200, 105.8300) ( 32.3610, 103.7310) ( 31.7050, 102.3060) ( 30.6300, 104.0800) ( 32.4480, 105.8300)
+( 32.4480, 105.8300) ( 30.6300, 103.6300) ( 31.4660, 102.0950) ( 31.0300, 102.4000) ( 32.5900, 103.6130)
+( 31.1000, 105.1000) ( 31.8700, 105.9800) ( 31.7700, 101.6150) ( 32.8510, 103.5700) ( 30.9600, 101.8700)
+( 32.7850, 102.5000) ( 33.0000, 104.6250) ( 32.9300, 103.4350) ( 30.4050, 104.5300) ( 30.3750, 104.5360)
+( 32.9010, 101.7060) ( 32.8000, 105.7800) ( 31.1430, 100.9300) ( 33.4300, 105.0100) ( 30.9500, 101.1630)
+( 25.3410, 100.4960) ( 25.4810, 100.5480) ( 25.6080, 103.2410) ( 25.6410, 101.9010) ( 25.7310, 101.3200)
+( 26.0010, 102.5310) ( 26.0500, 101.6810) ( 26.1050, 103.1650) ( 26.5030, 101.7480) ( 26.6200, 102.6100)
+( 26.6900, 102.2630) ( 26.6900, 101.8550) ( 26.9310, 102.9060) ( 27.0480, 101.9580) ( 27.1380, 100.9330)
+( 27.4200, 101.5130) ( 27.4530, 102.1880) ( 27.5400, 101.7100) ( 27.6560, 101.2380) ( 27.7480, 100.6530)
+( 27.8750, 102.2310) ( 28.3000, 102.4360) ( 28.5150, 102.1250) ( 28.9630, 101.5180) ( 29.2280, 103.2610)
+( 29.6000, 103.8000) ( 29.6880, 102.0800) ( 29.7900, 102.8160) ( 29.8460, 101.5580) ( 29.8480, 102.2900)
+( 30.0410, 103.8450) ( 30.0730, 101.7880) ( 30.0750, 101.4850) ( 30.1060, 101.0230) ( 30.2510, 102.8400)
+( 30.4150, 103.4100) ( 30.4950, 101.4960) ( 30.5000, 105.7800) ( 30.8000, 106.2000) ( 34.4030, 104.0730)
+( 32.8000, 106.1800) ( 32.8500, 107.1700) ( 33.1000, 106.3300) ( 33.1160, 106.6800) ( 33.1900, 106.5800)
+( 33.2280, 104.2250) ( 33.2760, 103.8880) ( 33.3400, 105.8050) ( 33.3400, 106.1550) ( 33.4000, 105.6280)
+( 33.4230, 104.8230) ( 33.5710, 102.9910) ( 33.6960, 105.5950) ( 33.6960, 105.5950) ( 33.7800, 105.2850)
+( 33.7860, 104.4010) ( 33.8910, 105.8150) ( 33.9150, 106.5080) ( 33.9360, 103.7260) ( 34.0000, 104.4200)
+( 34.0200, 105.3000) ( 34.0460, 104.3830) ( 34.1080, 105.3060) ( 34.1080, 103.1460) ( 34.2510, 105.8110)
+( 34.3600, 104.5000) ( 34.3600, 104.8300) ( 34.4660, 104.9150) ( 34.5000, 105.8600) ( 34.5510, 108.9130)
+( 34.5930, 105.6960) ( 34.7130, 104.9400) ( 34.7480, 106.1580) ( 34.8500, 104.4800) ( 34.8710, 105.6550)
+( 35.1410, 105.3780) ( 35.1730, 106.0110) ( 34.7930, 105.3680) ( 35.0380, 104.1050) ( 26.6760, 101.2450)
+( 27.3700, 102.5480) ( 33.9100, 106.2300) ( 26.6650, 100.7560) ( 29.2630, 102.4380) ( 28.9580, 103.8930)
+( 26.8260, 102.1000) ( 34.5150, 106.4000) ( 25.5760, 102.5050) ( 26.2110, 100.5960) ( 35.0050, 106.2060)
+( 29.6010, 103.4680) ( 25.7980, 102.9410) ( 30.1200, 103.1000) ( 35.0800, 105.7950) ( 34.4960, 108.2330)
+( 25.7960, 100.5600) ( 28.7700, 104.6000) ( 34.9460, 106.6780) ( 34.4330, 107.5800) ( 35.0580, 108.0860)
+( 28.9550, 102.7660) ( 34.0700, 107.6400) ( 34.9730, 108.9980) ( 35.0460, 104.5410) ( 29.9750, 103.0030)
+( 28.2500, 103.6400) ( 29.9750, 103.0030) ( 34.1100, 108.1560) ( 27.6930, 102.7900) ( 34.8950, 106.8210)
+( 28.6710, 102.5310) ( 26.4050, 103.2260) ( 27.9980, 102.8330) ( 33.8800, 109.9230) ( 25.0360, 100.5210)
+( 34.4260, 107.1430) ( 34.0880, 107.2950) ( 34.4710, 107.3780) ( 34.7200, 104.3800) ( 27.7700, 103.8910)
+( 33.6160, 106.9250) ( 34.3010, 108.1950) ( 34.3460, 109.9680) ( 31.0000, 107.1000) ( 34.9500, 109.9700)
+( 27.3560, 103.6860) ( 34.4600, 109.7060) ( 27.6830, 103.2680) ( 28.8430, 103.5260) ( 34.0500, 108.9080)
+( 28.6050, 103.9780) ( 29.3480, 102.6550) ( 33.5400, 107.9800) ( 33.5000, 109.2000) ( 28.3110, 103.1210)
+#
+#  1
+#  51
+#  (0.0, -100.0), (0.0, 400.0)0
+#
+#  2
+# 101     30.59521  31.92271
+# 101    103.49411 105.00661
+#===============================================================================
+# OUTPUTS
+# =======
+#
+# 1. select (1/0) output for los displacement (only for snapshots, see below),
+#    x, y, and z-cosines to the INSAR orbit: insar, xlos, ylos, zlos
+#
+#    if this option is selected, the snapshots will include additional data:
+#    LOS_Dsp = los displacement to the given satellite orbit.
+#
+# 2. select (1/0) output for Coulomb stress changes (only for snapshots, see
+#    below): icmb, friction, Skempton ratio, strike, dip, and rake angles [deg]
+#    describing the uniform regional master fault mechanism, the uniform regional
+#    principal stresses: sigma1, sigma2 and sigma3 [Pa] in arbitrary order (the
+#    orietation of the pre-stress field will be derived by assuming that the
+#    master fault is optimally oriented according to Coulomb failure criterion)
+#
+#    if this option is selected (icmb = 1), the snapshots will include additional
+#    data:
+#    CMB_Fix, Sig_Fix = Coulomb and normal stress changes on master fault;
+#    CMB_Op1/2, Sig_Op1/2 = Coulomb and normal stress changes on the two optimally
+#                       oriented faults;
+#    Str_Op1/2, Dip_Op1/2, Slp_Op1/2 = strike, dip and rake angles of the two
+#                       optimally oriented faults.
+#
+#    Note: the 1. optimally orieted fault is the one closest to the master fault.
+#
+# 3. output directory in char format: outdir
+#
+# 4. select outputs for displacement components (1/0 = yes/no): itout(i), i=1,3
+#
+# 5. the file names in char format for the x, y, and z components:
+#    toutfile(i), i=1,3
+#
+# 6. select outputs for stress components (1/0 = yes/no): itout(i), i=4,9
+#
+# 7. the file names in char format for the xx, yy, zz, xy, yz, and zx components:
+#    toutfile(i), i=4,9
+#
+# 8. select outputs for vertical NS and EW tilt components, block rotation, geoid
+#    and gravity changes (1/0 = yes/no): itout(i), i=10,14
+#
+# 9. the file names in char format for the NS tilt (positive if borehole top
+#    tilts to north), EW tilt (positive if borehole top tilts to east), block
+#    rotation (clockwise positive), geoid and gravity changes: toutfile(i), i=10,14
+#
+#    Note that all above outputs are time series with the time window as same
+#    as used for the Green's functions
+#
+#10. number of scenario outputs ("snapshots": spatial distribution of all above
+#    observables at given time points; <= NSCENMAX (see pscglob.h): nsc
+#
+#11. the time [day], and file name (in char format) for the 1. snapshot;
+#12. the time [day], and file name (in char format) for the 2. snapshot;
+#13. ...
+#
+#    Note that all file or directory names should not be longer than 80
+#    characters. Directories must be ended by / (unix) or \ (dos)!
+#===============================================================================
+ 1    0.0  0.0  -1.00      !displacement upward positive
+ 0     0.700  0.000  330.000   90.000  180.000    0.0E+00    0.0E+00    0.0E+00
+'.\'
+  0           0           0
+  'ux.dat'    'uy.dat'    'uz.dat'
+  0           0           0           0            0           0
+  'sxx.dat'   'syy.dat'   'szz.dat'   'sxy.dat'    'syz.dat'   'szx.dat'
+  0           0           0           0           0
+  'tx.dat'    'ty.dat'    'rot.dat'   'gd.dat'    'gr.dat'
+  1
+     0.00  'coseis-gps.dat'          |0 co-seismic
+#===============================================================================
+#
+# GREEN'S FUNCTION DATABASE
+# =========================
+# 1. directory where the Green's functions are stored: grndir
+#
+# 2. file names (without extensions!) for the 13 Green's functions:
+#    3 displacement komponents (uz, ur, ut): green(i), i=1,3
+#    6 stress components (szz, srr, stt, szr, srt, stz): green(i), i=4,9
+#    radial and tangential components measured by a borehole tiltmeter,
+#    rigid rotation around z-axis, geoid and gravity changes (tr, tt, rot, gd, gr):
+#    green(i), i=10,14
+#
+#    Note that all file or directory names should not be longer than 80
+#    characters. Directories must be ended by / (unix) or \ (dos)! The
+#    extensions of the file names will be automatically considered. They
+#    are ".ep", ".ss", ".ds" and ".cl" denoting the explosion (inflation)
+#    strike-slip, the dip-slip and the compensated linear vector dipole
+#    sources, respectively.
+#
+#===============================================================================
+ '..\wcpsgrnfcts\'
+ 'uz'  'ur'  'ut'
+ 'szz' 'srr' 'stt' 'szr' 'srt' 'stz'
+ 'tr'  'tt'  'rot' 'gd'  'gr'
+#===============================================================================
+# RECTANGULAR SUBFAULTS
+# =====================
+# 1. number of subfaults (<= NSMAX in pscglob.h), latitude [deg] and east
+#    longitude [deg] of the regional reference point as  origin of the Cartesian
+#    coordinate system: ns, lat0, lon0
+#
+# 2. parameters for the 1. rectangular subfault: geographic coordinates
+#    (O_lat, O_lon) [deg] and O_depth [km] of the local reference point on
+#    the present fault plane, length (along strike) [km] and width (along down
+#    dip) [km], strike [deg], dip [deg], number of equi-size fault
+#    patches along the strike (np_st) and along the dip (np_di) (total number of
+#    fault patches = np_st x np_di), and the start time of the rupture; the
+#    following data lines describe the slip distribution on the present sub-
+#    fault:
+#
+#    pos_s[km]  pos_d[km]  slip_along_strike[m]  slip_along_dip[m]  opening[m]
+#
+#    where (pos_s,pos_d) defines the position of the center of each patch in
+#    the local coordinate system with the origin at the reference point:
+#    pos_s = distance along the length (positive in the strike direction)
+#    pos_d = distance along the width (positive in the down-dip direction)
+#
+#
+# 3. ... for the 2. subfault ...
+# ...
+#                   N
+#                  /
+#                 /| strike
+#                +------------------------
+#                |\        p .            \ W
+#                :-\      i .              \ i
+#                |  \    l .                \ d
+#                :90 \  S .                  \ t
+#                |-dip\  .                    \ h
+#                :     \. | rake               \
+#                Z      -------------------------
+#                              L e n g t h
+#
+#    Note that a point inflation can be simulated by three point openning
+#    faults (each causes a third part of the volume of the point inflation)
+#    with orientation orthogonal to each other. the results obtained should
+#    be multiplied by a scaling factor 3(1-nu)/(1+nu), where nu is the Poisson
+#    ratio at the source. The scaling factor is the ratio of the seismic
+#    moment (energy) of an inflation source to that of a tensile source inducing
+#    a plate openning with the same volume change.
+#===============================================================================
+# n_faults (Slip model by Ji Chen, USGS)
+#-------------------------------------------------------------------------------
+  1
+#-------------------------------------------------------------------------------
+# n   O_lat   O_lon    O_depth length  width strike dip   np_st np_di start_time
+# [-] [deg]   [deg]    [km]    [km]     [km] [deg]  [deg] [-]   [-]   [day]
+#     pos_s   pos_d    slp_stk slp_dip open
+#     [km]    [km]     [m]     [m]     [m]
+#-------------------------------------------------------------------------------
+ 1    32.5224 105.4260 0.7411  315.00  40.00 229.00 33.00  21   8     0.00
+    7.50    2.50    0.00    0.00    0.00
+   22.50    2.50    0.57   -0.11    0.00
+   37.50    2.50    1.18   -0.38    0.00
+   52.50    2.50    0.85   -0.03    0.00
+   67.50    2.50   -0.03   -0.27    0.00
+   82.50    2.50   -0.54   -0.47    0.00
+   97.50    2.50   -0.37   -1.16    0.00
+  112.50    2.50    0.53   -1.68    0.00
+  127.50    2.50    0.50   -2.67    0.00
+  142.50    2.50    1.02   -2.57    0.00
+  157.50    2.50    0.21   -2.18    0.00
+  172.50    2.50   -0.82   -1.52    0.00
+  187.50    2.50   -1.47   -1.12    0.00
+  202.50    2.50   -2.24   -0.75    0.00
+  217.50    2.50   -2.58   -0.78    0.00
+  232.50    2.50   -2.00   -1.33    0.00
+  247.50    2.50   -1.01   -0.17    0.00
+  262.50    2.50    0.15   -0.15    0.00
+  277.50    2.50    0.48   -1.60    0.00
+  292.50    2.50    0.75   -1.34    0.00
+  307.50    2.50   -0.03   -0.04    0.00
+    7.50    7.50   -0.01    0.00    0.00
+   22.50    7.50    1.12   -0.07    0.00
+   37.50    7.50    1.06   -0.02    0.00
+   52.50    7.50    0.21   -0.25    0.00
+   67.50    7.50   -1.35   -0.60    0.00
+   82.50    7.50   -1.55   -1.04    0.00
+   97.50    7.50   -0.89   -2.67    0.00
+  112.50    7.50   -0.47   -3.92    0.00
+  127.50    7.50   -0.52   -5.11    0.00
+  142.50    7.50    0.33   -4.93    0.00
+  157.50    7.50   -0.03   -3.99    0.00
+  172.50    7.50   -1.25   -2.75    0.00
+  187.50    7.50   -2.64   -2.56    0.00
+  202.50    7.50   -4.41   -2.49    0.00
+  217.50    7.50   -5.55   -2.88    0.00
+  232.50    7.50   -4.60   -2.46    0.00
+  247.50    7.50   -2.85   -0.35    0.00
+  262.50    7.50   -0.42   -0.11    0.00
+  277.50    7.50    0.44   -0.89    0.00
+  292.50    7.50    0.61   -1.69    0.00
+  307.50    7.50    0.00   -0.08    0.00
+    7.50   12.50   -0.01   -0.04    0.00
+   22.50   12.50    0.63   -0.05    0.00
+   37.50   12.50   -0.06   -0.15    0.00
+   52.50   12.50   -2.17   -0.51    0.00
+   67.50   12.50   -4.15   -1.11    0.00
+   82.50   12.50   -3.91   -2.45    0.00
+   97.50   12.50   -2.88   -3.56    0.00
+  112.50   12.50   -2.43   -4.50    0.00
+  127.50   12.50   -2.12   -5.64    0.00
+  142.50   12.50   -1.30   -5.64    0.00
+  157.50   12.50   -2.00   -4.71    0.00
+  172.50   12.50   -3.09   -3.88    0.00
+  187.50   12.50   -3.95   -3.60    0.00
+  202.50   12.50   -6.18   -4.41    0.00
+  217.50   12.50   -7.78   -5.01    0.00
+  232.50   12.50   -6.52   -3.60    0.00
+  247.50   12.50   -4.02   -0.88    0.00
+  262.50   12.50   -1.76   -0.11    0.00
+  277.50   12.50   -0.84   -0.36    0.00
+  292.50   12.50   -0.47   -1.87    0.00
+  307.50   12.50    0.07   -0.06    0.00
+    7.50   17.50    0.02   -0.01    0.00
+   22.50   17.50    0.31   -0.15    0.00
+   37.50   17.50   -1.47   -0.02    0.00
+   52.50   17.50   -4.91   -0.14    0.00
+   67.50   17.50   -7.21   -1.25    0.00
+   82.50   17.50   -7.52   -1.88    0.00
+   97.50   17.50   -5.79   -1.79    0.00
+  112.50   17.50   -4.80   -2.79    0.00
+  127.50   17.50   -3.92   -3.21    0.00
+  142.50   17.50   -3.71   -4.62    0.00
+  157.50   17.50   -3.70   -4.03    0.00
+  172.50   17.50   -4.29   -2.87    0.00
+  187.50   17.50   -4.69   -2.63    0.00
+  202.50   17.50   -6.46   -4.26    0.00
+  217.50   17.50   -7.50   -5.18    0.00
+  232.50   17.50   -6.14   -4.57    0.00
+  247.50   17.50   -4.25   -2.55    0.00
+  262.50   17.50   -1.55   -1.43    0.00
+  277.50   17.50   -1.39   -1.01    0.00
+  292.50   17.50   -1.11   -2.57    0.00
+  307.50   17.50    0.00    0.00    0.00
+    7.50   22.50    0.02   -0.01    0.00
+   22.50   22.50    0.28   -0.04    0.00
+   37.50   22.50   -0.63   -0.06    0.00
+   52.50   22.50   -5.03   -0.24    0.00
+   67.50   22.50   -7.51   -2.18    0.00
+   82.50   22.50   -8.91   -2.61    0.00
+   97.50   22.50   -7.24   -1.05    0.00
+  112.50   22.50   -5.72   -1.23    0.00
+  127.50   22.50   -5.06   -1.17    0.00
+  142.50   22.50   -4.42   -2.54    0.00
+  157.50   22.50   -4.19   -2.16    0.00
+  172.50   22.50   -4.29   -0.89    0.00
+  187.50   22.50   -4.80   -1.75    0.00
+  202.50   22.50   -5.11   -3.93    0.00
+  217.50   22.50   -4.98   -5.16    0.00
+  232.50   22.50   -4.69   -5.14    0.00
+  247.50   22.50   -3.12   -3.77    0.00
+  262.50   22.50   -1.31   -2.97    0.00
+  277.50   22.50   -1.59   -2.25    0.00
+  292.50   22.50   -1.59   -3.28    0.00
+  307.50   22.50   -0.02   -0.01    0.00
+    7.50   27.50   -0.03   -0.03    0.00
+   22.50   27.50    0.09   -0.08    0.00
+   37.50   27.50   -0.66   -0.09    0.00
+   52.50   27.50   -4.27   -0.07    0.00
+   67.50   27.50   -6.66   -1.32    0.00
+   82.50   27.50   -7.54   -2.60    0.00
+   97.50   27.50   -5.74   -1.71    0.00
+  112.50   27.50   -4.45   -0.63    0.00
+  127.50   27.50   -3.16   -0.12    0.00
+  142.50   27.50   -3.12   -1.20    0.00
+  157.50   27.50   -2.97   -1.41    0.00
+  172.50   27.50   -2.52   -0.45    0.00
+  187.50   27.50   -2.48   -1.17    0.00
+  202.50   27.50   -2.36   -2.49    0.00
+  217.50   27.50   -2.47   -4.33    0.00
+  232.50   27.50   -2.12   -4.76    0.00
+  247.50   27.50   -0.98   -3.36    0.00
+  262.50   27.50   -0.45   -2.94    0.00
+  277.50   27.50   -0.74   -3.59    0.00
+  292.50   27.50   -2.06   -3.60    0.00
+  307.50   27.50   -0.03   -0.04    0.00
+    7.50   32.50    0.01   -0.04    0.00
+   22.50   32.50   -0.15    0.00    0.00
+   37.50   32.50    0.00   -0.01    0.00
+   52.50   32.50   -2.19   -0.01    0.00
+   67.50   32.50   -3.58   -0.08    0.00
+   82.50   32.50   -3.91   -1.27    0.00
+   97.50   32.50   -2.45   -0.84    0.00
+  112.50   32.50   -1.61   -0.39    0.00
+  127.50   32.50   -0.92   -0.12    0.00
+  142.50   32.50   -1.02   -0.17    0.00
+  157.50   32.50   -1.08   -0.55    0.00
+  172.50   32.50   -0.45   -0.24    0.00
+  187.50   32.50   -0.22   -0.05    0.00
+  202.50   32.50   -0.68   -0.76    0.00
+  217.50   32.50   -1.37   -2.36    0.00
+  232.50   32.50   -1.23   -2.35    0.00
+  247.50   32.50   -0.19   -1.36    0.00
+  262.50   32.50   -0.08   -1.38    0.00
+  277.50   32.50   -0.29   -2.50    0.00
+  292.50   32.50   -1.82   -2.41    0.00
+  307.50   32.50   -0.02   -0.05    0.00
+    7.50   37.50    0.06   -0.03    0.00
+   22.50   37.50    0.04   -0.03    0.00
+   37.50   37.50    0.01   -0.04    0.00
+   52.50   37.50   -0.02   -0.02    0.00
+   67.50   37.50   -0.03   -0.01    0.00
+   82.50   37.50   -0.01   -0.04    0.00
+   97.50   37.50   -0.02   -0.05    0.00
+  112.50   37.50   -0.05   -0.07    0.00
+  127.50   37.50    0.05   -0.01    0.00
+  142.50   37.50   -0.03   -0.01    0.00
+  157.50   37.50    0.03   -0.09    0.00
+  172.50   37.50    0.03   -0.02    0.00
+  187.50   37.50    0.06   -0.04    0.00
+  202.50   37.50    0.07   -0.03    0.00
+  217.50   37.50    0.00   -0.01    0.00
+  232.50   37.50   -0.04   -0.04    0.00
+  247.50   37.50    0.04   -0.08    0.00
+  262.50   37.50    0.00   -0.04    0.00
+  277.50   37.50   -0.03    0.00    0.00
+  292.50   37.50   -0.06   -0.02    0.00
+  307.50   37.50    0.09   -0.02    0.00
+#================================end of input===================================
--- a/examples/pscmp08_input.dat
+++ b/examples/pscmp08_input.dat
--- a/examples/psgrn08-wenchuan.inp
+++ b/examples/psgrn08-wenchuan.inp
@ -0,0 +1,145 @@
+#=============================================================================
+# This is input file of FORTRAN77 program "psgrn07a" for computing responses
+# (Green's functions) of a multi-layered viscoelastic halfspace to point
+# dislocation sources buried at different depths. All results will be stored in
+# the given directory and provide the necessary data base for the program
+# "pscmp07a" for computing time-dependent deformation, geoid and gravity changes
+# induced by an earthquake with extended fault planes via linear superposition.
+# For more details, please read the accompanying READ.ME file.
+#
+# written by Rongjiang Wang
+# GeoForschungsZentrum Potsdam
+# e-mail: wang@gfz-potsdam.de
+# phone +49 331 2881209
+# fax +49 331 2881204
+#
+# Last modified: Potsdam, Jan, 2008
+#
+#################################################################
+##                                                             ##
+## Cylindrical coordinates (Z positive downwards!) are used.   ##
+##                                                             ##
+## If not specified otherwise, SI Unit System is used overall! ##
+##                                                             ##
+#################################################################
+#
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR SOURCE-OBSERVATION CONFIGURATIONS
+#	================================================
+# 1. the uniform depth of the observation points [km], switch for oceanic (0)
+#    or continental(1) earthquakes;
+# 2. number of (horizontal) observation distances (> 1 and <= nrmax defined in
+#    psgglob.h), start and end distances [km], ratio (>= 1.0) between max. and
+#    min. sampling interval (1.0 for equidistant sampling);
+# 3. number of equidistant source depths (>= 1 and <= nzsmax defined in
+#    psgglob.h), start and end source depths [km];
+#
+#    r1,r2 = minimum and maximum horizontal source-observation
+#            distances (r2 > r1).
+#    zs1,zs2 = minimum and maximum source depths (zs2 >= zs1 > 0).
+#
+#    Note that the same sampling rates dr_min and dzs will be used later by the
+#    program "pscmp07a" for discretizing the finite source planes to a 2D grid
+#    of point sources.
+#------------------------------------------------------------------------------
+        0.0       1
+  101   0.0   1000.0  10.0
+   59   0.5    29.5
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR TIME SAMPLING
+#	============================
+# 1. number of time samples (<= ntmax def. in psgglob.h) and time window [days].
+#
+#    Note that nt (> 0) should be power of 2 (the fft-rule). If nt = 1, the
+#    coseismic (t = 0) changes will be computed; If nt = 2, the coseismic
+#    (t = 0) and steady-state (t -> infinity) changes will be computed;
+#    Otherwise, time series for the given time samples will be computed.
+#
+#------------------------------------------------------------------------------
+ 1   15330.0
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR WAVENUMBER INTEGRATION
+#	=====================================
+# 1. relative accuracy of the wave-number integration (suggested: 0.1 - 0.01)
+# 2. factor (> 0 and < 1) for including influence of earth's gravity on the
+#    deformation field (e.g. 0/1 = without / with 100% gravity effect).
+#------------------------------------------------------------------------------
+ 0.025
+ 0.00
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR OUTPUT FILES
+#	===========================
+#
+# 1. output directory
+# 2. file names for 3 displacement components (uz, ur, ut)
+# 3. file names for 6 stress components (szz, srr, stt, szr, srt, stz)
+# 4. file names for radial and tangential tilt components (as measured by a
+#    borehole tiltmeter), rigid rotation of horizontal plane, geoid and gravity
+#    changes (tr, tt, rot, gd, gr)
+#
+#    Note that all file or directory names should not be longer than 80
+#    characters. Directory and subdirectoy names must be separated and ended
+#    by / (unix) or \ (dos)! All file names should be given without extensions
+#    that will be appended automatically by ".ep" for the explosion (inflation)
+#    source, ".ss" for the strike-slip source, ".ds" for the dip-slip source,
+#    and ".cl" for the compensated linear vector dipole source)
+#
+#------------------------------------------------------------------------------
+ './'
+ 'uz'  'ur'  'ut'
+ 'szz' 'srr' 'stt' 'szr' 'srt' 'stz'
+ 'tr'  'tt'  'rot' 'gd'  'gr'
+#------------------------------------------------------------------------------
+#
+#	GLOBAL MODEL PARAMETERS
+#	=======================
+# 1. number of data lines of the layered model (<= lmax as defined in psgglob.h)
+#
+#    The surface and the upper boundary of the half-space as well as the
+#    interfaces at which the viscoelastic parameters are continuous, are all
+#    defined by a single data line; All other interfaces, at which the
+#    viscoelastic parameters are discontinuous, are all defined by two data
+#    lines (upper-side and lower-side values). This input format could also be
+#    used for a graphic plot of the layered model. Layers which have different
+#    parameter values at top and bottom, will be treated as layers with a
+#    constant gradient, and will be discretised to a number of homogeneous
+#    sublayers. Errors due to the discretisation are limited within about 5%
+#    (changeable, see psgglob.h).
+#
+# 2....	parameters of the multilayered model
+#
+#    Burgers rheology (a Kelvin-Voigt body and a Maxwell body in series
+#    connection) for relaxation of shear modulus is implemented. No relaxation
+#    of compressional modulus is considered.
+#
+#    eta1  = transient viscosity (dashpot of the Kelvin-Voigt body; <= 0 means
+#            infinity value)
+#    eta2  = steady-state viscosity (dashpot of the Maxwell body; <= 0 means
+#            infinity value)
+#    alpha = ratio between the effective and the unrelaxed shear modulus
+#            = mu1/(mu1+mu2) (> 0 and <= 1)
+#
+#    Special cases:
+#        (1) Elastic: eta1 and eta2 <= 0 (i.e. infinity); alpha meaningless
+#        (2) Maxwell body: eta1 <= 0 (i.e. eta1 = infinity)
+#                          or alpha = 1 (i.e. mu1 = infinity)
+#        (3) Standard-Linear-Solid: eta2 <= 0 (i.e. infinity)
+#------------------------------------------------------------------------------
+  9                               |int: no_model_lines;
+#------------------------------------------------------------------------------
+# no  depth[km]  vp[km/s]  vs[km/s]  rho[kg/m^3] eta1[Pa*s] eta2[Pa*s] alpha
+#------------------------------------------------------------------------------
+1       0.0       2.5000    1.2000    2100.0     0.0E+00    0.0E+00    1.000
+2       1.5       2.5000    1.2000    2100.0     0.0E+00    0.0E+00    1.000
+3       1.5       4.5000    2.6000    2500.0     0.0E+00    0.0E+00    1.000
+4       8.0       4.5000    2.6000    2500.0     0.0E+00    0.0E+00    1.000
+5       8.0       6.2000    3.6000    2800.0     0.0E+00    0.0E+00    1.000
+6      17.0       6.2000    3.6000    2800.0     0.0E+00    0.0E+00    1.000
+7      17.0       6.4000    3.6000    2850.0     0.0E+00    0.0E+00    1.000
+8      29.0       6.4000    3.6000    2850.0     0.0E+00    0.0E+00    1.000
+9      29.0       6.8000    3.8000    2950.0     0.0E+00    0.0E+00    1.000
+#=======================end of input===========================================
--- a/examples/psgrn08_input.dat
+++ b/examples/psgrn08_input.dat
@ -0,0 +1,159 @@
+#=============================================================================
+# This is input file of FORTRAN77 program "psgrn08" for computing responses
+# (Green's functions) of a multi-layered viscoelastic halfspace to point
+# dislocation sources buried at different depths. All results will be stored in
+# the given directory and provide the necessary data base for the program
+# "pscmp07a" for computing time-dependent deformation, geoid and gravity changes
+# induced by an earthquake with extended fault planes via linear superposition.
+# For more details, please read the accompanying READ.ME file.
+#
+# written by Rongjiang Wang
+# GeoForschungsZentrum Potsdam
+# e-mail: wang@gfz-potsdam.de
+# phone +49 331 2881209
+# fax +49 331 2881204
+#
+# Last modified: Potsdam, July, 2008
+#
+# References:
+#
+# (1) Wang, R., F. Lorenzo-Martín and F. Roth (2003), Computation of deformation
+#     induced by earthquakes in a multi-layered elastic crust - FORTRAN programs
+#     EDGRN/EDCMP, Computer and Geosciences, 29(2), 195-207.
+# (2) Wang, R., F. Lorenzo-Martin and F. Roth (2006), PSGRN/PSCMP - a new code for
+#     calculating co- and post-seismic deformation, geoid and gravity changes
+#     based on the viscoelastic-gravitational dislocation theory, Computers and
+#     Geosciences, 32, 527-541. DOI:10.1016/j.cageo.2005.08.006.
+# (3) Wang, R. (2005), The dislocation theory: a consistent way for including the
+#     gravity effect in (visco)elastic plane-earth models, Geophysical Journal
+#     International, 161, 191-196.
+#
+#################################################################
+##                                                             ##
+## Cylindrical coordinates (Z positive downwards!) are used.   ##
+##                                                             ##
+## If not specified otherwise, SI Unit System is used overall! ##
+##                                                             ##
+#################################################################
+#
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR SOURCE-OBSERVATION CONFIGURATIONS
+#	================================================
+# 1. the uniform depth of the observation points [km], switch for oceanic (0)
+#    or continental(1) earthquakes;
+# 2. number of (horizontal) observation distances (> 1 and <= nrmax defined in
+#    psgglob.h), start and end distances [km], ratio (>= 1.0) between max. and
+#    min. sampling interval (1.0 for equidistant sampling);
+# 3. number of equidistant source depths (>= 1 and <= nzsmax defined in
+#    psgglob.h), start and end source depths [km];
+#
+#    r1,r2 = minimum and maximum horizontal source-observation
+#            distances (r2 > r1).
+#    zs1,zs2 = minimum and maximum source depths (zs2 >= zs1 > 0).
+#
+#    Note that the same sampling rates dr_min and dzs will be used later by the
+#    program "pscmp08" for discretizing the finite source planes to a 2D grid
+#    of point sources.
+#------------------------------------------------------------------------------
+        0.0       0
+   73   0.0  2000.0  10.0
+   30   1.0    59.0
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR TIME SAMPLING
+#	============================
+# 1. number of time samples (<= ntmax def. in psgglob.h) and time window [days].
+#
+#    Note that nt (> 0) should be power of 2 (the fft-rule). If nt = 1, the
+#    coseismic (t = 0) changes will be computed; If nt = 2, the coseismic
+#    (t = 0) and steady-state (t -> infinity) changes will be computed;
+#    Otherwise, time series for the given time samples will be computed.
+#
+#------------------------------------------------------------------------------
+ 512   46628.75
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR WAVENUMBER INTEGRATION
+#	=====================================
+# 1. relative accuracy of the wave-number integration (suggested: 0.1 - 0.01)
+# 2. factor (> 0 and < 1) for including influence of earth's gravity on the
+#    deformation field (e.g. 0/1 = without / with 100% gravity effect).
+#------------------------------------------------------------------------------
+ 0.025
+ 1.00
+#------------------------------------------------------------------------------
+#
+#	PARAMETERS FOR OUTPUT FILES
+#	===========================
+#
+# 1. output directory
+# 2. file names for 3 displacement components (uz, ur, ut)
+# 3. file names for 6 stress components (szz, srr, stt, szr, srt, stz)
+# 4. file names for radial and tangential tilt components (as measured by a
+#    borehole tiltmeter), rigid rotation of horizontal plane, geoid and gravity
+#    changes (tr, tt, rot, gd, gr)
+#
+#    Note that all file or directory names should not be longer than 80
+#    characters. Directory and subdirectoy names must be separated and ended
+#    by / (unix) or \ (dos)! All file names should be given without extensions
+#    that will be appended automatically by ".ep" for the explosion (inflation)
+#    source, ".ss" for the strike-slip source, ".ds" for the dip-slip source,
+#    and ".cl" for the compensated linear vector dipole source)
+#
+#------------------------------------------------------------------------------
+ './'
+ 'uz'  'ur'  'ut'
+ 'szz' 'srr' 'stt' 'szr' 'srt' 'stz'
+ 'tr'  'tt'  'rot' 'gd'  'gr'
+#------------------------------------------------------------------------------
+#
+#	GLOBAL MODEL PARAMETERS
+#	=======================
+# 1. number of data lines of the layered model (<= lmax as defined in psgglob.h)
+#
+#    The surface and the upper boundary of the half-space as well as the
+#    interfaces at which the viscoelastic parameters are continuous, are all
+#    defined by a single data line; All other interfaces, at which the
+#    viscoelastic parameters are discontinuous, are all defined by two data
+#    lines (upper-side and lower-side values). This input format could also be
+#    used for a graphic plot of the layered model. Layers which have different
+#    parameter values at top and bottom, will be treated as layers with a
+#    constant gradient, and will be discretised to a number of homogeneous
+#    sublayers. Errors due to the discretisation are limited within about 5%
+#    (changeable, see psgglob.h).
+#
+# 2....	parameters of the multilayered model
+#
+#    Burgers rheology [a Kelvin-Voigt body (mu1, eta1) and a Maxwell body
+#    (mu2, eta2) in series connection] for relaxation of shear modulus is
+#    implemented. No relaxation of compressional modulus is considered.
+#
+#    eta1  = transient viscosity (dashpot of the Kelvin-Voigt body; <= 0 means
+#            infinity value)
+#    eta2  = steady-state viscosity (dashpot of the Maxwell body; <= 0 means
+#            infinity value)
+#    alpha = ratio between the effective and the unrelaxed shear modulus
+#            = mu1/(mu1+mu2) (> 0 and <= 1) (unrelaxed modulus mu2 is
+#            derived from S wave velocity and density)
+#
+#    Special cases:
+#        (1) Elastic: eta1 and eta2 <= 0 (i.e. infinity); alpha meaningless
+#        (2) Maxwell body: eta1 <= 0 (i.e. eta1 = infinity)
+#                          or alpha = 1 (i.e. mu1 = infinity)
+#        (3) Standard-Linear-Solid: eta2 <= 0 (i.e. infinity)
+#            fully relaxed modulus = alpha*unrelaxed_modulus
+#            characteristic relaxation time = eta1*alpha/unrelaxed_modulus
+#------------------------------------------------------------------------------
+  7                               |int: no_model_lines;
+#------------------------------------------------------------------------------
+# no  depth[km]  vp[km/s]  vs[km/s]  rho[kg/m^3] eta1[Pa*s] eta2[Pa*s] alpha
+#------------------------------------------------------------------------------
+  1      0.000   6.0000    3.4600    2600.0      0.0E+00    0.0E+00    1.000
+  2     16.000   6.0000    3.4600    2600.0      0.0E+00    0.0E+00    1.000
+  3     16.000   6.7000    3.8700    2800.0      0.0E+00    0.0E+00    1.000
+  4     30.000   6.7000    3.8700    2800.0      0.0E+00    0.0E+00    1.000
+  5     30.000   8.0000    4.6200    3400.0      0.0E+00    0.0E+00    1.000
+  6     60.000   8.0000    4.6200    3400.0      0.0E+00    0.0E+00    1.000
+  7     60.000   8.0000    4.6200    3400.0      0.0E+00    1.0E+19    1.000
+#=======================end of input===========================================
--- a/src/pscmp/Makefile.am
+++ b/src/pscmp/Makefile.am
@ -0,0 +1,2 @@
+bin_PROGRAMS = fomosto_pscmp2008a
+fomosto_pscmp2008a_SOURCES = cmbfix.f cmbopt.f dc3d.f disazi.f getdata.f mscorr.f prestress.f pscdisc.f pscglob.h pscgrn.f pscmain.f pscokada.f pscout.f roots3.f
--- a/src/pscmp/cmbfix.f
+++ b/src/pscmp/cmbfix.f
@ -0,0 +1,65 @@
+      subroutine cmbfix(sxx,syy,szz,sxy,syz,szx,
+     &                  p,f,cmb,sig,st,di,ra)
+      implicit none
+c
+c     calculate Coulomb stress
+c
+c     input:
+c     stress tensor, pore pressure, friction coefficient
+c     rupture orientation parameter (strike, dip and rake)
+c
+      double precision sxx,syy,szz,sxy,syz,szx,p,f,cmb,sig,st,di,ra
+c
+c     return:
+c     Coulomb stress (cmb) and normal stress (sig)
+c
+c     local memories:
+c
+      integer i,j
+      double precision pi,deg2rad,st0,di0,ra0,tau
+      double precision s(3,3),ns(3),ts(3),rst(3),rdi(3)
+c
+      pi=4.d0*datan(1.d0)
+      deg2rad=pi/180.d0
+      st0=st*deg2rad
+      di0=di*deg2rad
+      ra0=ra*deg2rad
+c
+      s(1,1)=sxx
+      s(1,2)=sxy
+      s(1,3)=szx
+      s(2,1)=sxy
+      s(2,2)=syy
+      s(2,3)=syz
+      s(3,1)=szx
+      s(3,2)=syz
+      s(3,3)=szz
+c
+      ns(1)=dsin(di0)*dcos(st0+0.5d0*pi)
+      ns(2)=dsin(di0)*dsin(st0+0.5d0*pi)
+      ns(3)=-dcos(di0)
+c
+      rst(1)=dcos(st0)
+      rst(2)=dsin(st0)
+      rst(3)=0.d0
+c
+      rdi(1)=dcos(di0)*dcos(st0+0.5d0*pi)
+      rdi(2)=dcos(di0)*dsin(st0+0.5d0*pi)
+      rdi(3)=dsin(di0)
+c
+      do i=1,3
+        ts(i)=rst(i)*dcos(ra0)-rdi(i)*dsin(ra0)
+      enddo
+c
+      sig=0.d0
+      tau=0.d0
+      do j=1,3
+        do i=1,3
+          sig=sig+ns(i)*s(i,j)*ns(j)
+          tau=tau+ts(i)*s(i,j)*ns(j)
+        enddo
+      enddo
+c
+      cmb=tau+f*(sig+p)
+      return
+      end
--- a/src/pscmp/cmbopt.f
+++ b/src/pscmp/cmbopt.f
@ -0,0 +1,230 @@
+      subroutine cmbopt(sxx,syy,szz,sxy,syz,szx,p,f,key,
+     &                  st0,di0,ra0,
+     &                  cmb,sig,st1,di1,ra1,st2,di2,ra2)
+      implicit none
+c
+c     Coulomb stress with the optimal orientation
+c
+c     input:
+c     stress tensor, pore pressure, friction coefficient
+c     key = 0: determine optimal Coulomb stress only;
+c           1: determine optimal Coulomb stress and orientations
+c
+      integer key
+      double precision sxx,syy,szz,sxy,syz,szx,p,f
+c
+c     output
+c     max. Coulomb stress at the two optimally oriented fault planes
+c     sig = normal stress
+c
+      double precision st0,di0,ra0,cmb,sig,st1,di1,ra1,st2,di2,ra2
+c
+c     local memories:
+c
+      integer i,j,j0,j1,j2,jmin,jmax
+      double precision pi,b,c,d,s1,s2,s3,snn,alpha,am,swap
+      double precision cmb1,cmb2,cmb3,det1,det2,det3,detmax,rmax
+      double precision s(3),r(3,3),ns(3,2),ts(3,2)
+      double precision mscorr
+c
+      pi=4.d0*datan(1.d0)
+c
+      if(sxy.eq.0.d0.and.syz.eq.0.d0.and.szx.eq.0.d0)then
+        s(1)=sxx
+        s(2)=syy
+        s(3)=szz
+      else
+        b=-(sxx+syy+szz)
+        c=sxx*syy+syy*szz+szz*sxx-sxy**2-syz**2-szx**2
+        d=sxx*syz**2+syy*szx**2+szz*sxy**2-2.d0*sxy*syz*szx-sxx*syy*szz
+        call roots3(b,c,d,s)
+      endif
+c
+      cmb1=0.5d0*dabs(s(2)-s(3))*dsqrt(1+f*f)+f*(0.5d0*(s(2)+s(3))+p)
+      cmb2=0.5d0*dabs(s(3)-s(1))*dsqrt(1+f*f)+f*(0.5d0*(s(3)+s(1))+p)
+      cmb3=0.5d0*dabs(s(1)-s(2))*dsqrt(1+f*f)+f*(0.5d0*(s(1)+s(2))+p)
+c
+      cmb=dmax1(cmb1,cmb2,cmb3)
+      st1=0.d0
+      di1=0.d0
+      ra1=0.d0
+      st2=0.d0
+      di2=0.d0
+      ra2=0.d0
+      if(key.eq.0.or.s(1).eq.s(2).and.s(2).eq.s(3))return
+c
+      if(cmb.eq.cmb1)then
+        s3=s(1)
+        s1=dmax1(s(2),s(3))
+        s2=dmin1(s(2),s(3))
+      else if(cmb.eq.cmb2)then
+        s1=dmax1(s(3),s(1))
+        s2=dmin1(s(3),s(1))
+        s3=s(2)
+      else
+        s1=dmax1(s(1),s(2))
+        s2=dmin1(s(1),s(2))
+        s3=s(3)
+      endif
+      sig=0.5d0*((s1-s2)*f/dsqrt(1+f*f)+s1+s2)
+      s(1)=s1
+      s(2)=s2
+      s(3)=s3
+c
+c     determine eigenvectors (the principal stress directions)
+c
+      j0=0
+      if(s(1).eq.s(2))then
+        j0=3
+        j1=1
+        j2=2
+      else if(s(2).eq.s(3))then
+        j0=1
+        j1=2
+        j2=3
+      else if(s(3).eq.s(1))then
+        j0=2
+        j1=1
+        j2=3
+      endif
+c
+      if(j0.eq.0)then
+        jmin=1
+        jmax=3
+      else
+        jmin=j0
+        jmax=j0
+        print *,' Warning: more than two optimal rupture orientations!'
+      endif
+c
+      do j=jmin,jmax
+        det1=syz*syz-(syy-s(j))*(szz-s(j))
+        det2=szx*szx-(sxx-s(j))*(szz-s(j))
+        det3=sxy*sxy-(sxx-s(j))*(syy-s(j))
+        detmax=dmax1(dabs(det1),dabs(det2),dabs(det3))
+        if(dabs(det1).eq.detmax)then
+          r(1,j)=det1
+          r(2,j)=(szz-s(j))*sxy-syz*szx
+          r(3,j)=(syy-s(j))*szx-syz*sxy
+        else if(dabs(det2).eq.detmax)then
+          r(1,j)=(szz-s(j))*sxy-szx*syz
+          r(2,j)=det2
+          r(3,j)=(sxx-s(j))*syz-szx*sxy
+        else
+          r(1,j)=(syy-s(j))*szx-sxy*syz
+          r(2,j)=(sxx-s(j))*syz-sxy*szx
+          r(3,j)=det3
+        endif
+      enddo
+c
+c     if any two eigenvalues are identical, their corresponding
+c     eigenvectors should be redetermined by orthogonalizing
+c     them to the 3. eigenvector as well as to each other
+c     
+      if(j0.gt.0)then
+        rmax=dmax1(dabs(r(1,j0)),dabs(r(2,j0)),dabs(r(3,j0)))
+        if(dabs(r(1,j0)).eq.rmax)then
+          r(1,j1)=-r(2,j0)
+          r(2,j1)=r(1,j0)
+          r(3,j1)=0.d0
+c
+          r(1,j2)=-r(3,j0)
+          r(2,j2)=0.d0
+          r(3,j2)=r(1,j0)
+          am=r(1,j1)*r(1,j2)/(r(1,j1)**2+r(2,j1)**2)
+          do i=1,3
+            r(i,j2)=r(i,j2)-am*r(i,j1)
+          enddo
+        else if(dabs(r(2,j0)).eq.rmax)then
+          r(1,j1)=r(2,j0)
+          r(2,j1)=-r(1,j0)
+          r(3,j1)=0.d0
+c
+          r(1,j2)=0.d0
+          r(2,j2)=-r(3,j0)
+          r(3,j2)=r(2,j0)
+          am=r(2,j1)*r(2,j2)/(r(1,j1)**2+r(2,j1)**2)
+          do i=1,3
+            r(i,j2)=r(i,j2)-am*r(i,j1)
+          enddo
+        else if(dabs(r(3,j0)).eq.rmax)then
+          r(1,j1)=r(3,j0)
+          r(2,j1)=0.d0
+          r(3,j1)=-r(1,j0)
+c
+          r(1,j2)=0.d0
+          r(2,j2)=r(3,j0)
+          r(3,j2)=-r(2,j0)
+          am=r(3,j1)*r(3,j2)/(r(1,j1)**2+r(3,j1)**2)
+          do i=1,3
+            r(i,j2)=r(i,j2)-am*r(i,j1)
+          enddo
+        endif
+      endif
+c
+      do j=1,3
+        am=dsqrt(r(1,j)**2+r(2,j)**2+r(3,j)**2)
+        do i=1,3
+          r(i,j)=r(i,j)/am
+        enddo
+      enddo
+c
+      alpha=0.5d0*datan2(1.d0,f)
+      snn=s(1)*dcos(alpha)**2+s(2)*dsin(alpha)**2
+c
+c     determine the two optimal fault-plane normals
+c
+      do i=1,3
+        ns(i,1)=r(i,1)*dcos(alpha)+r(i,2)*dsin(alpha)          
+        ns(i,2)=r(i,1)*dcos(alpha)-r(i,2)*dsin(alpha)
+      enddo
+c
+c     determine the direction of max. shear stress
+c
+      do j=1,2
+        am=dsqrt(ns(1,j)**2+ns(2,j)**2+ns(3,j)**2)
+        if(ns(3,j).gt.0.d0)am=-am
+        do i=1,3
+          ns(i,j)=ns(i,j)/am
+        enddo
+        ts(1,j)=(sxx-snn)*ns(1,j)+sxy*ns(2,j)+szx*ns(3,j)
+        ts(2,j)=sxy*ns(1,j)+(syy-snn)*ns(2,j)+syz*ns(3,j)
+        ts(3,j)=szx*ns(1,j)+syz*ns(2,j)+(szz-snn)*ns(3,j)
+        am=dsqrt(ts(1,j)**2+ts(2,j)**2+ts(3,j)**2)
+        do i=1,3
+          ts(i,j)=ts(i,j)/am
+        enddo
+      enddo
+c
+c     determine the two optimal focal mechanisms
+c    
+      st1=dmod(datan2(ns(2,1),ns(1,1))*180.d0/pi+270.d0,360.d0)
+	di1=dacos(-ns(3,1))*180.d0/pi
+      s1=dcos(st1*pi/180.d0)
+      s2=dsin(st1*pi/180.d0)
+      ra1=dacos(dmin1(dmax1(s1*ts(1,1)+s2*ts(2,1),-1.d0),1.d0))
+     &   *180.d0/pi
+      if(ts(3,1).gt.0.d0)ra1=-ra1
+c
+      st2=dmod(datan2(ns(2,2),ns(1,2))*180.d0/pi+270.d0,360.d0)
+	di2=dacos(-ns(3,2))*180.d0/pi
+      s1=dcos(st2*pi/180.d0)
+      s2=dsin(st2*pi/180.d0)
+      ra2=dacos(dmin1(dmax1(s1*ts(1,2)+s2*ts(2,2),-1.d0),1.d0))
+     &   *180.d0/pi
+      if(ts(3,2).gt.0.d0)ra2=-ra2
+c
+      if(mscorr(st0,di0,ra0,st1,di1,ra1).lt.
+     &   mscorr(st0,di0,ra0,st2,di2,ra2))then
+        swap=st1
+        st1=st2
+        st2=swap
+        swap=di1
+        di1=di2
+        di2=swap
+        swap=ra1
+        ra1=ra2
+        ra2=swap
+      endif
+      return
+      end
--- a/src/pscmp/dc3d.f
+++ b/src/pscmp/dc3d.f
@ -0,0 +1,665 @@
+      SUBROUTINE  DC3D(ALPHA,X,Y,Z,DEPTH,DIP,                           04610005
+     *              AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,                  04620005
+     *              UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET)  04630005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         04640005
+      REAL*4   ALPHA,X,Y,Z,DEPTH,DIP,AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3, 04650005
+     *         UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ             04660005
+C                                                                       04670005
+C********************************************************************   04680005
+C*****                                                          *****   04690005
+C*****    DISPLACEMENT AND STRAIN AT DEPTH                      *****   04700005
+C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****   04710005
+C*****              CODED BY  Y.OKADA ... SEP.1991              *****   04720005
+C*****              REVISED ... NOV.1991, APR.1992, MAY.1993,   *****   04730005
+C*****                          JUL.1993                        *****   04740005
+C********************************************************************   04750005
+C                                                                       04760005
+C***** INPUT                                                            04770005
+C*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)           04780005
+C*****   X,Y,Z : COORDINATE OF OBSERVING POINT                          04790005
+C*****   DEPTH : DEPTH OF REFERENCE POINT                               04800005
+C*****   DIP   : DIP-ANGLE (DEGREE)                                     04810005
+C*****   AL1,AL2   : FAULT LENGTH RANGE                                 04820005
+C*****   AW1,AW2   : FAULT WIDTH RANGE                                  04830005
+C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS              04840005
+C                                                                       04850005
+C***** OUTPUT                                                           04860005
+C*****   UX, UY, UZ  : DISPLACEMENT ( UNIT=(UNIT OF DISL)               04870005
+C*****   UXX,UYX,UZX : X-DERIVATIVE ( UNIT=(UNIT OF DISL) /             04880005
+C*****   UXY,UYY,UZY : Y-DERIVATIVE        (UNIT OF X,Y,Z,DEPTH,AL,AW) )04890005
+C*****   UXZ,UYZ,UZZ : Z-DERIVATIVE                                     04900005
+C*****   IRET        : RETURN CODE  ( =0....NORMAL,   =1....SINGULAR )  04910005
+C                                                                       04920005
+      COMMON /C0/DUMMY(5),SD,CD,dumm(5)                                 04930005
+      DIMENSION  XI(2),ET(2),KXI(2),KET(2)                              04940005
+      DIMENSION  U(12),DU(12),DUA(12),DUB(12),DUC(12)                   04950005
+      DATA  F0,EPS/ 0.D0, 1.D-06 /                                      04960005
+C-----                                                                  04970005
+      IF(Z.GT.0.) WRITE(*,'('' ** POSITIVE Z WAS GIVEN IN SUB-DC3D'')') 04980005
+      DO 111 I=1,12                                                     04990005
+        U  (I)=F0                                                       05000005
+        DUA(I)=F0                                                       05010005
+        DUB(I)=F0                                                       05020005
+        DUC(I)=F0                                                       05030005
+  111 CONTINUE                                                          05040005
+      AALPHA=ALPHA                                                      05050005
+      DDIP=DIP                                                          05060005
+      CALL DCCON0(AALPHA,DDIP)                                          05070005
+C-----                                                                  05080005
+      ZZ=Z                                                              05090005
+      DD1=DISL1                                                         05100005
+      DD2=DISL2                                                         05110005
+      DD3=DISL3                                                         05120005
+      XI(1)=X-AL1                                                       05130005
+      XI(2)=X-AL2                                                       05140005
+      IF(DABS(XI(1)).LT.EPS) XI(1)=F0                                   05150005
+      IF(DABS(XI(2)).LT.EPS) XI(2)=F0                                   05160005
+C======================================                                 05170005
+C=====  REAL-SOURCE CONTRIBUTION  =====                                 05180005
+C======================================                                 05190005
+      D=DEPTH+Z                                                         05200005
+      P=Y*CD+D*SD                                                       05210005
+      Q=Y*SD-D*CD                                                       05220005
+      ET(1)=P-AW1                                                       05230005
+      ET(2)=P-AW2                                                       05240005
+      IF(DABS(Q).LT.EPS)  Q=F0                                          05250005
+      IF(DABS(ET(1)).LT.EPS) ET(1)=F0                                   05260005
+      IF(DABS(ET(2)).LT.EPS) ET(2)=F0                                   05270005
+C--------------------------------                                       05280005
+C----- REJECT SINGULAR CASE -----                                       05290005
+C--------------------------------                                       05300005
+C----- ON FAULT EDGE                                                    05310014
+      IF(Q.EQ.F0                                                        05320014
+     *   .AND.(    (XI(1)*XI(2).LE.F0 .AND. ET(1)*ET(2).EQ.F0)          05330014
+     *         .OR.(ET(1)*ET(2).LE.F0 .AND. XI(1)*XI(2).EQ.F0) ))       05340014
+     *   GO TO 99                                                       05350005
+C----- ON NEGATIVE EXTENSION OF FAULT EDGE                              05360014
+      KXI(1)=0                                                          05370005
+      KXI(2)=0                                                          05380005
+      KET(1)=0                                                          05390005
+      KET(2)=0                                                          05400005
+      R12=DSQRT(XI(1)*XI(1)+ET(2)*ET(2)+Q*Q)                            05410005
+      R21=DSQRT(XI(2)*XI(2)+ET(1)*ET(1)+Q*Q)                            05420005
+      R22=DSQRT(XI(2)*XI(2)+ET(2)*ET(2)+Q*Q)                            05430005
+      IF(XI(1).LT.F0 .AND. R21+XI(2).LT.EPS) KXI(1)=1                   05440011
+      IF(XI(1).LT.F0 .AND. R22+XI(2).LT.EPS) KXI(2)=1                   05450011
+      IF(ET(1).LT.F0 .AND. R12+ET(2).LT.EPS) KET(1)=1                   05460011
+      IF(ET(1).LT.F0 .AND. R22+ET(2).LT.EPS) KET(2)=1                   05470011
+C=====                                                                  05480015
+      DO 223 K=1,2                                                      05490005
+      DO 222 J=1,2                                                      05500005
+        CALL DCCON2(XI(J),ET(K),Q,SD,CD,KXI(K),KET(J))                  05510014
+        CALL UA(XI(J),ET(K),Q,DD1,DD2,DD3,DUA)                          05520005
+C-----                                                                  05530005
+        DO 220 I=1,10,3                                                 05540005
+          DU(I)  =-DUA(I)                                               05550005
+          DU(I+1)=-DUA(I+1)*CD+DUA(I+2)*SD                              05560005
+          DU(I+2)=-DUA(I+1)*SD-DUA(I+2)*CD                              05570005
+          IF(I.LT.10) GO TO 220                                         05580005
+          DU(I)  =-DU(I)                                                05590005
+          DU(I+1)=-DU(I+1)                                              05600005
+          DU(I+2)=-DU(I+2)                                              05610005
+  220   CONTINUE                                                        05620005
+        DO 221 I=1,12                                                   05630005
+          IF(J+K.NE.3) U(I)=U(I)+DU(I)                                  05640005
+          IF(J+K.EQ.3) U(I)=U(I)-DU(I)                                  05650005
+  221   CONTINUE                                                        05660005
+C-----                                                                  05670005
+  222 CONTINUE                                                          05680005
+  223 CONTINUE                                                          05690005
+C=======================================                                05700005
+C=====  IMAGE-SOURCE CONTRIBUTION  =====                                05710005
+C=======================================                                05720005
+      D=DEPTH-Z                                                         05730005
+      P=Y*CD+D*SD                                                       05740005
+      Q=Y*SD-D*CD                                                       05750005
+      ET(1)=P-AW1                                                       05760005
+      ET(2)=P-AW2                                                       05770005
+      IF(DABS(Q).LT.EPS)  Q=F0                                          05780005
+      IF(DABS(ET(1)).LT.EPS) ET(1)=F0                                   05790005
+      IF(DABS(ET(2)).LT.EPS) ET(2)=F0                                   05800005
+C--------------------------------                                       05810005
+C----- REJECT SINGULAR CASE -----                                       05820005
+C--------------------------------                                       05830005
+C----- ON FAULT EDGE                                                    05840015
+      IF(Q.EQ.F0                                                        05850015
+     *   .AND.(    (XI(1)*XI(2).LE.F0 .AND. ET(1)*ET(2).EQ.F0)          05860015
+     *         .OR.(ET(1)*ET(2).LE.F0 .AND. XI(1)*XI(2).EQ.F0) ))       05870015
+     *   GO TO 99                                                       05880015
+C----- ON NEGATIVE EXTENSION OF FAULT EDGE                              05890015
+      KXI(1)=0                                                          05900005
+      KXI(2)=0                                                          05910005
+      KET(1)=0                                                          05920005
+      KET(2)=0                                                          05930005
+      R12=DSQRT(XI(1)*XI(1)+ET(2)*ET(2)+Q*Q)                            05940005
+      R21=DSQRT(XI(2)*XI(2)+ET(1)*ET(1)+Q*Q)                            05950005
+      R22=DSQRT(XI(2)*XI(2)+ET(2)*ET(2)+Q*Q)                            05960005
+      IF(XI(1).LT.F0 .AND. R21+XI(2).LT.EPS) KXI(1)=1                   05970011
+      IF(XI(1).LT.F0 .AND. R22+XI(2).LT.EPS) KXI(2)=1                   05980011
+      IF(ET(1).LT.F0 .AND. R12+ET(2).LT.EPS) KET(1)=1                   05990011
+      IF(ET(1).LT.F0 .AND. R22+ET(2).LT.EPS) KET(2)=1                   06000011
+C=====                                                                  06010015
+      DO 334 K=1,2                                                      06020005
+      DO 333 J=1,2                                                      06030005
+        CALL DCCON2(XI(J),ET(K),Q,SD,CD,KXI(K),KET(J))                  06040014
+        CALL UA(XI(J),ET(K),Q,DD1,DD2,DD3,DUA)                          06050005
+        CALL UB(XI(J),ET(K),Q,DD1,DD2,DD3,DUB)                          06060005
+        CALL UC(XI(J),ET(K),Q,ZZ,DD1,DD2,DD3,DUC)                       06070005
+C-----                                                                  06080005
+        DO 330 I=1,10,3                                                 06090005
+          DU(I)=DUA(I)+DUB(I)+Z*DUC(I)                                  06100005
+          DU(I+1)=(DUA(I+1)+DUB(I+1)+Z*DUC(I+1))*CD                     06110005
+     *           -(DUA(I+2)+DUB(I+2)+Z*DUC(I+2))*SD                     06120005
+          DU(I+2)=(DUA(I+1)+DUB(I+1)-Z*DUC(I+1))*SD                     06130005
+     *           +(DUA(I+2)+DUB(I+2)-Z*DUC(I+2))*CD                     06140005
+          IF(I.LT.10) GO TO 330                                         06150005
+          DU(10)=DU(10)+DUC(1)                                          06160005
+          DU(11)=DU(11)+DUC(2)*CD-DUC(3)*SD                             06170005
+          DU(12)=DU(12)-DUC(2)*SD-DUC(3)*CD                             06180005
+  330   CONTINUE                                                        06190005
+        DO 331 I=1,12                                                   06200005
+          IF(J+K.NE.3) U(I)=U(I)+DU(I)                                  06210005
+          IF(J+K.EQ.3) U(I)=U(I)-DU(I)                                  06220005
+  331   CONTINUE                                                        06230005
+C-----                                                                  06240005
+  333 CONTINUE                                                          06250005
+  334 CONTINUE                                                          06260005
+C=====                                                                  06270005
+      UX=U(1)                                                           06280005
+      UY=U(2)                                                           06290005
+      UZ=U(3)                                                           06300005
+      UXX=U(4)                                                          06310005
+      UYX=U(5)                                                          06320005
+      UZX=U(6)                                                          06330005
+      UXY=U(7)                                                          06340005
+      UYY=U(8)                                                          06350005
+      UZY=U(9)                                                          06360005
+      UXZ=U(10)                                                         06370005
+      UYZ=U(11)                                                         06380005
+      UZZ=U(12)                                                         06390005
+      IRET=0                                                            06400005
+      RETURN                                                            06410005
+C===========================================                            06420005
+C=====  IN CASE OF SINGULAR (ON EDGE)  =====                            06430005
+C===========================================                            06440005
+   99 UX=F0                                                             06450005
+      UY=F0                                                             06460005
+      UZ=F0                                                             06470005
+      UXX=F0                                                            06480005
+      UYX=F0                                                            06490005
+      UZX=F0                                                            06500005
+      UXY=F0                                                            06510005
+      UYY=F0                                                            06520005
+      UZY=F0                                                            06530005
+      UXZ=F0                                                            06540005
+      UYZ=F0                                                            06550005
+      UZZ=F0                                                            06560005
+      IRET=1                                                            06570005
+      RETURN                                                            06580005
+      END                                                               06590005
+      SUBROUTINE  UA(XI,ET,Q,DISL1,DISL2,DISL3,U)                       06600005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         06610005
+      DIMENSION U(12),DU(12)                                            06620005
+C                                                                       06630005
+C********************************************************************   06640005
+C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-A)             *****   06650005
+C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****   06660005
+C********************************************************************   06670005
+C                                                                       06680005
+C***** INPUT                                                            06690005
+C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM                  06700005
+C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS              06710005
+C***** OUTPUT                                                           06720005
+C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES                     06730005
+C                                                                       06740005
+      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D  06750005
+      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,  06760005
+     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ                                06770005
+      DATA F0,F2,PI2/0.D0,2.D0,6.283185307179586D0/                     06780005
+C-----                                                                  06790005
+      DO 111  I=1,12                                                    06800005
+  111 U(I)=F0                                                           06810005
+      XY=XI*Y11                                                         06820005
+      QX=Q *X11                                                         06830005
+      QY=Q *Y11                                                         06840005
+C======================================                                 06850005
+C=====  STRIKE-SLIP CONTRIBUTION  =====                                 06860005
+C======================================                                 06870005
+      IF(DISL1.NE.F0) THEN                                              06880005
+        DU( 1)=    TT/F2 +ALP2*XI*QY                                    06890005
+        DU( 2)=           ALP2*Q/R                                      06900005
+        DU( 3)= ALP1*ALE -ALP2*Q*QY                                     06910005
+        DU( 4)=-ALP1*QY  -ALP2*XI2*Q*Y32                                06920005
+        DU( 5)=          -ALP2*XI*Q/R3                                  06930005
+        DU( 6)= ALP1*XY  +ALP2*XI*Q2*Y32                                06940005
+        DU( 7)= ALP1*XY*SD        +ALP2*XI*FY+D/F2*X11                  06950005
+        DU( 8)=                    ALP2*EY                              06960005
+        DU( 9)= ALP1*(CD/R+QY*SD) -ALP2*Q*FY                            06970005
+        DU(10)= ALP1*XY*CD        +ALP2*XI*FZ+Y/F2*X11                  06980005
+        DU(11)=                    ALP2*EZ                              06990005
+        DU(12)=-ALP1*(SD/R-QY*CD) -ALP2*Q*FZ                            07000005
+        DO 222 I=1,12                                                   07010005
+  222   U(I)=U(I)+DISL1/PI2*DU(I)                                       07020005
+      ENDIF                                                             07030005
+C======================================                                 07040005
+C=====    DIP-SLIP CONTRIBUTION   =====                                 07050005
+C======================================                                 07060005
+      IF(DISL2.NE.F0) THEN                                              07070005
+        DU( 1)=           ALP2*Q/R                                      07080005
+        DU( 2)=    TT/F2 +ALP2*ET*QX                                    07090005
+        DU( 3)= ALP1*ALX -ALP2*Q*QX                                     07100005
+        DU( 4)=        -ALP2*XI*Q/R3                                    07110005
+        DU( 5)= -QY/F2 -ALP2*ET*Q/R3                                    07120005
+        DU( 6)= ALP1/R +ALP2*Q2/R3                                      07130005
+        DU( 7)=                      ALP2*EY                            07140005
+        DU( 8)= ALP1*D*X11+XY/F2*SD +ALP2*ET*GY                         07150005
+        DU( 9)= ALP1*Y*X11          -ALP2*Q*GY                          07160005
+        DU(10)=                      ALP2*EZ                            07170005
+        DU(11)= ALP1*Y*X11+XY/F2*CD +ALP2*ET*GZ                         07180005
+        DU(12)=-ALP1*D*X11          -ALP2*Q*GZ                          07190005
+        DO 333 I=1,12                                                   07200005
+  333   U(I)=U(I)+DISL2/PI2*DU(I)                                       07210005
+      ENDIF                                                             07220005
+C========================================                               07230005
+C=====  TENSILE-FAULT CONTRIBUTION  =====                               07240005
+C========================================                               07250005
+      IF(DISL3.NE.F0) THEN                                              07260005
+        DU( 1)=-ALP1*ALE -ALP2*Q*QY                                     07270005
+        DU( 2)=-ALP1*ALX -ALP2*Q*QX                                     07280005
+        DU( 3)=    TT/F2 -ALP2*(ET*QX+XI*QY)                            07290005
+        DU( 4)=-ALP1*XY  +ALP2*XI*Q2*Y32                                07300005
+        DU( 5)=-ALP1/R   +ALP2*Q2/R3                                    07310005
+        DU( 6)=-ALP1*QY  -ALP2*Q*Q2*Y32                                 07320005
+        DU( 7)=-ALP1*(CD/R+QY*SD)  -ALP2*Q*FY                           07330005
+        DU( 8)=-ALP1*Y*X11         -ALP2*Q*GY                           07340005
+        DU( 9)= ALP1*(D*X11+XY*SD) +ALP2*Q*HY                           07350005
+        DU(10)= ALP1*(SD/R-QY*CD)  -ALP2*Q*FZ                           07360005
+        DU(11)= ALP1*D*X11         -ALP2*Q*GZ                           07370005
+        DU(12)= ALP1*(Y*X11+XY*CD) +ALP2*Q*HZ                           07380005
+        DO 444 I=1,12                                                   07390005
+  444   U(I)=U(I)+DISL3/PI2*DU(I)                                       07400005
+      ENDIF                                                             07410005
+      RETURN                                                            07420005
+      END                                                               07430005
+      SUBROUTINE  UB(XI,ET,Q,DISL1,DISL2,DISL3,U)                       07440005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         07450005
+      DIMENSION U(12),DU(12)                                            07460005
+C                                                                       07470005
+C********************************************************************   07480005
+C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-B)             *****   07490005
+C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****   07500005
+C********************************************************************   07510005
+C                                                                       07520005
+C***** INPUT                                                            07530005
+C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM                  07540005
+C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS              07550005
+C***** OUTPUT                                                           07560005
+C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES                     07570005
+C                                                                       07580005
+      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D  07590005
+      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,  07600005
+     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ                                07610005
+      DATA  F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/            07620005
+C-----                                                                  07630005
+      RD=R+D                                                            07640005
+      D11=F1/(R*RD)                                                     07650005
+      AJ2=XI*Y/RD*D11                                                   07660005
+      AJ5=-(D+Y*Y/RD)*D11                                               07670005
+      IF(CD.NE.F0) THEN                                                 07680005
+        IF(XI.EQ.F0) THEN                                               07690005
+          AI4=F0                                                        07700005
+        ELSE                                                            07710005
+          X=DSQRT(XI2+Q2)                                               07720005
+          AI4=F1/CDCD*( XI/RD*SDCD                                      07730005
+     *       +F2*DATAN((ET*(X+Q*CD)+X*(R+X)*SD)/(XI*(R+X)*CD)) )        07740005
+        ENDIF                                                           07750005
+        AI3=(Y*CD/RD-ALE+SD*DLOG(RD))/CDCD                              07760005
+        AK1=XI*(D11-Y11*SD)/CD                                          07770005
+        AK3=(Q*Y11-Y*D11)/CD                                            07780005
+        AJ3=(AK1-AJ2*SD)/CD                                             07790005
+        AJ6=(AK3-AJ5*SD)/CD                                             07800005
+      ELSE                                                              07810005
+        RD2=RD*RD                                                       07820005
+        AI3=(ET/RD+Y*Q/RD2-ALE)/F2                                      07830005
+        AI4=XI*Y/RD2/F2                                                 07840005
+        AK1=XI*Q/RD*D11                                                 07850005
+        AK3=SD/RD*(XI2*D11-F1)                                          07860005
+        AJ3=-XI/RD2*(Q2*D11-F1/F2)                                      07870005
+        AJ6=-Y/RD2*(XI2*D11-F1/F2)                                      07880005
+      ENDIF                                                             07890005
+C-----                                                                  07900005
+      XY=XI*Y11                                                         07910005
+      AI1=-XI/RD*CD-AI4*SD                                              07920005
+      AI2= DLOG(RD)+AI3*SD                                              07930005
+      AK2= F1/R+AK3*SD                                                  07940005
+      AK4= XY*CD-AK1*SD                                                 07950005
+      AJ1= AJ5*CD-AJ6*SD                                                07960005
+      AJ4=-XY-AJ2*CD+AJ3*SD                                             07970005
+C=====                                                                  07980005
+      DO 111  I=1,12                                                    07990005
+  111 U(I)=F0                                                           08000005
+      QX=Q*X11                                                          08010005
+      QY=Q*Y11                                                          08020005
+C======================================                                 08030005
+C=====  STRIKE-SLIP CONTRIBUTION  =====                                 08040005
+C======================================                                 08050005
+      IF(DISL1.NE.F0) THEN                                              08060005
+        DU( 1)=-XI*QY-TT -ALP3*AI1*SD                                   08070005
+        DU( 2)=-Q/R      +ALP3*Y/RD*SD                                  08080005
+        DU( 3)= Q*QY     -ALP3*AI2*SD                                   08090005
+        DU( 4)= XI2*Q*Y32 -ALP3*AJ1*SD                                  08100005
+        DU( 5)= XI*Q/R3   -ALP3*AJ2*SD                                  08110005
+        DU( 6)=-XI*Q2*Y32 -ALP3*AJ3*SD                                  08120005
+        DU( 7)=-XI*FY-D*X11 +ALP3*(XY+AJ4)*SD                           08130005
+        DU( 8)=-EY          +ALP3*(F1/R+AJ5)*SD                         08140005
+        DU( 9)= Q*FY        -ALP3*(QY-AJ6)*SD                           08150005
+        DU(10)=-XI*FZ-Y*X11 +ALP3*AK1*SD                                08160005
+        DU(11)=-EZ          +ALP3*Y*D11*SD                              08170005
+        DU(12)= Q*FZ        +ALP3*AK2*SD                                08180005
+        DO 222 I=1,12                                                   08190005
+  222   U(I)=U(I)+DISL1/PI2*DU(I)                                       08200005
+      ENDIF                                                             08210005
+C======================================                                 08220005
+C=====    DIP-SLIP CONTRIBUTION   =====                                 08230005
+C======================================                                 08240005
+      IF(DISL2.NE.F0) THEN                                              08250005
+        DU( 1)=-Q/R      +ALP3*AI3*SDCD                                 08260005
+        DU( 2)=-ET*QX-TT -ALP3*XI/RD*SDCD                               08270005
+        DU( 3)= Q*QX     +ALP3*AI4*SDCD                                 08280005
+        DU( 4)= XI*Q/R3     +ALP3*AJ4*SDCD                              08290005
+        DU( 5)= ET*Q/R3+QY  +ALP3*AJ5*SDCD                              08300005
+        DU( 6)=-Q2/R3       +ALP3*AJ6*SDCD                              08310005
+        DU( 7)=-EY          +ALP3*AJ1*SDCD                              08320005
+        DU( 8)=-ET*GY-XY*SD +ALP3*AJ2*SDCD                              08330005
+        DU( 9)= Q*GY        +ALP3*AJ3*SDCD                              08340005
+        DU(10)=-EZ          -ALP3*AK3*SDCD                              08350005
+        DU(11)=-ET*GZ-XY*CD -ALP3*XI*D11*SDCD                           08360005
+        DU(12)= Q*GZ        -ALP3*AK4*SDCD                              08370005
+        DO 333 I=1,12                                                   08380005
+  333   U(I)=U(I)+DISL2/PI2*DU(I)                                       08390005
+      ENDIF                                                             08400005
+C========================================                               08410005
+C=====  TENSILE-FAULT CONTRIBUTION  =====                               08420005
+C========================================                               08430005
+      IF(DISL3.NE.F0) THEN                                              08440005
+        DU( 1)= Q*QY           -ALP3*AI3*SDSD                           08450005
+        DU( 2)= Q*QX           +ALP3*XI/RD*SDSD                         08460005
+        DU( 3)= ET*QX+XI*QY-TT -ALP3*AI4*SDSD                           08470005
+        DU( 4)=-XI*Q2*Y32 -ALP3*AJ4*SDSD                                08480005
+        DU( 5)=-Q2/R3     -ALP3*AJ5*SDSD                                08490005
+        DU( 6)= Q*Q2*Y32  -ALP3*AJ6*SDSD                                08500005
+        DU( 7)= Q*FY -ALP3*AJ1*SDSD                                     08510005
+        DU( 8)= Q*GY -ALP3*AJ2*SDSD                                     08520005
+        DU( 9)=-Q*HY -ALP3*AJ3*SDSD                                     08530005
+        DU(10)= Q*FZ +ALP3*AK3*SDSD                                     08540005
+        DU(11)= Q*GZ +ALP3*XI*D11*SDSD                                  08550005
+        DU(12)=-Q*HZ +ALP3*AK4*SDSD                                     08560005
+        DO 444 I=1,12                                                   08570005
+  444   U(I)=U(I)+DISL3/PI2*DU(I)                                       08580005
+      ENDIF                                                             08590005
+      RETURN                                                            08600005
+      END                                                               08610005
+      SUBROUTINE  UC(XI,ET,Q,Z,DISL1,DISL2,DISL3,U)                     08620005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         08630005
+      DIMENSION U(12),DU(12)                                            08640005
+C                                                                       08650005
+C********************************************************************   08660005
+C*****    DISPLACEMENT AND STRAIN AT DEPTH (PART-C)             *****   08670005
+C*****    DUE TO BURIED FINITE FAULT IN A SEMIINFINITE MEDIUM   *****   08680005
+C********************************************************************   08690005
+C                                                                       08700005
+C***** INPUT                                                            08710005
+C*****   XI,ET,Q,Z   : STATION COORDINATES IN FAULT SYSTEM              08720005
+C*****   DISL1-DISL3 : STRIKE-, DIP-, TENSILE-DISLOCATIONS              08730005
+C***** OUTPUT                                                           08740005
+C*****   U(12) : DISPLACEMENT AND THEIR DERIVATIVES                     08750005
+C                                                                       08760005
+      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D  08770005
+      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,  08780005
+     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ                                08790005
+      DATA F0,F1,F2,F3,PI2/0.D0,1.D0,2.D0,3.D0,6.283185307179586D0/     08800005
+C-----                                                                  08810005
+      C=D+Z                                                             08820005
+      X53=(8.D0*R2+9.D0*R*XI+F3*XI2)*X11*X11*X11/R2                     08830005
+      Y53=(8.D0*R2+9.D0*R*ET+F3*ET2)*Y11*Y11*Y11/R2                     08840005
+      H=Q*CD-Z                                                          08850005
+      Z32=SD/R3-H*Y32                                                   08860005
+      Z53=F3*SD/R5-H*Y53                                                08870005
+      Y0=Y11-XI2*Y32                                                    08880005
+      Z0=Z32-XI2*Z53                                                    08890005
+      PPY=CD/R3+Q*Y32*SD                                                08900005
+      PPZ=SD/R3-Q*Y32*CD                                                08910005
+      QQ=Z*Y32+Z32+Z0                                                   08920005
+      QQY=F3*C*D/R5-QQ*SD                                               08930005
+      QQZ=F3*C*Y/R5-QQ*CD+Q*Y32                                         08940005
+      XY=XI*Y11                                                         08950005
+      QX=Q*X11                                                          08960005
+      QY=Q*Y11                                                          08970005
+      QR=F3*Q/R5                                                        08980005
+      CQX=C*Q*X53                                                       08990005
+      CDR=(C+D)/R3                                                      09000005
+      YY0=Y/R3-Y0*CD                                                    09010005
+C=====                                                                  09020005
+      DO 111  I=1,12                                                    09030005
+  111 U(I)=F0                                                           09040005
+C======================================                                 09050005
+C=====  STRIKE-SLIP CONTRIBUTION  =====                                 09060005
+C======================================                                 09070005
+      IF(DISL1.NE.F0) THEN                                              09080005
+        DU( 1)= ALP4*XY*CD           -ALP5*XI*Q*Z32                     09090005
+        DU( 2)= ALP4*(CD/R+F2*QY*SD) -ALP5*C*Q/R3                       09100005
+        DU( 3)= ALP4*QY*CD           -ALP5*(C*ET/R3-Z*Y11+XI2*Z32)      09110005
+        DU( 4)= ALP4*Y0*CD                  -ALP5*Q*Z0                  09120005
+        DU( 5)=-ALP4*XI*(CD/R3+F2*Q*Y32*SD) +ALP5*C*XI*QR               09130005
+        DU( 6)=-ALP4*XI*Q*Y32*CD            +ALP5*XI*(F3*C*ET/R5-QQ)    09140005
+        DU( 7)=-ALP4*XI*PPY*CD    -ALP5*XI*QQY                          09150005
+        DU( 8)= ALP4*F2*(D/R3-Y0*SD)*SD-Y/R3*CD                         09160005
+     *                            -ALP5*(CDR*SD-ET/R3-C*Y*QR)           09170005
+        DU( 9)=-ALP4*Q/R3+YY0*SD  +ALP5*(CDR*CD+C*D*QR-(Y0*CD+Q*Z0)*SD) 09180005
+        DU(10)= ALP4*XI*PPZ*CD    -ALP5*XI*QQZ                          09190005
+        DU(11)= ALP4*F2*(Y/R3-Y0*CD)*SD+D/R3*CD -ALP5*(CDR*CD+C*D*QR)   09200005
+        DU(12)=         YY0*CD    -ALP5*(CDR*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD) 09210005
+        DO 222 I=1,12                                                   09220005
+  222   U(I)=U(I)+DISL1/PI2*DU(I)                                       09230005
+      ENDIF                                                             09240005
+C======================================                                 09250005
+C=====    DIP-SLIP CONTRIBUTION   =====                                 09260005
+C======================================                                 09270005
+      IF(DISL2.NE.F0) THEN                                              09280005
+        DU( 1)= ALP4*CD/R -QY*SD -ALP5*C*Q/R3                           09290005
+        DU( 2)= ALP4*Y*X11       -ALP5*C*ET*Q*X32                       09300005
+        DU( 3)=     -D*X11-XY*SD -ALP5*C*(X11-Q2*X32)                   09310005
+        DU( 4)=-ALP4*XI/R3*CD +ALP5*C*XI*QR +XI*Q*Y32*SD                09320005
+        DU( 5)=-ALP4*Y/R3     +ALP5*C*ET*QR                             09330005
+        DU( 6)=    D/R3-Y0*SD +ALP5*C/R3*(F1-F3*Q2/R2)                  09340005
+        DU( 7)=-ALP4*ET/R3+Y0*SDSD -ALP5*(CDR*SD-C*Y*QR)                09350005
+        DU( 8)= ALP4*(X11-Y*Y*X32) -ALP5*C*((D+F2*Q*CD)*X32-Y*ET*Q*X53) 09360005
+        DU( 9)=  XI*PPY*SD+Y*D*X32 +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53)   09370005
+        DU(10)=      -Q/R3+Y0*SDCD -ALP5*(CDR*CD+C*D*QR)                09380005
+        DU(11)= ALP4*Y*D*X32       -ALP5*C*((Y-F2*Q*SD)*X32+D*ET*Q*X53) 09390005
+        DU(12)=-XI*PPZ*SD+X11-D*D*X32-ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53) 09400005
+        DO 333 I=1,12                                                   09410005
+  333   U(I)=U(I)+DISL2/PI2*DU(I)                                       09420005
+      ENDIF                                                             09430005
+C========================================                               09440005
+C=====  TENSILE-FAULT CONTRIBUTION  =====                               09450005
+C========================================                               09460005
+      IF(DISL3.NE.F0) THEN                                              09470005
+        DU( 1)=-ALP4*(SD/R+QY*CD)   -ALP5*(Z*Y11-Q2*Z32)                09480005
+        DU( 2)= ALP4*F2*XY*SD+D*X11 -ALP5*C*(X11-Q2*X32)                09490005
+        DU( 3)= ALP4*(Y*X11+XY*CD)  +ALP5*Q*(C*ET*X32+XI*Z32)           09500005
+        DU( 4)= ALP4*XI/R3*SD+XI*Q*Y32*CD+ALP5*XI*(F3*C*ET/R5-F2*Z32-Z0)09510005
+        DU( 5)= ALP4*F2*Y0*SD-D/R3 +ALP5*C/R3*(F1-F3*Q2/R2)             09520005
+        DU( 6)=-ALP4*YY0           -ALP5*(C*ET*QR-Q*Z0)                 09530005
+        DU( 7)= ALP4*(Q/R3+Y0*SDCD)   +ALP5*(Z/R3*CD+C*D*QR-Q*Z0*SD)    09540005
+        DU( 8)=-ALP4*F2*XI*PPY*SD-Y*D*X32                               09550005
+     *                    +ALP5*C*((Y+F2*Q*SD)*X32-Y*Q2*X53)            09560005
+        DU( 9)=-ALP4*(XI*PPY*CD-X11+Y*Y*X32)                            09570005
+     *                    +ALP5*(C*((D+F2*Q*CD)*X32-Y*ET*Q*X53)+XI*QQY) 09580005
+        DU(10)=  -ET/R3+Y0*CDCD -ALP5*(Z/R3*SD-C*Y*QR-Y0*SDSD+Q*Z0*CD)  09590005
+        DU(11)= ALP4*F2*XI*PPZ*SD-X11+D*D*X32                           09600005
+     *                    -ALP5*C*((D-F2*Q*CD)*X32-D*Q2*X53)            09610005
+        DU(12)= ALP4*(XI*PPZ*CD+Y*D*X32)                                09620005
+     *                    +ALP5*(C*((Y-F2*Q*SD)*X32+D*ET*Q*X53)+XI*QQZ) 09630005
+        DO 444 I=1,12                                                   09640005
+  444   U(I)=U(I)+DISL3/PI2*DU(I)                                       09650005
+      ENDIF                                                             09660005
+      RETURN                                                            09670005
+      END                                                               09680005
+      SUBROUTINE  DCCON0(ALPHA,DIP)                                     09690005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         09700005
+C                                                                       09710005
+C*******************************************************************    09720005
+C*****   CALCULATE MEDIUM CONSTANTS AND FAULT-DIP CONSTANTS    *****    09730005
+C*******************************************************************    09740005
+C                                                                       09750005
+C***** INPUT                                                            09760005
+C*****   ALPHA : MEDIUM CONSTANT  (LAMBDA+MYU)/(LAMBDA+2*MYU)           09770005
+C*****   DIP   : DIP-ANGLE (DEGREE)                                     09780005
+C### CAUTION ### IF COS(DIP) IS SUFFICIENTLY SMALL, IT IS SET TO ZERO   09790005
+C                                                                       09800005
+      COMMON /C0/ALP1,ALP2,ALP3,ALP4,ALP5,SD,CD,SDSD,CDCD,SDCD,S2D,C2D  09810005
+      DATA F0,F1,F2,PI2/0.D0,1.D0,2.D0,6.283185307179586D0/             09820005
+      DATA EPS/1.D-6/                                                   09830005
+C-----                                                                  09840005
+      ALP1=(F1-ALPHA)/F2                                                09850005
+      ALP2= ALPHA/F2                                                    09860005
+      ALP3=(F1-ALPHA)/ALPHA                                             09870005
+      ALP4= F1-ALPHA                                                    09880005
+      ALP5= ALPHA                                                       09890005
+C-----                                                                  09900005
+      P18=PI2/360.D0                                                    09910005
+      SD=DSIN(DIP*P18)                                                  09920005
+      CD=DCOS(DIP*P18)                                                  09930005
+      IF(DABS(CD).LT.EPS) THEN                                          09940005
+        CD=F0                                                           09950005
+        IF(SD.GT.F0) SD= F1                                             09960005
+        IF(SD.LT.F0) SD=-F1                                             09970005
+      ENDIF                                                             09980005
+      SDSD=SD*SD                                                        09990005
+      CDCD=CD*CD                                                        10000005
+      SDCD=SD*CD                                                        10010005
+      S2D=F2*SDCD                                                       10020005
+      C2D=CDCD-SDSD                                                     10030005
+      RETURN                                                            10040005
+      END                                                               10050005
+      SUBROUTINE  DCCON1(X,Y,D)                                         10060005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         10070005
+C                                                                       10080005
+C********************************************************************** 10090005
+C*****   CALCULATE STATION GEOMETRY CONSTANTS FOR POINT SOURCE    ***** 10100005
+C********************************************************************** 10110005
+C                                                                       10120005
+C***** INPUT                                                            10130005
+C*****   X,Y,D : STATION COORDINATES IN FAULT SYSTEM                    10140005
+C### CAUTION ### IF X,Y,D ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZERO  10150005
+C                                                                       10160005
+      COMMON /C0/DUMMY(5),SD,CD,dumm(5)                                 10170005
+      COMMON /C1/P,Q,S,T,XY,X2,Y2,D2,R,R2,R3,R5,QR,QRX,A3,A5,B3,C3,     10180005
+     *           UY,VY,WY,UZ,VZ,WZ                                      10190005
+      DATA  F0,F1,F3,F5,EPS/0.D0,1.D0,3.D0,5.D0,1.D-6/                  10200005
+C-----                                                                  10210005
+      IF(DABS(X).LT.EPS) X=F0                                           10220005
+      IF(DABS(Y).LT.EPS) Y=F0                                           10230005
+      IF(DABS(D).LT.EPS) D=F0                                           10240005
+      P=Y*CD+D*SD                                                       10250005
+      Q=Y*SD-D*CD                                                       10260005
+      S=P*SD+Q*CD                                                       10270005
+      T=P*CD-Q*SD                                                       10280005
+      XY=X*Y                                                            10290005
+      X2=X*X                                                            10300005
+      Y2=Y*Y                                                            10310005
+      D2=D*D                                                            10320005
+      R2=X2+Y2+D2                                                       10330005
+      R =DSQRT(R2)                                                      10340005
+      IF(R.EQ.F0) RETURN                                                10350005
+      R3=R *R2                                                          10360005
+      R5=R3*R2                                                          10370005
+      R7=R5*R2                                                          10380005
+C-----                                                                  10390005
+      A3=F1-F3*X2/R2                                                    10400005
+      A5=F1-F5*X2/R2                                                    10410005
+      B3=F1-F3*Y2/R2                                                    10420005
+      C3=F1-F3*D2/R2                                                    10430005
+C-----                                                                  10440005
+      QR=F3*Q/R5                                                        10450005
+      QRX=F5*QR*X/R2                                                    10460005
+C-----                                                                  10470005
+      UY=SD-F5*Y*Q/R2                                                   10480005
+      UZ=CD+F5*D*Q/R2                                                   10490005
+      VY=S -F5*Y*P*Q/R2                                                 10500005
+      VZ=T +F5*D*P*Q/R2                                                 10510005
+      WY=UY+SD                                                          10520005
+      WZ=UZ+CD                                                          10530005
+      RETURN                                                            10540005
+      END                                                               10550005
+      SUBROUTINE  DCCON2(XI,ET,Q,SD,CD,KXI,KET)                         10560005
+      IMPLICIT REAL*8 (A-H,O-Z)                                         10570005
+C                                                                       10580005
+C********************************************************************** 10590005
+C*****   CALCULATE STATION GEOMETRY CONSTANTS FOR FINITE SOURCE   ***** 10600005
+C********************************************************************** 10610005
+C                                                                       10620005
+C***** INPUT                                                            10630005
+C*****   XI,ET,Q : STATION COORDINATES IN FAULT SYSTEM                  10640005
+C*****   SD,CD   : SIN, COS OF DIP-ANGLE                                10650005
+C*****   KXI,KET : KXI=1, KET=1 MEANS R+XI<EPS, R+ET<EPS, RESPECTIVELY  10660005
+C                                                                       10670005
+C### CAUTION ### IF XI,ET,Q ARE SUFFICIENTLY SMALL, THEY ARE SET TO ZER010680005
+C                                                                       10690005
+      COMMON /C2/XI2,ET2,Q2,R,R2,R3,R5,Y,D,TT,ALX,ALE,X11,Y11,X32,Y32,  10700005
+     *           EY,EZ,FY,FZ,GY,GZ,HY,HZ                                10710005
+      DATA  F0,F1,F2,EPS/0.D0,1.D0,2.D0,1.D-6/                          10720005
+C-----                                                                  10730005
+      IF(DABS(XI).LT.EPS) XI=F0                                         10740005
+      IF(DABS(ET).LT.EPS) ET=F0                                         10750005
+      IF(DABS( Q).LT.EPS)  Q=F0                                         10760005
+      XI2=XI*XI                                                         10770005
+      ET2=ET*ET                                                         10780005
+      Q2=Q*Q                                                            10790005
+      R2=XI2+ET2+Q2                                                     10800005
+      R =DSQRT(R2)                                                      10810005
+      IF(R.EQ.F0) RETURN                                                10820005
+      R3=R *R2                                                          10830005
+      R5=R3*R2                                                          10840005
+      Y =ET*CD+Q*SD                                                     10850005
+      D =ET*SD-Q*CD                                                     10860005
+C-----                                                                  10870005
+      IF(Q.EQ.F0) THEN                                                  10880005
+        TT=F0                                                           10890005
+      ELSE                                                              10900005
+        TT=DATAN(XI*ET/(Q*R))                                           10910005
+      ENDIF                                                             10920005
+C-----                                                                  10930005
+      IF(KXI.EQ.1) THEN                                                 10940005
+        ALX=-DLOG(R-XI)                                                 10950005
+        X11=F0                                                          10960005
+        X32=F0                                                          10970005
+      ELSE                                                              10980005
+        RXI=R+XI                                                        10990005
+        ALX=DLOG(RXI)                                                   11000005
+        X11=F1/(R*RXI)                                                  11010005
+        X32=(R+RXI)*X11*X11/R                                           11020005
+      ENDIF                                                             11030005
+C-----                                                                  11040005
+      IF(KET.EQ.1) THEN                                                 11050005
+        ALE=-DLOG(R-ET)                                                 11060005
+        Y11=F0                                                          11070005
+        Y32=F0                                                          11080005
+      ELSE                                                              11090005
+        RET=R+ET                                                        11100005
+        ALE=DLOG(RET)                                                   11110005
+        Y11=F1/(R*RET)                                                  11120005
+        Y32=(R+RET)*Y11*Y11/R                                           11130005
+      ENDIF                                                             11140005
+C-----                                                                  11150005
+      EY=SD/R-Y*Q/R3                                                    11160005
+      EZ=CD/R+D*Q/R3                                                    11170005
+      FY=D/R3+XI2*Y32*SD                                                11180005
+      FZ=Y/R3+XI2*Y32*CD                                                11190005
+      GY=F2*X11*SD-Y*Q*X32                                              11200005
+      GZ=F2*X11*CD+D*Q*X32                                              11210005
+      HY=D*Q*X32+XI*Q*Y32*SD                                            11220005
+      HZ=Y*Q*X32+XI*Q*Y32*CD                                            11230005
+      RETURN                                                            11240005
+      END                                                               11250005
--- a/src/pscmp/disazi.f
+++ b/src/pscmp/disazi.f
@ -0,0 +1,74 @@
+	subroutine disazi(rearth,lateq,loneq,latst,lonst,xnorth,yeast)
+	implicit none
+c
+	double precision rearth,lateq,loneq,latst,lonst,xnorth,yeast
+c
+	integer iangle
+	double precision latb,lonb,latc,lonc,angleb,anglec
+	double precision dis,a,b,c,s,aa
+c
+	double precision PI,PI2,DEGTORAD
+	data PI,PI2/3.14159265358979d0,6.28318530717959d0/
+	data DEGTORAD/1.745329252d-02/
+c
+c	A is north pole
+c
+	latb=lateq*DEGTORAD
+	lonb=loneq*DEGTORAD
+	latc=latst*DEGTORAD
+	lonc=lonst*DEGTORAD
+c        print *,(latc-latb)*rearth,(lonc-lonb)*rearth*dcos(lonb)
+c
+	if(lonb.lt.0.d0)lonb=lonb+PI2
+	if(lonc.lt.0.d0)lonc=lonc+PI2
+c
+	b=0.5d0*PI-latb
+	c=0.5d0*PI-latc
+	if(lonc.gt.lonb)then
+	  aa=lonc-lonb
+	  if(aa.le.PI)then
+	    iangle=1
+	  else
+	    aa=PI2-aa
+	    iangle=-1
+	  endif
+	else
+	  aa=lonb-lonc
+	  if(aa.le.PI)then
+	    iangle=-1
+	  else
+	    aa=PI2-aa
+	    iangle=1
+	  endif
+	endif
+	s=dcos(b)*dcos(c)+dsin(b)*dsin(c)*dcos(aa)
+	a=dacos(dsign(dmin1(dabs(s),1.d0),s))
+	dis=a*rearth
+	if(a*b*c.eq.0.d0)then
+	  angleb=0.d0
+	  anglec=0.d0
+	else
+	  s=0.5d0*(a+b+c)
+	  a=dmin1(a,s)
+	  b=dmin1(b,s)
+	  c=dmin1(c,s)
+	  anglec=2.d0*dasin(dmin1(1.d0,
+     &            dsqrt(dsin(s-a)*dsin(s-b)/(dsin(a)*dsin(b)))))
+	  angleb=2.d0*dasin(dmin1(1.d0,
+     &            dsqrt(dsin(s-a)*dsin(s-c)/(dsin(a)*dsin(c)))))
+	  if(iangle.eq.1)then
+	    angleb=PI2-angleb
+	  else
+	    anglec=PI2-anglec
+	  endif
+	endif
+c
+c	cartesian coordinates of the sation
+c
+	xnorth=dis*dcos(anglec)
+	yeast=dis*dsin(anglec)
+c        print *,xnorth,yeast
+c        pause
+c
+	return
+	end
--- a/src/pscmp/getdata.f
+++ b/src/pscmp/getdata.f
@ -0,0 +1,30 @@
+	subroutine getdata(unit,line)
+	implicit none
+c
+c	First implemented in Potsdam, Feb, 1999
+c	Last modified: Potsdam, Nov, 2001, by R. Wang
+c
+	integer unit
+	character line*180,char*1
+c
+	integer i