Browse Source

added package

modified README from fomosto_qssp
pull/1/head
hvasbath 3 years ago
parent
commit
fa43afe2fe
64 changed files with 11740 additions and 2 deletions
  1. +14
    -0
      .gitignore
  2. +3
    -0
      Makefile.am
  3. +36
    -2
      README.md
  4. +8
    -0
      configure.ac
  5. +302
    -0
      doc/examples/qssp2017-test.inp
  6. +9
    -0
      src/Makefile.am
  7. +44
    -0
      src/banddpasss.f
  8. +22
    -0
      src/brunewvlet.f
  9. +37
    -0
      src/butterworth.f
  10. +20
    -0
      src/caxcb.f
  11. +107
    -0
      src/cdsvd500.f
  12. +15
    -0
      src/cmemcpy.f
  13. +73
    -0
      src/disazi.f
  14. +77
    -0
      src/four1w.f
  15. +36
    -0
      src/fsimpson.f
  16. +41
    -0
      src/moments.f
  17. +119
    -0
      src/qpalloc.f
  18. BIN
      src/qpalloc.mod
  19. +91
    -0
      src/qpdifmat0.f
  20. +116
    -0
      src/qpdifmatl.f
  21. +75
    -0
      src/qpdifmats.f
  22. +224
    -0
      src/qpfftinv.f
  23. +672
    -0
      src/qpgetinp.f
  24. +729
    -0
      src/qpgrnspec.f
  25. +28
    -0
      src/qplegendre.f
  26. +77
    -0
      src/qpmain.f
  27. +257
    -0
      src/qppsvkern.f
  28. +124
    -0
      src/qppsvkerng.f
  29. +56
    -0
      src/qpqmodel.f
  30. +90
    -0
      src/qpshkern.f
  31. +284
    -0
      src/qpsmat.f
  32. +124
    -0
      src/qpsmat0.f
  33. +142
    -0
      src/qpsmatc.f
  34. +1747
    -0
      src/qpsprop.f
  35. +212
    -0
      src/qpsprop0.f
  36. +2030
    -0
      src/qpspropg.f
  37. +195
    -0
      src/qpspropg0.f
  38. +40
    -0
      src/qpstart0.f
  39. +46
    -0
      src/qpstart0a.f
  40. +28
    -0
      src/qpstart0g.f
  41. +36
    -0
      src/qpstart2t.f
  42. +87
    -0
      src/qpstart4.f
  43. +65
    -0
      src/qpstart4a.f
  44. +51
    -0
      src/qpstart4g.f
  45. +123
    -0
      src/qpstart6.f
  46. +103
    -0
      src/qpstart6a.f
  47. +29
    -0
      src/qpstart6g.f
  48. +714
    -0
      src/qpsublayer.f
  49. +99
    -0
      src/qptmat.f
  50. +195
    -0
      src/qptprop.f
  51. +1099
    -0
      src/qpwvint.f
  52. +113
    -0
      src/ruku.f
  53. +20
    -0
      src/skipdoc.f
  54. +48
    -0
      src/spbdphj.f
  55. +48
    -0
      src/spbdphy.f
  56. +41
    -0
      src/spbh12.f
  57. +162
    -0
      src/spbjh.f
  58. +44
    -0
      src/spbphj.f
  59. +44
    -0
      src/spbphy.f
  60. +44
    -0
      src/spbpsj.f
  61. +44
    -0
      src/spbpsy.f
  62. +29
    -0
      src/swavelet.f
  63. +31
    -0
      src/taper.f
  64. +121
    -0
      src/transfs2t.f

+ 14
- 0
.gitignore View File

@ -0,0 +1,14 @@
/aclocal.m4
/autom4te.cache
/configure
/install-sh
/missing
/py-compile
/Makefile.in
/src/Makefile.in
/Makefile
/src/Makefile
/config.log
/config.status
*.o
/src/fomosto_qssp2017

+ 3
- 0
Makefile.am View File

@ -0,0 +1,3 @@
SUBDIRS = src
EXTRA_DIST = doc

+ 36
- 2
README.md View File

@ -1,2 +1,36 @@
# fomosto-qssp2017
QSSP with dynamic memory allocation and many improvements
# QSSP2017 (packaged as Fomosto backend)
Code for calculating complete synthetic seismograms of a spherical earth using
normal mode theory.
QSSP2017 has been written by Rongjiang Wang.
This is the newest version of the legacy code of QSSP.
New features are:
- dynamic memory allocation
- todo ...
Packaging has been done by Hannes Vasyura-Bathke.
## References
- Gilbert, F. and Backus, G.: Elastic-gravitational vibrations of a radially
stratified sphere, in: Dynamics of Stratified Solids, edited by: Herrmann,
G., American Society of Mechanical Engineers, New York, 82–95, 1968.
Takeuchi, H., and M. Saito (1972). Seismic surface waves, in Methods in
Computational Physics, vol. 11, Bolt, B. A. , Editor Academic Press, New
York, 217- 295.
- Wang, R., (1997), Tidal response of the solid earth, in Tidal Phenomena,
edited by H. Wilhelm, W. Zürn and H.G. Wenzel, Lecture Notes in Earth
Sciences, Vol. 66, pp. 27-57, Springer-Verlag, Berlin/Heidelberg, Germany,
1997.
## Compile and install
```
autoreconf -i # only if 'configure' script is missing
./configure
make
sudo make install
```

+ 8
- 0
configure.ac View File

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

+ 302
- 0
doc/examples/qssp2017-test.inp View File

