ids: add all functions, make files, Readme and examples

.gitignore

@@ -1,131 +1,15 @@
-# ---> Python
-# Byte-compiled / optimized / DLL files
-__pycache__/
-*.py[cod]
-*$py.class
-
-# C extensions
-*.so
-
-# Distribution / packaging
-.Python
-build/
-develop-eggs/
-dist/
-downloads/
-eggs/
-.eggs/
-lib/
-lib64/
-parts/
-sdist/
-var/
-wheels/
-pip-wheel-metadata/
-share/python-wheels/
-*.egg-info/
-.installed.cfg
-*.egg
-MANIFEST
-
-# PyInstaller
-#  Usually these files are written by a python script from a template
-#  before PyInstaller builds the exe, so as to inject date/other infos into it.
-*.manifest
-*.spec
-
-# Installer logs
-pip-log.txt
-pip-delete-this-directory.txt
-
-# Unit test / coverage reports
-htmlcov/
-.tox/
-.nox/
-.coverage
-.coverage.*
-.cache
-nosetests.xml
-coverage.xml
-*.cover
-*.py,cover
-.hypothesis/
-.pytest_cache/
-
-# Translations
-*.mo
-*.pot
-
-# Django stuff:
-*.log
-local_settings.py
-db.sqlite3
-db.sqlite3-journal
-
-# Flask stuff:
-instance/
-.webassets-cache
-
-# Scrapy stuff:
-.scrapy
-
-# Sphinx documentation
-docs/_build/
-
-# PyBuilder
-target/
-
-# Jupyter Notebook
-.ipynb_checkpoints
-
-# IPython
-profile_default/
-ipython_config.py
-
-# pyenv
-.python-version
-
-# pipenv
-#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
-#   However, in case of collaboration, if having platform-specific dependencies or dependencies
-#   having no cross-platform support, pipenv may install dependencies that don't work, or not
-#   install all needed dependencies.
-#Pipfile.lock
-
-# PEP 582; used by e.g. github.com/David-OConnor/pyflow
-__pypackages__/
-
-# Celery stuff
-celerybeat-schedule
-celerybeat.pid
-
-# SageMath parsed files
-*.sage.py
-
-# Environments
-.env
-.venv
-env/
-venv/
-ENV/
-env.bak/
-venv.bak/
-
-# Spyder project settings
-.spyderproject
-.spyproject
-
-# Rope project settings
-.ropeproject
-
-# mkdocs documentation
-/site
-
-# mypy
-.mypy_cache/
-.dmypy.json
-dmypy.json
-
-# Pyre type checker
-.pyre/
-
+/aclocal.m4
+/autom4te.cache
+/configure
+/install-sh
+/missing
+/py-compile
+/Makefile.in
+/src/Makefile.in
+/Makefile
+/src/Makefile
+/config.log
+/config.status
+*.o
+*.mod
+/src/ids2020

Makefile.am

@@ -0,0 +1 @@
+SUBDIRS = src

README.md

@@ -1,7 +1,33 @@
-# ids
-
+# IDS2020 (packaged as pyrocko backend)
+
 Code to do kinematic earthquake source imaging based on the Iterative Deconvolution and Stacking (IDS) method using seismic waveform and static displacement jointly. 
 
 IDS has been written by Rongjiang Wang.
 
-Packaging has been done by Malte Metz.
+Packaging has been done by Malte Metz.
+
+
+## References
+
+- Zhang, Y., Wang, R., Chen, Y., Xu, L., Du, F., Jin, M., Tu, H. and Dahm, T. (2014), Kinematic
+  Rupture Model and Hypocenter Relocation of the 2013 Mw 6.6 Lushan Earthquake Constrained by
+  Strong‐Motion and Teleseismic Data, Seismological Research Letters, 85 (1), 15–22,
+  doi: https://doi.org/10.1785/0220130126.
+
+- Zhang, Y., Wang, R., and Chen, Y.‐T. (2015), Stability of rapid finite‐fault inversion for the
+  2014 Mw6.1 South Napa earthquake, Geophys. Res. Lett., 42, 10, 263–10, 272,
+  doi: https://doi.org/10.1002/2015GL066244. 
+
+- Diao, F. , Wang, R. , Aochi, H. , Walter, T.R., Zhang, Y. , Zheng, Y. and Xiong, X. (2016),
+  Rapid kinematic finite-fault inversion for an Mw 7+ scenario earthquake in the Marmara Sea: an
+  uncertainty study, Geophysical Journal International, 204, 2, 813–824,
+  https://doi.org/10.1093/gji/ggv459
+
+## Compile and install
+
+```
+autoreconf -i   # only if 'configure' script is missing
+./configure
+make
+sudo make install
+```

configure.ac

@@ -0,0 +1,8 @@
+AC_INIT([ids], [2020])
+AM_INIT_AUTOMAKE([-Wall -Werror foreign])
+AC_PROG_F77
+AC_CONFIG_FILES([
+                 Makefile
+                 src/Makefile
+                 ])
+AC_OUTPUT

examples/ids2020-illapel.inp

