C23456789012345678901234567890123456789012345678901234567890123456789012
C*******************Joachim Ungerer***********Stuttgart**02.10.1990*****
c
c Mit Korrekturen von Karl Koch, 1991, 1995
c Diese Version enthaelt die Besselfunktions-Subroutinen nach
c "Numerical Recipes in Fortran", nicht optimal da nicht 
c echt doppeltgenau, E.W. Maerz 2004
c Kleinere Anpassungen an den NAS FortranPlus Compiler, E.W. Maerz 04
c minor updates October 2009, E.W.
c
c *******************************************************************
c
c DEFINITIONS (translation of the German text below, by EW)
c The author's definitions are not quite clear. According to a few
c test cases, it appears that the following definitions are used:
c The moment tensor is defined as in Aki and Richard's book (first
c edition, p. 117, Box 4.4):
c strike: cw from North, fault plane dipping to the right
c dip: measured against the horizontal
c rake: direction of displacement of the upper block, measured in the
c fault plane from the horizontal upward (ccw as seen from above)
c
c Azimuth: from epicenter to station, cw from North
c Radial component: outward from the epicenter
c Transverse component: cw from radial
c Vertical: up
c
c *******************************************************************
c
c Bem. von EW, August 2004:
c Ungerers Definitionen des Momenten-Tensors, des Azimuts und der
c Seismogrammkomponenten koennten missverstanden werden.
c
c Anhand einiger Testfaelle halte ich folgende Def. fuer zutreffend:
c => Momenten-Tensor wie in Aki & Richards (1. Auflage) S. 117 Box 4.4
c    Die Formeln von A&R sind im Matlab-Skript "Momentens" programmiert.
c    Dieses Skript erzeugt auch ein Abstrahldiagramm fuer P-Wellen.
c => Stations-Azimut vom Epizentrum zur Station, von Nord ueber Ost (cw)
c => Radialkomponente weg vom Epizentrum
c => Transversalkomponente rechts (90 Grad cw) von der radialen
c => Vertikalkomponente nach oben (also alles wie ueblich)
c
c Definitionen der Herdwinkel in A&R und im Matlab-Skript:
c => strike analog zum Azimut, Orientierung so dass Flaeche nach
c           rechts einfaellt
c => dip    gegen die Horizontale, <= 90 Grad
c => rake   Verschiebungsrichtung des Hangenden (also des Blocks
c           rechts von der Bruchflaeche). Der Winkel wird in der 
c           Bruchflaeche von der Horizontalen ausgehend ccw gemessen
c           (also nach oben positiv)
c           0 Grad: Horizontalverschiebung, 90 Grad: Aufschiebung
c
c ********************************************************************
C                                                                      C
C Programm REFSEIS (Version 02.10.1990)                                C
C Developed for SUN 4/110 workstation, written in SUN-FORTRAN          C
C                                                                      C 
C Part of the thesis of Joachim Ungerer                                C
C Work performed at Institute of Geophysics, Stuttgart University      C
C                                                                      C
C Functions and limitations of REFSEIS.F:                              C
C The program REFSEIS.F generates complete (body and surface wave)     C
C synthetic near- and far-field seismograms for a generalized point    C
C source of Brustle-Muller model or an impulse excitation which is     C
C located within horizontal layered medium. Theoretical basis of the   C
C program is presented in the thesis [Ungerer, 1990].                  C
C Medium can be arbitrarily layered. As a special case, full space can C
C be modeled. Layer parameters include velocities of P and S waves,    C
C density, Q of P and S waves. In the framework of acausal damping     C
C model, strong damping effect can be satisfactorily taken into        C
C account which depends on the station distribution and layer          C
C parameters. Elastic case can be approximated by very high Q.         C
C REFSEIS.F calculates near- and far-field seismograms separately      C
C [Ungerer, 1990]. Either displacement (Typ=1) or velocity (Typ=2) or  C
C acceleration (Typ=3) is calculated. The source excitation function   C
C is implemented in the first part of the program. In addition to the  C
C seismograms, their corresponding spectra can also be obtained in the C
C output.                                                              C
C Seismograms are represented in a cylindrical coordinate system with  C
C origin at the epicenter as three components -- radial(r),            C
C azimuthal(phi) and vertical(z).                                      C
C The number of receivers can be arbitrary in principle. their         C
C locations with respect to the source is described by epicentral      C
C distance and azimuth.                                                C
C There is only one limitation on the location of the source. It can   C
C not be put at the surface, since the analytic expression on which    C
C the program is based would no longer be defined [Ungerer, 1990].     C
c ***** The source must not be vertically under a station *****        c
C The source mechanism and its orientation are described by the        C
C Cartician components of the source moment tensor.                    C
C Difficulties and numerical disturbances arise mainly in the slowness C
C integration. The integral is exact only when the upper bound is      C
C unlimited. Numerically we need to truncate it. Either a small        C
C integration range or a rough sampling can lead to considerable       C
C disturbances or the truncation phase in the resultant seismograms    C
C which are dependent on the source depth, the epicentral distance of  C
C the receiver and the Q factor. These disturbances can never be       C
C completely eliminated. But with proper selection of integration      C
C range, sampling interval, source depth and Q, useful seismograms can C
C be produced by REFSEIS which allows the investication of the near-   C
C and far-field seismograms from a particular event.                   C
C                                                                      C
C                                                                      C
C Brief description of the program:                                    C
C REFSEIS.F calculates the seismogram components of a generalized      C
C point source with strength M0 (seismic moment) at depth ZQ within a   C
C horizontally layered medium according to reflectivity method.        C
C Type and orientation of the source are defined by the components of  C
C its moment tensor. Excitation of the source is Brustle-Muller model  C
C or a delta function.                                                 C
C Seismogram components are calculated separately as:                  C
C      Far-field components   SFr,SFphi,SFz                            C
C      Near-field components  SNr,SNphi,SNz                            C
C They are represented in a cylindrical coordinate system with the     C
C origin at the epicenter.                                             C
C The seismogram can be reduced to receiver dependent first arrival    C
C time by reduction velocity Vred. They can also be shifted recerver-  C
C independently to the left by Tli amount or to the right by Tre       C
C amount.                                                              C
C Different seismogram types can be generated:                         C
C      -Displacement seimogram   (Typ=1)                               C
C      -Velocity seismogram   (Typ=2)                                  C
C      -Acceleration seismogram    (Typ=3)                             C
C                                                                      C
C NOTE!!!! IF THE SOURCE IS AN IMPULSE (Thd = 0), ONLY DISPLACEMENT IS C
C          GENERATED NO MATTER WHAT TYPE IS SET !!!!!!!!!!!!!!!!!!!    C
C                                                                      C
C Input is read from the file:   refseis.par                           C
C Protocol is written to the file:  refseis.out                        C
C Spectra are written to refspec.plt                                   C
C Seismograms are written to refseis.sig                               C
C                                                                      C
C The files refseis.spe and refseis.sig are optionally generated.       C
C                                                                      C
C The file refseis.out is a protocol. The kind of                      C
C information contained is determined through the switches.            C
C  (1) vertical slownesses of layers and reflection and transmission   C
C      coefficients of the interfaces.                                 C
C  (2) reflectivities and transmissivities of the source layer         C
C          Reflectivity matrix      RM(inus) = (RM11,RM12,RM21,RM22)   C
C                    "              RP(lus)  = (RP11,RP12,RP21,RP22)   C
C          Reflectivity scalar      rm(inus) = rm                      C
C                    "              rp(lus)  = rp                      C
C          Transmissivity matrix    TP(lus)  = (TP11,TP12,TP21,TP22)   C
C          Transmissivity scalar    tp(lus)  = tp                      C
C  (3) surface amplitutes         B0(0,1,2),D0(0,1,2),F0(1,2)          C
C  (4) normalized impules spectra (complex) VF,VN in [s*s/(g*cm)]*1E-20C
C  (5) spectrum moduli    BSF,BSN of source strength   M0              C
C  (7) far- and near-field seismograms  in [cm,cm/s,cm/(s*s)]          C
C                                                                      C
C                                                                      C
C subroutines:                                                         C
C       RTKC      calculates reflection and transmission coefficients  C
C       DFORK     fast Fourier transform                               C
C       MATINV    Matrix invertion                                     C
C       MAMULT    Matrix multiplication                                C
C All subroutines are at the end of the program.                       C
C                                                                      C
C Deviations from the FORTRAN standard:                                C
C special functions:                                                   C
C       d_j0,d_j1,d_jn  Bessel functions of 0-th, 1-th and n-th order  C
C                       (when called in SUN FORTRAN, should be         C
C                        declared as DOUBLE PRECISION)                 C
C       dtime(tarray)   determination of calculation time              C
C                                                                      C
C extended command and functions:                                      C
C      Die Berechnungen finden mit doppelter Genauigkeit statt.        C
C      Daher laeuft dieses Programm nur, wenn in FORTRAN folgende      C
C      Befehle implementiert sind:  DOUBLE COMPLEX bzw COMPLEX*16,     C
C                                   CDABS,CDEXP,CDSQRT,DCMPLX,DCONJG,  C
C                                   DREAL                              C
C      fuer NAS-Fortran statt cdabs usw.: abs, exp, sqrt               C
C                                                                      C
C Note1: As the Fourier transform of the displacement has a pole at    C
C        frequency zero, we define fr(1)=1.D-30 for displacement       C
C        calculation.                                                  C
C                                                                      C
C Note2: The occurrence of the truncation phase may be reduced by      C
C        the extention of the slowness integration range or the        C
C        increase of the sampling points Nu.                           C
C                                                                      C
C                                                                      C
C Description of variables:                                            C
C (The names of the variables are chosen mostly in accordance to the   C
C thesis)                                                              C
C                                                                      C
C    Eingabevariablen:                                                 C
C          Die Eingabe des Modells sollte nach der Vorlage des Ein-    C
C          gabefiles refseis.par erfolgen.                             C
C          Bis auf die Schaltervariable werden alle Eingabevariablen   C
C          formatiert eingelesen, sie sind deshalb rechtsbuendig in    C
C          ihre jeweiligen Eingabefelder zu setzen.                    C
C          text0         Kopfzeile fuer Text                           C
C      I.  Schalter (0=nein, 1=ja, 2=ja und Ende)                      C
C          Sch(1)        vert.Langsamkeiten,Ref.u.Transmissionskoeff.  C
C          Sch(2)        Reflektivitaeten und Transmissivitaeten       C
C          Sch(3)        Oberflaechenamplituden                        C
C          Sch(4)        normierte Impuls-Spektren(komplex)            C
C          Sch(5)        endgueltige Spektren(Betraege)                C
C          Sch(6)        Plot-File fuer die endgueltigen Spektren      C
C          Sch(7)        Endgueltige Nah- und Fernfeldseismogramme     C
C          Sch(8)        Plot-File fuer Seismogramme                   C
C          Sch(9)        noch frei                                     C
C          Sch(10)       noch frei                                     C
C      II. Schichtdaten                                                C
C          z(i)          Tiefe der Schichtdecke der i-ten Schicht [km] C
C                        (die "nullte" Schicht repraesentiert den obe- C
C                        ren Halbraum,z.B. Atmosphaere,ihre Tiefe ist  C
C                        definitionsgemaess Null)                      C
C          alpha(i)      P-Wellengeschw.[km/s] in Schicht i            C
C          beta(i)       S-   "                "                       C
C          rho(i)        Dichte[g/ccm] der Schicht i                   C
C          Qa(i)         Q-Faktor fuer P-Wellen                        C
C          Qb(i)            "       "  S-Wellen                        C
C          Bem: Grosse Schichtdicken koennen zum Ueberlauf durch Ex-   C
C          ponentialterme fuehren. Abhilfe durch Einfuegen virtueller  C
C          Trennflaechen.                                              C
C      III.Herddaten                                                   C
C          ZQ            Tiefe[km] der ErdbebenQuelle                  C
C          M0            Staerke des Momententensors in [dyn*cm]       C
C          Mxx,Myy,Mzz   normierter Momententensor  M=(ni*fj+nj*fi)    C
C          Mxy,Mxz,Myz   (Bem:seismischer Momententensor MS=M0*M)      C
C          Bem1: Ein grosser Abstand des Herdes zur Schichtdecke kann  C
C          zum Ueberlauf durch die Exponentialterme ea,eb fuehren. Ab- C
C          hilfe durch Einfuegen einer virtuellen Trennflaeche nahe    C
C          ueber dem Herd.                                             C
C          Bem2: Die Umrechnung von Streichen(S),Fallen(F) und Dislo-  C
C          kationswinkel(D) beim Double-Couple in die Vektoren f und n C
C          erfolgt durch: f=(cosD*cosS+cosF*sinD*sinS,cosD*sinS-cosF*  C
C          sinD*cosS,-sinD*sinF); n=(sinF*sinS,-sinF*cosS,cosF)        C
C          [siehe hierzu auch Aki/Richards S.106 ff]                   C
C      IV. Daten fuer Normierte Anregungsfunktion                      C
C          typ      1=Verschiebung,2=Geschwindigkeit,3=Beschleunigung  C
C          The      Einsatzzeit des Herdes                             C
C          Thd      Dauer des Zeitverlaufs von The bis Endwert Null    C
C                   (Momentenanstiegszeit)                             C
C          Bem: Impulsanregung fuer The=Thd=0                          C
C      V.  Reduktionsdaten                                             C
C          Vred     Reduktionsgeschwindigkeit[km/s] (bei Vred=0 keine  C
C                   empfaengerabhaengige Reduktion),die zugehoerige    C
C                   Reduktionszeit ist hypozentralbezogen.             C
C          Tli      Zeitverschiebung[s] des Seismogramms nach rechts   C
C          Tre      Zeitverschiebung[s] des Seismogramms nach links    C
C      VI. Empfaengerdaten                                             C
C          r(E)      Radialentf.[km] des Empfaengers E vom Epizentrum  C
C          phi(E)    Winkel[Grad] des Empfaengers E zum Epizentrum     C
C     VII. Langsamkeitsdaten                                           C
C          umin,umax  Langsamkeitsbereich [s/km]                       C
C          uwil,uwir  Langsamkeits-Window links,rechts [s/km]          C
C                     (Ueblich sind umin=uwil=0, uwir sollte reziprok  C
C                     der kleinsten Wellengeschwindigkeit der Herd-    C
C                     schicht oder groesser sein)                      C
C          Nu         Zahl der gleichabstaendigen Langsamkeit.(NU muss C
C                     mindestens 2 sein)                               C
C          Bem:Langsamkeitsbereich und Integrationsgenauigkeit Nu sind C
C          die masgebenden Groessen um moeglichst stoerungsfreie Seis- C
C          mogramme zu erzeugen. Das rechteckigen Langsamkeitswindow   C
C          [uwil,uwir] wird beidseitig mit einem Cos-Taper bis zu den  C
C          jeweiligen Grenzen [umin,umax] auf Null heruntergezogen ->  C
C          evt. Trennung von Raum und Oberflaechenwellen durch Wahl    C
C          geeigneter Langsamkeitsbereiche.                            C
C     VIII.Frequenz- und Zeitdaten                                     C
C          fmin,fmax     Frequenzbereich [Hz](physikalisch)            C
C          fwil,fwir     Frequenz-Window links,rechts [Hz]             C 
C          Dt            zeitliches Abtastintervall [sec]              C
C          SL            Seismogrammlaenge (Zweierpotenz)=Anzahl der   C
C                        Stuetzstellen des Seismogrammes(maximal=MSL)  C
C          Bem: das rechteckige Frequenzwindow [fwil,fwir] wird beid-  C
C          seitig mit einem Cos-Taper bis zu den Grenzen [fmin,fmax]   C
C          auf Null heruntergezogen -> Frequenzfilter                  C
C                                                                      C
C    Laufvariablen                                                     C
C          E       fuer Empfaenger                                     C
C          f       fuer Frequenzen, teilweise auch fuer Zeiten         C
C          i       fuer Schichten, und als Zaehlindex                  C
C          k       Zaehlindex                                          C
C          l       fuer Langsamkeiten                                  C
C                                                                      C
C    Parameter(koennen nach Belieben im Programm eingestellt werden)   C
C          IME     Imaginaere Einheit                                  C
C          ME      Max.Zahl der Empfaenger                             C
C          Mf          "        physikalischen Frequenzen(=MSL/2+1)    C
C                               bis zur Nyquistfrequenz Fny            C
C          MS          "        Schichten                              C
C          MSL         "        Stuetzstellen (Max.Seismogrammlaenge)  C
C          hin,rueck  -1.,+1. zur Steuerung der Fouriertransformation  C
C          pi      =3.14159265358979                                   C
C                                                                      C
C    Charaktervariablen                                                C
C          typ0    ='Verschiebung','Geschwindigkeit' oder              C
C                   'Beschleunigung', je nach Schalter(9)              C
C          impuls0 ='Impulsanregung'                                   C
C                                                                      C
C    Variablen fuer Rechenzeit                                         C
C          dtime    Funktion zur Bestimmung der Rechenzeit             C
C          tarray   ihre zweidim. Variable(User-Zeit,System-Zeit)      C
C          time     Aufnahmevariable fuer Rechenzeit                   C
C                                                                      C
C    Hilfsvariablen zur Zwischenspeicherung                            C
C          help[reell],hilf,gew                                        C
C          H11,H12,H21,H22                                             C
C          I11,I12,I21,I22                                             C
C                                                                      C
C    Rechenvariablen(fest)                                             C
C          ap   =Wurzel(SL)*Dt Aperiodizitaetsfaktor fuer Fourierkoeff.C
C          alphC(i)  komplexe P-Wellengeschwindigkeit                  C
C          betaC(i)     "     S-Wellen     "                           C
C          d(i)      Dicke der Schicht i in [km]                       C
C          Df        Schrittweite der Frequenz,Frequenzintervall       C
C          Du        Schrittweite der Langsamkeit u                    C
C          DDu       getaperte und gewichtete Integrationsschrittweite C
C          FL        Frequenzlaenge des vollstaendigen Spektrums       C
C          fmi,fma   Index der Frequenzen fmin,fmax                    C
C          Fny       Nyquistfrequenz (fmax muss < Fny sein)            C
C          fr(f)     f-te Frequenz in [Hz]                             C
C          fwl,fwr   Index der Frequenzen fwil,fwir                    C
C          h         Schicht in der der Herd sich befindet             C
C          M0dim     =M0/1.D20 zur Dimensionierung der Seismogramme    C
C          n         Zahl der Schichten                                C
C          NE        Zahl der Empfaenger                               C
C          Nf        Zahl der phys. Frequenzen im Bereich [fmax,fmin]  C
C          Nff        "    "    "       "      im Fenster [fwil,fwir]  C
C          Ny        Index der Nyquistfrequenz Fny (max.phys.Freqenz)  C
C          phiB(E)   Winkel des Empfaengers E in Bogenmass             C
C          t(f)      f-te Zeit in [s], d.h. Zeitwerte des Seismogramms C
C          TL        Zeitlaenge des Seismogramms (Periode)             C
C          tred(f)   reduzierte Zeitskala  / modified 8-2-91 K.Koch /  C
C                                                                      C
C          Variablen der normierten Herdfunktion                       C
C          hfkt(tvar,T1,T2) Definiton der Herdfunktion                 C
C          tvar           ihre Zeitvariable                            C
C          T1,T2          ihre Begrenzungsparameter(T1=The,T2=The+Thd) C
C          g(f)           Feld in dem hfkt je Frequenz abgelegt wird   C
C          Die Dimension der Herdfunktion richtet sich nach ihrer In-  C
C          terpretation:als Verschiebung [ ]  ->  Bodenverschiebung    C
C                       als 1.Ableitung [1/s] ->  Bodengeschwindigkeit C
C                       als 2.Ableitung [1/(s*s)]->  "  beschleunigung C
C                                                                      C
C    Rechenvariablen(veraenderlich)                                    C
C          a(i)     vert.Langsamkeit der P-Wellen in Schicht i [s/km]  C
C          b(i)      "        "      der S-Wellen       "              C
C          E11(i),E22(i)     Phasenmatrix der Schicht i (Diagonalel.)  C
C          ea         =exp[IME*w*a(h)*(ZQ-z(h))]                       C
C          eb         =   "      b(h)    "                             C
C          earg       Argument fuer Exponentialfunktionen              C
C          ftap(f)    Fensterfunktion der Frequenz mit Cos-Taper zur   C
C                     Verwendung als Frequenzfilter,das Rechteckwindow C
C                     [fwil,fwir] wird beidseitig mit einem Cos-Taper  C
C                     bis zu den Frequenzen fmax bzw. fmin auf Null    C
C                     heruntergezogen.                                 C
C          red        Reduktionsfaktor:red=exp[i*w*(R(E)/Vred+Tli-Tre] C
C                     mit R(E)=Wurzel[ZQ*ZQ+r(E)*r(E)]=Hypozentralentf.C
C          u          horizontale Langsamkeit in [s/km]                C
C          uQ         ihr Quadrat                                      C
C          utap       Fensterfunktion der Langsamkeit mit Cos-Taper    C
C                     zum Glaetten der Stoerphasen                     C
C          w(f)       Kreisfrequenz Omega [1/s]                        C
C                                                                      C
C          Reflektions und Transmissionskoeffizienten d.Trennflaeche i C
C          Rppd(i),Rpsd(i),Tppd(i),Tpsd(i)  fuer  P-Wellen von oben    C
C          Rssd(i),Rspd(i),Tssd(i),Tspd(i)  fuer SV-Wellen      "      C
C          rd(i),td(i)                      fuer SH-Wellen      "      C
C          Rppu(i),Rpsu(i),Tppu(i),Tpsu(i)  fuer  P-Wellen von unten   C
C          Rssu(i),Rspu(i),Tssu(i),Tspu(i)  fuer SV-Wellen      "      C
C          ru(i),tu(i)                      fuer SH-Wellen      "      C
C                                                                      C
C          Reflektivitaets- und Transmissivitaetsmatrizen              C
C          RTD11(i),RTD12(i) Reflektivitaetsmatrix f.Einfall von oben, C
C        & RTD21(i),RTD22(i) fuer unteren Stapel ab Decke Schicht i    C 
C          RBU11(i),RBU12(i) Reflektivitaetsmatrix f.Einfall v. unten, C
C        & RBU21(i),RBU21(i) fuer oberen Stapel ab Boden Schicht i     C
C          RTU11(i),RTU12(i) Reflektivitaetsmatrix f.Einfall v. unten, C
C        & RTU21(i),RTU22(i) fuer oberen Stapel ab Decke Schicht i     C
C          TTUo11(i),TTUo12(i) Transmissivitaetsmatrix f.Einf.v.unten, C
C        & TTUo21(i),TTUo22(i) ab Decke Schicht i bis Decke Schicht 0  C
C          S11(i),S12(i)     Schichtmatrix von Schicht i,zur Berech-   C
C        & S21(i),S22(i)     nung der Transmissivitaetsmatrix TTUo     C
C                                                                      C
C          Reflektivitaets- und Transmissivitaetsskalare               C
C          rtd(i)   SH-Reflektivitaet ab Decke Schicht i abwaerts      C
C          rbu(i)   SH-Reflektivitaet ab Boden Schicht i aufwaerts     C
C          rtu(i)   SH-Reflektivitaet ab Decke Schicht i    "          C
C          ttuo(i)  SH-Transmissivitaet ab Decke Schicht i  "          C
C          s(i)     Schichtskalar von Schicht i,fuer Skalar TTU        C  
C                                                                      C
C          Reflektivitaeten und Transmissivitaeten der Herdschicht     C
C          RM11,RM12,RM21,RM22   Reflektivitaetsmatrix RM(inus)        C
C          RP11,RP12,RP21,RP22   Reflektivitaetsmatrix RP(lus)         C
C          TP11,TP12,TP21,TP22   Transmissivitaetsmatrix TP(lus)       C
C          rm                    SH-Reflektivitaet rm(inus)            C
C          rp                    SH-Reflektivitaet rp(lus)             C
C          tp                    SH-Transmissivitaet tp(lus)           C
C                                                                      C
C          Amplituden der Wellenfelder                                 C
C          AQ(0,1,2)  Quell-Amplituden des abwaertslaufenden           C
C          CQ(0,1,2)  Wellenfeldes der Quelle                          C
C          EQ(1,2)         "                                           C
C          BQ(0,1,2)  Quell-Amplituden des aufwaertslaufenden          C
C          DQ(0,1,2)  Wellenfeldes der Quelle                          C
C          FQ(1,1)         "                                           C
C          B0(0,1,2)  Oberflaechenamplituden                           C
C          D0(0,1,2)           "                                       C
C          F0(1,2)             "                                       C
C                                                                      C
C          Empfaengerabhaengige Variablen                              C
C          K0(E),K1(E)  Faktoren des Winkelanteils,enthalten auch      C
C          K2(E),K3(E)  die Komponenten des Momententensors            C
C          L1(E),L2(E)                                                 C
C          jarg         Argument fuer Besselfunktionen                 C
C          J0,J1,J2     zur Aufnahme der Werte der Besselfunktionen    C
C                       (d_j0,d_j1,d_jn, wobei n=2)                    C
C                                                                      C
C          Integrationsvariablen                                       C
C          In1...In6  Hilfsvariablen zur Berechnung der Integranden    C
C          FFI        =w*w/(4pi*rho(h))  Fernfeld-Integrandenfaktor    C
C          NFI        =w/(r*4pi*rho(h))  Nahfeld-    "                 C
C          Integrandenglieder je Langsamkeit fuer f-te Frequenz        C
C          des E-ten Empfaengers:  /modified 8-2-91 k.koch /           C
C          IDFr       Int.glied der radialen    Fernfeldkomponente     C
C          IDNr          "             "        Nahfeldkomponente      C
C          IDFphi        "        transversalen Fernfeldkomponente     C
C          IDNphi        "             "        Nahfeldkomponente      C
C          IDFz          "        vertikalen    Fernfeldkomponente     C
C          IDNz          "             "        Nahfeldkomponente      C
C                                                                      C
C          normierte Teil-Impulsspektren fuer f-te Frequenz des E-ten  C
C          Empfaengers, fuer Verschiebung in [s*s/(g*cm)]* 1E-20:      C
C          VFr(f,E)   Radiales      Fernfeld-Impulsspektrum            C
C          VNr(f,E)      "           Nahfeld-      "                   C
C          VFphi(f,E) Transversales Fernfeld-      "                   C
C          VNphi(f,E)    "           Nahfeld-      "                   C
C          VFz(f,E)   Vertikales    Fernfeld-      "                   C
C          VNz(f,E)      "           Nahfeld-      "                   C
C                                                                      C
C          End-Spektren(Betraege) fuer Staerke M0:                     C
C          BSFr(f),BSFphi(f),BSFz(f)       / modified 8-2-91 K.Koch /  C
C          BSNr(f),BSNphi(f),BSNz(f)       / modified 8-2-91 K.Koch /  C
C          in [cm*s],[cm],[cm/s]                                       C
C                                                                      C
C                                                                      C
C    Ergebnissvariablen                                                C
C          Vollstaendige Seismogramme  bzw. Spektren                   C
C          SFr(f),SFphi(f),SFz(f)   Fernfeld  /modified 8-2-91 K.Koch/ C
C          SNr(f),SNphi(f),SNz(f)   Nahfeld   /modified 8-2-91 K.Koch/ C
C          fuer Bodenverschiebung in     [cm]         bzw. [cm*s]      C
C          fuer Bodengeschwindigkeit in  [cm/s]       bzw. [cm]        C
C          fuer Bodenbeschleunigung in   [cm/(s*s)]   bzw. [cm/s]      C
C**********************************************************************C



      PARAMETER  (ME=50)
      PARAMETER  (MSL=2**12) 
      PARAMETER  (Mf=MSL/2+1)
      PARAMETER  (MS=50)
      REAL*8      hin,rueck,pi
      PARAMETER  (hin=-1.D0,rueck=1.D0,pi=3.14159265358979)
      COMPLEX*16  IME
      PARAMETER  (IME=(0.D0,1.D0))


      CHARACTER impuls0*15,text0*70,typ0*15,nix*5,zeile*70

      INTEGER   E,f,fmi,fma,fwl,fwr,h,i,k,l,n,NE,Nf,Nff,
     &          Nu,Ny,SL,typ,Sch(10),estim

      REAL*4    dtime,tarray(2),time

C Dimension of w() set to MSL to prevent overwriting of other
C arrays                          / modified  8-20-91 K.Koch /
C
      REAL*8    alpha(0:MS),beta(0:MS),d(0:MS),rho(0:MS),z(0:MS),
     &          Du,u,umin,umax,uwil,uwir,utap,uQ,gew,ZQ,
     &          Df,Dt,fmin,fmax,fwil,fwir,FL,Fny,TL,
     &          fr(MSL),ftap(Mf),t(MSL),w(MSL),
     &          phi(ME),phiB(ME),r(ME),FFI,NFI,
     &          ap,hfkt,Thd,The,T1,T2,tvar,
     &          tred(MSL),Tli,Tre,Vred,help,
     &          M0,M0dim,Mxx,Myy,Mzz,Mxy,Mxz,Myz,
     &          K0(ME),K1(ME),K2(ME),K3(ME),L1(ME),L2(ME),
     &          BSFr(MSL),BSFphi(MSL),BSFz(MSL),
     &          BSNr(MSL),BSNphi(MSL),BSNz(MSL),Qa(0:MS),Qb(0:MS)

C  Change the Q parameters from integer to real number.  D. Yang 8-8-95

      COMPLEX*16  a(0:MS),b(0:MS),alphC(0:MS),betaC(0:MS),
     &            Rppd(MS),Rpsd(MS),Rssd(MS),Rspd(MS),rd(MS),
     &            Tppd(MS),Tpsd(MS),Tssd(MS),Tspd(MS),td(MS),
     &            Rppu(MS),Rpsu(MS),Rssu(MS),Rspu(MS),ru(MS),
     &            Tppu(MS),Tpsu(MS),Tssu(MS),Tspu(MS),tu(MS),
     &            RTD11(0:MS),RTD12(0:MS),RTD21(0:MS),RTD22(0:MS),
     &            RBU11(0:MS),RBU12(0:MS),RBU21(0:MS),RBU22(0:MS),
     &            RTU11(0:MS),RTU12(0:MS),RTU21(0:MS),RTU22(0:MS),
     &            TTUo11(0:MS),TTUo12(0:MS),TTUo21(0:MS),TTUo22(0:MS),
     &            rtd(0:MS),rbu(0:MS),rtu(0:MS),ttuo(0:MS),
     &            S11(0:MS),S12(0:MS),S21(0:MS),S22(0:MS),s(0:MS),
     &            RM11,RM12,RM21,RM22,rm,RP11,RP12,RP21,RP22,rp,
     &            TP11,TP12,TP21,TP22,tp

      COMPLEX*16  H11,H12,H21,H22,I11,I12,I21,I22,hilf,
     &            DDu,ea,eb,earg,red,g(MSL),
     &            In1,In2,In3,In4,In5,In6,
     &            AQ(0:2),BQ(0:2),CQ(0:2),DQ(0:2),EQ(2),FQ(2),
     &            B0(0:2),D0(0:2),F0(2),E11(0:MS),E22(0:MS)

      COMPLEX*16  IDFr,IDFphi,IDFz,IDNr,IDNphi,IDNz

      COMPLEX*16  VFr(Mf,ME),VFphi(Mf,ME),VFz(Mf,ME),
     &            VNr(Mf,ME),VNphi(Mf,ME),VNz(Mf,ME)

      COMPLEX*16  SFr(MSL),SFphi(MSL),SFz(MSL),
     &            SNr(MSL),SNphi(MSL),SNz(MSL)

      DOUBLE PRECISION  jarg,d_j0,d_j1,d_jn,J0,J1,J2



C**********************************************************************C
C 1. Herdfunktion (als Funktionsdefinition)                            C
C                                                                      C
C Die zugrunde liegende Momentenfunktion ist entnommen aus:            C
C Bruestle,Mueller:Physics of the earth 32,83 und lautet:              C
C M(t)=Mo*9/16* {[1-cos(pi*t/T)] + 1/9*[cos(3*pi*t/T)-1]}              C
C Wir gehen aus von ihrer 1.Ableitung:3*pi/4T*{[sin(pi*t/T)]**3} ,um   C
C daraus die Bodengeschwindigkeit direkt, die Bodenverschiebung und    C
C Bodenbeschleunigung indirekt ueber die Fouriertransformation zu bere-C
C chnen(siehe Programmteil Nr.9). Die Polstelle der Verschiebungstrans-C
C formierten an der Stelle fr(1)=0 wird abgefangen durch fr(1)=1.D-30).C
C Bei Einbau neuer Momentenfunktionen ist lediglich zu beachten, dass  C
C sie als Geschwindigkeitsfunktionen behandelt  werden. Ihre Anfangs-  C
C und Endwerte muessen Null sein.                                      C
C**********************************************************************C
C     Geschwindigkeitsanregung(1.Ableitung d. Bruestle-Mueller-Funktion)
      hfkt(tvar,T1,T2)=0.75*pi/(T2-T1)*dsin(pi*(tvar-T1)/(T2-T1))**3
 