@ -0,0 +1,302 @@
# This is the input file of FORTRAN77 program "qssp2017" for calculating
# synthetic seismograms of a self-gravitating, spherically symmetric,
# isotropic and viscoelastic earth.
#
# by
# Rongjiang Wang <wang@gfz-potsdam.de>
# Helmholtz-Centre Potsdam
# GFZ German Reseach Centre for Geosciences
# Telegrafenberg, D-14473 Potsdam, Germany
#
# Last modified: Potsdam, October 2017
#
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# If not specified, SI Unit System is used overall!
#
# Coordinate systems:
# spherical (r,t,p) with r = radial, t = co-latitude, and p = east longitude.
# local cartesian (e,n,z) with e = east, n = north, and z = vertical (upwards positve).
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
#
# UNIFORM RECEIVER DEPTH
# ======================
# 1. uniform receiver depth [km]
#--------------------------------------------------------------------------------------------------------
660.00
#--------------------------------------------------------------------------------------------------------
#
# SPACE-TIME SAMPLING PARAMETERS
# ==============================
# 1. time window [sec], sampling interval [sec]
# 2. max. frequency [Hz] of Green's functions
# 3. max. slowness [s/km] of Green's functions
# Note: if the near-field static displacement is desired, the maximum slowness should not
# be smaller than the S wave slowness in the receiver layer
# 4. anti-aliasing factor (> 0 & < 1), if it is <= 0 or >= 1/e (~ 0.4), then default value of
# 1/e is used (e.g., 0.1 = alias phases will be suppressed to 10% of their original
# amplitude)
# 5. switch (1/0 = yes/no) of turning-point filter, the range (d1, d2) of max. penetration
# depth [km] (d1 is meaningless if it is smaller than the receiver/source depth, and
# d2 is meaningless if it is equal to or larger than the earth radius)
#
# Note: The turning-point filter (Line 5) works only for the extended QSSP code (e.g.,
# qssp2016). if this filter is selected, all phases with the turning point
# shallower than d1 or deeper than d2 will be filtered.
#
# 6. Earth radius [km], switch of free-surface-reflection filter (1/0 = with/without free
# surface reflection)
#
# Note: The free-surface-reflection filter (Line 6) works only for the extended QSSP
# code (e.g., qssp2016). if this filter is selected, all phases with the turning
# point shallower than d1 or deeper than d2 will be filtered.
#--------------------------------------------------------------------------------------------------------
8191.0 1.00
0.25
0.40
0.01
0 2891.5 6371.0
6371.0 1
#--------------------------------------------------------------------------------------------------------
#
# SELF-GRAVITATING EFFECT
# =======================
# 1. the critical frequency [Hz] and the critical harmonic degree, below which
# the self-gravitating effect should be included
#--------------------------------------------------------------------------------------------------------
0.05 500
#--------------------------------------------------------------------------------------------------------
#
# WAVE TYPES
# ==========
# 1. selection (1/0 = yes/no) of spheroidal modes (P-SV waves), selection of toroidal modes
# (SH waves), minimum and maximum cutoff harmonic degrees
# Note: if the near-field static displacement is desired, the minimum cutoff harmonic
# degree should not be smaller than, e.g., 2000.
#--------------------------------------------------------------------------------------------------------
1 1 2000 10000
#--------------------------------------------------------------------------------------------------------
# GREEN'S FUNCTION FILES
# ======================
# 1. number of discrete source depths, estimated radius of each source patch [km] and
# directory for Green's functions
# 2. list of the source depths [km], the respective file names of the Green's functions
# (spectra) and the switch number (0/1) (0 = do not calculate this Green's function because
# it exists already, 1 = calculate or update this Green's function).
# Note: Green's functions need to be recalculated if any of the above parameters is changed.
#--------------------------------------------------------------------------------------------------------
1 10.0 './GreenFunctions/'
0.00 'Green_0km' 0
#--------------------------------------------------------------------------------------------------------
#
# MULTI-EVENT SOURCE PARAMETERS
# =============================
# 1. number of discrete point sources and selection of the source data format
# (1, 2 or 3)
# 2. list of the multi-event sources
#
# Format 1 (full moment tensor):
# Unit Mrr Mtt Mpp Mrt Mrp Mtp Lat Lon Depth T_origin T_rise
# [Nm] [deg] [deg] [km] [sec] [sec]
#
# Format 2 (double couple):
# Unit Strike Dip Rake Lat Lon Depth T_origin T_rise
# [Nm] [deg] [deg] [deg] [deg] [deg] [km] [sec] [sec]
#
# Format 3 (single force):
# Unit Feast Fnorth Fvertical Lat Lon Depth T_origin T_rise
# [N] [deg] [deg] [km] [sec] [sec]
#
# Note: for each point source, the default moment (force) rate time function is used, defined by a
# squared half-period (T_rise) sinusoid starting at T_origin.
#--------------------------------------------------------------------------------------------------------
1 3
1.0E+20 0.0 0.0 1.0 0.0 0.0 0.0 0.0 25.0
#--------------------------------------------------------------------------------------------------------
#
# RECEIVER PARAMETERS
# ===================
# 1. select output observables (1/0 = yes/no)
# Note: the gravity change defined here is space based, i.e., the effect due to free-air
# gradient and inertial are not included. the vertical component is positve upwards.
# 2. output file name
# 3. output time window [sec] (<= Green's function time window)
# 4. selection of order of Butterworth bandpass filter (if <= 0, then no filtering), lower
# and upper corner frequencies (<= cut-off frequency defined above)
# 5. lower and upper slowness cut-off [s/km] (slowness band-pass filter)
# 6. number of receiver
# 7. list of the station parameters
# Format:
# Lat Lon Name Time_reduction
# [deg] [deg] [sec]
# (Note: Time_reduction = start time of the time window)
#--------------------------------------------------------------------------------------------------------
# disp | velo | acce | strain | strain_rate | stress | stress_rate | rotation | rotation_rate | gravity
#--------------------------------------------------------------------------------------------------------
1 1 1 0 0 0 0 0 0 1
'test'
8000.0
0 0.05 0.25
0.00 0.40
5
25.0 25.0 'Loc1' 0.0
0.0 35.0 'Loc2' 0.0
35.0 0.0 'Loc3' 0.0
0.0 0.0 'Loc4' 0.0
0.0 180.0 'Loc5' 0.0
#--------------------------------------------------------------------------------------------------------
#
# MULTI-LAYERED EARTH MODEL (AK135)
# =================================
# 1. number of data lines of the layered model and selection for including
# the physical dispersion according to Kamamori & Anderson (1977)
#--------------------------------------------------------------------------------------------------------
142 0
#--------------------------------------------------------------------------------------------------------
#
# MODEL PARAMETERS
# ================
# no depth[km] vp[km/s] vs[km/s] ro[g/cm^3] qp qs
#--------------------------------------------------------------------------------------------------------
1 0.000 5.8000 3.2000 2.6000 1478.30 599.99
2 3.300 5.8000 3.2000 2.6000 1478.30 599.99
3 3.300 5.8000 3.2000 2.6000 1478.30 599.99
4 10.000 5.8000 3.2000 2.6000 1478.30 599.99
5 10.000 6.8000 3.9000 2.9200 1368.02 599.99
6 18.000 6.8000 3.9000 2.9200 1368.02 599.99
7 18.000 8.0355 4.4839 3.6410 950.50 394.62
8 43.000 8.0379 4.4856 3.5801 972.77 403.93
9 80.000 8.0400 4.4800 3.5020 1008.71 417.59
10 80.000 8.0450 4.4900 3.5020 182.03 75.60
11 120.000 8.0505 4.5000 3.4268 182.57 76.06
12 165.000 8.1750 4.5090 3.3711 188.72 76.55
13 210.000 8.3007 4.5184 3.3243 200.97 79.40
14 210.000 8.3007 4.5184 3.3243 338.47 133.72
15 260.000 8.4822 4.6094 3.3663 346.37 136.38
16 310.000 8.6650 4.6964 3.4110 355.85 139.38
17 360.000 8.8476 4.7832 3.4577 366.34 142.76
18 410.000 9.0302 4.8702 3.5068 377.93 146.57
19 410.000 9.3601 5.0806 3.9317 413.66 162.50
20 460.000 9.5280 5.1864 3.9273 417.32 164.87
21 510.000 9.6962 5.2922 3.9233 419.94 166.80
22 560.000 9.8640 5.3989 3.9218 422.55 168.78
23 610.000 10.0320 5.5047 3.9206 425.51 170.82
24 660.000 10.2000 5.6104 3.9201 428.69 172.93
25 660.000 10.7909 5.9607 4.2387 1350.54 549.45
26 710.000 10.9222 6.0898 4.2986 1311.17 543.48
27 760.000 11.0553 6.2100 4.3565 1277.93 537.63
28 809.500 11.1355 6.2424 4.4118 1269.44 531.91
29 859.000 11.2228 6.2799 4.4650 1260.68 526.32
30 908.500 11.3068 6.3164 4.5162 1251.69 520.83
31 958.000 11.3897 6.3519 4.5654 1243.02 515.46
32 1007.500 11.4704 6.3860 4.5926 1234.54 510.20
33 1057.000 11.5493 6.4182 4.6198 1226.52 505.05
34 1106.500 11.6265 6.4514 4.6467 1217.91 500.00
35 1156.000 11.7020 6.4822 4.6735 1210.02 495.05
36 1205.500 11.7768 6.5131 4.7001 1202.04 490.20
37 1255.000 11.8491 6.5431 4.7266 1193.99 485.44
38 1304.500 11.9208 6.5728 4.7528 1186.06 480.77
39 1354.000 11.9891 6.6009 4.7790 1178.19 476.19
40 1403.500 12.0571 6.6285 4.8050 1170.53 471.70
41 1453.000 12.1247 6.6554 4.8307 1163.16 467.29
42 1502.500 12.1912 6.6813 4.8562 1156.04 462.96
43 1552.000 12.2558 6.7070 4.8817 1148.76 458.72
44 1601.500 12.3181 6.7323 4.9069 1141.32 454.55
45 1651.000 12.3813 6.7579 4.9321 1134.01 450.45
46 1700.500 12.4427 6.7820 4.9570 1127.02 446.43
47 1750.000 12.5030 6.8056 4.9817 1120.09 442.48
48 1799.500 12.5638 6.8289 5.0062 1108.58 436.68
49 1849.000 12.6226 6.8517 5.0306 1097.16 431.03
50 1898.500 12.6807 6.8743 5.0548 1085.97 425.53
51 1948.000 12.7384 6.8972 5.0789 1070.38 418.41
52 1997.500 12.7956 6.9194 5.1027 1064.23 414.94
53 2047.000 12.8524 6.9416 5.1264 1058.03 411.52
54 2096.500 12.9093 6.9625 5.1499 1048.09 406.50
55 2146.000 12.9663 6.9852 5.1732 1042.07 403.23
56 2195.500 13.0226 7.0069 5.1963 1032.14 398.41
57 2245.000 13.0786 7.0286 5.2192 1018.38 392.16
58 2294.500 13.1337 7.0504 5.2420 1008.79 387.60
59 2344.000 13.1895 7.0722 5.2646 999.44 383.14
60 2393.500 13.2465 7.0932 5.2870 990.77 378.79
61 2443.000 13.3017 7.1144 5.3092 985.63 375.94
62 2492.500 13.3584 7.1368 5.3313 976.81 371.75
63 2542.000 13.4156 7.1584 5.3531 968.46 367.65
64 2591.500 13.4741 7.1804 5.3748 960.36 363.64
65 2640.000 13.5311 7.2031 5.3962 952.00 359.71
66 2690.000 13.5899 7.2253 5.4176 940.88 354.61
67 2740.000 13.6498 7.2485 5.4387 933.21 350.88
68 2740.000 13.6498 7.2485 5.6934 722.73 271.74
69 2789.670 13.6533 7.2593 5.7196 726.87 273.97
70 2839.330 13.6570 7.2700 5.7458 725.11 273.97
71 2891.500 13.6601 7.2817 5.7721 723.12 273.97
72 2891.500 8.0000 0.0000 9.9145 57822.00 0.00
73 2939.330 8.0382 0.0000 9.9942 57822.00 0.00
74 2989.660 8.1283 0.0000 10.0722 57822.00 0.00
75 3039.990 8.2213 0.0000 10.1485 57822.00 0.00
76 3090.320 8.3122 0.0000 10.2233 57822.00 0.00
77 3140.660 8.4001 0.0000 10.2964 57822.00 0.00
78 3190.990 8.4861 0.0000 10.3679 57822.00 0.00
79 3241.320 8.5692 0.0000 10.4378 57822.00 0.00
80 3291.650 8.6496 0.0000 10.5062 57822.00 0.00
81 3341.980 8.7283 0.0000 10.5731 57822.00 0.00
82 3392.310 8.8036 0.0000 10.6385 57822.00 0.00
83 3442.640 8.8761 0.0000 10.7023 57822.00 0.00
84 3492.970 8.9461 0.0000 10.7647 57822.00 0.00
85 3543.300 9.0138 0.0000 10.8257 57822.00 0.00
86 3593.640 9.0792 0.0000 10.8852 57822.00 0.00
87 3643.970 9.1426 0.0000 10.9434 57822.00 0.00
88 3694.300 9.2042 0.0000 11.0001 57822.00 0.00
89 3744.630 9.2634 0.0000 11.0555 57822.00 0.00
90 3794.960 9.3205 0.0000 11.1095 57822.00 0.00
91 3845.290 9.3760 0.0000 11.1623 57822.00 0.00
92 3895.620 9.4297 0.0000 11.2137 57822.00 0.00
93 3945.950 9.4814 0.0000 11.2639 57822.00 0.00
94 3996.280 9.5306 0.0000 11.3127 57822.00 0.00
95 4046.620 9.5777 0.0000 11.3604 57822.00 0.00
96 4096.950 9.6232 0.0000 11.4069 57822.00 0.00
97 4147.280 9.6673 0.0000 11.4521 57822.00 0.00
98 4197.610 9.7100 0.0000 11.4962 57822.00 0.00
99 4247.940 9.7513 0.0000 11.5391 57822.00 0.00
100 4298.270 9.7914 0.0000 11.5809 57822.00 0.00
101 4348.600 9.8304 0.0000 11.6216 57822.00 0.00
102 4398.930 9.8682 0.0000 11.6612 57822.00 0.00
103 4449.260 9.9051 0.0000 11.6998 57822.00 0.00
104 4499.600 9.9410 0.0000 11.7373 57822.00 0.00
105 4549.930 9.9761 0.0000 11.7737 57822.00 0.00
106 4600.260 10.0103 0.0000 11.8092 57822.00 0.00
107 4650.590 10.0439 0.0000 11.8437 57822.00 0.00
108 4700.920 10.0768 0.0000 11.8772 57822.00 0.00
109 4751.250 10.1095 0.0000 11.9098 57822.00 0.00
110 4801.580 10.1415 0.0000 11.9414 57822.00 0.00
111 4851.910 10.1739 0.0000 11.9722 57822.00 0.00
112 4902.240 10.2049 0.0000 12.0001 57822.00 0.00
113 4952.580 10.2329 0.0000 12.0311 57822.00 0.00
114 5002.910 10.2565 0.0000 12.0593 57822.00 0.00
115 5053.240 10.2745 0.0000 12.0867 57822.00 0.00
116 5103.570 10.2854 0.0000 12.1133 57822.00 0.00
117 5153.500 10.2890 0.0000 12.1391 57822.00 0.00
118 5153.500 11.0427 3.5043 12.7037 633.26 85.03
119 5204.610 11.0585 3.5187 12.7289 629.89 85.03
120 5255.320 11.0718 3.5314 12.7530 626.87 85.03
121 5306.040 11.0850 3.5435 12.7760 624.08 85.03
122 5356.750 11.0983 3.5551 12.7980 621.50 85.03
123 5407.460 11.1166 3.5661 12.8188 619.71 85.03
124 5458.170 11.1316 3.5765 12.8387 617.78 85.03
125 5508.890 11.1457 3.5864 12.8574 615.93 85.03
126 5559.600 11.1590 3.5957 12.8751 614.21 85.03
127 5610.310 11.1715 3.6044 12.8917 612.62 85.03
128 5661.020 11.1832 3.6126 12.9072 611.12 85.03
129 5711.740 11.1941 3.6202 12.9217 609.74 85.03
130 5762.450 11.2041 3.6272 12.9351 608.48 85.03
131 5813.160 11.2134 3.6337 12.9474 607.31 85.03
132 5863.870 11.2219 3.6396 12.9586 606.26 85.03
133 5914.590 11.2295 3.6450 12.9688 605.28 85.03
134 5965.300 11.2364 3.6498 12.9779 604.44 85.03
135 6016.010 11.2424 3.6540 12.9859 603.69 85.03
136 6066.720 11.2477 3.6577 12.9929 603.04 85.03
137 6117.440 11.2521 3.6608 12.9988 602.49 85.03
138 6168.150 11.2557 3.6633 13.0036 602.05 85.03
139 6218.860 11.2586 3.6653 13.0074 601.70 85.03
140 6269.570 11.2606 3.6667 13.0100 601.46 85.03
141 6320.290 11.2618 3.6675 13.0117 601.32 85.03
142 6371.000 11.2622 3.6678 13.0122 601.27 85.03
#---------------------------------end of all inputs------------------------------------------------------