@@ -0,0 +1,235 @@
+# Example for the input file of FORTRAN77 program "ids2020" for kinematic earthquake
+# source imaging based on the Iterative Deconvolution and Stacking (IDS) Method using
+# seismic waveform and static displacement data jointly.
+#
+# written by Rongjiang Wang
+# Helmholtz Centre Potsdam
+# GFZ German Research Centre for Geosciences
+# e-mail: wang@gfz-potsdam.de
+# Phone +49 331 2881209
+# Fax +49 331 2881204
+#
+# Last modified: Potsdam, April, 2020
+#
+# Reference:
+# Zhang and Wang et al. (2014). Automatic imaging of earthquake rupture processes
+# by iterative deconvolution and stacking of high-rate GNSS and strong motion
+# seismograms, JGR Solid Earth, 199(7), 5633-5650.
+#
+#################################################################
+##                                                             ##
+## If not specified otherwise, SI Unit System is used overall! ##
+##                                                             ##
+##                                                             ##
+## Each comment line should start with "#"                     ##
+## and no blank line is allowed!                               ##
+##                                                             ##
+#################################################################
+#
+#================================================================================
+# EARTHQUAKE LOCATION PARAMETERS
+# ------------------------------
+# 1. origin time [s] (year, month, day, hour, minute [positive integer number]
+#    and second [positive or negative real number]), which will be also used as
+#    the reference time of the output data
+# 2. hypocentre (lat[deg], lon[deg], depth[km])
+# 3. preliminary estimate of moment magnitude
+#    Note: the hypocentre location and origin time (at least to minute) should be
+#          consistent with those given in the "StationInfo.dat" file of the
+#          waveform data (s. below)
+#--------------------------------------------------------------------------------
+  2015 9 16 22 54 32.0
+  -31.573  -72.674   22.40
+  8.3
+#================================================================================
+# WAVEFORM DATA (HIGH-RATE-GNSS/STRONG-MOTION SEISMOGRAMS)
+# --------------------------------------------------------------------
+# 1. folder name of the data
+#    Note: this folder should include a file 'StationInfo.dat' giving all
+#          necessary information about the local monitoring network.
+# 2. lower and upper thresholds of the epicentral distances [km] (stations
+#    outside the threshold distances will not be used)
+# 3. weight of the static offset data (obtained automatically from the velocity
+#    waveforms after an empirical baseline correction) relative to the waveform
+#    data (suggested value 0.0-1.0)
+#--------------------------------------------------------------------------------
+   './WaveformData/'
+  0.0  8000.0
+  0.0
+#================================================================================
+# GREEN'S FUNCTION DATABASE
+# -------------------------
+# 1. Green's function database (folder name)
+# 2. Parameters of Butterworth bandpass filter applied to the synthetics:
+#    order, lower and upper cutoff frequencies [Hz]
+#    computation [Hz]
+#-------------------------------------------------------------------------------- 
+ '../GreenChileSM/GreenFunctions/'
+  3  0.000  0.50
+#================================================================================
+# GNSS STATIC OFFSETS
+# -------------------
+# 1. number of GNSS sites, file name of the data
+#    Note: if the number of GNSS sites <= 0, then no GNSS data are available.
+#    the GNSS data file includes one header line and 8 culomns,
+#    Lat(deg)  Lon(deg)  East(m)  North(m)  Up(m)  wf_E  wf_N   wf_Z
+#    where wf_E/N/Z are within-dataset relative weighting factor proportional to
+#    inverse error variance
+# 2. weight of the whole GNSS dataset relative to the waveform dataset (suggested
+#    value 0.0-1.0)
+#--------------------------------------------------------------------------------
+ 0  './GPSData/GPS_Pisagua.dat'
+ 0.0
+#================================================================================
+# InSAR LOS DISPLACEMENTS
+# -----------------------
+# 1. number of InSAR grids, file name of the data
+#    Note: if the number of InSAR grids <= 0, then no InSAR data are used.
+#          the InSAR data file includes one header line and 8 culomns,
+#          Lat(deg)  Lon(deg) LOS(m)    wf_los   incidence(deg)   azimuth(deg)
+#          where wf_los is the within-dataset relative weighting factor
+#          proportional to inverse error variance
+# 2. weight of the whole InSAR dataset in the misfit variance relative to the
+#    waveform dataset (suggested value 0.0-1.0)
+#--------------------------------------------------------------------------------
+ 0  './InSAR/NoInSAR.dat'
+ 0.00
+#================================================================================
+# Geodetic GREEN'S FUNCTION DATABASE
+# ----------------------------------
+# 1. geodetic Green's function database (folder name)
+# 2. 3 file names of the geodetic Green's functions
+#-------------------------------------------------------------------------------- 
+ '../GreenChileGD/DGreen/'
+ 'uz'  'ur'  'ut'
+#================================================================================
+# FAULT PARAMETERS
+# ----------------
+# 1. ipatch: switch (0/1) for subfaults input form,
+#    0 = read discrete patches from files,
+#    1 = automatic discretization of a single rectangular fault
+#-------------------------------------------------------------------------------- 
+#   idisc
+#-------------------------------------------------------------------------------- 
+    1
+#-------------------------------------------------------------------------------- 
+# IF(idisc = 0)THEN
+#
+#   subfault paramters
+#   ==================
+#   file_name:         data file for discrete patches
+#   number of sub-faults
+#   iref               selection of the patch reference location
+#
+#   Note: the file includes one header line followed by 8 columns of data:
+#
+#   latitude[deg], longitude[deg], depth[km], length[km], width[km],
+#   strike[deg], dip[deg], rake[deg]
+#
+#   where (lat, lon, dep) is the reference location of the patch depending
+#   on iref: 1 = upper-left corner, 2 = upper-right corner,
+#            3 = lower-left corner, 4 = lower-right corner,
+#            5 = central point (s. Fig.)
+#
+#        definitions for a rectangular fault patch
+#        =========================================
+#
+#                 north(x)
+#                  /
+#                 / | strike
+#                1-----------------------2
+#                |\        p .            \ W
+#                :-\      i .              \ i
+#                |  \    l .     5          \ d
+#                :90 \  S .                  \ t
+#                |-dip\  .                    \ h
+#                :     \. ) rake               \ 
+#          downward(z)  3-----------------------4
+#                         L   e   n   g   t   h
+#--------------------------------------------------------------------------------
+#   file_name               nsf   iref
+#--------------------------------------------------------------------------------
+#   './fault_patches.dat'    288   5
+#--------------------------------------------------------------------------------
+# ELSE IF(idisc = 1)THEN
+#
+# 2. focal mechanism: strike, dip and rake angles [deg]
+#--------------------------------------------------------------------------------
+  7.0     19.0    109.0
+#================================================================================
+# IDS PARAMETERS
+# --------------
+# 1. maximum numbers of iterations.
+#    Note: if max. iteration = 0, forward modelling for any dummy data set is
+#          made using the given rupture model and the filter parameters
+#          if max. iteration < 0, same as iteration = 0, but without filtering.
+#--------------------------------------------------------------------------------
+	100
+#================================================================================
+# OUTPUT
+# ------
+# 1. folder name for the output files (should exist already)
+# 2. sub-folder name for comparisons of observed and modelled waveform data
+#    (should exist already)
+# 3. sub-folder name for results from empirical baseline correction of
+#    waveform data (this folder should exist already)
+# 4. file name of major earthquake parameters
+# 5. file name of earthquake source-time-function (moment rate history)
+# 6. file name of sub-fault source-time-function (moment rate history)
+# 7. file name of rupture propagation
+# 8. file name of final slip model
+#    Note: if forward modelling is selected (s. above), the above two files
+#          should provide data for the given kinematic rupture model and will
+#          be overwritten after the IDS inversion
+# 9. file name of normalized residual waveform variance
+#10. file name of predicted static offsets at seismic waveform data sites
+#11. file name of comparison between observed and modelled GNSS data
+#12. file name of comparison between observed and modelled InSAR data
+#13. number of snapshots of the rupture process
+#14. file name of the 1. snapshot for the slip released within the given
+#    time window (t1[s], t2[s])
+#--------------------------------------------------------------------------------
+   './Output'
+   'O2S'
+   'EBC'
+   'EQ_MP.dat'
+   'EQ_STF.dat'
+   'PT_STF.dat'
+   'Rup_Front.dat'
+   'SlipModel_Tp_Mean.dat'
+   'WFMisfitvar.dat'
+   'SMStaticOffset.dat'
+   'GNSSDatafit.dat'
+   'InSARDatafit.dat'
+   30
+   'Snapshot_010s.dat'   0.0  10.0
+   'Snapshot_015s.dat'   0.0  15.0
+   'Snapshot_020s.dat'   0.0  20.0
+   'Snapshot_025s.dat'   0.0  25.0
+   'Snapshot_030s.dat'   0.0  30.0
+   'Snapshot_035s.dat'   0.0  35.0
+   'Snapshot_040s.dat'   0.0  40.0
+   'Snapshot_045s.dat'   0.0  45.0
+   'Snapshot_050s.dat'   0.0  50.0
+   'Snapshot_055s.dat'   0.0  55.0
+   'Snapshot_060s.dat'   0.0  60.0
+   'Snapshot_065s.dat'   0.0  65.0
+   'Snapshot_070s.dat'   0.0  70.0
+   'Snapshot_075s.dat'   0.0  75.0
+   'Snapshot_080s.dat'   0.0  80.0
+   'SnapshotR_010s.dat'   5.0  10.0
+   'SnapshotR_015s.dat'  10.0  15.0
+   'SnapshotR_020s.dat'  15.0  20.0
+   'SnapshotR_025s.dat'  20.0  25.0
+   'SnapshotR_030s.dat'  25.0  30.0
+   'SnapshotR_035s.dat'  30.0  35.0
+   'SnapshotR_040s.dat'  35.0  40.0
+   'SnapshotR_045s.dat'  40.0  45.0
+   'SnapshotR_050s.dat'  45.0  50.0
+   'SnapshotR_055s.dat'  50.0  55.0
+   'SnapshotR_060s.dat'  55.0  60.0
+   'SnapshotR_065s.dat'  60.0  65.0
+   'SnapshotR_070s.dat'  65.0  70.0
+   'SnapshotR_075s.dat'  70.0  75.0
+   'SnapshotR_080s.dat'  75.0  80.0
+#================================end of input====================================