C**********************************************************************C
C 2. Dateieroeffnung und Start des Zeitaufrufes                        C
C**********************************************************************C
      OPEN(UNIT=1,FILE='refseis.par')
      OPEN(UNIT=2,FILE='refseis.out')
C     Die optionalen Dateien refseis.spe und refseis.sig werden bei
C     Schalter 6 und 8 eroeffnet

C     Start des Zeitaufrufes
c      time=dtime(tarray)
      call CPU_TIME(tarray(1))
c
C**********************************************************************C
C 3. Eingabe der Daten (Freiformat)            C 
C    und Bestimmung der Anzahl von Schichten und Empfangsstationen     C
C**********************************************************************C

C     Kopfzeile einlesen
      READ(1,100) text0 
  100 FORMAT(A)
      write(*,*)
      write(*,*) "REFSEIS ", text0
      write(*,*) "Maximum seismogram length MSL is ",MSL," samples"
      write(*,*) "This version ignores all frequeny parameters in the in
     &put file"

C     Ueberlesen von 3 Zeilen
      READ(1,101)
  101 FORMAT(//)

C     Schalter einlesen(formatfrei)
      READ(1,*)Sch(1),Sch(2),Sch(3),Sch(4),Sch(5),Sch(6),Sch(7),
     &         Sch(8),Sch(9),Sch(10)

C     Ueberlesen von 3 Zeilen
      READ(1,101)

C     Schichtdaten einlesen
      i=-1
  120 i=i+1
      read(1,'(a70)') zeile
      if (index(zeile,'.').gt.0) then
        READ(zeile,*) nix,z(i),alpha(i),beta(i),rho(i),Qa(i),Qb(i)
C 130   FORMAT(T5,F10.4,3F11.4,2F11.1)
        goto 120
      endif

C     Anzahl der Schichten
      n=i-1

C     Uberlesen von 2 Zeilen
      READ(1,131)
  131 FORMAT(/)

C     Herddaten einlesen
      READ(1,*) ZQ,M0,Mxx,Myy,Mzz,Mxy,Mxz,Myz
C 140 FORMAT(F10.4,D11.3,6F8.3)

C     Ueberlesen von 3 Zeilen
      READ(1,101)

C     Daten der Anregungsfunktion einlesen
      READ(1,*) typ,The,Thd
C 145 FORMAT(I4,2F11.4)     

C     Ueberlesen von 3 Zeilen
      READ(1,101)

C     Reduktionsdaten einlesen
      READ(1,*) Vred,Tli,Tre
C 147 FORMAT(3F11.4)     

C     Ueberlesen von 3 Zeilen
      READ(1,101)

C     Empfaengerdaten einlesen
      E=0
  150 E=E+1
      read(1,'(a70)') zeile
      if(index(zeile,'.').gt.0) then
      READ(zeile,*) r(E),phi(E)
      if(phi(E).lt.0.) phi(E)=phi(E)+360.
C 160 FORMAT(F10.4,F11.4)
      GOTO 150
      endif

C     Anzahl der Stationen 
      NE=E-1

C     Ueberlesen von 2 Zeilen
      READ(1,131)

C     Langsamkeitsdaten einlesen
      READ (1,*) umin,uwil,uwir,umax,Nu
C 170 FORMAT(4F11.6,I7)

C     Ueberlesen von 3 Zeilen
      READ(1,101)

C     Frequenz- und Zeitdaten einlesen
      READ(1,*) fmin,fwil,fwir,fmax,Dt,SL
C 180 FORMAT(4F11.4,F10.4,I15)
C *** frequency parameters will be ignored and replaced by 0 .. 2.5/Thd
C     The cutoff is at the first zero of the source spectrum
      Thd=max(Thd,5*dt)
      write(*,182) Thd   
  182 format(' Bandwidth is determined by rise time (',f6.2,
     &  ' sec) or Nyqvist frequency')
      fmin=0.                  
      fwil=fmin
      fwir=2.5/Thd
      fmax=fwir

C**********************************************************************C
C 4. Initialisierungen von Feldern und Variablen                       C
C**********************************************************************C
C     Initialisierung des Frequenzfilters ftap(f)
      DO 185 f=1,SL/2+1
         ftap(f)=0.
  185 CONTINUE

C     Initialisierung der normierten Teilspektren bis zur Nyquistfrequenz
      DO 190 E=1,NE
         DO 184 f=1,SL/2+1
            VFr(f,E)  =(0.,0.)
            VFphi(f,E)=(0.,0.)
            VFz(f,E)  =(0.,0.)
            VNr(f,E)  =(0.,0.)
            VNphi(f,E)=(0.,0.)
            VNz(f,E)  =(0.,0.)
  184    CONTINUE
  190 CONTINUE





C**********************************************************************C
C 5. Vorberechnungen fuer Langsamkeit                                  C
C     -Schrittweite der Langsamkeit Du                                 C
C     -Korrektur von uwil und uwir auf die Schrittweitenskala          C
C**********************************************************************C
C     Schrittweite fuer u
      Du=(umax-umin)/DBLE(Nu-1)
C     Korrektur von uwil und uwir auf die Schrittweitenskala
      uwil=DBLE( IDINT((uwil-umin)/Du) )*Du+umin
      uwir=DBLE( IDINT((uwir-umin)/Du) )*Du+umin




C**********************************************************************C
C 6. Vorberechnungen fuer Frequenz und Zeit                            C
C     -Zeitlaenge TL des Seismogramms und Schrittweite Df der Frequenz C
C     -Berechnung der Frequenzen fr(f) und Zeiten t(f)                 C
C      und Umrechnung in Kreisfrequenz w                               C
C     -Korrektur von fmin,fwil,fwir,fmax auf die Schrittweitenskala    C
C      der Frequenz ,wegen Anpassung an die Fouriertransformation      C
C     -weitere Frequenz und Zeitdaten: Nf,Nff,FL,Ny,Fny                C
C     -Fensterfunktion ftap(f) der Frequenz mit Cos-Taper (verwendbar  C
C      als Frequenzfilter).                                            C
C**********************************************************************C

C     Zeitlaenge des Seismogramms und Schrittweite der Frequenz
      TL=SL*Dt
      Df=1.D0/TL

C     Berechnung der Frequenzen und Zeiten
      DO 200 f=1,SL
         fr(f)=DBLE(f-1)*Df
          t(f)=DBLE(f-1)*Dt
C        Umrechnung in Kreisfrequenz w
         w(f)=2.D0*pi*fr(f)
  200 CONTINUE  

C     Korrektur von fmin,fwil,fwir,famx auf die Schrittweitenskala der
C     Frequenz ,wegen Anpassung an die Fouriertransformation
      fmi=IDINT(fmin/Df)+1
      fwl=IDINT(fwil/Df)+1
      fwr=IDINT(fwir/Df)+1
      fma=IDINT(fmax/Df)+1
      fmin=fr(fmi)
      fwil=fr(fwl)
      fwir=fr(fwr)
      fmax=fr(fma)

C     weitere Frequenz- und Zeitdaten:Nf,Nff,FL,Ny,Fny
      Nf=fma-fmi+1   
      Nff=fwr-fwl+1
      FL=1.D0/Dt
      Ny=SL/2+1
      Fny=fr(Ny)

C     Fensterfunktion ftap(f) der Frequenz mit Cos-Taper,
C     zur Verwendung als Frequenzfilter: beidseitig des rechteckigen
C     Frequenzwindows [fwil,fwir] wird das Frequenzfenster mit einem
C     Cos-Taper jeweils bis zu den Grenzfrequenzen fmin und fmax auf 
C     Null heruntergezogen.
      DO 210 f=fmi,fma
         IF(fmin.LE.fr(f) .and. fr(f).LT.fwil) THEN
            ftap(f)=0.5D0*(1.D0-dcos(pi*(fr(f)-fmin)/(fwil-fmin)))
         ELSE IF(fwil.LE.fr(f) .and. fr(f).LE.fwir) THEN
            ftap(f)=1.D0 
         ELSE IF(fwir.LT.fr(f) .and. fr(f).LE.fmax) THEN
            ftap(f)=0.5D0*(1.D0+dcos(pi*(fr(f)-fwir)/(fmax-fwir)))
         ENDIF
  210 CONTINUE




C**********************************************************************C
C 7. Vorberechnungen fuer Empfaenger                                   C
C     -Umrechnung der Winkel phi(E) in Bogenmass phiB(E)               C
C     -Winkel- und Tensorparameter K0(E),K1(E),K2(E),K3(E),L1(E),L2(E) C
C**********************************************************************C

      DO 220 E=1,NE
C      Umrechnung der Winkel in Bogenmass
       phiB(E)=2.D0*pi*phi(E)/360.D0

C      Winkel- und Tensorparameter
       K0(E)=Mzz
       K1(E)=Mxz*dcos(phiB(E))+Myz*dsin(phiB(E))
       K2(E)=0.5D0*(Myy-Mxx)*dcos(2.D0*phiB(E))-Mxy*dsin(2.D0*phiB(E))
       help=Mxx*dcos(phiB(E))**2 + Myy*dsin(phiB(E))**2
       K3(E)=help + Mxy*dsin(2.D0*phiB(E))
       L1(E)=Mxz*dsin(phiB(E))-Myz*dcos(phiB(E))
       L2(E)=0.5D0*(Myy-Mxx)*dsin(2.D0*phiB(E))+Mxy*dcos(2.D0*phiB(E))
  220 CONTINUE





C**********************************************************************C
C 8. Vorberechnungen fuer Schichtmodell                                C
C     -Dimensionierung der Seismogramme und Spektren mit der Dimensio- C
C      nierungsvariablen M0dim (M0 -> M0dim) auf:                      C
C      Spektren:Verschiebungsdichte[cm*s],Geschwindigkeitsdichte[cm],  C
C               Beschleunigungsdichte[cm/s]                            C
C      Seismogramme:Verschiebung[cm],Geschwindigkeit[cm/s],Beschleu-   C
C               nigung[cm/(s*s)]                                       C
C     -Schichtdicken d(i)                                              C
C     -Herdschicht h                                                   C
C     -komplexe Wellengeschwindigkeiten                                C
C**********************************************************************C

C     Dimensionierung
      M0dim=M0/1.D20

C     Berechnung der Schichtdicken der Schichten 1 bis n-1
      DO 230 i=0,n-1
         d(i)=z(i+1)-z(i)
  230 CONTINUE
C     Die Dicke des unteren Halbraumes (Schicht n) ist keine Rechen-
C     groesse ,fuer den Kontrollausdruck wird sie auf Null gesetzt
      d(n)=0. 

C     Bestimmung der Herdschicht h (liegt der Herd auf der Trenn-
C     flaeche z(i), gehoert er der Schicht i an) 
      IF(ZQ.LT.z(n)) THEN
        i=0
  240   i=i+1
        IF(z(i)-ZQ.LE.0.) GOTO 240
        h=i-1
      ELSE IF(ZQ.GE.z(n))THEN
        h=n
      ENDIF

C     Komplexe Wellengeschwindigkeiten(hier gehen die Q-Faktoren ein)
      DO 250 i=0,n
         alphC(i)=alpha(i)*(1.D0+ 0.5D0*IME/Qa(i))
         betaC(i)= beta(i)*(1.D0+ 0.5D0*IME/Qb(i))
  250 CONTINUE





C**********************************************************************C
C 9. Normierte Anregungsfunktion g(f) des Herdes                       C
C  - Aperidiozitaetsfaktor ap                                          C   
C  - Impuls-Seismogramm fuer The=Thd=0 (frequenzbegrenzt durch ftap)   C
C  - Abspeichern der normierten Herdfunktion im Zeitbereich auf Feld g C
C  - Fouriertransformation in den Frequenzbereich(Aufruf von DFORK)    C
C  - Typ der Anregungsfunktion(typ=1 Verschiebung,typ=2 GeschwindigkeitC
C    typ=3 Beschleunigung)                                             C
C  - Frequenzfilterung mit ftap(f)                                     C
C**********************************************************************C
C     Aperidiozitaetsfaktor(nur relevant fuer die Spektren) 
c      ap=1.D0
      ap=DSQRT(DBLE(SL))*Dt

C     Impuls-Anregung fuer The=Thd=0 (frequenzbegrenzt)
      IF(The.EQ.0 .and. Thd.EQ.0) THEN
        DO 260 f=1,Ny
           g(f)=DCMPLX(1.)*ftap(f)
  260   CONTINUE

C     sonst Berechnung der Herdfunktion und Fouriertransformation
      ELSE
C       Herdfunktion im Zeitbereich
        DO 270 f=1,SL
           IF(t(f).LT.The) THEN
             g(f)=DCMPLX(0.)
           ELSE IF(t(f).GE.The .and. t(f).LE.The+Thd) THEN
             g(f)=DCMPLX(hfkt(t(f),The,The+Thd))
           ELSE IF(t(f).GT.The+Thd) THEN
             g(f)=DCMPLX(0.)
           ENDIF
  270   CONTINUE

C       Schnelle Fouriertransformation:Aufruf von DFORK
        CALL DFORK(SL,g,hin)

C       spektrale Anregungsfunktionen fuer Verschiebung,Geschwindigkeit
C       oder Beschleunigung (mit Beruecksichtigung der Aperidiozitaet)
C       fuer das Teilspektrum bis zur Nyqistfrequenz,zusaetzlich wird 
C       der Frequenzfilter ftap(f) angelegt.
        DO 280 f=1,Ny
C          Anregung fuer Bodenverschiebung(spektrale Integration)
C          mit Zuweisung fuer Frequenz und Kreisfrequenz 0 
           IF(typ.EQ.1) THEN
             fr(1)=1.D-30
             w(1)=2.D0*pi*fr(1)
             g(f)=ap*g(f)/(IME*w(f)) *ftap(f)
C          Anregung fuer Bodengeschwindigkeit
           ELSE IF(typ.EQ.2) THEN
             g(f)=ap*g(f) *ftap(f)
C          Anregung fuer Bodenbeschleunigung(spektrale Differentiation)
           ELSE IF(typ.EQ.3) THEN
             g(f)=ap*g(f)*IME*w(f) *ftap(f)
           ENDIF
  280   CONTINUE
      ENDIF






C**********************************************************************C
C 10. Ausgabe des aufbereiteten Erdmodells auf refseis.out             C 
C**********************************************************************C
	
C     Seismogramminfo
      WRITE(2,290)
  290 FORMAT(1x,'*** SEISMOGRAM INFORMATION ***********************')
C     Kopfzeile ausdrucken
      WRITE(2,291) text0
  291 FORMAT(1X,A)
C     Seismogrammtyp  ausdrucken
      impuls0=' '
      IF(The.EQ.0 .and. Thd.EQ.0) impuls0='-Delta excitation'
      IF(typ.EQ.1) typ0='Displacement'
      IF(typ.EQ.2) typ0='Velocity'       
      IF(typ.EQ.3) typ0='Acceleration'  
      WRITE(2,292) typ0,impuls0
  292 FORMAT(1X,'Seismogram type: ',2A)
C     Reduktionsdaten ausdrucken
      WRITE(2,295) Vred,Tre,Tli
  295 FORMAT(1X,'Reduction with:  Vred=',F8.4,'    Tre=',F6.2,
     &          '     Tli=',F6.2/)

C     Schichtmodell ausdrucken
      WRITE(2,299)
  299 FORMAT(1X,'*** VELOCITY MODEL *******************************')
      WRITE(2,300)
  300 FORMAT(1X,'Layers  Depth[km] Thick[km] Alpha[km/s] Beta[km/s]',
     &          ' Rho[g/ccm] Q-alpha  Q-beta')
      WRITE(2,310)(i,z(i),d(i),alpha(i),beta(i),rho(i),
     &             Qa(i),Qb(i), i=0,n)
  310 FORMAT(1X,I7,2F10.4,F12.4,2F11.4,2f10.0)  

C     Herdmodell und Momententensor ausdrucken
      WRITE(2,320)
  320 FORMAT(/,1X,'*** SOURCE MODEL **********************',
     &          '*** MOMENT TENSOR ***')
      WRITE(2,330) ZQ,Mxx
  330 FORMAT(1X,'Depth[km] :',F8.4,T45,' Mxx=',F7.3)
      WRITE(2,340) h,Myy
  340 FORMAT(1X,'in  Layer :',I3,T45,' Myy=',F7.3)    
      WRITE(2,350) M0,Mzz
  350 FORMAT(1X, 'Strength Mo[dyn*cm]:',E11.4,T45,' Mzz=',F7.3)
      WRITE(2,360) typ0,impuls0,Mxy 
  360 FORMAT(1X,'Excitation: ',2A,T45,' Mxy=',F7.3)
      WRITE(2,365) The,Mxz
  365 FORMAT(1X, 'Onset    The[s]:',F10.4,T45,' Mxz=',F7.3)
      WRITE(2,370) Thd,Myz
  370 FORMAT(1X, 'Duration Thd[s]:',F10.4,T45,' Myz=',F7.3/)

C     Empfaengermodell ausdrucken
      WRITE(2,375)
  375 FORMAT(1X,'*** RECEIVER PARAMETERS **************************')
      WRITE(2,380)
  380 FORMAT(1X,'Station  Epicentral dis.[km]  Azimuth[Grad]')
      WRITE(2,390)(E,r(E),phi(E), E=1,NE)
  390 FORMAT(1X,I7,F21.4,F15.4)

C     Langsamkeitsdaten ausdrucken
      WRITE(2,400)
  400 FORMAT(1X,/' *** SLOWNESS PARAMETERS **************************')
      WRITE(2,410)
  410 FORMAT(1X,'umin[s/km]  uwil[s/km]  uwir[s/km]  umax[s/km]')
      WRITE(2,420) umin,uwil,uwir,umax
  420 FORMAT(1X,F10.6,3F12.6)
      WRITE(2,425) Nu
  425 FORMAT(1X,'No. evenly-spaced slownesses in [umin,umax] Nu:',I5/)

C     Kontrollausgabe auf Bildschirm
      WRITE(*,400)
      WRITE(*,410)
      WRITE(*,420) umin,uwil,uwir,umax
      WRITE(*,425) Nu

C     Seismogrammdaten ausdrucken
      WRITE(2,460)
  460 FORMAT(1X,'*** SEISMOGRAM DATA ******************************')
      WRITE(2,470) SL,Fny
  470 FORMAT(1X,'Seismogram length SL:',I10,T37,
     &          'Nyquist frequency Fny[Hz] :',F10.5)
      WRITE(2,480) Dt,Df
  480 FORMAT(1X,'Sampling intv. Dt[s]:',F10.5,T37,
     &          'Frequency interval   Df[Hz]:',F10.5)
      WRITE(2,485) TL,fmax
  485 FORMAT(1X,'Time length TL[s]   :',F10.5,T37,
     &          'Cutoff freq.          [Hz]:',F10.5/)
      write(2,488) n,Nu,Nf
  488 format(i6,' Layers,',i6,' Slownesses,',i6,' Frequencies'/)
C     Kontrollausgabe auf Bildschirm
      WRITE(*,490)
  490 FORMAT(1X,'*** SEISMOGRAM DATA ******************************')
      WRITE(*,492) SL,Fny
  492 FORMAT(1X,'Seismogram length SL:',I10,T37,
     &          'Nyquist frequency Fny[Hz] :',F10.5)
      WRITE(*,494) Dt,Df
  494 FORMAT(1X,'Sampling intv. Dt[s]:',F10.5,T37,
     &          'Frequency interval  Df[Hz]:',F10.5)
      WRITE(*,496) TL,fmax
  496 FORMAT(1X,'Time length TL[s]   :',F10.5,T37,
     &          'Cutoff freq.          [Hz]:',F10.5/)
      write(*,488) n,Nu,Nf
      if(fmax.gt.fny) stop "fmax exceeds nyquist frequency - stop"

C**********************************************************************C
C                                                                      C
C Hauptteil des Programms:Berechnung des Seismogrammes                 C
C                                                                      C
C***********************************************************************

C~~~~~Grosse Langsamkeitsschleife ueber u ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C
      DO 4000 l=1,Nu
C----------------------------------------------------------------------C
C     Abschaetzung der Rechendauer: auf dem Bildschirm wird angezeigt  C 
C     die wievielte Langsamkeit berechnet wird(kann kommentiert werden)C
      if(l/500.ne.(l-1)/500.or.l.eq.100) then
      call CPU_TIME(tarray(2))
      estim=int((nu-l)*(tarray(2)-tarray(1))/l)+1
      WRITE(*,500) l, estim
  500 FORMAT(1X,I5,'-th slowness. Remaining time: ',i5,' sec')
      endif
C----------------------------------------------------------------------C


C**********************************************************************C
C 11.Vorberechnungen von Ausdruecken derselben Langsamkeit u:          C
C     -Langsamkeitsparameter u, und sein Quadrat uQ                    C
C     -vertikale Langsamkeiten aller Schichten                         C
C     -Reflektions- und Transmissionskoeffizienten aller Schichten     C
C     -Fensterfunktion utap der Langsamkeit mit Cos-Taper              C
C     -Gewichtungsfaktor gew und Multiplikationsfaktor DDu fuer Trapez-C
C      regel                                                           C
C**********************************************************************C

C     Langsamkeitsparameter u
      u=umin+DBLE(l-1)*Du
      uQ = u*u

C     Berechnung der vert. Langsamkeiten aller Schichten i:
C     Der Realteil soll fuer w>0 positiv, der Imaginaerteil negativ 
C     sein.Die Wurzelfunktion CSQRT arbeitet so, dass der Realteil der
C     gezogenen Wurzel stets positiv ist und der Imaginaerteil mit dem
C     Vorzeichen versehen wird,das der Imaginaerteil des Arguments hat.
C     Weil der Imaginaerteil des Arguments zur Bildung der vertikalen
C     Langsamkeiten negativ ist, wird obige Forderung erfuellt.
      DO 590 i=0,n
         H11= 1.D0/(alphC(i)*alphC(i))
         H22= 1.D0/(betaC(i)*betaC(i))
         a(i) = sqrt(H11-uQ)
         b(i) = sqrt(H22-uQ)
  590 CONTINUE

C     Berechnung er Reflektions und Transmissionskoeffizienten fuer
C     die Trennflaeche i+1 (Aufruf von RTKC)
      DO 600 i=0,n-1
         CALL RTKC (u,a(i),b(i),betaC(i),rho(i),
     &              a(i+1),b(i+1),betaC(i+1),rho(i+1),
     &              Rppd(i+1),Rpsd(i+1),Tppd(i+1),Tpsd(i+1),Rssd(i+1),
     &              Rspd(i+1),Tssd(i+1),Tspd(i+1),rd(i+1),td(i+1),
     &              Rppu(i+1),Rpsu(i+1),Tppu(i+1),Tpsu(i+1),Rssu(i+1),
     &              Rspu(i+1),Tssu(i+1),Tspu(i+1),ru(i+1),tu(i+1))
  600 CONTINUE           


C     Fensterfunktion utap der Langsamkeit mit Cos-Taper
      IF (umin.LE.u .and. u.LT.uwil) THEN
         utap=0.5D0*(1.D0-dcos(pi*(u-umin)/(uwil-umin)))
      ELSE IF(uwil.LE.u .and. u.LE.uwir) THEN
         utap=1.D0
      ELSE IF (uwir.LT.u .and. u.LE.umax) THEN
         utap=0.5D0*(1.D0+dcos(pi*(u-uwir)/(umax-uwir)))
      ENDIF

C     Gewichtungsfaktor fuer Trapezregel.Der erste und letzte Funk-
C     tionswert werden nur halb gewichtet(gew=0.5)
      gew=1.D0
      IF ( l.EQ.1  .OR. l.EQ.Nu ) gew=0.5D0
C     getaperte und gewichtete Integrationsschrittweite fuer Trapezregel
      DDu=DCMPLX(Du*gew*utap)



C-----Kopfzeile der Ausgabe fuer Schalter 1,2,3------------------------C
      IF(Sch(1).NE.0 .OR. Sch(2).NE.0 .OR. Sch(3).NE.0) THEN
      WRITE(2,700) u
  700 FORMAT(/////1X,'horizontal slowness u=',F10.6,'*************',
     &               '*********************************'/' ')
      ENDIF
C----------------------------------------------------------------------C


C-----Schalter 1 ------------------------------------------------------C
      IF(Sch(1).NE.0) THEN
C     Testausgabe der vertikalen Langsamkeiten ,sowie der Reflektions- C
C     und Transmissionskoeffizienten derTrennflaechen.                 C
C     Schleife ueber die Schichten 
      DO 790 i=0,n-1
C        vert. Langsamkeiten der Schichten i und i+1
         WRITE(2,710) i,i+1
  710    FORMAT(' vertical slowness of layers',I3,'  and',I3)            
         WRITE(2,720)i,a(i),i+1,a(i+1),i,b(I),i+1,b(i+1)
  720    FORMAT(' a',I3,'=',2F10.6,'   a',I3,'=',2F10.6/
     &          ' b',I3,'=',2F10.6,'   b',I3,'=',2F10.6/' ')
C        Reflektions- u.Transmissionskoeffizienten der Trennflaechen                 
         WRITE(2,730) i+1
  730    FORMAT(' Reflection and Transmission coefficients of',    
     &          ' Interface',I3)   
         WRITE(2,740)Rppd(i+1),Rssd(i+1),rd(i+1),Rpsd(i+1),Rspd(i+1),
     &                Tppd(i+1),Tssd(i+1),td(i+1),Tpsd(i+1),Tspd(i+1)
  740    FORMAT(' Rppd=',2F10.6,'  Rssd=',2F10.6,'  rd=',2F10.6/
     &          ' Rpsd=',2F10.6,'  Rspd=',2F10.6/ 
     &          ' Tppd=',2F10.6,'  Tssd=',2F10.6,'  td=',2F10.6/
     &          ' Tpsd=',2F10.6,'  Tspd=',2F10.6/' ')   
         WRITE(2,750)Rppu(i+1),Rssu(i+1),ru(i+1),Rpsu(i+1),Rspu(i+1),
     &                Tppu(i+1),Tssu(i+1),tu(i+1),Tpsu(i+1),Tspu(i+1)
  750    FORMAT(' Rppu=',2F10.6,'  Rssu=',2F10.6,'  ru=',2F10.6/
     &          ' Rpsu=',2F10.6,'  Rspu=',2F10.6/ 
     &          ' Tppu=',2F10.6,'  Tssu=',2F10.6,'  tu=',2F10.6/
     &          ' Tpsu=',2F10.6,'  Tspu=',2F10.6/' '/' '/' ')   
  790 CONTINUE
C     Ende der Schichtschleife
      ENDIF
      IF(Sch(1).EQ.2) GOTO 4000
C----------------------------------------------------------------------C






C~~~~~Grosse Frequenzschleife ueber w ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C
      DO 3000 f=fmi,fma


C**********************************************************************C
C 12. Vorberechnung von Ausdruecken derselben Frequenz w:              C
C     -Phasenmatrizen                                                  C
C**********************************************************************C

C     Phasenmatrizen der Schichten 0 bis n-1 
      DO 900 i=0,n-1
         hilf=-IME*d(i)*w(f)
         earg=a(i)*hilf
         E11(i)=exp(earg)
         earg=b(i)*hilf
         E22(i)=exp(earg)
  900 CONTINUE




C**********************************************************************C
C 13. Reflektivitaeten RM(inus) und rm(inus) fuer Welleneinfall von    C
C     oben (pro Langsamkeit):                                          C
C     Berechnet wird im Rekursionsverfahren die Reflektivitaetsmatrix  C
C     RTD(h),sowie die SH-Reflektivitaet rtd(h) fuer die Herdschicht h,C
C     das sind aber gerade RM(inus) und rm(inus).                      C
C     Die Reflektivitaeten bis zu den jeweiligen Schichten i werden    C
C     fuer die aktuelle Langsamkeit zwischengespeichert, so dass auch  C
C     alle RTD(i) und rt(i) bekannt sind (dadurch auch erweiterbar auf C
C     mehrere Herde).                                                  C
C**********************************************************************C

C     Startbedingung fuer Rekursion
C     fuer RTD 
      RTD11(n)=(0.,0.)
      RTD12(n)=(0.,0.)
      RTD21(n)=(0.,0.)
      RTD22(n)=(0.,0.)
C     fuer rtd
      rtd(n)  =(0.,0.)

C     Rekursion von Schicht n-1 bis h  
      DO 1000 i=n-1,h,-1

C        Rekursion fuer RTD
         CALL MAMULT(RTD11(i+1),RTD12(i+1),RTD21(i+1),RTD22(i+1),
     &               Rppu(i+1) ,Rspu(i+1) ,Rpsu(i+1) ,Rssu(i+1) ,
     &               H11       ,H12       ,H21       ,H22       )
         H11=1.D0-H11
         H12=  -H12
         H21=  -H21
         H22=1.D0-H22
         CALL MATINV(H11,H12,H21,H22, I11,I12,I21,I22)
         CALL MAMULT(Tppu(i+1),Tspu(i+1),Tpsu(i+1),Tssu(i+1),
     &               I11      ,I12      ,I21      ,I22      ,
     &               H11      ,H12      ,H21      ,H22      )         
         CALL MAMULT(H11       ,H12       ,H21       ,H22       ,
     &               RTD11(i+1),RTD12(i+1),RTD21(i+1),RTD22(i+1),
     &               I11       ,I12       ,I21       ,I22       )
         CALL MAMULT(I11      ,I12      ,I21      ,I22      ,
     &               Tppd(i+1),Tspd(i+1),Tpsd(i+1),Tssd(i+1),
     &               H11      ,H12      ,H21      ,H22      )
         H11=Rppd(i+1)+H11
         H12=Rspd(i+1)+H12
         H21=Rpsd(i+1)+H21
         H22=Rssd(i+1)+H22
         RTD11(i)=E11(i)*E11(i)*H11
         RTD12(i)=E11(i)*E22(i)*H12
         RTD21(i)=E11(i)*E22(i)*H21
         RTD22(i)=E22(i)*E22(i)*H22
C        Ende Rekursion RTD

C        Rekursion fuer rtd 
         hilf =tu(i+1)*rtd(i+1)*td(i+1)/(1.D0-rtd(i+1)*ru(I+1))
         rtd(i)=(rd(i+1)+hilf)*E22(i)*E22(i)

 1000 CONTINUE
C     Ende Rekursion fuer RTD und rtd 


C     Reflektivitaeten RM(inus) und rm(inus) der Herdschicht
      RM11=RTD11(h)
      RM12=RTD12(h)
      RM21=RTD21(h)
      RM22=RTD22(h)
      rm  =rtd(h)




C**********************************************************************C
C 14. Reflektivitaeten RP(lus) und rp(lus) sowie Transmissivitaeten    C
C     TP(lus) und tp(plus),fuer Welleneinfall von unten:               C
C     Berechnet werden im Rekursionsverfahren die Reflektivitaeten     C
C     RTU(h), rtu(h), sowie die Transmissivitaeten TTUo(h),ttuo(h) des C
C     oberen Stapelbereichs (Schicht 0 bis h-1); das sind aber gerade  C
C     RP,rp und TP,tp.                                                 C
C     Die Reflektivitaeten und Transmissivitaeten ab Decke Schicht 0   C
C     (Bezugsnivaeu des oberen Halbraumes) bis zu den Schichten i      C
C     werden zwischengespeichert.                                      C
C**********************************************************************C

C     Startbedingung fuer Rekursion
      RBU11(0)=(0.,0.)
      RBU12(0)=(0.,0.)
      RBU21(0)=(0.,0.)
      RBU22(0)=(0.,0.)
      rbu(0)=(0.,0.)

      TTUo11(0)=(1.,0.)
      TTUo12(0)=(0.,0.)
      TTUo21(0)=(0.,0.)
      TTUo22(0)=(1.,0.)
      ttuo(0)=(1.,0.)


C     Rekursion von Schicht 0 bis Schicht h-1
      DO 1100 i=0,h-1

C        Rekursion fuer TTUo(i) zuerst
         CALL MAMULT (Rppd(i+1),Rspd(i+1),Rpsd(i+1),Rssd(i+1),
     &                RBU11(i) ,RBU12(i) ,RBU21(i) ,RBU22(i) ,
     &                H11      ,H12      ,H21      ,H22      )
         H11 =1.D0-H11
         H12 =  -H12
         H21 =  -H21
         H22 =1.D0-H22
         CALL MATINV (H11,H12,H21,H22,I11,I12,I21,I22)
         H11 =E11(i)*I11
         H12 =E11(i)*I12
         H21 =E22(i)*I21
         H22 =E22(i)*I22
         CALL MAMULT(H11      ,H12      ,H21      ,H22      ,
     &               Tppu(i+1),Tspu(i+1),Tpsu(i+1),Tssu(i+1),
     &               S11(i)   ,S12(i)   ,S21(i)   ,S22(i)   )
         CALL MAMULT(TTUo11(i),TTUo12(i) ,TTUo21(i)   ,TTUo22(i)  ,
     &               S11(i)   ,S12(i)    ,S21(i)      ,S22(i)     ,
     &             TTUo11(i+1),TTUo12(i+1),TTUo21(i+1),TTUo22(i+1)) 
C        Ende Rekursion TTUo 


C        Rekursion fuer RTU(i) [liefert auch RBU(i) fuer TTUo(i)]
         CALL MAMULT(RBU11(i) ,RBU12(i) ,RBU21(i) ,RBU22(i) ,
     &               Rppd(i+1),Rspd(i+1),Rpsd(i+1),Rssd(i+1),
     &               H11      ,H12      ,H21      ,H22      )
         H11=1.D0-H11
         H12=  -H12
         H21=  -H21
         H22=1.D0-H22
         CALL MATINV(H11,H12,H21,H22,I11,I12,I21,I22)
         CALL MAMULT(Tppd(i+1),Tspd(i+1),Tpsd(i+1),Tssd(i+1),
     &               I11,I12,I21,I22,       H11,H12,H21,H22 ) 
         CALL MAMULT(H11     ,H12     ,H21     ,H22     ,
     &               RBU11(i),RBU12(i),RBU21(i),RBU22(i),
     &               I11     ,I12     ,I21     ,I22     )
         CALL MAMULT(I11      ,I12      ,I21      ,I22      ,
     &               Tppu(i+1),Tspu(i+1),Tpsu(i+1),Tssu(i+1),
     &               H11      ,H12      ,H21      ,H22      )
         RTU11(i+1)=Rppu(i+1)+H11
         RTU12(i+1)=Rspu(i+1)+H12
         RTU21(i+1)=Rpsu(i+1)+H21
         RTU22(i+1)=Rssu(i+1)+H22

C        neue RBU werden auch gebraucht zur Rekursion von TTUo
         RBU11(i+1)=E11(i+1)*E11(i+1)*RTU11(i+1)
         RBU12(i+1)=E11(i+1)*E22(i+1)*RTU12(i+1)
         RBU21(i+1)=E11(i+1)*E22(i+1)*RTU21(i+1)
         RBU22(i+1)=E22(i+1)*E22(i+1)*RTU22(i+1) 
C        Ende Rekursion fuer  RTU
 
 
C        Rekursion fuer ttuo
         s(i)=tu(i+1)/(1.D0-rd(i+1)*rbu(i)) *E22(i)
         ttuo(i+1)=ttuo(i)*s(i)

C        Rekursion fuer rtu
         rtu(i+1)=ru(i+1)+td(i+1)*tu(i+1)*rbu(i)/(1.D0-rd(i+1)*rbu(i))
         rbu(i+1)=rtu(i+1)*E22(i+1)*E22(i+1)

 1100 CONTINUE
C     Ende der Rekursion fuer RTU,TTUO,rtu,ttuo



C     Reflektivitaeten RP(lus),rp(lus) der Herdschicht
      RP11 = RTU11(h)
      RP12 = RTU12(h)
      RP21 = RTU21(h)
      RP22 = RTU22(h)
        rp = rtu(h)
   
C     Transmissivitaeten TP(lus),tp(lus) der Herdschicht 
      TP11 = TTUo11(h)
      TP12 = TTUo12(h)
      TP21 = TTUo21(h)
      TP22 = TTUo22(h)
        tp = ttuo(h)



C-----Schalter 2 ------------------------------------------------------C  
C     Testausgabe der Reflektivitaeten und Transmissivitaeten          C
C     fuer die Herdschicht                                             C
      IF(Sch(2).NE.0) THEN
      WRITE(2,1200) fr(f),w(f)
 1200 FORMAT(//' Reflectivities and Transmissivities of the source'/
     &      'for frequency fr=',F10.5,' or circular frequency w=',F10.5)

C     Reflektivitaeten RM(inus),rm(inus)       
      WRITE(2,1210) RM11,RM12,RM21,RM22,rm
 1210 FORMAT(' RM11=',2F10.6,'    RM12=',2F10.6/
     &       ' RM21=',2F10.6,'    RM22=',2F10.6/
     &       ' rm  =',2F10.6/' ')  

C     Reflektivitaeten RP(lus),rp(lus)
      WRITE(2,1220) RP11,RP12,RP21,RP22,rp
 1220 FORMAT(' RP11=',2F10.6,'    RP12=',2F10.6/
     &       ' RP21=',2F10.6,'    RP22=',2F10.6/
     &       ' rp  =',2F10.6/' ')

C     Transmissivitaeten TP(lus),tp(lus)
      WRITE(2,1230) TP11,TP12,TP21,TP22,tp
 1230 FORMAT(' TP11=',2F10.6,'    TP12=',2F10.6/
     &       ' TP21=',2F10.6,'    TP22=',2F10.6/
     &       ' tp  =',2F10.6/' ')
      ENDIF
      IF(Sch(2).EQ.2) GOTO 3000
C----------------------------------------------------------------------C






C**********************************************************************C
C 15. Berechnung der Oberflaechenamplituden B0(i),D0(i) mit i=0,1,2    C
C     und F0(k) mit k=1,2.                                             C
C     Zuerst werden die Quellamplituden AQ(i),CQ(i),EQ(k) und BQ(i),   C
C     DQ(i),FQ(k) bestimmt und daraus dann die Oberflaechenamplituden  C
C     mit Hilfe der Reflektivitaeten und Transmissivitaeten des Herdes C
C     berechnet .                                                      C
C**********************************************************************C

C     Bestimmung der Amplitudenexponenten
      hilf=IME*w(f)*(ZQ-z(h))
      earg=a(h)*hilf
      ea=exp(earg)
      earg=b(h)*hilf
      eb=exp(earg)

C     Bestimmung der Quell-Amplituden
      AQ(0)=IME*a(h)*ea
      AQ(1)=2.D0*u*ea
      AQ(2)=IME*u*u/a(h)*ea 
      CQ(0)=IME*u*eb
      CQ(1)=(u*u/b(h)-b(h))*eb
      CQ(2)=-IME*u*eb
      EQ(1)=eb/(betaC(h)*betaC(h))
      EQ(2)=EQ(1)*IME*u/b(h)

      BQ(0)=IME*a(h)/ea
      BQ(1)=-2.D0*u/ea
      BQ(2)=IME*u*u/(a(h)*ea)
      DQ(0)=-IME*u/eb
      DQ(1)=(u*u/b(h)-b(h))/eb
      DQ(2)=IME*u/eb
      FQ(1)=-1.D0/(betaC(h)*betaC(h)*eb)
      FQ(2)=-FQ(1)*IME*u/b(h)


C     Berechnung der Oberflaechenamplituden B0(i),D0(i):
C     - zuerst Berechnung der inversen Matrix der unendlichen
C     Matrizenreihe
      CALL MAMULT(RM11,RM12,RM21,RM22,RP11,RP12,RP21,RP22,
     &            H11,H12,H21,H22)
      H11=1.D0-H11
      H12=  -H12
      H21=  -H21
      H22=1.D0-H22
      CALL MATINV(H11,H12,H21,H22,I11,I12,I21,I22)
C     - dann Multiplikation mit TP(lus) nach [H11,H12,H21,H22]
      CALL MAMULT(TP11,TP12,TP21,TP22,I11,I12,I21,I22,
     &            H11,H12,H21,H22)
C     - schliesslich Berechnung von B0(i),D0(i)
      DO 1300 i=0,2
C        Multiplikation von RM(inus) mit Vektor [AQ(i),CQ(i)] und
C        Addition zu Vektor [BQ(i),DQ(i)]; abgelegt auf [I11,I22].
         I11=BQ(i)+RM11*AQ(i)+RM12*CQ(i)
         I22=DQ(i)+RM21*AQ(i)+RM22*CQ(i)
C        Multiplikation der Matrix [H11,H12,H21,H22] mit [I11,I22]
         B0(i)=H11*I11+H12*I22
         D0(i)=H21*I11+H22*I22
 1300 CONTINUE


C     Berechnung der Oberflaechenamplitude F0(k)
      DO 1400 k=1,2
         F0(k)=tp/(1.D0-rm*rp)*(FQ(k)+rm*EQ(k))
 1400 CONTINUE
   


C-----Schalter 3 ------------------------------------------------------C
C     Testausgabe der Oberflaechenamplituden                           C
      IF(Sch(3).NE.0) THEN 

      WRITE(2,1500) fr(f),w(f)
 1500 FORMAT(1X,'Surface amplitudes of the wavefield for'/     
     &          'frequency fr=',F10.5,' or circular frequency w=',F10.5)
      WRITE(2,1510) B0(0),D0(0),F0(1),B0(1),D0(1),F0(2),B0(2),D0(2)
 1510 FORMAT(' B0(0)=',2F9.6,' D0(0)=',2F9.6,' F0(1)=',2F9.6/
     &       ' B0(1)=',2F9.6,' D0(1)=',2F9.6,' F0(2)=',2F9.6/
     &       ' B0(2)=',2F9.6,' D0(2)=',2F9.6/' '/)
      ENDIF
      IF(Sch(3).EQ.2) GOTO 3000
C----------------------------------------------------------------------C 


 

C**********************************************************************C
C 16. Vorberechnungen fuer die Berechnung der Integrandenglieder       C
C**********************************************************************C
      In1=u*B0(0)+b(0)*D0(0)
      In2=u*B0(2)+b(0)*D0(2)
      In3=u*B0(1)+b(0)*D0(1)
      In4=a(0)*B0(0)-u*D0(0)
      In5=a(0)*B0(2)-u*D0(2)
      In6=a(0)*B0(1)-u*D0(1)
      FFI=w(f)*w(f)/(4.D0*pi*rho(h))



C~~~~~Grosse Empfaengerschleife ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C
      DO 2000 E=1,NE


C**********************************************************************C
C 17. Berechnung der einzelnen Integrandenglieder der spektralen Ver-  C
C     schiebungskomponenten fuer eine Langsamkeit,aber fuer alle Fre-  C
C     quenzen und Empfaenger:                                          C
C       -zuerst Aufruf der Besselfunktionen                            C
C       -Fernfeldintegrandenglieder IDFr,IDFphi,IDFz                   C
C       -Nahfeldintegrandenglieder  IDNr,IDNphi,IDNz                   C
C***********************************************************************

C     Besselfunktionen aufrufen   
      jarg=u*w(f)*r(E)
      J0=d_j0(jarg)  
      J1=d_j1(jarg)
      J2=d_jn(2,jarg) 

C     Berechnung der Fernfeld-Integrandenglieder /modified 8-2-91 K.Koch/
      hilf=J1*(K0(E)*In1+K3(E)*In2)
      IDFr  = FFI*u*(K1(E)*In3*J0-hilf)
      IDFphi=-FFI*u*(F0(1)*L1(E)*J0+F0(2)*L2(E)*J1)
      hilf=J0*(K0(E)*In4+K3(E)*In5)
      IDFz  = FFI*IME*u*(hilf+K1(E)*In6*J1)

C     Berechnung der Nahfeld-Integrandenglieder /modified 8-2-91 K.Koch/
      NFI=w(f)/(r(E)*4.D0*pi*rho(h))
      H11=  (F0(1)-In3)*J1
      H22=2.D0*(F0(2)-In2)*J2
      IDNr  =NFI*(K1(E)*H11+K2(E)*H22)
      IDNphi=NFI*(L1(E)*H11+L2(E)*H22)
      IDNz  =NFI*IME*2.D0*K2(E)*In5*J1





C**********************************************************************C
C 18. Berechnung der Komponenten des normierten Impulsspektrums:       C
C     (Berechnung der Langsamkeitsintegrale)                           C
C       -Fernfeldkomponenten VFr(f,E),VFphi(f,E),VFz(f,E)              C
C       -Nahfeldkomponenten  VNr(f,E),VNphi(f,E),VNz(f,E)              C
C     Fuer jeden Empfaenger E und jede Frequenz f werden die Komponen- C
C     ten des Impulsspektrums aus ihren Integrandengliedern mittels    C
C     der Trapezregel "aufintegriert" (d.h.aufsummiert):               C
C     Int(x1,xN)[f(x)dx] = dx*[0.5*f(x1)+f(x2)+....+f(xN-1)+0.5*f(xN)] C
C     Die Summation fuer jeden Empfaenger und jede Frequenz erfolgt    C
C     ueber die Langsamkeit.Die vollstaendigen Spektren liegen vor     C
C     wenn die Langsamkeitsschleife abgearbeitet ist.                  C
C       -zusaetzlich wird das Langsamkeitsfenster utap draufmultipli-  C
C        ziert (DDu=Du*gew*utap siehe Nr.11)                           C
C**********************************************************************C

C     Aufsummierung zu Fernfeldspektren    / modified 8-2-91 K.Koch /
      VFr(f,E)  =VFr(f,E)  +  IDFr*DDu
      VFphi(f,E)=VFphi(f,E)+IDFphi*DDu
      VFz(f,E)  =VFz(f,E)  +  IDFz*DDu

C     Aufsummierung zu Nahfeldspektren     / modified 8-2-91 K.Koch /
      VNr(f,E)  =VNr(f,E)  +  IDNr*DDu
      VNphi(f,E)=VNphi(f,E)+IDNphi*DDu
      VNz(f,E)  =VNz(f,E)  +  IDNz*DDu

 

 2000 CONTINUE
C~~~~~Ende der grossen Empfaengerschleife~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C


 3000 CONTINUE
C~~~~~Ende der grossen Frequenzschleife~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~C


 4000 CONTINUE
C~~~~~Ende der grossen Langsamkeitsschleife~~~~~~~~~~~~~~~~~~~~~~~~~~~~C
  
      IF (Sch(1).EQ.2) GOTO 99999
      IF (Sch(2).EQ.2) GOTO 99999
      IF (Sch(3).EQ.2) GOTO 99999


C-----Schalter 4 ------------------------------------------------------C
      IF(Sch(4).NE.0) THEN
C Testausgabe des normierten Teil-Impulsspektrums                      C
      WRITE(2,4110)
 4110 FORMAT(1X,////,'*** normalized impulse spectra in [s*s/(g*cm)]',
     &               '*1E-20  ********************') 
      DO 4150 E=1,NE
         WRITE(2,4120) E
 4120    FORMAT(1X,//'normalized impulse spectrum at recerver',I4)        
         DO 4140 f=fmi,fma
            WRITE(2,4130) fr(f),VFr(f,E),VFphi(f,E),VFz(f,E),
     &                           VNr(f,E),VNphi(f,E),VNz(f,E)
 4130       FORMAT(1X,'Frequency fr [Hz]=',F12.6/
     &                ' VFr=',2F9.5,' VFphi=',2F9.5,' VFz=',2F9.5/
     &                ' VNr=',2F9.5,' VNphi=',2F9.5,' VNz=',2F9.5)
 4140    CONTINUE
 4150 CONTINUE
      ENDIF
      IF(Sch(4).EQ.2) GOTO 99999
C----------------------------------------------------------------------C





C**********************************************************************C
C 19. vollstaendiges Spektrum fuer die Quelle der Staerke M0:          C
C     -Spektren SFr,SFphi,SFz und SNr,SNphi,SNz durch Multiplikation   C
C      mit der spektralen und gefilterten Anregungsfunktion g(f).      C
C     -Die Dimensionierung wird durch M0dim fuer die Verschiebungs-    C
C      dichte auf [cm] festgelegt.                                     C
C     -Reduktion der Spektren durch Multiplikation mit dem Reduktions- C
C      faktor red (bewirkt spaeter die Zeitverschiebung der Seismo-    C
C      gramme),zusaetzlich wird die Reduktionszeit tred(f) pro Em-     C
C      pfaenger ermittelt, um die Zeitskalierung anzupassen.           C
C     -Anschliessend Ergaenzung zu vollstaendigen komplexen Spektren   C
C     -dann Betraege der vollstaendigen komplexen Spektren:            C
C         Fernfeld:BSFr,BSFphi,BSFz                                    C 
C         Nahfeld :BSNr,BSNphi,BSNz                                    C
C**********************************************************************C

C   The section between Comment-Box (19.) and (21.) was redone to avoid
C   two-dimensional arrays. It saves a total of several MBytes of memory
C   (in the case of MSL=2048, ME=20,  the value is 8 MB !!!!!!!!)
C                                          / modified 8-2-91  K.Koch /


C-----Schalter 5 ------------------------------------------------------C
C Ausgabe der endgueltigen spektralen Komponenten (Betraege)           C
      IF(Sch(5).NE.0) THEN
      WRITE(2,6000)
 6000 FORMAT(1X,////,'*** Spectrum of strength MO', 
     &               ' in [cm*s],[cm],[cm/s)] *******************') 
      ENDIF

C-----Schalter 6: Plot-File fuer Spektren -----------------------------C
      IF(Sch(6).NE.0) THEN
      OPEN(UNIT=3,FILE='refspec.plt')
C Ausgabe als Plot-File                                                C
      WRITE(3,6320) text0
 6320 FORMAT(1X,A70)
      WRITE(3,6330) SL,NE,fr(1),fr(SL),Vred,typ
 6330 FORMAT(1x,'SL=',I5,' NE=',I3,' f(1)=',F10.4,' f(SL)=',F10.4,
     &          ' Vred=',F10.4,' typ=',I1)
      ENDIF

C-----Schalter 7 ------------------------------------------------------C
C Ausgabe der berechneten Seismogrammkomponenten                       C
      IF(Sch(7).NE.0) THEN
      WRITE(2,7100)
 7100 FORMAT(1X,////,'*** Seismogram of the source of strength M0 in',
     &               ' [cm,cm/s,cm/(s*s)] ******') 
      ENDIF

C-----Schalter 8: Plot-File fuer Seismogramme -------------------------C
      IF(Sch(8).NE.0) THEN
      OPEN(UNIT=4,FILE='refseis.sig')
C Ausgabe als Plot-File                                                C
      WRITE(4,7520) text0
 7520 FORMAT(1X,A70)
C     Reduktion der Spektren nach Verschiebungssatz
C        Bei Vred=0 keine Empfaengerreduktion
         IF(Vred.EQ.0.)THEN
           help=0.
         ELSE
C          Reduktion bezueglich Hypozentralenfernung
           help=DSQRT(ZQ*ZQ+r(E)*r(E))/Vred
         ENDIF
      thilf=help+Tli-Tre
      WRITE(4,7530) SL,NE,t(1)+thilf,t(SL)+thilf,Vred,typ,ZQ
 7530 FORMAT(1x,'SL=',I5,' NE=',I3,' t(1)=',F10.4,' t(SL)=',F10.4,
     &          ' Vred=',F10.4,' typ=',I1,' ZQ=',F10.4)
      ENDIF


C   Here we use only one loop {8000} instead of several loops over the 
C   receivers as done in the original program
C                                          / modified 8-2-91  K.Koch /

      DO 8000 E=1,NE

         DO 188 f=1,SL/2+1
            SFr(f)  =(0.,0.)
            SFphi(f)=(0.,0.)
            SFz(f)  =(0.,0.)
            SNr(f)  =(0.,0.)
            SNphi(f)=(0.,0.)
            SNz(f)  =(0.,0.)
  188    CONTINUE
         DO 5000 f=fmi,fma
C           Fernfeld-Spektren 
            SFr(f)  =M0dim*VFr(f,E)*g(f)
            SFphi(f)=M0dim*VFphi(f,E)*g(f)
            SFz(f)  =M0dim*VFz(f,E)*g(f)
C           Nahfeld-Spektren
            SNr(f)  =M0dim*VNr(f,E)*g(f)
            SNphi(f)=M0dim*VNphi(f,E)*g(f)
            SNz(f)  =M0dim*VNz(f,E)*g(f)
 5000    CONTINUE

C     Reduktion der Spektren nach Verschiebungssatz
C        Bei Vred=0 keine Empfaengerreduktion
         IF(Vred.EQ.0.)THEN
           help=0.
         ELSE
C          Reduktion bezueglich Hypozentralenfernung
           help=DSQRT(ZQ*ZQ+r(E)*r(E))/Vred
         ENDIF
         DO 5110 f=1,SL
C           Berechnung der Reduktionszeit tred(f)
            tred(f)=t(f)+help+Tli-Tre
C           und des Reduktionsfaktors red
            red=exp(IME*w(f)*(help+Tli-Tre))
C           Reduktion der Fernfeld-Spektren
            SFr(f)  =SFr(f)*red
            SFphi(f)=SFphi(f)*red
            SFz(f)  =SFz(f)*red
C           Reduktion der Nahfeld-Spektren
            SNr(f)  =SNr(f)*red
            SNphi(f)=SNphi(f)*red
            SNz(f)  =SNz(f)*red
 5110    CONTINUE

C     Ergaenzung zu vollstaendigen komplexen Spektren(Symmetrisierung
C     zur Nyquistfrequenz)
         DO 5200 f=0,SL/2-2
            SFr(SL-f)  =DCONJG(SFr(f+2))
            SFphi(SL-f)=DCONJG(SFphi(f+2))
            SFz(SL-f)  =DCONJG(SFz(f+2))
            SNr(SL-f)  =DCONJG(SNr(f+2))
            SNphi(SL-f)=DCONJG(SNphi(f+2))
            SNz(SL-f)  =DCONJG(SNz(f+2))
 5200    CONTINUE

C     Spektren (Betraege der komplexen Spektren)
         DO 5400 f=1,SL
            BSFr(f)  =abs(SFr(f))
            BSFphi(f)=abs(SFphi(f))
            BSFz(f)  =abs(SFz(f))
            BSNr(f)  =abs(SNr(f))
            BSNphi(f)=abs(SNphi(f))
            BSNz(f)  =abs(SNz(f))
 5400    CONTINUE


C-----Schalter 5 ------------------------------------------------------C
C Ausgabe der endgueltigen spektralen Komponenten (Betraege)           C
      IF(Sch(5).NE.0) THEN
         WRITE(2,6020) E
 6020    FORMAT(1X,//' Modulus of the spectrum at',
     &               ' receiver',I4)         
         WRITE(2,6030)
 6030    FORMAT(1X,'    fr[Hz]       BSFr     BSFphi       BSFz',
     &             '       BSNr     BSNphi       BSNz')
         DO 6100 f=1,SL 
            WRITE(2,6040) fr(f),BSFr(f),BSFphi(f),BSFz(f),
     &                           BSNr(f),BSNphi(f),BSNz(f)
 6040       FORMAT(1X,F10.6,6F11.6)
 6100    CONTINUE
      ENDIF
      IF(Sch(5).EQ.2) GOTO 8000
C----------------------------------------------------------------------C



C-----Schalter 6: Plot-File fuer Spektren -----------------------------C
      IF(Sch(6).NE.0) THEN
         WRITE(3,6340) r(E),phi(E)
 6340    FORMAT(1X,'r=',F10.4,' phi=',F6.2)
         DO 6390 f=1,SL 
            WRITE(3,6350)fr(f),SNGL(BSFr(f)),SNGL(BSFphi(f)),
     &                         SNGL(BSFz(f)),SNGL(BSNr(f)),
     &                         SNGL(BSNphi(f)),SNGL(BSNz(f))
 6350       FORMAT(7(1X,E10.4))
 6390    CONTINUE
      ENDIF
      IF(Sch(6).EQ.2) GOTO 8000
C----------------------------------------------------------------------C




C**********************************************************************C
C 20. Fouriertransformation des vollstaendigen Spektrums in den        C
C     Zeitbereich (Ruecktransformation) zum Seismogramm; jeweils       C
C     pro Empfaenger.                                                  C
C      -zur Fouriertransformation in den Zeitbereich werden die        C
C       spektralen Komponenten SFr(f)....SNz(f) mit dem Aufruf von     C
C       DFORK transformiert, die nun als Seismogrammkomponenten gelten C
C      -dabei Beruecksichtigung der Aperidiozitaet durch den Faktor ap.C
C      -die z-Komponenten werden negiert, damit sich die Seismogramm-  C
C       ausschlaege auf die Oberflaechennormale(nach oben) beziehen    C
C***********************************************************************


C        Ruecktransformation in den Zeitbereich:Seismogramme
         CALL DFORK(SL,SFr,rueck)
         CALL DFORK(SL,SFphi,rueck)
         CALL DFORK(SL,SFz,rueck)
         CALL DFORK(SL,SNr,rueck)
         CALL DFORK(SL,SNphi,rueck)
         CALL DFORK(SL,SNz,rueck)                

C        Seismogrammkomponenten S..(f) mit Beruecksichtigung der
C        Aperidiozitaet durch ap
C        Bem: die z-Komponenten werden negiert,damit die Ausschlaege
C        sich auf die Oberflaechennormale(nach oben) beziehen
         DO 6600 f=1,SL
            SFr(f)  =SFr(f)*1.D0/ap
            SFphi(f)=SFphi(f)*1.D0/ap
            SFz(f)  =-SFz(f)*1.D0/ap
            SNr(f)  =SNr(f)*1.D0/ap
            SNphi(f)=SNphi(f)*1.D0/ap
            SNz(f)  =-SNz(f)*1.D0/ap
 6600    CONTINUE




C-----Schalter 7 ------------------------------------------------------C
C Ausgabe der berechneten Seismogrammkomponenten                       C
      IF(Sch(7).NE.0) THEN

         WRITE(2,7120) E
 7120    FORMAT(1X,//' Far-field Seismogram at Station    ',I4)
         WRITE(2,7130)
 7130    FORMAT(1X,'   tred[s]    Re(SFr)    Im(SFr)',
     &             '  Re(SFphi)  Im(SFphi)    Re(SFz)    Im(SFz)')
         DO 7200 f=1,SL 
            WRITE(2,7140) tred(f),SFr(f),SFphi(f),SFz(f)
 7140       FORMAT(1X,F10.6,6F11.6)
 7200    CONTINUE

         WRITE(2,7220) E 
 7220    FORMAT(1X,//' Near-field Seismogram at Station  ',I4)
         WRITE(2,7230)
 7230    FORMAT(1X,'   tred[s]    Re(SNr)    Im(SNr)',
     &             '  Re(SNphi)  Im(SNphi)    Re(SNz)    Im(SNz)')
         DO 7300 f=1,SL 
            WRITE(2,7240) tred(f),SNr(f),SNphi(f),SNz(f)
 7240       FORMAT(1X,F10.6,6F10.6)
 7300    CONTINUE

      ENDIF
      IF(Sch(7).EQ.2) GOTO 8000
C----------------------------------------------------------------------C




C-----Schalter 8: Plot-File fuer Seismogramme -------------------------C
      IF(Sch(8).NE.0) THEN
         WRITE(4,7540) r(E),phi(E)
 7540    FORMAT(1X,' r=',F10.4,' phi=',F6.2)
         DO 7600 f=1,SL 
            WRITE(4,7550) tred(f),SNGL(SFr(f)),SNGL(SFphi(f)),
     &                              SNGL(SFz(f)),SNGL(SNr(f)),
     &                              SNGL(SNphi(f)),SNGL(SNz(f))

C   Data format in the file 'refseis.sig' is changed here. D. Yang 95 98

 7550       FORMAT(1X,F8.3,6(1X,E10.4))
C 7550       FORMAT(7(1X,E10.4))
 7600    CONTINUE
      ENDIF
      IF(Sch(8).EQ.2) GOTO 8000
C----------------------------------------------------------------------C

 8000 CONTINUE


C-----Schalter 6: Plot-File fuer Spektren -----------------------------C
      IF(Sch(6).NE.0) CLOSE(UNIT=3)

C-----Schalter 8: Plot-File fuer Seismogramme -------------------------C
      IF(Sch(8).NE.0) CLOSE(UNIT=4)


C**********************************************************************C
C 21. Beendigung des Programms und Rechenzeit                          C
C**********************************************************************C 
c      time=dtime(tarray)
      call CPU_TIME(tarray(2))
      WRITE(2,97777) tarray(2)-tarray(1)
97777 FORMAT(1X,//'Computation Time [s]=',F10.2)
      WRITE(*,98888) tarray(2)-tarray(1)
98888 FORMAT(1X,//'Computation Time [s]=',F10.2)
99999 CLOSE(UNIT=1)
      CLOSE(UNIT=2)

      END
 
C *** Ende REFSEIS ****************************************************C





C23456789012345678901234567890123456789012345678901234567890123456789012
C******* Joachim Ungerer ********** Stuttgart *****17.3.90*************C
C                                                                      C
C Subroutine RTKC                                                      C
C                                                                      C
C Berechnet die Reflektions- und Transmissionskoeffizienten einer      C
C ebenen Trennflaeche zwischen zwei homogenen Halbraeumen; fuer        C
C komplexe Wellengeschwindigkeiten.                                    C
C                                                                      C
C Berechnung mit doppelter Genauigkeit:                                C
C   DOUBLE COMPLEX muss implementiert sein                             C
C                                                                      C
C Variablendokumentation:                                              C
C                                                                      C
C         Uebergabevariablen:                                          C
C         u       horizontale Langsamkeit                              C
C         a1,b1   vert.Langsamkeiten des oberen Halbraumes             C
C         betaC1  S-Wellengeschw. des      "        "     (komplex)    C
C         rho1    Dichte des               "        "                  C
C         a2,b2   vert.Langsamkeiten des unteren Halbraumes            C
C         betaC2  S-Wellengeschw. des      "        "     (komplex)    C
C         rho2    Dichte des               "        "                  C
C                                                                      C
C         Rueckgabevariablen:                                          C
C         Rppd,Rpsd,Tppd,Tpsd    fuer  P-Welle von oben                C 
C         Rssd,Rspd,Tssd,Tspd    fuer SV-Welle      "                  C
C         rd,td                  fuer SH-Welle      "                  C
C         Rppu,Rpsu,Tppu,Tpsu    fuer  P-Welle von unten               C
C         Rssu,Rspu,Tssu,Tspu    fuer SV-Welle      "                  C
C         ru,tu                  fuer SH-Welle      "                  C
C                                                                      C
C         Rechenvariablen:                                             C
C          ........        gehen aus der Deklaration und               C
C                          ihrer Berechnung hervor                     C
C                                                                      C
C**********************************************************************C


      SUBROUTINE RTKC(u,a1,b1,betaC1,rho1,a2,b2,betaC2,rho2,
     &                Rppd,Rpsd,Tppd,Tpsd,Rssd,Rspd,Tssd,Tspd,rd,td,
     &                Rppu,Rpsu,Tppu,Tpsu,Rssu,Rspu,Tssu,Tspu,ru,tu)


C Deklaration der Variablen ********************************************

C     Variablen fuer beide Faelle
      REAL*8      u,uQ,rho1,rho2,rho12

      COMPLEX*16  Rppd,Rpsd,Tppd,Tpsd,Rssd,Rspd,Tssd,Tspd,rd,td,
     &            Rppu,Rpsu,Tppu,Tpsu,Rssu,Rspu,Tssu,Tspu,ru,tu,
     &            betaC1,betaC2,a1,a2,b1,b2,a1b1,a2b2,a1b2,a2b1,ab,
     &            abm,rhoabm,C,C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,
     &            mue1,mue2,mue1b1,mue2b2,muebN,muebP

C     nur fuer Einfall von oben
      COMPLEX*16  D1D,D2D,DD,D1,D2,D3

C     nur fuer Einfall von unten
      COMPLEX*16  D1U,D2U,DU,U1,U2,U3          



C Berechnungen ********************************************************C

C Berechnungen der Rechenvariablen fuer beide Faelle
      uQ    = u*u
      mue1  = rho1*betaC1*betaC1
      mue2  = rho2*betaC2*betaC2
      mue1b1= mue1*b1
      mue2b2= mue2*b2
      muebN = mue1b1-mue2b2
      muebP = mue1b1+mue2b2
      rho12 = rho1*rho2
      a1b1  = a1*b1
      a2b2  = a2*b2
      a1b2  = a1*b2
      a2b1  = a2*b1
      ab    = a1b1*a2b2
      abm   = a1b2-a2b1
      rhoabm= 2.D0*rho12*abm
      C     = 2.D0*(mue1-mue2)
      C0    = C*uQ
      C1    = C0-rho1
      C2    = C0+rho2
      C3    = C1+rho2
      C4    = C2*a1-C1*a2
      C5    = C2*b1-C1*b2
      C6    = C3*C3*uQ
      C7    = C*C0*ab
      C8    = C1*C1*a2b2+rho12*a2b1
      C9    = C2*C2*a1b1+rho12*a1b2
      C10   = C3+C*a2b1
      C11   = C3+C*a1b2

C Koeffizienten fuer Einfall von oben
C     Rechenvariablen
      D1D   = C6+C8
      D2D   = C7+C9
      DD    = D1D+D2D
      D1    = 2.D0*a1/DD
      D2    = 2.D0*b1/DD
      D3    = C3*C2+C*C1*a2b2
C     Koeffizienten
      Rppd  = (D2D-D1D)/DD
      Rpsd  =-u*D1*D3
      Tppd  = rho1*D1*C5
      Tpsd  =-u*rho1*D1*C10
      Rssd  = Rppd-rhoabm/DD
      Rspd  = u*D2*D3
      Tssd  = rho1*D2*C4
      Tspd  = u*rho1*D2*C11
      rd    = muebN/muebP
      td    = rd + 1.D0

C Koeffizienten fuer Einfall von unten     
C     Rechenvariablen
      D1U   = C6+C9
      D2U   = C7+C8
      DU    = D1U+D2U
      U1    = 2.D0*a2/DU
      U2    = 2.D0*b2/DU
      U3    = C3*C1+C*C2*a1b1
C     Koeffizienten
      Rppu  = (D2U-D1U)/DU
      Rpsu  = u*U1*U3
      Tppu  = rho2*U1*C5
      Tpsu  =-u*rho2*U1*C11
      Rssu  = Rppu+rhoabm/DU
      Rspu  =-u*U2*U3
      Tssu  = rho2*U2*C4
      Tspu  = u*rho2*U2*C10
      ru    =-muebN/muebP
      tu    = ru + 1.D0

      RETURN
      END

C***Ende von RTKC *****************************************************C




C23456789012345678901234567890123456789012345678901234567890123456789012
C*** Subroutine fuer Fouriertransformation *****************************
C Die von Gerherd Mueller verwendetet schnelle Fouriertransformation   C
C FORK wurde umgeschrieben fuer DOUBLE COMPLEX                         C
C Es muessen implementiert sein: DOUBLE COMPLEX,DCMPLX,exp           C
C                                                                      C
C Zum Verfahren der schnellen Fouriertransformation(FFT) und zur Ar-   C
C beitsweise von FORK siehe G.Mueller: Digitale Signalverarbeitung I,  C
C Vorlesungsmanuskript.                                                C
C                                                                      C
C Variablen:                                                           C
C    LX       Seismogrammlaenge bzw. Anzahl der Stuetzstellen,Abtast-  C
C             werte des Seismogramms/Spektrums.Muss eine Zeier-Potenz  C
C             sein.                                                    C
C    cx(LX)   Feld auf dessen Realteil die Funktionswerte der Zeit-    C
C             funktion stehen und nach Transformation ihre Fourierko-  C
C             effizienten.                                             C
C    SIGNI    SIGNI=-1.D0 bedeutet Berechnung der Fourierkoeffizienten C
C             SIGNI=+1.D0 bedeutet Ruecktransformation                 C
C**********************************************************************C     

      SUBROUTINE DFORK(LX,CX,SIGNI)
      INTEGER     I,ISTEP,J,L,LX,M
      REAL*8      SC,PI,SIGNI
      COMPLEX*16  CX(LX),CARG,CW,CTEMP

      PI=3.14159265358979
      J=1
      SC=1.D0/DBLE(LX)
      SC=DSQRT(SC)
      DO 5  I=1,LX
      IF(I.GT.J) GOTO 2
      CTEMP=CX(J)*SC
      CX(J)=CX(I)*SC
      CX(I)=CTEMP
2     M=LX/2
3     IF(J.LE.M) GOTO 5
      J=J-M
      M=M/2
      IF(M.GE.1) GOTO3
5     J=J+M
      L=1
6     ISTEP=2*L
      DO 8  M=1,L
      CARG=DCMPLX(0.,1.)*(PI*SIGNI*DBLE(M-1))/DBLE(L)
      CW=exp(CARG)
      DO 8 I=M,LX,ISTEP
      CTEMP=CW*CX(I+L)
      CX(I+L)=CX(I)-CTEMP
8     CX(I)=CX(I)+CTEMP
      L=ISTEP
      IF(L.LT.LX) GOTO 6
      RETURN
      END





C23456789012345678901234567890123456789012345678901234567890123456789012
C Subroutine fuer Matrizeninversion
C Bem: DOUBLE COMPLEX muss implementiert sein
 
      SUBROUTINE MATINV(A11,A12,A21,A22,B11,B12,B21,B22)
      COMPLEX*16   A11,A12,A21,A22,B11,B12,B21,B22,D
      D=1.D0/(A11*A22-A12*A21)
      B11= A22*D
      B12=-A12*D
      B21=-A21*D
      B22=A11*D
      RETURN
      END



C23456789012345678901234567890123456789012345678901234567890123456789012
C Subroutine fuer Matrizenmultiplikation
C Bem: DOUBLE COMPLEX muss implementiert sein

      SUBROUTINE MAMULT(A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,C21,C22)
      COMPLEX*16   A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,C21,C22
      C11=A11*B11+A12*B21
      C12=A11*B12+A12*B22
      C21=A21*B11+A22*B21
      C22=A21*B12+A22*B22
      RETURN
      END


      double precision function d_jn(n,x)
      implicit none
      INTEGER n,IACC
      REAL*8 bessj,x,BIGNO,BIGNI,d_j0,d_j1
      PARAMETER (IACC=40,BIGNO=1.d10,BIGNI=1.d-10)
CU    USES bessj0,bessj1
c     bessj0, bessj1 renamed to d_j0, d_j1
      INTEGER j,jsum,m
      REAL*8 ax,bj,bjm,bjp,sum,tox
      if(n.lt.2) pause 'bad argument n in bessj'
      ax=dabs(x)
      if(ax.eq.0.d0)then
        bessj=0.d0
      else if(ax.gt.float(n))then
        tox=2.d0/ax
        bjm=d_j0(ax)
        bj=d_j1(ax)
        do 11 j=1,n-1
          bjp=j*tox*bj-bjm
          bjm=bj
          bj=bjp
11      continue
        bessj=bj
      else
        tox=2.d0/ax
        m=2*((n+int(sqrt(float(IACC*n))))/2)
        bessj=0.d0
        jsum=0
        sum=0.d0
        bjp=0.d0
        bj=1.d0
        do 12 j=m,1,-1
          bjm=j*tox*bj-bjp
          bjp=bj
          bj=bjm
          if(dabs(bj).gt.BIGNO)then
            bj=bj*BIGNI
            bjp=bjp*BIGNI
            bessj=bessj*BIGNI
            sum=sum*BIGNI
          endif
          if(jsum.ne.0)sum=sum+bj
          jsum=1-jsum
          if(j.eq.n)bessj=bjp
12      continue
        sum=2.d0*sum-bj
        bessj=bessj/sum
      endif
      if(x.lt.0.d0.and.mod(n,2).eq.1)bessj=-bessj
      d_jn=bessj
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software "*3.


      double precision function d_j0(x)
      implicit none
      REAL*8 bessj0,x
      REAL*8 ax,xx,z
      DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
     *s1,s2,s3,s4,s5,s6,y
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6
      DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
     *-.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
     *.1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
      DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
     *651619640.7d0,-11214424.18d0,77392.33017d0,-184.9052456d0/,s1,s2,
     *s3,s4,s5,s6/57568490411.d0,1029532985.d0,9494680.718d0,
     *59272.64853d0,267.8532712d0,1.d0/
      if(abs(x).lt.8.d0)then
        y=x**2
        bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+y*
     *(s4+y*(s5+y*s6)))))
      else
        ax=dabs(x)
        z=8.d0/ax
        y=z**2
        xx=ax-.785398164
        bessj0=dsqrt(.636619772/ax)*(dcos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))-z*dsin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
      endif
      d_j0=bessj0
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software "*3.


      double precision function d_j1(x)
      implicit none
      REAL*8 bessj1,x
      REAL*8 ax,xx,z
      DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
     *s1,s2,s3,s4,s5,s6,y
      SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,
     *s5,s6
      DATA r1,r2,r3,r4,r5,r6/72362614232.d0,-7895059235.d0,
     *242396853.1d0,-2972611.439d0,15704.48260d0,-30.16036606d0/,s1,s2,
     *s3,s4,s5,s6/144725228442.d0,2300535178.d0,18583304.74d0,
     *99447.43394d0,376.9991397d0,1.d0/
      DATA p1,p2,p3,p4,p5/1.d0,.183105d-2,-.3516396496d-4,
     *.2457520174d-5,-.240337019d-6/, q1,q2,q3,q4,q5/.04687499995d0,
     *-.2002690873d-3,.8449199096d-5,-.88228987d-6,.105787412d-6/
      if(abs(x).lt.8.d0)then
        y=x**2
        bessj1=x*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))/(s1+y*(s2+y*(s3+
     *y*(s4+y*(s5+y*s6)))))
      else
        ax=dabs(x)
        z=8.d0/ax
        y=z**2
        xx=ax-2.356194491d0
        bessj1=dsqrt(.636619772/ax)*(dcos(xx)*(p1+y*(p2+y*(p3+y*(p4+y*
     *p5))))-z*dsin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))*dsign(1.d0,x)
      endif
      d_j1=bessj1
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software "*3.