+ 9
- 0
src/Makefile.am View File

@ -0,0 +1,9 @@
bin_PROGRAMS = fomosto_qssp2017
fomosto_qssp2017_SOURCES = banddpasss.f brunewvlet.f butterworth.f caxcb.f cdsvd500.f cmemcpy.f disazi.f \
four1w.f fsimpson.f moments.f qpalloc.f qpdifmat0.f qpdifmatl.f qpdifmats.f \
qpfftinv.f qpgetinp.f qpgrnspec.f qplegendre.f qpmain.f qppsvkern.f qppsvkerng.f \
qpqmodel.f qpshkern.f qpsmat0.f qpsmatc.f qpsmat.f qpsprop0.f qpsprop.f qpspropg0.f \
qpspropg.f qpstart0a.f qpstart0.f qpstart0g.f qpstart2t.f qpstart4a.f qpstart4.f \
qpstart4g.f qpstart6a.f qpstart6.f qpstart6g.f qpsublayer.f qptmat.f qptprop.f \
qpwvint.f ruku.f skipdoc.f spbdphj.f spbdphy.f spbh12.f spbjh.f spbphj.f spbphy.f \
spbpsj.f spbpsy.f swavelet.f taper.f transfs2t.f

+ 44
- 0
src/banddpasss.f View File

@ -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

+ 22
- 0
src/brunewvlet.f View File

@ -0,0 +1,22 @@
subroutine brunewvlet(tau,df,mm2,wvf)
implicit none
integer*4 mm2
real*8 tau,df,fi
complex*16 wvf(mm2)
c
integer*4 l
real*8 f
c
real*8 pi2
complex*16 c1
data pi2/6.28318530717959d0/
data c1/(1.d0,0.d0)/
c
c for wavelet: Brune's source time function
c
do l=1,mm2
f=df*dble(l-1)
wvf(l)=c1/dcmplx(1.d0,pi2*f*tau)**2
enddo
return
end

+ 37
- 0
src/butterworth.f View File

@ -0,0 +1,37 @@
subroutine butterworth(n,fc,df,nf,lpf)
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

+ 20
- 0
src/caxcb.f View File

@ -0,0 +1,20 @@
subroutine caxcb(a,b,n,l,m,c)
implicit none
c
c calculate c=a*b
c
integer*4 n,l,m
complex*16 a(n,l),b(l,m),c(n,m)
c
integer*4 i,j,k
c
do j=1,m
do i=1,n
c(i,j)=dcmplx(0.d0,0.d0)
do k=1,l
c(i,j)=c(i,j)+a(i,k)*b(k,j)
enddo
enddo
enddo
return
end

+ 107
- 0
src/cdsvd500.f View File