examples/ids2020-iquique.inp

@@ -0,0 +1,236 @@
+# Example for the input file of FORTRAN77 program "ids2020" for kinematic earthquake
+# source imaging based on the Iterative Deconvolution and Stacking (IDS) Method using
+# seismic waveform and static displacement data jointly.
+#
+# written by Rongjiang Wang
+# Helmholtz Centre Potsdam
+# GFZ German Research Centre for Geosciences
+# e-mail: wang@gfz-potsdam.de
+# Phone +49 331 2881209
+# Fax +49 331 2881204
+#
+# Last modified: Potsdam, April, 2020
+#
+# Reference:
+# Zhang and Wang et al. (2014). Automatic imaging of earthquake rupture processes
+# by iterative deconvolution and stacking of high-rate GNSS and strong motion
+# seismograms, JGR Solid Earth, 199(7), 5633-5650.
+#
+#################################################################
+##                                                             ##
+## If not specified otherwise, SI Unit System is used overall! ##
+##                                                             ##
+##                                                             ##
+## Each comment line should start with "#"                     ##
+## and no blank line is allowed!                               ##
+##                                                             ##
+#################################################################
+#
+#================================================================================
+# EARTHQUAKE LOCATION PARAMETERS
+# ------------------------------
+# 1. origin time [s] (year, month, day, hour, minute [positive integer number]
+#    and second [positive or negative real number]), which will be also used as
+#    the reference time of the output data
+# 2. hypocentre (lat[deg], lon[deg], depth[km])
+# 3. preliminary estimate of moment magnitude
+#    Note: the hypocentre location and origin time (at least to minute) should be
+#          consistent with those given in the "StationInfo.dat" file of the
+#          waveform data (s. below)
+#--------------------------------------------------------------------------------
+  2014 4 3 2 42 72.66
+  -20.5984  -70.6320  27.33
+  7.70
+#================================================================================
+# WAVEFORM DATA (HIGH-RATE-GNSS/STRONG-MOTION/TELE-SEISMIC SEISMOGRAMS)
+# --------------------------------------------------------------------
+# 1. folder name of the data
+#    Note: this folder should include a file 'StationInfo.dat' giving all
+#          necessary information about the local monitoring network.
+# 2. lower and upper thresholds of the epicentral distances [km] (stations
+#    outside the threshold distances will not be used)
+# 3. weight of the static offset data (obtained automatically from the velocity
+#    waveforms after an empirical baseline correction) relative to the waveform
+#    data (suggested value 0.0-1.0)
+#--------------------------------------------------------------------------------
+   './WaveformData/'
+  0.0  500.0
+  0.1
+#================================================================================
+# GREEN'S FUNCTION DATABASE
+# -------------------------
+# 1. Green's function database (folder name)
+# 2. Parameters of Butterworth bandpass filter applied to the synthetics:
+#    order, lower and upper cutoff frequencies [Hz]
+#    computation [Hz]
+#-------------------------------------------------------------------------------- 
+ '../GreenChileSM/GreenFunctions/'
+  3  0.00  0.50
+#================================================================================
+# GNSS STATIC OFFSETS
+# -------------------
+# 1. number of GNSS sites, file name of the data
+#    Note: if the number of GNSS sites <= 0, then no GNSS data are available.
+#    the GNSS data file includes one header line and 8 culomns,
+#    Lat(deg)  Lon(deg)  East(m)  North(m)  Up(m)  wf_E  wf_N   wf_Z
+#    where wf_E/N/Z are within-dataset relative weighting factor proportional to
+#    inverse error variance
+# 2. weight of the whole GNSS dataset relative to the waveform dataset (suggested
+#    value 0.0-1.0)
+#--------------------------------------------------------------------------------
+ 21  './GPSData/GPS_Iquique.dat'
+ 1.0
+#================================================================================
+# InSAR LOS DISPLACEMENTS
+# -----------------------
+# 1. number of InSAR grids, file name of the data
+#    Note: if the number of InSAR grids <= 0, then no InSAR data are used.
+#          the InSAR data file includes one header line and 8 culomns,
+#          Lat(deg)  Lon(deg) LOS(m)    wf_los   incidence(deg)   azimuth(deg)
+#          where wf_los is the within-dataset relative weighting factor
+#          proportional to inverse error variance
+# 2. weight of the whole InSAR dataset in the misfit variance relative to the
+#    waveform dataset (suggested value 0.0-1.0)
+#--------------------------------------------------------------------------------
+ 0  './InSAR/NoInSAR.dat'
+ 0.00
+#================================================================================
+# Geodetic GREEN'S FUNCTION DATABASE
+# ----------------------------------
+# 1. geodetic Green's function database (folder name)
+# 2. 3 file names of the geodetic Green's functions
+#-------------------------------------------------------------------------------- 
+ '../GreenChileGD/DGreen/'
+ 'uz'  'ur'  'ut'
+#================================================================================
+# FAULT PARAMETERS
+# ----------------
+# 1. ipatch: switch (0/1) for subfaults input form,
+#    0 = read discrete patches from files,
+#    1 = automatic discretization of a single rectangular fault
+#-------------------------------------------------------------------------------- 
+#   idisc
+#-------------------------------------------------------------------------------- 
+    1
+#-------------------------------------------------------------------------------- 
+# IF(idisc = 0)THEN
+#
+#   subfault paramters
+#   ==================
+#   file_name:         data file for discrete patches
+#   number of sub-faults
+#   iref               selection of the patch reference location
+#
+#   Note: the file includes one header line followed by 8 columns of data:
+#
+#   latitude[deg], longitude[deg], depth[km], length[km], width[km],
+#   strike[deg], dip[deg], rake[deg]
+#
+#   where (lat, lon, dep) is the reference location of the patch depending
+#   on iref: 1 = upper-left corner, 2 = upper-right corner,
+#            3 = lower-left corner, 4 = lower-right corner,
+#            5 = central point (s. Fig.)
+#
+#        definitions for a rectangular fault patch
+#        =========================================
+#
+#                 north(x)
+#                  /
+#                 / | strike
+#                1-----------------------2
+#                |\        p .            \ W
+#                :-\      i .              \ i
+#                |  \    l .     5          \ d
+#                :90 \  S .                  \ t
+#                |-dip\  .                    \ h
+#                :     \. ) rake               \ 
+#          downward(z)  3-----------------------4
+#                         L   e   n   g   t   h
+#--------------------------------------------------------------------------------
+#   file_name               nsf   iref
+#--------------------------------------------------------------------------------
+#   './fault_patches.dat'    288   5
+#--------------------------------------------------------------------------------
+#
+# ELSE IF(idisc = 1)THEN
+#
+# 2. focal mechanism: strike, dip and rake angles [deg]
+#--------------------------------------------------------------------------------
+   347.0     25.00    87.00
+#================================================================================
+# IDS PARAMETERS
+# --------------
+# 1. maximum numbers of iterations.
+#    Note: if max. iteration = 0, forward modelling for any dummy data set is
+#          made using the given rupture model and the filter parameters
+#          if max. iteration < 0, same as iteration = 0, but without filtering.
+#--------------------------------------------------------------------------------
+	100
+#================================================================================
+# OUTPUT
+# ------
+# 1. folder name for the output files (should exist already)
+# 2. sub-folder name for comparisons of observed and modelled waveform data
+#    (should exist already)
+# 3. sub-folder name for results from empirical baseline correction of
+#    waveform data (this folder should exist already)
+# 4. file name of major earthquake parameters
+# 5. file name of earthquake source-time-function (moment rate history)
+# 6. file name of sub-fault source-time-function (moment rate history)
+# 7. file name of rupture propagation
+# 8. file name of final slip model
+#    Note: if forward modelling is selected (s. above), the above two files
+#          should provide data for the given kinematic rupture model and will
+#          be overwritten after the IDS inversion
+# 9. file name of normalized residual waveform variance
+#10. file name of predicted static offsets at seismic waveform data sites
+#11. file name of comparison between observed and modelled GNSS data
+#12. file name of comparison between observed and modelled InSAR data
+#13. number of snapshots of the rupture process
+#14. file name of the 1. snapshot for the slip released within the given
+#    time window (t1[s], t2[s])
+#--------------------------------------------------------------------------------
+   './Output'
+   'O2S'
+   'EBC'
+   'EQ_MP.dat'
+   'EQ_STF.dat'
+   'PT_STF.dat'
+   'Rup_Front.dat'
+   'SlipModel.dat'
+   'WFMisfitvar.dat'
+   'SMStaticOffset.dat'
+   'GNSSDatafit.dat'
+   'InSARDatafit.dat'
+   30
+   'Snapshot_010s.dat'   0.0  10.0
+   'Snapshot_015s.dat'   0.0  15.0
+   'Snapshot_020s.dat'   0.0  20.0
+   'Snapshot_025s.dat'   0.0  25.0
+   'Snapshot_030s.dat'   0.0  30.0
+   'Snapshot_035s.dat'   0.0  35.0
+   'Snapshot_040s.dat'   0.0  40.0
+   'Snapshot_045s.dat'   0.0  45.0
+   'Snapshot_050s.dat'   0.0  50.0
+   'Snapshot_055s.dat'   0.0  55.0
+   'Snapshot_060s.dat'   0.0  60.0
+   'Snapshot_065s.dat'   0.0  65.0
+   'Snapshot_070s.dat'   0.0  70.0
+   'Snapshot_075s.dat'   0.0  75.0
+   'Snapshot_080s.dat'   0.0  80.0
+   'SnapshotR_010s.dat'   5.0  10.0
+   'SnapshotR_015s.dat'  10.0  15.0
+   'SnapshotR_020s.dat'  15.0  20.0
+   'SnapshotR_025s.dat'  20.0  25.0
+   'SnapshotR_030s.dat'  25.0  30.0
+   'SnapshotR_035s.dat'  30.0  35.0
+   'SnapshotR_040s.dat'  35.0  40.0
+   'SnapshotR_045s.dat'  40.0  45.0
+   'SnapshotR_050s.dat'  45.0  50.0
+   'SnapshotR_055s.dat'  50.0  55.0
+   'SnapshotR_060s.dat'  55.0  60.0
+   'SnapshotR_065s.dat'  60.0  65.0
+   'SnapshotR_070s.dat'  65.0  70.0
+   'SnapshotR_075s.dat'  70.0  75.0
+   'SnapshotR_080s.dat'  75.0  80.0
+#================================end of input====================================

examples/ids2020-pisagua.inp

@@ -0,0 +1,236 @@
+# Example for the input file of FORTRAN77 program "ids2020" for kinematic earthquake
+# source imaging based on the Iterative Deconvolution and Stacking (IDS) Method using
+# seismic waveform and static displacement data jointly.
+#
+# written by Rongjiang Wang
+# Helmholtz Centre Potsdam
+# GFZ German Research Centre for Geosciences
+# e-mail: wang@gfz-potsdam.de
+# Phone +49 331 2881209
+# Fax +49 331 2881204
+#
+# Last modified: Potsdam, April, 2020
+#
+# Reference:
+# Zhang and Wang et al. (2014). Automatic imaging of earthquake rupture processes
+# by iterative deconvolution and stacking of high-rate GNSS and strong motion
+# seismograms, JGR Solid Earth, 199(7), 5633-5650.
+#
+#################################################################
+##                                                             ##
+## If not specified otherwise, SI Unit System is used overall! ##
+##                                                             ##
+##                                                             ##
+## Each comment line should start with "#"                     ##
+## and no blank line is allowed!                               ##
+##                                                             ##
+#################################################################
+#
+#================================================================================
+# EARTHQUAKE LOCATION PARAMETERS
+# ------------------------------
+# 1. origin time [s] (year, month, day, hour, minute [positive integer number]
+#    and second [positive or negative real number]), which will be also used as
+#    the reference time of the output data
+# 2. hypocentre (lat[deg], lon[deg], depth[km])
+# 3. preliminary estimate of moment magnitude
+#    Note: the hypocentre location and origin time (at least to minute) should be
+#          consistent with those given in the "StationInfo.dat" file of the
+#          waveform data (s. below)
+#--------------------------------------------------------------------------------
+  2014 4 1 23 45  105.55
+  -19.5764    -70.9037    28.54
+  8.10
+#================================================================================
+# WAVEFORM DATA (HIGH-RATE-GNSS/STRONG-MOTION SEISMOGRAMS)
+# --------------------------------------------------------------------
+# 1. folder name of the data
+#    Note: this folder should include a file 'StationInfo.dat' giving all
+#          necessary information about the local monitoring network.
+# 2. lower and upper thresholds of the epicentral distances [km] (stations
+#    outside the threshold distances will not be used)
+# 3. weight of the static offset data (obtained automatically from the velocity
+#    waveforms after an empirical baseline correction) relative to the waveform
+#    data (suggested value 0.0-1.0)
+#--------------------------------------------------------------------------------
+   './WaveformData/'
+  0.0  800.0
+  1.0
+#================================================================================
+# GREEN'S FUNCTION DATABASE
+# -------------------------
+# 1. Green's function database (folder name)
+# 2. Parameters of Butterworth bandpass filter applied to the synthetics:
+#    order, lower and upper cutoff frequencies [Hz]
+#    computation [Hz]
+#-------------------------------------------------------------------------------- 
+ '../GreenChileSM/GreenFunctions/'
+  3  0.00  0.50
+#================================================================================
+# GNSS STATIC OFFSETS
+# -------------------
+# 1. number of GNSS sites, file name of the data
+#    Note: if the number of GNSS sites <= 0, then no GNSS data are available.
+#    the GNSS data file includes one header line and 8 culomns,
+#    Lat(deg)  Lon(deg)  East(m)  North(m)  Up(m)  wf_E  wf_N   wf_Z
+#    where wf_E/N/Z are within-dataset relative weighting factor proportional to
+#    inverse error variance
+# 2. weight of the whole GNSS dataset relative to the waveform dataset (suggested
+#    value 0.0-1.0)
+#--------------------------------------------------------------------------------
+ 21  './GPSData/GPS_Pisagua.dat'
+ 1.0
+#================================================================================
+# InSAR LOS DISPLACEMENTS
+# -----------------------
+# 1. number of InSAR grids, file name of the data
+#    Note: if the number of InSAR grids <= 0, then no InSAR data are used.
+#          the InSAR data file includes one header line and 8 culomns,
+#          Lat(deg)  Lon(deg) LOS(m)    wf_los   incidence(deg)   azimuth(deg)
+#          where wf_los is the within-dataset relative weighting factor
+#          proportional to inverse error variance
+# 2. weight of the whole InSAR dataset in the misfit variance relative to the
+#    waveform dataset (suggested value 0.0-1.0)
+#--------------------------------------------------------------------------------
+ 0  './InSAR/NoInSAR.dat'
+ 0.00
+#================================================================================
+# Geodetic GREEN'S FUNCTION DATABASE
+# ----------------------------------
+# 1. geodetic Green's function database (folder name)
+# 2. 3 file names of the geodetic Green's functions
+#-------------------------------------------------------------------------------- 
+ '../GreenChileGD/DGreen/'
+ 'uz'  'ur'  'ut'
+#================================================================================
+# FAULT PARAMETERS
+# ----------------
+# 1. ipatch: switch (0/1) for subfaults input form,
+#    0 = read discrete patches from files,
+#    1 = automatic discretization of a single rectangular fault
+#-------------------------------------------------------------------------------- 
+#   idisc
+#-------------------------------------------------------------------------------- 
+    1
+#-------------------------------------------------------------------------------- 
+# IF(idisc = 0)THEN
+#
+#   subfault paramters
+#   ==================
+#   file_name:         data file for discrete patches
+#   number of sub-faults
+#   iref               selection of the patch reference location
+#
+#   Note: the file includes one header line followed by 8 columns of data:
+#
+#   latitude[deg], longitude[deg], depth[km], length[km], width[km],
+#   strike[deg], dip[deg], rake[deg]
+#
+#   where (lat, lon, dep) is the reference location of the patch depending
+#   on iref: 1 = upper-left corner, 2 = upper-right corner,
+#            3 = lower-left corner, 4 = lower-right corner,
+#            5 = central point (s. Fig.)
+#
+#        definitions for a rectangular fault patch
+#        =========================================
+#
+#                 north(x)
+#                  /
+#                 / | strike
+#                1-----------------------2
+#                |\        p .            \ W
+#                :-\      i .              \ i
+#                |  \    l .     5          \ d
+#                :90 \  S .                  \ t
+#                |-dip\  .                    \ h
+#                :     \. ) rake               \ 
+#          downward(z)  3-----------------------4
+#                         L   e   n   g   t   h
+#--------------------------------------------------------------------------------
+#   file_name               nsf   iref
+#--------------------------------------------------------------------------------
+#   './fault_patches.dat'    288   5
+#--------------------------------------------------------------------------------
+#
+# ELSE IF(idisc = 1)THEN
+#
+# 2. focal mechanism: strike, dip and rake angles [deg]
+#--------------------------------------------------------------------------------
+   346.0     13.00     84.00
+#================================================================================
+# IDS PARAMETERS
+# --------------
+# 1. maximum numbers of iterations.
+#    Note: if max. iteration = 0, forward modelling for any dummy data set is
+#          made using the given rupture model and the filter parameters
+#          if max. iteration < 0, same as iteration = 0, but without filtering.
+#--------------------------------------------------------------------------------
+	100
+#================================================================================
+# OUTPUT
+# ------
+# 1. folder name for the output files (should exist already)
+# 2. sub-folder name for comparisons of observed and modelled waveform data
+#    (should exist already)
+# 3. sub-folder name for results from empirical baseline correction of
+#    waveform data (this folder should exist already)
+# 4. file name of major earthquake parameters
+# 5. file name of earthquake source-time-function (moment rate history)
+# 6. file name of sub-fault source-time-function (moment rate history)
+# 7. file name of rupture propagation
+# 8. file name of final slip model
+#    Note: if forward modelling is selected (s. above), the above two files
+#          should provide data for the given kinematic rupture model and will
+#          be overwritten after the IDS inversion
+# 9. file name of normalized residual waveform variance
+#10. file name of predicted static offsets at seismic waveform data sites
+#11. file name of comparison between observed and modelled GNSS data
+#12. file name of comparison between observed and modelled InSAR data
+#13. number of snapshots of the rupture process
+#14. file name of the 1. snapshot for the slip released within the given
+#    time window (t1[s], t2[s])
+#--------------------------------------------------------------------------------
+   './Output'
+   'O2S'
+   'EBC'
+   'EQ_MP.dat'
+   'EQ_STF.dat'
+   'PT_STF.dat'
+   'Rup_Front.dat'
+   'SlipModel.dat'
+   'WFMisfitvar.dat'
+   'SMStaticOffset.dat'
+   'GNSSDatafit.dat'
+   'InSARDatafit.dat'
+   30
+   'Snapshot_010s.dat'   0.0  10.0
+   'Snapshot_015s.dat'   0.0  15.0
+   'Snapshot_020s.dat'   0.0  20.0
+   'Snapshot_025s.dat'   0.0  25.0
+   'Snapshot_030s.dat'   0.0  30.0
+   'Snapshot_035s.dat'   0.0  35.0
+   'Snapshot_040s.dat'   0.0  40.0
+   'Snapshot_045s.dat'   0.0  45.0
+   'Snapshot_050s.dat'   0.0  50.0
+   'Snapshot_055s.dat'   0.0  55.0
+   'Snapshot_060s.dat'   0.0  60.0
+   'Snapshot_065s.dat'   0.0  65.0
+   'Snapshot_070s.dat'   0.0  70.0
+   'Snapshot_075s.dat'   0.0  75.0
+   'Snapshot_080s.dat'   0.0  80.0
+   'SnapshotR_010s.dat'   5.0  10.0
+   'SnapshotR_015s.dat'  10.0  15.0
+   'SnapshotR_020s.dat'  15.0  20.0
+   'SnapshotR_025s.dat'  20.0  25.0
+   'SnapshotR_030s.dat'  25.0  30.0
+   'SnapshotR_035s.dat'  30.0  35.0
+   'SnapshotR_040s.dat'  35.0  40.0
+   'SnapshotR_045s.dat'  40.0  45.0
+   'SnapshotR_050s.dat'  45.0  50.0
+   'SnapshotR_055s.dat'  50.0  55.0
+   'SnapshotR_060s.dat'  55.0  60.0
+   'SnapshotR_065s.dat'  60.0  65.0
+   'SnapshotR_070s.dat'  65.0  70.0
+   'SnapshotR_075s.dat'  70.0  75.0
+   'SnapshotR_080s.dat'  75.0  80.0
+#================================end of input====================================