@ -0,0 +1,107 @@
subroutine cdsvd500(ca,cb,n,m,eps,key)
implicit none
c-------------------------------------------------------------------------------
c Solve complex linear equation system by single-value decomposiztion i
c method (modified from ludcmp and lubksb in the <Numerical Recipies> i
c ca: coefficient matrix(n,n); i
c cb: right-hand matrix(n,m) by input, i
c solution matrix(n,m) by return; i
c cunit: unit of the culomn vectors i
c eps: control constant; i
c key: if the main term of a column is i
c smaller than eps, key=0: anormal return, i
c else key=1: normal return. i
c i
c Note: n <= 500 will be NOT CHECKED! i
c-------------------------------------------------------------------------------
integer*4 n,m,key
real*8 eps
complex*16 ca(n,n),cb(n,m)
c
integer*4 NMAX
parameter (NMAX=500)
integer*4 i,ii,imax,j,k,ll
integer*4 indx(NMAX)
real*8 aamax,dum
real*8 vv(NMAX)
complex*16 cdum,csum
c
do i=1,n
aamax=0.d0
do j=1,n
aamax=dmax1(aamax,cdabs(ca(i,j)))
enddo
if(aamax.le.eps)then
key=0
return
endif
vv(i)=1.d0/aamax
enddo
do j=1,n
do i=1,j-1
csum=ca(i,j)
do k=1,i-1
csum=csum-ca(i,k)*ca(k,j)
enddo
ca(i,j)=csum
enddo
aamax=0.d0
do i=j,n
csum=ca(i,j)
do k=1,j-1
csum=csum-ca(i,k)*ca(k,j)
enddo
ca(i,j)=csum
dum=vv(i)*cdabs(csum)
if(dum.ge.aamax) then
imax=i
aamax=dum
endif
enddo
if(j.ne.imax) then
do k=1,n
cdum=ca(imax,k)
ca(imax,k)=ca(j,k)
ca(j,k)=cdum
enddo
vv(imax)=vv(j)
endif
indx(j)=imax
if(cdabs(ca(j,j)).le.eps)then
key=0
return
endif
if(j.ne.n) then
cdum=(1.d0,0.d0)/ca(j,j)
do i=j+1,n
ca(i,j)=ca(i,j)*cdum
enddo
endif
enddo
c
do k=1,m
ii=0
do i=1,n
ll=indx(i)
csum=cb(ll,k)
cb(ll,k)=cb(i,k)
if(ii.ne.0) then
do j=ii,i-1
csum=csum-ca(i,j)*cb(j,k)
enddo
else if(cdabs(csum).ne.0.d0)then
ii=i
endif
cb(i,k)=csum
enddo
do i=n,1,-1
csum=cb(i,k)
do j=i+1,n
csum=csum-ca(i,j)*cb(j,k)
enddo
cb(i,k)=csum/ca(i,i)
enddo
enddo
key=1
return
end

+ 15
- 0
src/cmemcpy.f View File

@ -0,0 +1,15 @@
subroutine cmemcpy(a,b,n)
implicit none
c
c copy complex array a to b
c
integer*4 n
complex*16 a(n),b(n)
c
integer*4 i
c
do i=1,n
b(i)=a(i)
enddo
return
end

+ 73
- 0
src/disazi.f View File

@ -0,0 +1,73 @@
subroutine disazi(rearth,lateq,loneq,latst,lonst,xnorth,yeast)
implicit none
c
real*8 rearth,lateq,loneq,latst,lonst,xnorth,yeast
c
integer*4 iangle
real*8 latb,lonb,latc,lonc,angleb,anglec
real*8 dis,a,b,c,s,aa
c
real*8 PI,PI2,DEGTORAD
data PI,PI2/3.14159265358979d0,6.28318530717959d0/
data DEGTORAD/1.745329251994328d-02/
c
c A is north pole
c B is station
c C is epicentre (local cartesian coordinate origin)
c
latb=lateq*DEGTORAD
lonb=loneq*DEGTORAD
latc=latst*DEGTORAD
lonc=lonst*DEGTORAD
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 station
c
xnorth=dis*dcos(anglec)
yeast=dis*dsin(anglec)
c
return
end

+ 77
- 0
src/four1w.f View File

@ -0,0 +1,77 @@
subroutine four1w(cdat,ddat,nn,isign)
implicit none
integer*4 nn,isign
complex*16 cdat(nn)
real*8 ddat(2*nn)
c
c fast fourier transform (fft)
c convention: f(t)=\int F(f)e^{-i2\pi ft} df
c replace ddat by its discrete fourier transform, if isign is input
c as 1; or replace data by nn times its inverse discrete fourier
c transform, if isign is input as -1. data is a double complex array of
c length nn or, equivalently, a double precision array of length 2*nn.
c nn must be an integer power of 2 (this is not checked!)
c
c note for convention: f(t)=\int F(f)e^{i2\pi f t} df, t-domain ==>
c f-domain, if isign=-1, and f-domain ==> t-domain, if isign=1.
c
integer*4 i,j,n,m,mmax,istep
real*8 tempr,tempi
real*8 wr,wi,wpr,wpi,wtemp,theta
c
do i=1,nn
ddat(2*i-1)=dreal(cdat(i))
ddat(2*i)=dimag(cdat(i))
enddo
c
n=2*nn
j=1
do i=1,n,2
if(j.gt.i)then
tempr=ddat(j)
tempi=ddat(j+1)
ddat(j)=ddat(i)
ddat(j+1)=ddat(i+1)
ddat(i)=tempr
ddat(i+1)=tempi
endif
m=n/2
1 if((m.ge.2).and.(j.gt.m))then
j=j-m
m=m/2
goto 1
endif
j=j+m
enddo
mmax=2
2 if(n.gt.mmax)then
istep=2*mmax
theta=6.28318530717959d0/dble(isign*mmax)
wpr=-2.d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.d0
wi=0.d0
do m=1,mmax,2
do i=m,n,istep
j=i+mmax
tempr=wr*ddat(j)-wi*ddat(j+1)
tempi=wr*ddat(j+1)+wi*ddat(j)
ddat(j)=ddat(i)-tempr
ddat(j+1)=ddat(i+1)-tempi
ddat(i)=ddat(i)+tempr
ddat(i+1)=ddat(i+1)+tempi
enddo
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
enddo
mmax=istep
goto 2
endif
c
do i=1,nn
cdat(i)=dcmplx(ddat(2*i-1),ddat(2*i))
enddo
c
return
end

+ 36
- 0
src/fsimpson.f View File

@ -0,0 +1,36 @@
double precision function fsimpson(xa,xb,c)
implicit none
real*8 xa,xb,c
c
c integrate dsqrt(c^2-x^2)/x from xa to xb
c
integer*4 i,n
real*8 x,dx,y,f1,f2
c
integer*4 nmax
real*8 eps
data nmax/2048/
data eps/1.0d-03/
c
dx=(xb-xa)/dble(n)
c
n=1
dx=xb-xa
y=0.5d0*(dsqrt(c*c-xa*xa)/xa+dsqrt(c*c-xb*xb)/xb)
f1=0.d0
f2=y*dx
10 n=2*n
dx=0.5d0*dx
do i=1,n-1,2
x=xa+dble(i)*dx
y=y+dsqrt(c*c-x*x)/x
enddo
f2=y*dx
if(dabs(f2-f1).gt.eps*f2.and.n.lt.nmax)then
f1=f2
goto 10
endif
fsimpson=f2
c
return
end

+ 41
- 0
src/moments.f View File

@ -0,0 +1,41 @@
subroutine moments(munit,strike,dip,rake,
& mtt,mpp,mrr,mtp,mpr,mrt)
implicit none
real*8 munit,strike,dip,rake,
& mtt,mpp,mrr,mtp,mpr,mrt
c
real*8 DEG2RAD
parameter(DEG2RAD=1.745329251994328d-02)
c
real*8 st,di,ra
real*8 sss,sss2,ss2s,css,css2,cs2s
real*8 ssd,ss2d,csd,cs2d,ssr,csr
c
st=strike*DEG2RAD
di=dip*DEG2RAD
ra=rake*DEG2RAD
c
sss=dsin(st)
sss2=sss*sss
ss2s=dsin(2.d0*st)
css=dcos(st)
css2=css*css
cs2s=dcos(2.d0*st)
c
ssd=dsin(di)
ss2d=dsin(2.d0*di)
csd=dcos(di)
cs2d=dcos(2.d0*di)
c
ssr=dsin(ra)
csr=dcos(ra)
c
mtt=munit*(-ssd*csr*ss2s-ss2d*ssr*sss2)
mpp=munit*(ssd*csr*ss2s-ss2d*ssr*css2)
mrr=-(mtt+mpp)
mtp=munit*(-ssd*csr*cs2s-0.5d0*ss2d*ssr*ss2s)
mpr=munit*(csd*csr*sss-cs2d*ssr*css)
mrt=munit*(-csd*csr*css-cs2d*ssr*sss)
c
return
end

+ 119
- 0
src/qpalloc.f View File

@ -0,0 +1,119 @@
module qpalloc
c
c CONSTANTS
c =========
integer*4 ndmax
parameter(ndmax=2)
real*8 PI,PI2
parameter(PI=3.14159265358979d0,PI2=6.28318530717959d0)
real*8 DEG2RAD,KM2M
parameter(DEG2RAD=1.745329251994328d-02,KM2M=1.0d+03)
real*8 BIGG
parameter(BIGG=6.6732d-11)
real*8 RESOLUT
parameter(RESOLUT=0.01d0)
real*8 FLTAPER
parameter(FLTAPER=0.2d0)
real*8 FSBUP,FSBLW,FSBREF
parameter(FSBUP=1.0d+02,FSBLW=2.5d-04,FSBREF=1.0d+00)
real*8 fbvatm,fbvocean,fbvcore
parameter(fbvatm=1.0d-06,fbvocean=0.0d+00,fbvcore=0.0d+00)
c
logical*2 selpsv,selsh,setzh12a,setzh12t
logical*2 nogravity,freesurf
logical*2 dispersion
c
integer*4 ngrn,nt,ntcut,ntcutout,nf,nfcut,nlpf
integer*4 lyadd,ipatha,ipathb,ldeggr,ldegmin,ldegcut,ldegmax
integer*4 nr,igfirst,iglast
integer*4 ns,l0,lymax
integer*4 lys,lyr,lylwa,lylwb,lyos,lyob,lycm,lycc,ly0
integer*4 icmp(11)
c
real*8 dt,df,fi,fcut,fgr,rratmos,depatmos
real*8 rearth,rr0,minpath,maxpath
real*8 slwmax,slwlwcut,slwupcut,f1corner,f2corner
real*8 dpr,freeairgrd
real*8 qsmin,togsmin
c
complex*16 comi,comi2
c
character*80 dispout(3),veloout(3),acceout(3),
& rotaout(3),rotarateout(3),
& strainout(6),strainrateout(6),
& stressout(6),stressrateout(6),
& gravout(3),grmout
c
logical*2,allocatable:: ksmallp(:),ksmalls(:),ksmallt(:)
c
integer*4,allocatable:: lygrn(:),grnsel(:),ldegpsv(:),ldegsh(:)
integer*4,allocatable:: nruku(:,:)
integer*4,allocatable:: isg1(:),isg2(:),nsg(:)
integer*4,allocatable:: idr(:,:)
integer*4,allocatable:: lylwp(:),lylws(:),lylwt(:),
& lyupp(:),lyups(:),lyupt(:)
integer*4,allocatable:: ldegtap(:,:,:)
c
real*8,allocatable:: latr(:),lonr(:),tred(:),grndep(:)
real*8,allocatable:: sfr(:),sft(:),sfp(:),
& mtt(:),mpp(:),mrr(:),mtp(:),mpr(:),mrt(:),
& lats(:),lons(:),deps(:),togs(:),trss(:)
real*8,allocatable:: dis(:,:),plm(:,:)
real*8,allocatable:: dp0(:),dp0up(:),dp0lw(:),
& vp0(:),vp0up(:),vp0lw(:),
& vs0(:),vs0up(:),vs0lw(:),
& ro0(:),ro0up(:),ro0lw(:),
& qp0(:),qp0up(:),qp0lw(:),
& qs0(:),qs0up(:),qs0lw(:)
real*8,allocatable:: rrup(:),rrlw(:),vpup(:),vplw(:),
& vsup(:),vslw(:),roup(:),rolw(:),
& qpup(:),qplw(:),qsup(:),qslw(:)
real*8,allocatable:: xp(:),xs(:),xt(:)
real*8,allocatable:: tap(:),disk(:)
c
complex*16,allocatable:: ul0(:,:),vl0(:,:),wl0(:,:),
& el0(:,:),fl0(:,:),gl0(:,:),
& pl0(:,:),ql0(:,:)
complex*16,allocatable:: urlm(:,:,:),utlm(:,:,:),uplm(:,:,:),
& errlm(:,:,:),ertlm(:,:,:),erplm(:,:,:),etrlm(:,:,:),
& ett0lm(:,:,:),ettalm(:,:,:),ettblm(:,:,:),etp0lm(:,:,:),
& etpalm(:,:,:),etpblm(:,:,:),eprlm(:,:,:),ept0lm(:,:,:),
& eptalm(:,:,:),eptblm(:,:,:),epp0lm(:,:,:),eppalm(:,:,:),
& eppblm(:,:,:),grlm(:,:,:),gtlm(:,:,:),gplm(:,:,:)
complex*16,allocatable:: ue(:,:),un(:,:),uz(:,:),
& roe(:,:),ron(:,:),roz(:,:),
& uee(:,:),uen(:,:),uez(:,:),
& unn(:,:),unz(:,:),uzz(:,:),
& see(:,:),sen(:,:),sez(:,:),
& snn(:,:),snz(:,:),szz(:,:),
& ge(:,:),gn(:,:),gz(:,:),gm(:,:)
complex*16,allocatable:: ssa(:,:),ss2a(:,:),csa(:,:),cs2a(:,:),
& ssb(:,:),csb(:,:),ssd(:,:),csd(:,:),ssf(:,:)
complex*16,allocatable:: crrup(:),crrlw(:),
& cro(:),croup(:),crolw(:),cla(:),claup(:),clalw(:),
& cmu(:),cmuup(:),cmulw(:),cvp(:),cvpup(:),cvplw(:),
& cvs(:),cvsup(:),cvslw(:),cgr(:),cgrup(:),cgrlw(:),
& cga(:),cgaup(:),cgalw(:)
complex*16,allocatable:: cps(:,:),cpt(:,:),kp(:),ks(:),kt(:)
complex*16,allocatable:: mat2x2up(:,:,:),mat2x2lw(:,:,:),
& mat2x2inv(:,:,:),mas3x3up(:,:,:),mas3x3lw(:,:,:),
& mas3x3inv(:,:,:),mas4x4up(:,:,:),mas4x4lw(:,:,:),
& mas4x4inv(:,:,:),mas6x6up(:,:,:),mas6x6lw(:,:,:),
& mas6x6inv(:,:,:)
complex*16,allocatable:: zjup(:,:,:),zjlw(:,:,:),
& zhup(:,:,:),zhlw(:,:,:),wj(:,:,:),wh(:,:,:),
& zh12p(:,:,:),zh12sv(:,:,:),zh12sh(:,:,:)
complex*16,allocatable:: zjupg(:),zjlwg(:),
& zhupg(:),zhlwg(:),wjg(:),whg(:)
complex*16,allocatable:: cua(:,:),cypnorm(:,:)
complex*16,allocatable:: wvf(:,:)
complex*16,allocatable:: sf1(:),sf2(:),sf3(:),
& expl(:),clvd(:),ss12(:),
& ss11(:),ds31(:),ds23(:)
c
character*80,allocatable:: specfile(:),uspecfile(:),vspecfile(:),
& wspecfile(:),especfile(:),fspecfile(:),gspecfile(:),
& pspecfile(:),qspecfile(:)
character*10,allocatable:: rname(:)
c
end module

BIN
src/qpalloc.mod View File


+ 91
- 0
src/qpdifmat0.f View File

@ -0,0 +1,91 @@
subroutine qpdifmat0(ly,ldeg,rr,mat)
use qpalloc
implicit none
integer ly,ldeg
real*8 rr
complex*16 mat(6,6)
c
c 3x3 coefficient matrix for spheroidal mode l = 0
c
real*8 mass,rorr,beta,dro,ro1
complex*16 cup,clw,crr,crorr,cxirr,clarr,cmurr,cgrrr,cvprr
complex*16 cgarr,gamma
c
complex*16 c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
crr=dcmplx(rr,0.d0)
cup=(crr-crrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
clarr=cup*claup(ly)+clw*clalw(ly)
cmurr=cup*cmuup(ly)+clw*cmulw(ly)
cvprr=cup*cvpup(ly)+clw*cvplw(ly)
cxirr=clarr+c2*cmurr
beta=dlog(roup(ly)/rolw(ly))/(rrup(ly)-rrlw(ly))
c
if(rr.le.rrlw(ly))then
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else if(ly.ge.lyos)then
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rolw(ly)-dro*rrlw(ly)
mass=PI*(rr-rrlw(ly))*((4.d0/3.d0)*ro1
& *(rr**2+rr*rrlw(ly)+rrlw(ly)**2)
& +dro*(rr**3+rr**2*rrlw(ly)
& +rr*rrlw(ly)**2+rrlw(ly)**3))
else
rorr=rolw(ly)*dexp(dlog(roup(ly)/rolw(ly))
& *(rr-rrlw(ly))/(rrup(ly)-rrlw(ly)))
crorr=dcmplx(rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
mass=2.d0*PI2/beta**3
& *(rorr*(beta*rr*(beta*rr-2.d0)+2.d0)
& -rolw(ly)*(beta*rrlw(ly)*(beta*rrlw(ly)-2.d0)+2.d0))
endif
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
c
c stabilization at the static limit
c
if(ly.lt.lyos)then
if(cdabs(comi)/PI2.lt.fbvatm)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvatm))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else if(ly.ge.lyos.and.ly.lt.lyob)then
if(cdabs(comi)/PI2.lt.fbvocean)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvocean))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else if(ly.ge.lycm.and.ly.lt.lycc)then
if(cdabs(comi)/PI2.lt.fbvcore)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvcore))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
endif
c
mat(1,1)=(c1-c2*clarr/cxirr)/crr
mat(1,2)=c1/(crorr*cvprr**2)/crr
mat(1,3)=(0.d0,0.d0)
c
mat(2,1)=c4*(cmurr*(c1+c2*clarr/cxirr)/crr-crorr*cgrrr)
& -crorr*comi2*crr
mat(2,2)=c2*clarr/cxirr/crr
mat(2,3)=(0.d0,0.d0)
c
mat(3,1)=cgarr/crr
mat(3,2)=(0.d0,0.d0)
mat(3,3)=(0.d0,0.d0)
c
return
end