src/Makefile.am

@@ -0,0 +1,9 @@
+bin_PROGRAMS = ids2020
+ids2020_SOURCES = banddpasss.f bscbilin.f bscmono.f butterworth.f d2dfit.f dc3d.f deepphase.f \
+                   depthp.f depthphase.f disazi.f four1w.f idsalloc.f idsdifgrn.f idsforward.f \
+                   idsgetinp.f idsgnsdat.f idsgnsgrn.f idsinverse.f idskernel.f idsmain.f \
+                   idsmodel.f idsosmgrn.f idsoutput.f idspstt.f idssardat.f idssargrn.f \
+                   idsscale.f idssmdat.f idssmgrn.f idssmwei.f moments.f rampfit.f skipdoc.f \
+                   stepfit.f swavelet.f tpssub.f
+clean-local:
+	-rm -f *.mod

src/banddpasss.f

@@ -0,0 +1,44 @@
+      subroutine bandpass(n,f1,f2,df,nf,h)
+      implicit none
+      integer*4 n,nf
+      real*8 f1,f2,df
+      complex*16 h(nf)
+c
+      integer*4 l,k
+      real*8 w,w1,w2,delt,arg
+c
+      real*8 PI
+      data PI/3.14159265358979d0/
+c
+      if(n.le.0.or.f1.ge.f2)then
+        do l=1,nf
+          h(l)=(1.d0,0.d0)
+        enddo
+      else if(f1.le.0.d0)then
+        w2=2.d0*PI*f2
+        do l=1,nf
+          w=2.d0*PI*dble(l-1)*df
+          delt=w2-w1
+          h(l)=(1.d0,0.d0)
+          do k=1,n
+            arg=PI*dble(2*k-1)/dble(2*n)
+            h(l)=h(l)*dcmplx(0.d0,-delt)
+     &          /dcmplx(w-delt*dcos(arg),-delt*dsin(arg))
+          enddo
+        enddo
+      else
+        w1=2.d0*PI*f1
+        w2=2.d0*PI*f2
+        do l=1,nf
+          w=2.d0*PI*dble(l-1)*df
+          delt=w*(w2-w1)
+          h(l)=(1.d0,0.d0)
+          do k=1,n
+            arg=PI*dble(2*k-1)/dble(2*n)
+            h(l)=h(l)*dcmplx(0.d0,-delt)
+     &          /dcmplx(w*w-w1*w2-delt*dcos(arg),-delt*dsin(arg))
+          enddo
+        enddo
+      endif
+      return
+      end