+ 116
- 0
src/qpdifmatl.f View File

@ -0,0 +1,116 @@
subroutine qpdifmatl(ly,ldeg,rr,mat)
use qpalloc
implicit none
integer*4 ly,ldeg
real*8 rr
complex*16 mat(6,6)
c
c 4x4 coefficient matrix for spheroidal mode l > 0 in liquid media
c
real*8 mass,rorr,beta,dro,ro1
complex*16 cldeg,clp1,cll1
complex*16 cup,clw,crr,crorr,cgrrr,cvprr
complex*16 cgarr,gamma,cn2rr
logical*2 bv
c
complex*16 c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
cldeg=dcmplx(dble(ldeg),0.d0)
clp1=cldeg+c1
cll1=cldeg*clp1
c
crr=dcmplx(rr,0.d0)
cup=(crr-crrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
cvprr=cup*cvpup(ly)+clw*cvplw(ly)
beta=dlog(roup(ly)/rolw(ly))/(rrup(ly)-rrlw(ly))
c
if(rr.le.rrlw(ly))then
mass=0.d0
crorr=crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
else if(ly.ge.lyos)then
crorr=cup*croup(ly)+clw*crolw(ly)
rorr=dreal(crorr)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
dro=(rorr-rolw(ly))/(rr-rrlw(ly))
ro1=rorr-dro*rrlw(ly)
mass=PI*(rr-rrlw(ly))*((4.d0/3.d0)*ro1
& *(rr**2+rr*rrlw(ly)+rrlw(ly)**2)
& +dro*(rr**3+rr**2*rrlw(ly)
& +rr*rrlw(ly)**2+rrlw(ly)**3))
else
rorr=rolw(ly)*dexp(dlog(roup(ly)/rolw(ly))
& *(rr-rrlw(ly))/(rrup(ly)-rrlw(ly)))
crorr=dcmplx(rorr,0.d0)
cgarr=dcmplx(2.d0*PI2*BIGG*rorr,0.d0)
c
mass=2.d0*PI2/beta**3
& *(rorr*(beta*rr*(beta*rr-2.d0)+2.d0)
& -rolw(ly)*(beta*rrlw(ly)*(beta*rrlw(ly)-2.d0)+2.d0))
endif
cgrrr=cgrlw(ly)*(crrlw(ly)/crr)**2
& +dcmplx(BIGG*mass/rr**2,0.d0)
c
c stabilization at the static limit
c
if(ly.lt.lyos)then
bv=fbvatm.gt.0.d0
if(cdabs(comi)/PI2.lt.fbvatm)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvatm))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else if(ly.ge.lyos.and.ly.lt.lyob)then
bv=fbvocean.gt.0.d0
if(cdabs(comi)/PI2.lt.fbvocean)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvocean))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else if(ly.ge.lycm.and.ly.lt.lycc)then
bv=fbvcore.gt.0.d0
if(cdabs(comi)/PI2.lt.fbvcore)then
gamma=dcmplx((dsin(0.25d0*cdabs(comi)/fbvcore))**3,0.d0)
cvprr=cvprr/cdsqrt(gamma
& +(gamma-c1)*dcmplx(beta,0.d0)*cvprr**2/cgrrr)
endif
else
stop 'Error in qpdifmatl: unknown fiquid layer!'
endif
c
c Brunt-Väisälä frequency
c
cn2rr=-cgrrr*(dcmplx(beta,0.d0)+cgrrr/cvprr**2)
c
mat(1,1)=cgrrr/cvprr**2-c1/crr
mat(1,2)=cll1/crr-comi2*crr/cvprr**2
mat(1,3)=-crr/cvprr**2
mat(1,4)=(0.d0,0.d0)
c
if(bv.and.cdabs(comi).gt.0.d0)then
mat(2,1)=(c1-cn2rr/comi2)/crr
mat(2,2)=cn2rr/cgrrr
mat(2,3)=cn2rr/(comi2*cgrrr)
mat(2,4)=(0.d0,0.d0)
else
mat(2,1)=c1/crr
mat(2,2)=(0.d0,0.d0)
mat(2,3)=(0.d0,0.d0)
mat(2,4)=(0.d0,0.d0)
endif
c
mat(3,1)=cgarr/crr
mat(3,2)=(0.d0,0.d0)
mat(3,3)=-clp1/crr
mat(3,4)=c1/crr
c
mat(4,1)=cgarr*clp1/crr
mat(4,2)=-cgarr*cll1/crr
mat(4,3)=(0.d0,0.d0)
mat(4,4)=cldeg/crr
return
end

+ 75
- 0
src/qpdifmats.f View File

@ -0,0 +1,75 @@
subroutine qpdifmats(ly,ldeg,rr,mat)
use qpalloc
implicit none
integer*4 ly,ldeg
real*8 rr
complex*16 mat(6,6)
c
c 6x6 coefficient matrix for spheroidal mode l > 0 in solid media
c
complex*16 cldeg,clp1,cll1,cup,clw
complex*16 crr,drr,crorr,clarr,cmurr,cxirr,cgarr,cgrrr
c
complex*16 c1,c2,c4
data c1,c2,c4/(1.d0,0.d0),(2.d0,0.d0),(4.d0,0.d0)/
c
cldeg=dcmplx(dble(ldeg),0.d0)
clp1=cldeg+c1
cll1=cldeg*clp1
crr=dcmplx(rr,0.d0)
cup=(crr-rrlw(ly))/(crrup(ly)-crrlw(ly))
clw=c1-cup
crorr=cup*croup(ly)+clw*crolw(ly)
clarr=cup*claup(ly)+clw*clalw(ly)
cmurr=cup*cmuup(ly)+clw*cmulw(ly)
cgarr=cup*cgaup(ly)+clw*cgalw(ly)
cgrrr=cup*cgrup(ly)+clw*cgrlw(ly)
cxirr=clarr+c2*cmurr
c
mat(1,1)=(c1-c2*clarr/cxirr)/crr
mat(1,2)=c1/cxirr/crr
mat(1,3)=cll1*clarr/cxirr/crr
mat(1,4)=(0.d0,0.d0)
mat(1,5)=(0.d0,0.d0)
mat(1,6)=(0.d0,0.d0)
c
mat(2,1)=c4*(cmurr*(c1+c2*clarr/cxirr)/crr-crorr*cgrrr)
& -crorr*comi2*crr
mat(2,2)=c2*clarr/cxirr/crr
mat(2,3)=cll1*(crorr*cgrrr
& -c2*cmurr*(c1+c2*clarr/cxirr)/crr)
mat(2,4)=cll1/crr
mat(2,5)=crorr*clp1*crr
mat(2,6)=-crorr*crr
c
mat(3,1)=-c1/crr
mat(3,2)=(0.d0,0.d0)
mat(3,3)=c2/crr
mat(3,4)=c1/cmurr/crr
mat(3,5)=(0.d0,0.d0)
mat(3,6)=(0.d0,0.d0)
c
mat(4,1)=crorr*cgrrr-c2*cmurr*(c1+c2*clarr/cxirr)/crr
mat(4,2)=-clarr/cxirr/crr
mat(4,3)=c2*cmurr*(c2*cll1*(c1-cmurr/cxirr)-c1)/crr
& -crorr*comi2*crr
mat(4,4)=-c1/crr
mat(4,5)=-crorr*crr
mat(4,6)=(0.d0,0.d0)
c
mat(5,1)=cgarr/crr
mat(5,2)=(0.d0,0.d0)
mat(5,3)=(0.d0,0.d0)
mat(5,4)=(0.d0,0.d0)
mat(5,5)=-clp1/crr
mat(5,6)=c1/crr
c
mat(6,1)=cgarr*clp1/crr
mat(6,2)=(0.d0,0.d0)
mat(6,3)=-cgarr*cll1/crr
mat(6,4)=(0.d0,0.d0)
mat(6,5)=(0.d0,0.d0)
mat(6,6)=cldeg/crr
c
return
end

+ 224
- 0
src/qpfftinv.f View File