src/bscbilin.f

@@ -0,0 +1,127 @@
+      subroutine bscbilin(n,ip,vel,bserr,dt,offset,sigma)
+      implicit none
+c
+      integer*4 n,ip
+      real*8 dt,offset,sigma
+      real*8 vel(n),bserr(n)
+c
+      integer i,j,k,i0,i1,i2,iw,ia,ib,id,ierr
+      real*8 a,b,a0,a1,delta,gamma,smin
+c
+      integer*4 ndeg
+      real*8 PI
+      data ndeg,PI/1,3.14159265358979d0/
+c
+      real*8, allocatable:: poly(:,:),bat(:),swap(:)
+c
+      allocate (swap(n),stat=ierr)
+      if(ierr.ne.0)stop ' Error in bscmono: swap not allocated!'
+      allocate (poly(n,0:ndeg),stat=ierr)
+      if(ierr.ne.0)stop ' Error in bscmono: poly not allocated!'
+      allocate (bat(0:ndeg),stat=ierr)
+      if(ierr.ne.0)stop ' Error in bscmono: bat not allocated!'
+c
+c     remove pre-event trend
+c
+      call d2dfit(ip,vel,poly,bat,1,bserr)
+      a0=bserr(1)
+      a1=(bserr(ip)-bserr(1))/dble(ip-1)
+c
+      do i=1,n
+        vel(i)=vel(i)-(a0+a1*dble(i-1))
+      enddo
+c
+c     estimate signal duration
+c
+      swap(1)=0.d0
+      do i=2,n
+        swap(i)=swap(i-1)+dabs(vel(i)-vel(i-1))
+      enddo
+      iw=2
+      do i=2,n-1
+        if(swap(i).le.0.85d0*swap(n))iw=i
+      enddo
+c
+c     empirical bi-linear baseline correction
+c
+      call d2dfit(1+n-iw,vel(iw),poly,bat,1,bserr(iw))
+c
+      a=bserr(iw)
+      b=(bserr(n)-bserr(iw))/dble(n-iw)
+      do i=1,iw-1
+        bserr(i)=a+b*dble(i-iw)
+      enddo
+c
+      swap(1)=0.d0
+      do i=2,iw-1
+        swap(i)=swap(i-1)+vel(i)
+      enddo
+      do i=iw,n
+        swap(i)=swap(i-1)+vel(i)-bserr(i)
+      enddo
+      call stepfit(n,swap,ip,iw,i0,delta,smin)
+      ia=iw-1
+      ib=iw
+c
+      id=max0(1,(iw-ip)/50)
+c
+      do k=ip,iw-id,id
+        do j=max0(k+id,(ip+iw)/2),iw,id
+          do i=1,ip-1
+            swap(i)=0.d0
+          enddo
+          do i=ip,k
+            swap(i)=swap(i-1)+vel(i)
+          enddo
+          do i=k+1,j-1
+            swap(i)=swap(i-1)+vel(i)-bserr(j)*dble(i-k)/dble(j-k)
+          enddo
+          do i=j,n
+            swap(i)=swap(i-1)+vel(i)-bserr(i)
+          enddo
+c
+          call stepfit(n,swap,ip,iw,i0,delta,gamma)
+c
+          if(gamma.le.smin)then
+            ia=k
+            ib=j
+            smin=gamma
+          endif
+        enddo
+      enddo
+c
+      do i=1,ia
+        bserr(i)=0.d0
+      enddo
+      a=bserr(ib)
+      do i=ia+1,ib-1
+        bserr(i)=a*dble(i-ia)/dble(ib-ia)
+      enddo
+c
+      do i=1,n
+        vel(i)=vel(i)-bserr(i)
+      enddo
+c
+      swap(1)=0.d0
+      do i=2,n
+        swap(i)=swap(i-1)+vel(i)*dt
+      enddo
+      call rampfit(n,swap,ip,iw,i1,i2,offset,sigma,poly)
+      sigma=0.d0
+      do i=i2,iw+1
+        sigma=sigma+dabs(swap(i)-offset)
+      enddo
+      sigma=sigma/dble(1+iw-i2)
+      delta=0.d0
+      do i=ip,i2
+        delta=delta+bserr(i)*dt
+      enddo
+      sigma=dmax1(sigma,dabs(delta))
+c
+      do i=1,n
+        bserr(i)=bserr(i)+a0+a1*dble(i-1)
+      enddo
+c
+      deallocate(swap,poly,bat)
+      return
+      end

src/bscmono.f