@ -0,0 +1,224 @@
subroutine qpfftinv(ierr)
use qpalloc
implicit none
integer*4 ierr
c
integer*4 lf,mf,ir
real*8 f
c
real*8,allocatable:: y0(:),y(:),dswap(:)
complex*16,allocatable:: cy(:,:),lpf(:),cswap(:)
c
allocate(y0(nr),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: y0 not allocated!'
allocate(y(nr),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: y not allocated!'
allocate(cy(nf,nr),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: cy not allocated!'
allocate(lpf(nf),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: lpf not allocated!'
allocate(dswap(4*nf),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: dswap not allocated!'
allocate(cswap(2*nf),stat=ierr)
if(ierr.ne.0)stop 'Error in qpfftinv: cswap not allocated!'
c
do lf=1,nf
f=dble(lf-1)*df
comi=dcmplx(PI2*f,PI2*fi)
comi2=comi**2
do ir=1,nr
gm(lf,ir)=-gz(lf,ir)+(dcmplx(freeairgrd,0.d0)-comi2)*uz(lf,ir)
enddo
enddo
c
c low-pass filter
c
if(nlpf.gt.0)then
if(f1corner.le.0.d0)then
call butterworth(nlpf,f2corner,df,nf,lpf)
else
call bandpass(nlpf,f1corner,f2corner,df,nf,lpf)
endif
do lf=1,nf
do ir=1,nr
ue(lf,ir)=ue(lf,ir)*lpf(lf)
un(lf,ir)=un(lf,ir)*lpf(lf)
uz(lf,ir)=uz(lf,ir)*lpf(lf)
c
ge(lf,ir)=ge(lf,ir)*lpf(lf)
gn(lf,ir)=gn(lf,ir)*lpf(lf)
gz(lf,ir)=gz(lf,ir)*lpf(lf)
c
roe(lf,ir)=roe(lf,ir)*lpf(lf)
ron(lf,ir)=ron(lf,ir)*lpf(lf)
roz(lf,ir)=roz(lf,ir)*lpf(lf)
c
uee(lf,ir)=uee(lf,ir)*lpf(lf)
uen(lf,ir)=uen(lf,ir)*lpf(lf)
uez(lf,ir)=uez(lf,ir)*lpf(lf)
unn(lf,ir)=unn(lf,ir)*lpf(lf)
unz(lf,ir)=unz(lf,ir)*lpf(lf)
uzz(lf,ir)=uzz(lf,ir)*lpf(lf)
c
see(lf,ir)=see(lf,ir)*lpf(lf)
sen(lf,ir)=sen(lf,ir)*lpf(lf)
sez(lf,ir)=sez(lf,ir)*lpf(lf)
snn(lf,ir)=snn(lf,ir)*lpf(lf)
snz(lf,ir)=snz(lf,ir)*lpf(lf)
szz(lf,ir)=szz(lf,ir)*lpf(lf)
c
gm(lf,ir)=gm(lf,ir)*lpf(lf)
enddo
enddo
endif
c
if(icmp(1).eq.1)then
c
c output seismograms: displacement vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ue,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,dispout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,un,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,dispout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,dispout(3))
endif
if(icmp(2).eq.1)then
c
c output seismograms: velocity vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ue,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,veloout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,un,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,veloout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,veloout(3))
endif
if(icmp(3).eq.1)then
c
c output seismograms: acceleration vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ue,cy,cswap,dswap,
& 2,ntcutout,rname,y0,y,acceout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,un,cy,cswap,dswap,
& 2,ntcutout,rname,y0,y,acceout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uz,cy,cswap,dswap,
& 2,ntcutout,rname,y0,y,acceout(3))
endif
c
if(icmp(4).eq.1)then
c
c output seismograms: strain tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uee,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uen,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uez,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unn,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uzz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,strainout(6))
endif
c
if(icmp(5).eq.1)then
c
c output seismograms: strain rate tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uee,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uen,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uez,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unn,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,unz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,uzz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,strainrateout(6))
endif
if(icmp(6).eq.1)then
c
c output seismograms: stress tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,see,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sen,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sez,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snn,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,szz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,stressout(6))
endif
c
if(icmp(7).eq.1)then
c
c output seismograms: stress rate tensor
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,see,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sen,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,sez,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(3))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snn,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(4))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,snz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(5))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,szz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,stressrateout(6))
endif
c
if(icmp(8).eq.1)then
c
c output seismograms: rotation vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roe,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,rotaout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ron,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,rotaout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,rotaout(3))
endif
c
if(icmp(9).eq.1)then
c
c output seismograms: rotation rate vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roe,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,rotarateout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ron,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,rotarateout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,roz,cy,cswap,dswap,
& 1,ntcutout,rname,y0,y,rotarateout(3))
endif
if(icmp(10).eq.1)then
c
c output seismograms: local gravitational acceleration vector
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,ge,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,gravout(1))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,gn,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,gravout(2))
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,gz,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,gravout(3))
endif
c
if(icmp(11).eq.1)then
c
c output seismograms: gravimeter response
c
call transfs2t(nf,nr,df,togsmin,dt,fi,tred,gm,cy,cswap,dswap,
& 0,ntcutout,rname,y0,y,grmout)
endif
return
end

+ 672
- 0
src/qpgetinp.f View File