@@ -0,0 +1,313 @@
+      subroutine bscmono(nw,ipp,vel,bse,dt,offset,sigma)
+      implicit none
+c
+      integer*4 nw,ipp
+      real*8 dt,offset,sigma
+      real*8 vel(nw),bse(nw)
+c
+      integer*4 i,j,i0,i1,i2,i3,i4,i5,i6,iw,ierr
+      integer*4 k,k1,k2,k3,k4,k5,k6,ks,id,ip,npeak,ipga
+      real*8 a0,a1,delta,smin,gamma,pga,pgv
+      real*8 d1,d2,d3,d4,d5,d6,d1opt,d2opt,d3opt,d4opt,d5opt
+c
+      integer*4 ndeg
+      parameter(ndeg=2)
+c
+      real*8 pi,sdw
+      data pi,sdw/3.14159265358979d+00,0.75d+00/
+c
+      real*8, allocatable:: poly(:,:),bat(:),swp(:)
+      allocate (poly(nw,0:ndeg),stat=ierr)
+      if(ierr.ne.0)stop ' Error in bscmono: poly not allocated!'
+      allocate (bat(0:ndeg),stat=ierr)
+      if(ierr.ne.0)stop ' Error in bscmono: bat not allocated!'
+      allocate (swp(nw),stat=ierr)
+      if(ierr.ne.0)stop ' Error in bscmono: swp not allocated!'
+c
+c     remove pre-event trend
+c
+      ip=ipp
+      call d2dfit(ip,vel,poly,bat,1,bse)
+c
+      a0=bse(1)
+      a1=(bse(ip)-bse(1))/dble(ip-1)
+c
+      do i=1,nw
+        vel(i)=vel(i)-(a0+a1*dble(i-1))
+      enddo
+c
+c     uncorrected acceleration
+c
+      do i=2,nw
+        swp(i)=(vel(i)-vel(i-1))/dt
+      enddo
+      swp(1)=swp(2)
+c
+c     estimate pga time
+c
+      ipga=ip
+      pga=dabs(swp(ip))
+      do i=ip+1,nw
+        if(pga.lt.dabs(swp(i)))then
+          ipga=i
+          pga=dabs(swp(i))
+        endif
+      enddo
+c
+      swp(1)=0.d0
+      do i=2,nw
+        swp(i)=swp(i-1)+swp(i)**2
+      enddo
+c
+c     estimate energy signal window
+c
+      iw=ip
+      do i=ip+1,nw
+        if(swp(i).ge.sdw*swp(nw))then
+          iw=i
+          goto 100
+        endif
+      enddo
+100   continue
+c
+      iw=min0((2*nw+ip)/3,max0(iw,2*ipga-ip))
+c
+c     estimate earliest start time of baseline shift
+c
+      swp(1)=0.d0
+      do i=2,nw
+        swp(i)=swp(i-1)+vel(i)*dt
+      enddo
+c
+      i0=ip
+      do i=1+(2*ip+ipga)/3,ip+1,-1
+        if(swp(i)*swp(i-1).le.0.d0)then
+          i0=i
+          goto 200
+        endif
+      enddo
+200   continue
+      ip=i0
+c
+c     main smoothing procedure
+c
+      i0=1+(2*nw+iw)/3
+      call d2dfit(2+nw-i0,swp(i0-1),poly,bat,2,bse(i0-1))
+      do i=i0,nw
+        swp(i)=(bse(i)-bse(i-1))/dt
+      enddo
+      d1=(swp(nw)-swp(i0))/dble(nw-i0)
+c
+      do i=1,ip
+        bse(i)=0.d0
+      enddo
+      do i=ip+1,i0
+        bse(i)=vel(i)
+      enddo
+      do i=i0+1,nw
+        bse(i)=vel(i)*dcos(0.5d0*pi*dble(i-i0)/dble(nw-i0))**2
+     &        +swp(i)*dsin(0.5d0*pi*dble(i-i0)/dble(nw-i0))**2
+      enddo
+c
+      k=0
+400   k=k+1
+      do i=ip,nw
+        swp(i)=bse(i)
+      enddo
+      do i=ip+1,nw-1
+        bse(i)=(swp(i-1)+swp(i)+swp(i+1))/3.d0
+      enddo
+c
+      npeak=0
+      do i=iw+1,nw-1
+        if(bse(i).lt.bse(i-1)+d1.and.
+     &     bse(i).lt.bse(i+1)-d1.or.
+     &     bse(i).gt.bse(i-1)+d1.and.
+     &     bse(i).gt.bse(i+1)-d1)npeak=npeak+1
+      enddo
+      if(npeak.gt.1)goto 400
+c
+      gamma=0.d0
+      do i=iw,nw
+        gamma=dmax1(gamma,dabs(vel(i)-bse(i)))
+      enddo
+c
+      pgv=gamma
+      do i=1,iw-1
+        pgv=dmax1(pgv,dabs(vel(i)-bse(i)))
+      enddo
+c
+      i0=iw
+      npeak=0
+      do i=iw-1,ip+1,-1
+        gamma=dmax1(gamma,dabs(vel(i)-bse(i)))
+        if(bse(i).lt.bse(i-1)+d1.and.
+     &     bse(i).lt.bse(i+1)-d1.or.
+     &     bse(i).gt.bse(i-1)+d1.and.
+     &     bse(i).gt.bse(i+1)-d1)npeak=npeak+1
+        if(npeak.gt.0.or.
+     &     gamma.ge.0.25d0*pgv.and.dabs(bse(i)).le.0.5d0*pgv)then
+          goto 500
+        else
+          i0=i
+        endif
+      enddo
+500   continue
+      if(npeak.gt.0)then
+        iw=(2*i0+iw)/3
+      else
+        iw=i0
+      endif
+c
+      call d2dfit(1+nw-iw,bse(iw),poly,bat,1,swp(iw))
+c
+      d1=swp(iw+1)-swp(iw)
+      d2=swp(iw)/dble(iw-ip)
+      if(d2*d1.ge.0.d0.and.dabs(d2).lt.dabs(d1))then
+        ip=(ip+iw-idint(swp(iw)/d1))/2
+      endif
+c
+      do i=1,iw
+        bse(i)=0.d0
+      enddo
+c
+c     estimate best-fit monotonic trend within signal window by grid search
+c     assuming that baseline shift increases faster than ~sqrt(t)
+c
+      id=1+(iw-ip)/6
+      if(ip+6*id.ge.nw)id=id-1
+      i1=ip+id
+      i2=i1+id
+      i3=i2+id
+      i4=i3+id
+      i5=i4+id
+      i6=i5+id
+      smin=0.d0
+c
+      k6=18
+      d6=bse(i6)
+      delta=d6/18.d0
+c
+      do k5=0,16
+        d5=dble(k5)*delta
+        do k4=0,min0(15,k5)
+          d4=dble(k4)*delta
+          do k3=0,min0(13,k4)
+            d3=dble(k3)*delta
+            do k2=0,min0(10,k3)
+              d2=dble(k2)*delta
+              do k1=0,min0(7,k2)
+                d1=dble(k1)*delta
+                sigma=0.d0
+                do i=ip+1,i1
+                  sigma=sigma+(vel(i)
+     &                 -d1*(dble(i-ip)/dble(id))**2)**2
+                enddo
+                do i=i1+1,i2
+                  sigma=sigma+(vel(i)
+     &                 -(d1+(d2-d1)*dble(i-i1)/dble(id)))**2
+                enddo
+                do i=i2+1,i3
+                  sigma=sigma+(vel(i)
+     &                 -(d2+(d3-d2)*dble(i-i2)/dble(id)))**2
+                enddo
+                do i=i3+1,i4
+                  sigma=sigma+(vel(i)
+     &                 -(d3+(d4-d3)*dble(i-i3)/dble(id)))**2
+                enddo
+                do i=i4+1,i5
+                  sigma=sigma+(vel(i)
+     &                 -(d4+(d5-d4)*dble(i-i4)/dble(id)))**2
+                enddo
+                do i=i5+1,i6
+                  sigma=sigma+(vel(i)
+     &                 -(d5+(d6-d5)*dble(i-i5)/dble(id)))**2
+                enddo
+                if(smin.le.0.d0.or.sigma.le.smin)then
+                  smin=sigma
+                  d1opt=d1
+                  d2opt=d2
+                  d3opt=d3
+                  d4opt=d4
+                  d5opt=d5
+                endif
+              enddo
+            enddo
+          enddo
+        enddo
+      enddo
+c
+      do i=1,ip
+        bse(i)=0.d0
+      enddo
+      do i=ip+1,i1
+        bse(i)=d1opt*(dble(i-ip)/dble(id))**2
+      enddo
+      do i=i1+1,i2
+        bse(i)=d1opt+(d2opt-d1opt)*dble(i-i1)/dble(id)
+      enddo
+      do i=i2+1,i3
+        bse(i)=d2opt+(d3opt-d2opt)*dble(i-i2)/dble(id)
+      enddo
+      do i=i3+1,i4
+        bse(i)=d3opt+(d4opt-d3opt)*dble(i-i3)/dble(id)
+      enddo
+      do i=i4+1,i5
+        bse(i)=d4opt+(d5opt-d4opt)*dble(i-i4)/dble(id)
+      enddo
+      do i=i5+1,i6-1
+        bse(i)=d5opt+(d6-d5opt)*dble(i-i5)/dble(id)
+      enddo
+c
+600   continue
+c
+c     re-smoothing
+c
+      d1=bse(i6)/dble(i6-ip)
+      ks=1+(i6-ip)/2
+c
+      k=0
+800   k=k+1
+      do i=ipp,min0(nw,i6+k)
+        swp(i)=bse(i)
+      enddo
+      do i=ipp+1,min0(nw,i6+k)-1
+        bse(i)=(swp(i-1)+swp(i)+swp(i+1))/3.d0
+      enddo
+c
+      npeak=0
+      do i=ipp+1,i6-1
+        if(bse(i).lt.bse(i-1)+d1.and.
+     &     bse(i).lt.bse(i+1)-d1.or.
+     &     bse(i).gt.bse(i-1)+d1.and.
+     &     bse(i).gt.bse(i+1)-d1)npeak=npeak+1
+      enddo
+      if(k.lt.ks.or.npeak.gt.1)goto 800
+c
+      do i=1,nw
+        vel(i)=vel(i)-bse(i)
+      enddo
+c
+      swp(1)=0.d0
+      do i=2,nw
+        swp(i)=swp(i-1)+vel(i)*dt
+      enddo
+      call rampfit(nw,swp,ipp,iw,i1,i2,offset,sigma,poly)
+      sigma=0.d0
+      do i=i2,iw+1
+        sigma=sigma+dabs(swp(i)-offset)
+      enddo
+      sigma=sigma/dble(1+iw-i2)
+      delta=0.d0
+      do i=ip,iw
+        delta=delta+bse(i)*dt
+      enddo
+      sigma=dmax1(sigma,dabs(delta))
+c
+      do i=1,nw
+        bse(i)=bse(i)+a0+a1*dble(i-1)
+      enddo
+c
+      deallocate(poly,bat,swp)
+      return
+      end

src/butterworth.f

@@ -0,0 +1,46 @@
+      subroutine butterworth(n,fc,df,nf,lpf)
+c
+c     butterworth lowpass filter
+c
+c     n: order
+c     fc: cutoff frequency
+c     df: frequency sampling interval
+c     nf: number of freqquency samples
+c     lpf: complex lowpass filter (return values)
+c
+      implicit none
+      integer*4 n,nf
+      real*8 fc,df
+      complex*16 lpf(nf)
+c
+      integer*4 l,k
+      complex*16 s,cs
+c
+      real*8 PI
+      data PI/3.14159265358979d0/
+c
+      if(n.le.0.or.fc.le.0.d0)then
+        do l=1,nf
+          lpf(l)=(1.d0,0.d0)
+        enddo
+      else if(n.eq.2*(n/2))then
+        do l=1,nf
+          s=dcmplx(0.d0,df*dble(l-1)/fc)
+          lpf(l)=(1.d0,0.d0)
+          do k=1,n/2
+            cs=dcmplx(2.d0*dcos(PI*dble(2*k+n-1)/dble(2*n)),0.d0)
+            lpf(l)=lpf(l)/(s*s-s*cs+(1.d0,0.d0))
+          enddo
+        enddo
+      else
+        do l=1,nf
+          s=dcmplx(0.d0,df*dble(l-1)/fc)
+          lpf(l)=(1.d0,0.d0)/(s+(1.d0,0.d0))
+          do k=1,(n-1)/2
+            cs=dcmplx(2.d0*dcos(PI*dble(2*k+n-1)/dble(2*n)),0.d0)
+            lpf(l)=lpf(l)/(s*s-s*cs+(1.d0,0.d0))
+          enddo
+        enddo
+      endif
+      return
+      end

src/d2dfit.f

@@ -0,0 +1,38 @@
+      subroutine d2dfit(n,dis,p,b,ndeg,disfit)
+      implicit none
+      integer*4 n,ndeg
+      real*8 dis(n),disfit(n),p(n,0:ndeg),b(0:ndeg)
+c
+      integer*4 i,ideg
+      real*8 x,dx
+c
+      dx=2.d0/dble(n-1)
+      do i=1,n
+        x=-1.d0+dble(i-1)*dx
+        p(i,0)=1.d0
+        if(ndeg.gt.0)then
+          p(i,1)=x
+          do ideg=2,ndeg
+            p(i,ideg)=(dble(2*ideg-1)*x*p(i,ideg-1)
+     &                -dble(ideg-1)*p(i,ideg-2))/dble(ideg)
+          enddo
+        endif
+      enddo
+c
+      do ideg=0,ndeg
+        b(ideg)=0.5d0*(dis(1)*p(1,ideg)+dis(n)*p(n,ideg))
+        do i=2,n-1
+          b(ideg)=b(ideg)+dis(i)*p(i,ideg)
+        enddo
+        b(ideg)=b(ideg)*dx*0.5d0*dble(2*ideg+1)
+      enddo
+c
+      do i=1,n
+        x=-1.+dble(i-1)*dx
+        disfit(i)=0.d0
+        do ideg=0,ndeg
+          disfit(i)=disfit(i)+b(ideg)*p(i,ideg)
+        enddo
+      enddo
+      return
+      end

src/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