@ -0,0 +1,672 @@
subroutine qpgetinp(unit)
use qpalloc
implicit none
integer unit
c
c work space
c
integer*4 i,j,l,ir,ig,isg,is,is1,flen,iswap,nhypo,sdfsel,ierr
integer*4 ipath,isurf
real*8 twindow,twinout,suppress,munit,sfe,sfn,sfz
real*8 strike,dip,rake,depdif,dswap(14)
character*80 grndir,outfile,fswap
c
c uniform receiver depth
c ======================
c
lyadd=1
c
call skipdoc(unit)
read(unit,*)dpr
dpr=KM2M*dpr
c
c time (frequency) sampling
c =========================
c
call skipdoc(unit)
read(unit,*)twindow,dt
ntcut=1+idnint(twindow/dt)
nt=2
100 nt=2*nt
if(nt.lt.ntcut)goto 100
nf=nt/2
df=1.d0/(dble(nt)*dt)
c
call skipdoc(unit)
read(unit,*)fcut
nfcut=min0(nf,1+idnint(fcut/df))
fcut=dble(nfcut-1)*df
call skipdoc(unit)
read(unit,*)slwmax
if(slwmax.le.0.d0)then
stop ' Error in qpgetinp: bad selection of max. slowness!'
else
slwmax=slwmax/KM2M
endif
c
call skipdoc(unit)
read(unit,*)suppress
if(suppress.le.0.d0.or.suppress.ge.1.d0)then
suppress=dexp(-1.d0)
endif
fi=dlog(suppress)*df/PI2
call skipdoc(unit)
read(unit,*)ipath,minpath,maxpath
minpath=minpath*KM2M
maxpath=maxpath*KM2M
call skipdoc(unit)
read(unit,*)rearth,isurf
rearth=rearth*KM2M
freesurf=isurf.eq.1
if(ipath.ne.1)then
ipatha=0
ipathb=0
minpath=0.d0
maxpath=rearth
else
if(minpath.lt.dpr.or.minpath.ge.dmax1(rearth,maxpath))then
ipatha=0
minpath=0.d0
else
ipatha=1
lyadd=lyadd+1
endif
if(maxpath.ge.rearth.or.maxpath.le.minpath)then
ipathb=0
maxpath=rearth
else
ipathb=1
lyadd=lyadd+1
endif
endif
c
c cutoffs of spectra
c ==================
c
call skipdoc(unit)
read(unit,*)fgr,ldeggr
if(fgr.lt.0.d0)fgr=0.d0
if(ldeggr.lt.0)ldeggr=0
if(fgr.gt.0.d0.and.ldeggr.le.0.or.
& fgr.le.0.d0.and.ldeggr.gt.0)then
stop ' Bad fgr and ldeggr combination!'
endif
nogravity=fgr*dble(ldeggr).le.0.d0
c
call skipdoc(unit)
read(unit,*)i,j,ldegmin,ldegcut
selpsv=i.eq.1
selsh=j.eq.1
if(.not.(selpsv.or.selsh))then
stop ' Error in qpgetinp: none of PSV and SH is selected!'
endif
if(ldegcut.lt.ldegmin)ldegcut=ldegmin
ldegmin=min0(max0(1+ndmax,ldegmin),ldegcut)
ldegmax=ldegcut+ndmax+1
c
allocate(plm(0:ldegmax,0:2),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: plm not allocated!'
c
allocate(ul0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ul0 not allocated!'
allocate(vl0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: vl0 not allocated!'
allocate(wl0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: wl0 not allocated!'
allocate(el0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: el0 not allocated!'
allocate(fl0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: fl0 not allocated!'
allocate(gl0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: gl0 not allocated!'
c
allocate(pl0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: pl0 not allocated!'
allocate(ql0(0:ldegmax,6),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ql0 not allocated!'
c
allocate(urlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: urlm not allocated!'
allocate(utlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: utlm not allocated!'
allocate(uplm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: uplm not allocated!'
c
allocate(grlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: grlm not allocated!'
allocate(gtlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: gtlm not allocated!'
allocate(gplm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: gplm not allocated!'
c
allocate(errlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: errlm not allocated!'
allocate(ertlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ertlm not allocated!'
allocate(erplm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: erplm not allocated!'
allocate(etrlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: etrlm not allocated!'
allocate(ett0lm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ett0lm not allocated!'
allocate(ettalm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ettalm not allocated!'
allocate(ettblm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ettblm not allocated!'
allocate(etp0lm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: etp0lm not allocated!'
allocate(etpalm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: etpalm not allocated!'
allocate(etpblm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: etpblm not allocated!'
allocate(eprlm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: eprlm not allocated!'
allocate(ept0lm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: ept0lm not allocated!'
allocate(eptalm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: eptalm not allocated!'
allocate(eptblm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: eptblm not allocated!'
allocate(epp0lm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: epp0lm not allocated!'
allocate(eppalm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: eppalm not allocated!'
allocate(eppblm(0:ldegmax,6,0:ndmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgetinp: eppblm not allocated!'
c
allocate(lyupp(0:ldegmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lyupp not allocated!'
allocate(lyups(0:ldegmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lyups not allocated!'
allocate(lyupt(0:ldegmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lyupt not allocated!'
allocate(lylwp(0:ldegmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lylwp not allocated!'
allocate(lylws(0:ldegmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lylws not allocated!'
allocate(lylwt(0:ldegmax),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lylwt not allocated!'
c
c Green's function files
c ======================
c
call skipdoc(unit)
read(unit,*)ngrn,rr0,grndir
if(ngrn.le.0)then
stop ' bad number of source depths!'
endif
rr0=rr0*KM2M
lyadd=lyadd+ngrn
c
allocate(grndep(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: grndep not allocated!'
allocate(grnsel(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: grnsel not allocated!'
allocate(lygrn(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lygrn not allocated!'
c
allocate(specfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: specfile not allocated!'
allocate(uspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: uspecfile not allocated!'
allocate(vspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: vspecfile not allocated!'
allocate(wspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: wspecfile not allocated!'
allocate(especfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: especfile not allocated!'
allocate(fspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: fspecfile not allocated!'
allocate(gspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: gspecfile not allocated!'
allocate(pspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: pspecfile not allocated!'
allocate(qspecfile(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: qspecfile not allocated!'
c
allocate(isg1(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: isg1 not allocated!'
allocate(isg2(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: isg2 not allocated!'
allocate(nsg(ngrn),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: nsg not allocated!'
c
do ig=1,ngrn
call skipdoc(unit)
read(unit,*)grndep(ig),specfile(ig),grnsel(ig)
if(grnsel(ig).lt.0.or.grnsel(ig).gt.1)then
stop ' bad Green function selection!'
endif
grndep(ig)=grndep(ig)*KM2M
enddo
c
c sort green function files by source depth
c
do i=1,ngrn
do j=i+1,ngrn
if(grndep(j).lt.grndep(i))then
dswap(1)=grndep(i)
fswap=specfile(i)
iswap=grnsel(i)
c
grndep(i)=grndep(j)
specfile(i)=specfile(j)
grnsel(i)=grnsel(j)
c
grndep(j)=dswap(1)
specfile(j)=fswap
grnsel(j)=iswap
endif
enddo
enddo
c
do flen=80,1,-1
if(grndir(flen:flen).ne.' ')goto 200
enddo
200 continue
do ig=1,ngrn
uspecfile(ig)=grndir(1:flen)//'U_'//specfile(ig)
vspecfile(ig)=grndir(1:flen)//'V_'//specfile(ig)
wspecfile(ig)=grndir(1:flen)//'W_'//specfile(ig)
especfile(ig)=grndir(1:flen)//'E_'//specfile(ig)
fspecfile(ig)=grndir(1:flen)//'F_'//specfile(ig)
gspecfile(ig)=grndir(1:flen)//'G_'//specfile(ig)
pspecfile(ig)=grndir(1:flen)//'P_'//specfile(ig)
qspecfile(ig)=grndir(1:flen)//'Q_'//specfile(ig)
enddo
c
c multi-event source parameters
c =============================
c
call skipdoc(unit)
read(unit,*)ns,sdfsel
c
allocate(sfr(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: sfr not allocated!'
allocate(sft(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: sft not allocated!'
allocate(sfp(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: sfp not allocated!'
c
allocate(mrr(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: mrr not allocated!'
allocate(mtt(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: mtt not allocated!'
allocate(mpp(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: mpp not allocated!'
allocate(mrt(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: mrt not allocated!'
allocate(mpr(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: mpr not allocated!'
allocate(mtp(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: mtp not allocated!'
allocate(lats(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lats not allocated!'
allocate(lons(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: lons not allocated!'
allocate(deps(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: deps not allocated!'
allocate(togs(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: togs not allocated!'
allocate(trss(ns),stat=ierr)
if(ierr.ne.0)stop ' Error in qpgettinp: trss not allocated!'
c
if(sdfsel.eq.1)then
do is=1,ns
call skipdoc(unit)
c
c the six moment-tensor elements: Mrr, Mtt, Mpp, Mrt, Mrp, Mtp
c
read(unit,*)munit,mrr(is),mtt(is),mpp(is),
& mrt(is),mpr(is),mtp(is),
& lats(is),lons(is),deps(is),togs(is),trss(is)
mtt(is)=mtt(is)*munit
mpp(is)=mpp(is)*munit
mrr(is)=mrr(is)*munit
mtp(is)=mtp(is)*munit
mpr(is)=mpr(is)*munit
mrt(is)=mrt(is)*munit
deps(is)=deps(is)*KM2M
sfp(is)=0.d0
sft(is)=0.d0
sfr(is)=0.d0
enddo
else if(sdfsel.eq.2)then
do is=1,ns
call skipdoc(unit)
read(unit,*)munit,strike,dip,rake,
& lats(is),lons(is),deps(is),togs(is),trss(is)
call moments(munit,strike,dip,rake,
& mtt(is),mpp(is),mrr(is),
& mtp(is),mpr(is),mrt(is))
deps(is)=deps(is)*KM2M
sfp(is)=0.d0
sft(is)=0.d0
sfr(is)=0.d0
enddo
else if(sdfsel.eq.3)then
do is=1,ns
call skipdoc(unit)
read(unit,*)munit,sfe,sfn,sfz,
& lats(is),lons(is),deps(is),togs(is),trss(is)
sfp(is)=sfe*munit
sft(is)=-sfn*munit
sfr(is)=sfz*munit
deps(is)=deps(is)*KM2M
mtt(is)=0.d0
mpp(is)=0.d0
mrr(is)=0.d0
mtp(is)=0.d0
mpr(is)=0.d0
mrt(is)=0.d0
enddo
else
stop ' bad selection for the source data format!'
endif
c
togsmin=togs(1)
do is=1,ns
togsmin=dmin1(togsmin,togs(is))
enddo
do is=1,ns
togs(is)=togs(is)-togsmin
enddo
c
c sort sub-events by depth
c
do i=1,ns
do j=i+1,ns
if(deps(j).lt.deps(i))then
dswap(1)=lats(i)
dswap(2)=lons(i)
dswap(3)=deps(i)
dswap(4)=mtt(i)
dswap(5)=mpp(i)
dswap(6)=mrr(i)
dswap(7)=mtp(i)
dswap(8)=mpr(i)
dswap(9)=mrt(i)
dswap(10)=sft(i)
dswap(11)=sfp(i)
dswap(12)=sfr(i)
dswap(13)=togs(i)
dswap(14)=trss(i)
c
lats(i)=lats(j)
lons(i)=lons(j)
deps(i)=deps(j)
mtt(i)=mtt(j)
mpp(i)=mpp(j)
mrr(i)=mrr(j)
mtp(i)=mtp(j)
mpr(i)=mpr(j)
mrt(i)=mrt(j)
sft(i)=sft(j)
sfp(i)=sfp(j)
sfr(i)=sfr(j)
togs(i)=togs(j)
trss(i)=trss(j)
c
lats(j)=dswap(1)
lons(j)=dswap(2)
deps(j)=dswap(3)
mtt(j)=dswap(4)
mpp(j)=dswap(5)
mrr(j)=dswap(6)
mtp(j)=dswap(7)
mpr(j)=dswap(8)
mrt(j)=dswap(9)
sft(j)=dswap(10)
sfp(j)=dswap(11)
sfr(j)=dswap(12)
togs(j)=dswap(13)
trss(j)=dswap(14)
endif
enddo
enddo
c
isg1(1)=1
is1=1
do ig=1,ngrn-1
depdif=0.5d0*(grndep(ig)+grndep(ig+1))
isg2(ig)=isg1(ig)-1
do is=is1,ns
if(deps(is).lt.depdif)then
isg2(ig)=is
endif
enddo
isg1(ig+1)=isg2(ig)+1
is1=isg1(ig+1)
enddo
isg2(ngrn)=ns
c
do ig=1,ngrn
nsg(ig)=max0(0,1+isg2(ig)-isg1(ig))
enddo
c
c receiver parameters
c ===================
c
call skipdoc(unit)
read(unit,*)(icmp(i),i=1,11)
call skipdoc(unit)
read(unit,*)outfile
call skipdoc(unit)
read(unit,*)twinout
ntcutout=min0(nt,1+idnint(twinout/dt))
call skipdoc(unit)
read(unit,*)nlpf,f1corner,f2corner
call skipdoc(unit)
read(unit,*)slwlwcut,slwupcut
slwlwcut=slwlwcut/KM2M
slwupcut=slwupcut/KM2M