C *** FOR BUGS OR COMMENTS, PLEASE E-MAIL to "hisada@cc.kogakuin.ac.jp", 
C     OR MAIL TO "Yoshiaki Hisada, Dept. of Architecture, 
C     Kogakuin Univeristy, Nishi-Shinjuku 1-24-2, Tokyo 163-91, Japan" 
c 
c *** This code compute phase velocities with small amplitude of the secular 
c     functions, which are used for computing Green's functions.  
c 
c *** Note that some of them are not the phase velocities of the surface waves, 
c     but the body waves.  
c 
c *** This code is open only for the academic use. 
c 
c     (originally made in 1992, and modified in Feb, 1996 from "normal.f") 
C 
      PROGRAM phs3x
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER(MOM=16400)
      COMMON /OUTPUT/ISEC,IVCT
C
      WRITE(*,1000) 
 1000 FORMAT(//' "phs3.f" computes the normal mode solutions'/ 
     *  '   (the phase and group velocities, and eigen vectors of'/ 
     *  '    the Love and Rayleigh waves).')  
      WRITE(*,*)' For trouble, please check "fort.10".' 
      WRITE(*,*)' This code is open for the academic use.' 
      WRITE(*,*) 
     *' For questions, e-mail to "hisada@cc.kogakuin.ac.jp".' 
      WRITE(*,*) 
     *' (Originally made in 1992, and modified in 2005)'
C 
      CALL INDAT(MOM) 
C 
      CALL PHASE(MOM)
      CALL ENGUU(MOM) 
C 
      WRITE(*,*) 
      WRITE(*,*)'Data "mdl_dat, mdr_dat, phsl_dat, and phsr_dat"' 
      WRITE(*,*)'are used for "grflt4.f" as data of poles'
 
      WRITE(*,*) 
      WRITE(*,*)'You can plot the dispersion curves using "phsl.m"' 
      WRITE(*,*)'and "phsr.m" for "M-file" of "Matlab"' 
      WRITE(*,*)
      write(*,*)'The followings are output data for Excel format:'
      WRITE(*,*)'  Ldisper.csv and Rdisper.csv: the dipersion curves'
c
      if(ISEC.eq.1)then
      WRITE(*,*)'  sec-Lov.csv: Secular Functions of Love Wave'
      WRITE(*,*)'  sec-Ray.csv: Secular Functions of Rayleigh Wave'
      WRITE(*,*)'  Ldis.csv: Eigen Displacement Vector of Love Wave'
      WRITE(*,*)'  Lstr.csv: Eigen Stress Vector of Love Wave'
      WRITE(*,*)'  Rdis.csv: Eigen Displacement Vector of Rayleigh Wave'
      WRITE(*,*)'  Rstr.csv: Eigen Stress Vector of Rayleigh Wave'
      end if
C
      STOP 
      END 
C 
C 
C *** DATA INPUT *** 
C 
      SUBROUTINE INDAT(MOM) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      character infile*20 
      COMMON /FREQ/DOM,IOM,NOM 
     *    /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
     *    /MINMAX/VSMIN,VSMAX,CCMIN,CCMAX,CMIN,CMAX,TOL,NITER,NCC
     *    /OUTPUT/ISEC,IVCT  /NINTG/ZINT,NZINT 
C 
C *** PARAMETERS IN COMMOM STATEMENT *** 
C 
      NML =20 
C 
C *** OMEGA DATA *** 
C 
      write(*,*) 
      write(*,*)'Enter Input File Name for "grflto2x-in.csv"'
      read(*,'(a20)')infile
c      infile='grflto2x-in.csv'
      OPEN (UNIT=7,FILE=infile,STATUS='OLD') 
C 
      read(7,*)
      read(7,*)
      read(7,*)dt,nt
      a=nt
 1    a=a/2 
      if(a.gt.1.d0)then 
        goto 1 
      else 
        if(a.ne.1.d0)then 
          write(*,*) 
          write(*,*)'NT is not a power of 2 !!' 
          write(10,*)'NT is not a power of 2 !!' 
          stop 
        end if 
      end if 
      read(7,*)
      read(7,*)period
      duration=dt*nt
      pi2=datan(1.d0)*8.d0 
      dom=pi2/duration 
      iom=1 
      nom=int(duration/period)
c 
      write(*,*) 
      write(*,*)'Maximum frequency (Hz) =',1.d0/period
      write(*,*)'Check Secular Functions & Eigen Vectors? (N=0 or Y=1)'
      read(*,*)ichk 
      if(ichk.eq.1)then 
        nom=1 
        write(*,*)'Enter the frequency for checking (Hz)'
        read(*,*)fchk 
        dom=fchk*pi2 
      end if 
c 
      WRITE(10,*)'FROM indat:' 
      WRITE(10,*) 
      WRITE(10,*)'dt,nt,period',dt,nt,period 
      WRITE(10,*)'IOM,NOM,DOM=',IOM,NOM,DOM 
C 
      IF(NOM.GT.MOM)THEN 
        WRITE(10,*)' NOM>MOM !!!!!!!!!!!!!' 
        STOP 
      END IF 
C 
C *** MEDIUM DATA *** 
C 
      read(7,*)
      read(7,*)
      READ(7,*)NL
      read(7,*)
      DO 100 IL=1,NL
        READ(7,*)I,DNS(IL),VP(IL),QP,FP,VS(IL),QS,FS,THK(IL)
  100 CONTINUE
      Qp=Qp
      Qs=Qs
C
      WRITE(10,2200)NL 
      DO 150 IL=1,NL 
        WRITE(10,2210)IL,DNS(IL),VP(IL),VS(IL),THK(IL) 
  150 CONTINUE 
 2200 FORMAT(//' Total Number of Layers =',I5/ 
     *       /' Layer, Density ,   Vp    ' 
     *       ,',   Vs   ,Thickness(m)')
 2210 FORMAT(I5,3F10.3,F12.3) 
C 
      IF(NL.GT.NML)THEN 
        WRITE(10,*)' NL>NML !!!!!!!!!!!!!' 
        STOP 
      END IF 
C 
      CLOSE(7) 
C 
C *** MIN & MAX VS AND CC *** 
C 
      VSMIN=VS(NL) 
      VSMAX=0.0D0 
      VPMAX=0.0D0  
      MINIL=1 
      MAXIL=1 
      DO 200 IL=1,NL 
        IF(VS(IL).LE.VSMIN)THEN 
          VSMIN=VS(IL) 
          MINIL=IL 
        END IF 
        IF(VS(IL).GE.VSMAX)THEN 
          VSMAX=VS(IL) 
          MAXIL=IL 
        END IF 
        VPMAX=MAX(VPMAX,VP(IL)) 
  200 CONTINUE 
C 
      IF(VSMAX.NE.VS(NL))THEN 
        WRITE(*,*)' FROM SUBROUTINE INDAT :' 
        WRITE(*,*)' VS(NL).NE.VSMAX, THUS STOP !!!' 
        STOP 
      END IF 
      IF(VPMAX.NE.VP(NL))THEN 
        WRITE(*,*)' FROM SUBROUTINE INDAT :' 
        WRITE(*,*)' VP(NL).NE.VPMAX, THUS STOP !!!' 
        STOP 
      END IF 
C 
      WRITE(*,1602) 
      READ(*,*)TOL 
c 
 1602 FORMAT(/' Enter Tolerance for Secular Function (ex., 0.01):')
C 
      WRITE(*,1603) 
      READ(*,*)NITER 
 1603 FORMAT(/' Enter Max Number of Iterations (ex., 200):')
C  
      CCMIN=FRHF(VP(MINIL),VS(MINIL)) 
      CCMAX=FRHF(VP(MAXIL),VS(MAXIL)) 
C 
      WRITE(*,1600)VSMAX,VPMAX 
      WRITE(*,1610)CCMIN,CCMAX 
 1600 FORMAT(/' VSMAX, VPMAX =',2F9.3) 
 1610 FORMAT( ' CCMIN, CCMAX =',2F9.3) 
c
      NCC=1 
c      NSC=1
      IF(NL.EQ.1)GOTO 25 
c      WRITE(*,1620) 
c      READ(*,*)IVC 
c      IF(IVC.EQ.1)THEN 
c        WRITE(*,1630) 
c        READ(*,*)CMIN,CMAX 
c      ELSE 
        CMIN=CCMIN-CCMIN/50.D0 
        CMAX=VSMAX+VSMAX/100.D0 
c      END IF 
c 1620 FORMAT( ' Change CMIN and CMAX? (yes=1):') 
c 1630 FORMAT( ' Enter CMIN and CMAX:') 
c
      WRITE(*,1640)CMIN,CMAX  
 1640 FORMAT(/' Min and Max Velocities (CMIN and CMAX: m/s) =',2F9.3) 
c
      WRITE(*,1650)
 1650 FORMAT( ' Enter Partition Number from CMIN to CMAX (ex., 1000):')
      READ(*,*)NCC 
 
      DCC=(CMAX-CMIN)/NCC 
      WRITE(*,1660)DCC 
 1660 FORMAT( ' Increment of Phase Velocity (DCC: m/s) =',F10.5) 
c 
      ISEC=0 
      NZINT=1 
      IF((NOM-IOM).ne.0)return 
      IF(IOM.EQ.0)RETURN 
c 
c      WRITE(*,1680) 
c 1680 FORMAT( ' Output Secular Functions? (yes=1)                  :')
c      READ(*,*)ISEC 
      ISEC=1 
      IF(ISEC.EQ.1)THEN 
        WRITE(*,*) 
        WRITE(*,*)' You can plot the secular functions using' 
        WRITE(*,*)' "secl.m and secr.m" for M-files of "Matlab"' 
      END IF 
C 
   25 NZINT=1 
      IF((NOM-IOM).ne.0)return 
c      WRITE(*,1800) 
c      READ(*,*)IVCT 
      IVCT=1 
      IF(IVCT.NE.1)RETURN 
 
      WRITE(*,1810) 
      READ(*,*)ZINT 
      IF(ZINT.LT.THK(1))THEN 
        WRITE(*,*) 
        WRITE(*,*)' Depth is smaller than the thichness of 1st layer !' 
      END IF 
      WRITE(*,1820) 
      READ(*,*)NZINT 
      WRITE(*,*) 
      IF(IVCT.EQ.1)THEN 
        WRITE(*,*) 
        WRITE(*,*)' You can plot the eigen vectors using "vld.m, vrd.m' 
        WRITE(*,*)' vls.m, and vrs.m" for "M-file" of "Matlab"' 
      END IF 
c 1800 FORMAT(/' Output Eigen Vectors? (yes=1)              :')
 1810 FORMAT(/' Enter Depth for plotting the Eigen Vectors (m):')
 1820 FORMAT(/' Enter Partition Number along the Depth (ex., 200):')
C 
      RETURN 
      END 
C 
C 
C *** PHASE VELOCITY OF THE RAYLEIGH WAVE IN THE HOMOGENEOUS HALF-SPACE *** 
C 
      FUNCTION FRHF(VP,VS) 
      IMPLICIT REAL*8(A-H,O-Z) 
C 
      DZ=0.0001d0 
C 
      SP=VS/VP 
      IF(SP.GE.0.0D0) Z= 0.91263D0 
      IF(SP.GE.0.2D0) Z= 0.90767D0 
      IF(SP.GE.0.4D0) Z= 0.88899D0 
      IF(SP.GE.0.6D0) Z= 0.83575D0 
      IF(SP.GE.0.7D0) Z= 0.77062D0 
      IF(SP.GE.0.8D0) Z= 0.63338D0 
C 
      F0=FRHF1(SP,Z) 
  100 Z=Z-DZ 
      F1=FRHF1(SP,Z) 
      AF=F0/F1 
      F00=F0 
      F0 =F1 
      IF(AF.GE.0.D0)GOTO 100 
      FZ=Z+DZ*F0/(F0-F00) 
C 
      FRHF=VS*DSQRT(FZ) 
C 
      RETURN 
      END 
C 
C 
      FUNCTION FRHF1(SP,Z) 
      IMPLICIT REAL*8(A-H,O-Z) 
C 
      Z2=Z*Z 
      Z3=Z2*Z 
C 
      F1=16.D0*(Z-1.D0)*SP*SP 
      FRHF1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0 
C 
      RETURN 
      END 
C 
c 
C *** COMPUTATION OF PHASE VELOCITY BY THE BI-SECTION METHOD *** 
C 
      SUBROUTINE PHASE(MOM) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 GRU0L,GTDL,GRDL,GTDR,GRDR,SCL,SCR,AA(2,2) 
     *          ,SC1,SC2,SC3,CMS,CMP,ER(4,4),DDS,DDP 
      COMMON /FREQ/DOM,IOM,NOM  
     *  /MINMAX/VSMIN,VSMAX,CCMIN,CCMAX,CMIN,CMAX,TOL,NITER,NCC
     *   /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
     *   /GRTLV/GRU0L,GTDL(20),GRDL(20) 
     *   /GRTRY/GTDR(2,2,20),GRDR(2,2,20) 
     *   /CCMD/CCL(0:16401,0:998),CCR(0:16401,0:998)
     *         ,MDL(0:16401),MDR(0:16401)
     *   /OUTPUT/ISEC,IVCT   
C 
C *** CHECK OF DIMENSION *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM PHASE: NL.GT.20 !!!!!!!!!!!!!!!' 
         WRITE(10,*)' FROM PHASE: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
      NC=16401
      IF(MOM.GT.NC)THEN 
         WRITE(*,*)' FROM PHASE: MOM.GT.16401 !!!!!!!!!!!!!!!'
         WRITE(10,*)' FROM PHASE: MOM.GT.16401 !!!!!!!!!!!!!!!'
         STOP 
      END IF 
C 
C *** PHASE VELOCITY *** 
C 
      OPEN(20,FILE='om_dat',STATUS='unknown') 
C 
      DCC=(CMAX-CMIN)/NCC 
      write(*,*) 
      write(10,*) 
      DO 100 IO=IOM,NOM 
C 
        OM=DOM*DFLOAT(IO) 
        write(20,1100)om 
        write(*,*)
        write(*,1110)om,om/6.283185307d0 
        write(10,1110)om,om/6.283185307d0 
 1100   FORMAT(e15.7) 
 1110   FORMAT(' omega and freq.=',f8.3,f8.4) 
c 
        IF(IO.EQ.0)THEN 
          MDL(IO)=0 
          CCL(IO,0)=VSMAX 
          GOTO 20 
        END IF 
C 
        IF(ISEC.EQ.1)THEN 
          OPEN(11,FILE='sclc_dat',STATUS='unknown') 
          OPEN(12,FILE='sclr_dat',STATUS='unknown') 
          OPEN(13,FILE='scli_dat',STATUS='unknown') 
          OPEN(14,FILE='scrc_dat',STATUS='unknown') 
          OPEN(15,FILE='scrr_dat',STATUS='unknown') 
          OPEN(16,FILE='scri_dat',STATUS='unknown')
          open(17,file='sec-Lov.csv',STATUS='unknown')
          open(18,file='sec-Ray.csv',STATUS='unknown')
          write(17,*)'Phase Velocity (m/s), Real Part, Imag Part
     *,Absolute Value'
          write(18,*)'Phase Velocity (m/s), Real Part, Imag Part
     *,Absolute Value'
        END IF
C
        IF(NL.EQ.1)GOTO 20 
C 
C    ** SEARCH OF LOVE POLES BY THE BI-SECTION METHOD ** 
C 
        MD=0 
c
        C1=VSMIN-DCC
        WN=OM/C1 
        CALL GNRTLV(OM,WN) 
        SC1=1.D0-GRU0L*GRDL(1) 
        AS1=CDABS(SC1) 
C 
        IF(ISEC.EQ.1)then 
          WRITE(11,1100)C1 
          WRITE(12,1100)dreal(sc1) 
          WRITE(13,1100)dimag(sc1)
          write(17,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1)
        end if 
 1101   FORMAT(e15.7,3(',',e15.7))
C
        C2=C1+DCC 
        WN=OM/C2 
        CALL GNRTLV(OM,WN) 
        SC2=1.D0-GRU0L*GRDL(1) 
        AS2=CDABS(SC2) 
C 
   10   C3=C2+DCC 
        if(C3.gt.(CMAX+DCC))goto 15
        WN=OM/C3 
        CALL GNRTLV(OM,WN) 
        SC3=1.D0-GRU0L*GRDL(1) 
        AS3=CDABS(SC3) 
C  
        C1N=C1 
        C2N=C2 
        C3N=C3 
        AS2N=AS2  
c
c        write(2,9999)c1n,as1,c2n,as2,c3n,as3
c 9999   format(3(f10.4,f10.6))

        IF((AS2.LT.AS1).AND.(AS2.LT.AS3))THEN 
C 
c         write(2,*)0,c1n,as1
c         write(2,*)0,c2n,as2
c         write(2,*)0,c3n,as3
c
          do 110 iter=1,NITER
c
            C12=(C1N+C2N)/2.D0 
            WN=OM/C12 
            CALL GNRTLV(OM,WN) 
            SCL=1.D0-GRU0L*GRDL(1) 
            AS12=CDABS(SCL) 
c         write(2,*)iter,c12,as12
c
            if(as12.eq.as2n)goto 12
c
            C23=(C2N+C3N)/2.D0 
            WN=OM/C23 

            CALL GNRTLV(OM,WN) 
            SCL=1.D0-GRU0L*GRDL(1) 
            AS23=CDABS(SCL) 
c         write(2,*)iter,c23,as23
c         write(2,*)100,c12,as12
c         write(2,*)200,c2n,as2n
c         write(2,*)300,c23,as23
c
            IF(AS12.LT.AS2N)THEN 
              C3N =C2N 
              C2N =C12 
              AS2N=AS12 
            ELSE 
              IF(AS23.LT.AS2N)THEN 
                C1N =C2N 
                C2N =C23 
                AS2N=AS23 
              ELSE 
                C1N =C12 
                C3N =C23 
              END IF 
            END IF 
c
            IF(AS2N.LT.TOL)THEN 
              CCL(IO,MD)=C2N 
              MD=MD+1 
              goto 14
            END IF 
C 
  110     CONTINUE 
C 
   12     write(*,*)'Love Wave: Did not converge!'
          write(10,*)'Love Wave: Did not converge!'
          write(*,1200)MD,iter,C2N,AS2N
          write(10,1200)MD,iter,C2N,AS2N
 1200     format(' Mode,Iter,Phase Vel,Amp=',i3,i5,f9.2,f9.3)
c
   14     C1=C2N 
          WN=OM/C1 
          CALL GNRTLV(OM,WN) 
          SC1=1.D0-GRU0L*GRDL(1) 
          AS1=CDABS(SC1) 
c
          IF(ISEC.EQ.1)THEN 
            WRITE(11,1100)C1 
            WRITE(12,1100)dreal(sc1) 
            WRITE(13,1100)dimag(sc1) 
            write(17,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1)
          END IF
C 
          C2=C1+DCC 
          WN=OM/C2 
          CALL GNRTLV(OM,WN) 
          SC2=1.D0-GRU0L*GRDL(1) 
          AS2=CDABS(SC2) 
C 
        ELSE
c
          IF(ISEC.EQ.1)THEN 
            WRITE(11,1100)C2 
            WRITE(12,1100)dreal(sc2) 
            WRITE(13,1100)dimag(sc2) 
            write(17,1101)C2,dreal(sc2),dimag(sc2),cdabs(sc2)
          END IF
C      
          C1=C2 
          SC1=SC2
          AS1=AS2 
C 
          C2=C3 
          SC2=SC3 
          AS2=AS3 
C
        END IF 

        goto 10 
C 
   15   MDL(IO)=MD-1 
        NC=998
        IF(MDL(IO).GT.NC)THEN 
           WRITE(*,*)' FROM PHASE: LOVE MODE.GT.998 !!!!!!!!!!!!!!'
           STOP 
        END IF 
C 
C *** SEARCH OF RAYLEIGH POLES BY THE BI-SECTION METHOD *** 
C 
   20   MD=0 
C 
        write(*,*)'1;R start,om',om

        IF(NL.EQ.1)THEN 
          CCR(IO,MD)=CCMIN 
          GO TO 100 
        END IF 
c 
        IF(IO.EQ.0)THEN 
          MDR(IO)=0 
          CCR(IO,0)=CCMAX 
          GOTO 100 
        END IF 
C 
        SM2=DNS(1)*VS(1)*VS(1)*DNS(1)*VS(1)*VS(1) 
C 
        C1=CMIN 
        IF(dabs(C1-VS(1)).lt.1.0d-7)C1=C1+C1/1000.D0 
        IF(dabs(C1-VP(1)).lt.1.0d-7)C1=C1+C1/1000.D0 
        WN=OM/C1 
        CALL GNRTRY(OM,WN) 
        CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
        DDS=CDEXP(-CMS*THK(1)) 
        DDP=CDEXP(-CMP*THK(1)) 
        AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1)+ER(3,4)*DDS*GRDR(2,1,1) 
        AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1)+ER(3,4)*DDS*GRDR(2,2,1) 
        AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1)+ER(4,4)*DDS*GRDR(2,1,1) 
        AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1)+ER(4,4)*DDS*GRDR(2,2,1) 
        SC1=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
        AS1=CDABS(SC1) 
        AMAX=AS1 
C 
        IF(ISEC.EQ.1)then 
          WRITE(14,1100)C1 
          WRITE(15,1100)dreal(sc1) 
          WRITE(16,1100)dimag(sc1) 
          write(18,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1)
        end if
C 
        C2=C1+DCC 
        IF(dabs(C2-VS(1)).lt.1.0d-7)C2=C2+C2/1000.D0 
        IF(dabs(C2-VP(1)).lt.1.0d-7)C2=C2+C2/1000.D0 
        WN=OM/C2 
        CALL GNRTRY(OM,WN) 
        CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
        DDS=CDEXP(-CMS*THK(1)) 
        DDP=CDEXP(-CMP*THK(1)) 
        AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1)+ER(3,4)*DDS*GRDR(2,1,1) 
        AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1)+ER(3,4)*DDS*GRDR(2,2,1) 
        AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1)+ER(4,4)*DDS*GRDR(2,1,1) 
        AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1)+ER(4,4)*DDS*GRDR(2,2,1) 
        SC2=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
        AS2=CDABS(SC2) 
        IF(AS2.GT.AMAX)AMAX=AS2 
C 
   30   C3=C2+DCC 
        if(C3.gt.(CMAX+DCC))goto 35
        IF(dabs(C3-VS(1)).lt.1.0d-7)C3=C3+C3/1000.D0 
        IF(dabs(C3-VP(1)).lt.1.0d-7)C3=C3+C3/1000.D0 
        WN=OM/C3 
        CALL GNRTRY(OM,WN) 
        CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
        DDS=CDEXP(-CMS*THK(1)) 
        DDP=CDEXP(-CMP*THK(1)) 
        AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) 
     *                 +ER(3,4)*DDS*GRDR(2,1,1) 
        AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) 
     *                 +ER(3,4)*DDS*GRDR(2,2,1) 
        AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) 
     *                 +ER(4,4)*DDS*GRDR(2,1,1) 
        AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) 
     *                 +ER(4,4)*DDS*GRDR(2,2,1) 
        SC3=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
        AS3=CDABS(SC3) 
        IF(AS3.GT.AMAX)AMAX=AS3 
C 
        C1N=C1 
        C2N=C2 
        C3N=C3  
        AS2N=AS2 
c
c        write(2,9999)c1n,as1,c2n,as2,c3n,as3
c 9999   format(3(f10.4,f10.6))
C  
        IF((AS2.LT.AS1).AND.(AS2.LT.AS3))THEN 
c
c         write(2,*)0,c1n,as1
c         write(2,*)0,c2n,as2
c         write(2,*)0,c3n,as3
c
          do 310 iter=1,NITER       
C 
            C12=(C1N+C2N)/2.D0 
            IF(dabs(C12-VS(1)).lt.1.0d-7)C12=C12+C12/1000.D0 
            IF(dabs(C12-VP(1)).lt.1.0d-7)C12=C12+C12/1000.D0 
            WN=OM/C12 
            CALL GNRTRY(OM,WN) 
            CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
            DDS=CDEXP(-CMS*THK(1)) 
            DDP=CDEXP(-CMP*THK(1)) 
            AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) 
     *                     +ER(3,4)*DDS*GRDR(2,1,1) 
            AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) 
     *                     +ER(3,4)*DDS*GRDR(2,2,1) 
            AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) 
     *                     +ER(4,4)*DDS*GRDR(2,1,1) 
            AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) 
     *                     +ER(4,4)*DDS*GRDR(2,2,1) 
            SCR=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
            AS12=CDABS(SCR) 
c         write(2,*)iter,c12,as12
c
            if(as12.eq.as2n)goto 32
C 
            C23=(C2N+C3N)/2.D0 
            IF(dabs(C23-VS(1)).lt.1.0d-7)C23=C23+C23/1000.D0 
            IF(dabs(C23-VP(1)).lt.1.0d-7)C23=C23+C23/1000.D0 
            WN=OM/C23 
            CALL GNRTRY(OM,WN) 
            CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
            DDS=CDEXP(-CMS*THK(1)) 
            DDP=CDEXP(-CMP*THK(1)) 
            AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) 
     *                     +ER(3,4)*DDS*GRDR(2,1,1) 
            AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) 
     *                     +ER(3,4)*DDS*GRDR(2,2,1) 
            AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) 
     *                     +ER(4,4)*DDS*GRDR(2,1,1) 
            AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) 
     *                     +ER(4,4)*DDS*GRDR(2,2,1) 
            SCR=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
            AS23=CDABS(SCR) 
c
c         write(2,*)iter,c23,as23
c         write(2,*)100,c12,as12
c         write(2,*)200,c2n,as2n
c         write(2,*)300,c23,as23
C 
            IF(AS12.LT.AS2N)THEN 
              C3N =C2N 
              C2N =C12 
              AS2N=AS12 
            ELSE 
              IF(AS23.LT.AS2N)THEN 
                C1N =C2N 
                C2N =C23 
                AS2N=AS23 
              ELSE 
                C1N =C12 
                C3N =C23 
              END IF 
            END IF 
C 
            IF(AS2N.LT.TOL)THEN 
              CCR(IO,MD)=C2N 
              MD=MD+1 
              GOTO 34 
            END IF 
C 
  310     continue 
C 
   32     write(*,*)'Rayleigh Wave: Did not converge!'
          write(10,*)'Rayleigh Wave: Did not converge!'
          write(*,3200)MD,iter,C2N,AS2N
          write(10,3200)MD,iter,C2N,AS2N
 3200     format(' Mode,Iter,Phase Vel,Amp=',i3,i5,f9.2,f9.3)
C 
   34     C1=C2N
          WN=OM/C1 
          IF(C1.EQ.VS(1))C1=C1+C1/1000.D0 
          IF(C1.EQ.VP(1))C1=C1+C1/1000.D0 
          CALL GNRTRY(OM,WN) 
          CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
          DDS=CDEXP(-CMS*THK(1)) 
          DDP=CDEXP(-CMP*THK(1)) 
          AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) 
     *                   +ER(3,4)*DDS*GRDR(2,1,1) 
          AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) 
     *                   +ER(3,4)*DDS*GRDR(2,2,1) 
          AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) 
     *                   +ER(4,4)*DDS*GRDR(2,1,1) 
          AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) 
     *                   +ER(4,4)*DDS*GRDR(2,2,1) 
          SC1=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
          AS1=CDABS(SC1) 

          IF(ISEC.EQ.1)THEN 
            WRITE(14,1100)C1 
            WRITE(15,1100)dreal(sc1) 
            WRITE(16,1100)dimag(sc1) 
            write(18,1101)C1,dreal(sc1),dimag(sc1),cdabs(sc1)
          END IF
C 
          C2=C1+DCC
          WN=OM/C2 
          IF(C2.EQ.VS(1))C2=C2+C2/1000.D0 
          IF(C2.EQ.VP(1))C2=C2+C2/1000.D0 
          CALL GNRTRY(OM,WN) 
          CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
          DDS=CDEXP(-CMS*THK(1)) 
          DDP=CDEXP(-CMP*THK(1)) 
          AA(1,1)=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) 
     *                   +ER(3,4)*DDS*GRDR(2,1,1) 
          AA(1,2)=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) 
     *                   +ER(3,4)*DDS*GRDR(2,2,1) 
          AA(2,1)=ER(4,1)+ER(4,3)*DDP*GRDR(1,1,1) 
     *                   +ER(4,4)*DDS*GRDR(2,1,1) 
          AA(2,2)=ER(4,2)+ER(4,3)*DDP*GRDR(1,2,1) 
     *                   +ER(4,4)*DDS*GRDR(2,2,1) 
          SC2=(AA(1,1)*AA(2,2)-AA(1,2)*AA(2,1))/(CMS*CMP*SM2) 
          AS2=CDABS(SC2) 
c
        ELSE 
C 
          IF(ISEC.EQ.1)then 
            WRITE(14,1100)C2 
            WRITE(15,1100)dreal(sc2) 
            WRITE(16,1100)dimag(sc2)
            write(18,1101)C2,dreal(sc2),dimag(sc2),cdabs(sc2)
          end if
C 
          C1=C2 
          SC1=SC2 
          AS1=AS2 
C 
          C2=C3 
          SC2=SC3 
          AS2=AS3 
C 
        END IF 
C      
        GOTO 30 
C 
   35   MDR(IO)=MD-1 
        NC=998
        IF(MDR(IO).GT.NC)THEN 
           WRITE(*,*)' FROM PHASE: RAYLEIGH MODE.GT.998 !!!!!!!!!!!!!!!'
           STOP 
        END IF 
C 
      write(*,*) 
     * '  The highest modes for L- and R-waves',mdl(io),mdr(io) 
C 
  100 CONTINUE 
C 
      CLOSE(11) 
      CLOSE(12) 
      CLOSE(13) 
      CLOSE(14) 
      CLOSE(15) 
      CLOSE(16) 
      CLOSE(20) 
c
      OPEN(22,FILE='mdl_dat',STATUS='unknown')
      OPEN(23,FILE='mdr_dat',STATUS='unknown')
        DO 300 IO=IOM,NOM
          WRITE(22,1000)MDL(IO)
          WRITE(23,1000)MDR(IO)
  300   CONTINUE
 1000 FORMAT(I5)
      CLOSE(22)
      CLOSE(23)
c
      OPEN(25,FILE='phsl1_dat',STATUS='unknown')
      OPEN(26,FILE='phsr1_dat',STATUS='unknown')
        DO 210 IO=IOM,NOM
          DO 200 MD=0,MDL(IO)
            WRITE(25,1100)CCL(IO,MD)
  200     CONTINUE
          DO 205 MD=0,MDR(IO)
            WRITE(26,1100)CCR(IO,MD)
  205   CONTINUE
  210   CONTINUE
      CLOSE(25)
      CLOSE(26)
c
      RETURN
      END 
C 
C 
C *** COMPUTATION OF EIGEN-VECTORS, ENERGY INTEGRALS AND GROUP VELOCITIES *** 
C 
      SUBROUTINE ENGUU(MOM) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 GRU0L,GTDL,GRDL,GTDR,GRDR 
     *          ,EL(2,2),CDS,CUS,DDS,DUS,CMS,CDPU,VR20 
     *          ,ER(4,4),CDP,CUP,DDP,DUP,CMP,DR1,DR2,DR14,DR23 
     *          ,A1,A2,A3,A4,A5,A6,B1,B2,B3,B4,B5,B6 
      dimension res(0:998)
      COMMON /FREQ/DOM,IOM,NOM
     *   /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
     *   /OUTPUT/ISEC,IVCT  /NINTG/ZINT,NZINT 
     *   /MINMAX/VSMIN,VSMAX,CCMIN,CCMAX,CMIN,CMAX,TOL,NITER,NCC
     *   /GRTLV/GRU0L,GTDL(20),GRDL(20)
     *   /GRTRY/GTDR(2,2,20),GRDR(2,2,20) 
     *   /CCMD/CCL(0:16401,0:998),CCR(0:16401,0:998)
     *        ,MDL(0:16401),MDR(0:16401)
     *   /VCTL/VL1(16401,0:998),VL2(16401,0:998)
     *   /VCLR/VR1(16401,0:998),VR2(16401,0:998)
     *        ,VR3(16401,0:998),VR4(16401,0:998)
     *   /UEGL/UUL(0:16401,0:998),E1L(0:16401,0:998)
     *   /UEGR/UUR(0:16401,0:998),E1R(0:16401,0:998)
C 
C *** CHECK OF DIMENSION *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM ENGUU: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
      NC=16401
      IF((NZINT+2).GT.NC)THEN 
         WRITE(*,*)' FROM ENGUU: NZINT+2.GT.16401 !!!!!!!!!!!!!!!'
         STOP 
      END IF 
C 
      NC=16401
      IF(MOM.GT.NC)THEN 
         WRITE(*,*)' FROM ENGUU: MOM.GT.16401 !!!!!!!!!!!!!!!'
         STOP 
      END IF 
C 
C *** EIGEN VECTORS FOR LOVE WAVE *** 
C 
      IF(IVCT.EQ.1)THEN 
        OPEN(9,FILE='vzn_dat',STATUS='unknown')
          write(9,1000)zint
          write(9,1100)nzint
          write(9,1100)nl
          do 90 l=1,nl
            write(9,1000)thk(l)
 90       continue
        close(9)
        OPEN(10,FILE='vl1_dat',STATUS='unknown')
        OPEN(11,FILE='vl2_dat',STATUS='unknown')
        OPEN(20,FILE='vr1_dat',STATUS='unknown')
        OPEN(21,FILE='vr2_dat',STATUS='unknown')
        OPEN(22,FILE='vr3_dat',STATUS='unknown')
        OPEN(23,FILE='vr4_dat',STATUS='unknown')
        OPEN(24,FILE='Ldis.csv',STATUS='unknown')
        OPEN(25,FILE='Lstr.csv',STATUS='unknown')
        OPEN(26,FILE='Rdis.csv',STATUS='unknown')
        OPEN(27,FILE='Rstr.csv',STATUS='unknown')
        WRITE(24,2024)(MD,MD=MDL(NOM),0,-1)
        WRITE(25,2024)(MD,MD=MDL(NOM),0,-1)
        WRITE(26,2025)(MD,MD,MD=MDR(NOM),0,-1)
        WRITE(27,2025)(MD,MD,MD=MDR(NOM),0,-1)
      END IF
 2024 FORMAT('depth (km)',400(',',I3,' mode'))
 2025 FORMAT('depth (km)',400(',',I3,' mode (h),',I3,' mode (v)'))
 1000 FORMAT(e15.7)
 1100 format(i5) 
C 
      DZ=ZINT/NZINT 
      DO 100 IO=IOM,NOM 
C 
        OM=DOM*DFLOAT(IO) 
        OM2=OM*OM 
C 
        IF(NL.EQ.1)GOTO 20 
c 
        IF(IO.EQ.0)THEN 
          UUL(IO,0)=VSMAX 
          GOTO 20 
        END IF 
C 
        DO 110 MD=0,MDL(IO) 
C 
          WN=OM/CCL(IO,MD) 
          WN2=WN*WN 
          CALL GNRTLV(OM,WN) 
C 
          IF(IVCT.NE.1)GOTO 10 
C 
          LAY=1 
          UDEP=0.D0 
          BDEP=THK(1) 
C 
          VS2=VS(1)*VS(1) 
          WS2=OM2/VS2 
          CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
          IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
          SM=DNS(1)*VS2 
          EL(2,1)=-SM*CMS 
          EL(2,2)= SM*CMS 
C 
          CDS=(0.5D0,0.D0) 
          CUS=GRDL(1)*0.5D0 
C 
          DO 120 IZ=1,NZINT+1 
C 
            Z=DZ*(IZ-1) 
            IF(Z.GT.BDEP)THEN 
C 
              LAY=LAY+1 
              UDEP=BDEP 
              IF(LAY.EQ.NL)THEN 
                BDEP=BDEP+ZINT 
              ELSE 
                BDEP=BDEP+THK(LAY) 
              END IF 
              VS2=VS(LAY)*VS(LAY) 
              WS2=OM2/VS2 
              CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
              IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
              SM=DNS(LAY)*VS2 
              EL(2,1)=-SM*CMS 
              EL(2,2)= SM*CMS 
C 
              CDS=GTDL(LAY-1)*CDS 
              CUS=GRDL(LAY)*CDS 
C 
            END IF 
C 
            DDS=CDEXP(-CMS*(Z-UDEP)) 
            IF(LAY.EQ.NL)THEN 
              VL1(IZ,MD)=dreal(DDS*CDS)
              VL2(IZ,MD)=dreal(EL(2,1)*DDS*CDS)
            ELSE 
              DUS=CDEXP(-CMS*(BDEP-Z)) 
              VL1(IZ,MD)=dreal(DDS*CDS+DUS*CUS)
              VL2(IZ,MD)=dreal(EL(2,1)*DDS*CDS+EL(2,2)*DUS*CUS)
            END IF 
C      
  120     CONTINUE 
C 
          DO 140 I=1,NZINT+1 
          WRITE(10,1000)VL1(I,MD) 
          WRITE(11,1000)VL2(I,MD) 
  140     CONTINUE 
C 
C *** ENERGY INTEGRALS AND GROUP VELOCITIES FOR LOVE WAVE *** 
C 
   10     E1L(IO,MD)=0.D0 
          E2L       =0.D0 
C 
          DO 160 LAY=1,NL 
C 
            IF(LAY.EQ.1)THEN 
              CDS=(0.5D0,0.D0) 
              CUS=GRDL(1)*0.5D0 
            ELSE 
              CDS=GTDL(LAY-1)*CDS 
              CUS=GRDL(LAY)*CDS 
            END IF 
 
            VS2=VS(LAY)*VS(LAY) 
            WS2=OM2/VS2 
            CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
            IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
C 
            IF(LAY.EQ.NL)THEN 
              DR1=CDS*CDS/(2.D0*CMS) 
            ELSE 
              DDS=CDEXP(-2.D0*CMS*THK(LAY)) 
              DR1=(CDS*CDS+CUS*CUS)*(1.D0-DDS)/(2.D0*CMS) 
              DDS=CDEXP(-CMS*THK(LAY)) 
              DR1=DR1+2.D0*CDS*CUS*THK(LAY)*DDS 
            END IF 
C      
            SM=DNS(LAY)*VS2 
            E1L(IO,MD)=dreal(E1L(IO,MD)+DNS(LAY)*DR1)
            E2L       =dreal(E2L       +SM*DR1)
C 
  160     CONTINUE 
C 
          E1L(IO,MD)=E1L(IO,MD)/2.D0 
          E2L       =E2L       /2.D0 
          UUL(IO,MD)=E2L/(E1L(IO,MD)*CCL(IO,MD)) 
C          
  110   CONTINUE 
c 
        do MD=0,MDL(IO) 
          VDIS=0.001d0
          VSTR=0.001d0
          do IZ=1,NZINT+1 
            if(dabs(VL1(IZ,MD)).gt.VDIS)VDIS=dabs(VL1(IZ,MD)) 
            if(dabs(VL2(IZ,MD)).gt.VSTR)VSTR=dabs(VL2(IZ,MD)) 
          end do 
          do IZ=1,NZINT+1 
            VL1(IZ,MD)=VL1(IZ,MD)/VDIS 
            VL2(IZ,MD)=VL2(IZ,MD)/VSTR
          end do
        end do 
c 
      IF(IVCT.EQ.1)THEN
        do IZ=1,NZINT+1
          Z=DZ*(IZ-1) 
          write(24,2500)z/1000.d0,(VL1(IZ,MD)/2.d0+MD,MD=MDL(IO),0,-1)
          write(25,2500)z/1000.d0,(VL2(IZ,MD)/2.d0+MD,MD=MDL(IO),0,-1)
        end do
      end if
C 
C *** EIGEN VECTORS FOR RAYLEIGH WAVE *** 
C 
   20     IF(IO.EQ.0)THEN 
            UUR(IO,0)=CCMAX 
            GOTO 100 
          END IF 
c 
          DO 210 MD=0,MDR(IO) 
C 
          WN=OM/CCR(IO,MD) 
          CALL GNRTRY(OM,WN) 
C 
          IF(IVCT.NE.1)GOTO 30 
C 
          LAY=1 
          UDEP=0.D0 
          CALL MTXRY(OM,WN,1,CMS,CMP,ER) 
C 
C     *   HOMOGENEOUS HALF-SPACE SOLUTION * 
C 
          IF(NL.EQ.1)THEN 
C 
            B1=-ER(3,1)/ER(3,2) 
            B2= ER(2,1)+ER(2,2)*B1 
            CDP= 1.D0/B2 
            CDS= B1*CDP 
C 
            DO 215 IZ=1,NZINT+1 
C 
              Z=DZ*(IZ-1) 
              DDS=CDEXP(-CMS*Z) 
              DDP=CDEXP(-CMP*Z) 
              A1=DDP*CDP 
              A2=DDS*CDS 
              VR1(IZ,MD)=dreal(ER(1,1)*A1+ER(1,2)*A2)
              VR2(IZ,MD)=dreal(ER(2,1)*A1+ER(2,2)*A2)
              VR3(IZ,MD)=dreal(ER(3,1)*A1+ER(3,2)*A2)
              VR4(IZ,MD)=dreal(ER(4,1)*A1+ER(4,2)*A2)
C      
  215       CONTINUE 
C 
C     *   LAYERED HALF-SPACE SOLUTION * 
C 
          ELSE 
C 
            DDS=CDEXP(-CMS*THK(1)) 
            DDP=CDEXP(-CMP*THK(1)) 
            A1=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1)+ER(3,4)*DDS*GRDR(2,1,1) 
            A2=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1)+ER(3,4)*DDS*GRDR(2,2,1) 
            B1=-A1/A2 
C 
            CUP=GRDR(1,1,1)+GRDR(1,2,1)*B1 
            CUS=GRDR(2,1,1)+GRDR(2,2,1)*B1 
            A3=DDP*CUP 
            A4=DDS*CUS 
            VR20=ER(2,1)+ER(2,2)*B1+ER(2,3)*A3+ER(2,4)*A4 
C 
            CDP=1.D0/VR20 
            CDS=B1*CDP 
            CUP=GRDR(1,1,1)*CDP+GRDR(1,2,1)*CDS 
            CUS=GRDR(2,1,1)*CDP+GRDR(2,2,1)*CDS 
C 
            BDEP=THK(1) 
            DO 220 IZ=1,NZINT+1 
C 
              Z=DZ*(IZ-1) 
              IF(Z.GT.BDEP)THEN 
C 
                LAY=LAY+1 
                UDEP=BDEP 
                IF(LAY.EQ.NL)THEN 
                  BDEP=BDEP+ZINT 
                ELSE 
                  BDEP=BDEP+THK(LAY) 
                END IF 
                CALL MTXRY(OM,WN,LAY,CMS,CMP,ER) 
C 
                CDPU=CDP 
                CDP=GTDR(1,1,LAY-1)*CDPU+GTDR(1,2,LAY-1)*CDS 
                CDS=GTDR(2,1,LAY-1)*CDPU+GTDR(2,2,LAY-1)*CDS 
                CUP=GRDR(1,1,LAY)*CDP+GRDR(1,2,LAY)*CDS 
                CUS=GRDR(2,1,LAY)*CDP+GRDR(2,2,LAY)*CDS 
C 
              END IF 
C 
              DDS=CDEXP(-CMS*(Z-UDEP)) 
              DDP=CDEXP(-CMP*(Z-UDEP)) 
              A1=DDP*CDP 
              A2=DDS*CDS 
              IF(LAY.EQ.NL)THEN 
                VR1(IZ,MD)=dreal(ER(1,1)*A1+ER(1,2)*A2)
                VR2(IZ,MD)=dreal(ER(2,1)*A1+ER(2,2)*A2)
                VR3(IZ,MD)=dreal(ER(3,1)*A1+ER(3,2)*A2)
                VR4(IZ,MD)=dreal(ER(4,1)*A1+ER(4,2)*A2)
              ELSE 
                DUS=CDEXP(-CMS*(BDEP-Z)) 
                DUP=CDEXP(-CMP*(BDEP-Z)) 
                A3=DUP*CUP 
                A4=DUS*CUS 
        VR1(IZ,MD)=dreal(ER(1,1)*A1+ER(1,2)*A2+ER(1,3)*A3+ER(1,4)*A4)
        VR2(IZ,MD)=dreal(ER(2,1)*A1+ER(2,2)*A2+ER(2,3)*A3+ER(2,4)*A4)
        VR3(IZ,MD)=dreal(ER(3,1)*A1+ER(3,2)*A2+ER(3,3)*A3+ER(3,4)*A4)
        VR4(IZ,MD)=dreal(ER(4,1)*A1+ER(4,2)*A2+ER(4,3)*A3+ER(4,4)*A4)
              END IF 
C      
  220       CONTINUE 
C 
          END IF 
C 
          DO 240 I=1,NZINT+1 
          WRITE(20,1000)VR1(I,MD) 
          WRITE(21,1000)VR2(I,MD) 
          WRITE(22,1000)VR3(I,MD) 
          WRITE(23,1000)VR4(I,MD) 
  240     CONTINUE 
C 
C *** ENERGY INTEGRAL AND GROUP VELOCITY FOR RAYLEIGH WAVE *** 
C 
   30     E1R(IO,MD)=0.D0 
          E2R       =0.D0 
          E3R       =0.D0 
C 
          DO 260 LAY=1,NL 
C 
            CALL MTXRY(OM,WN,LAY,CMS,CMP,ER) 
C 
            IF(NL.EQ.1)THEN 
              B1=-ER(3,1)/ER(3,2) 
              B2= ER(2,1)+ER(2,2)*B1 
              CDP= 1.D0/B2 
              CDS= B1*CDP 
            ELSE 
              IF(LAY.EQ.1)THEN 
                DDS=CDEXP(-CMS*THK(1)) 
                DDP=CDEXP(-CMP*THK(1)) 
                A1=ER(3,1)+ER(3,3)*DDP*GRDR(1,1,1) 
     *                     +ER(3,4)*DDS*GRDR(2,1,1) 
                A2=ER(3,2)+ER(3,3)*DDP*GRDR(1,2,1) 
     *                     +ER(3,4)*DDS*GRDR(2,2,1) 
                B1=-A1/A2 
                CUP=GRDR(1,1,1)+GRDR(1,2,1)*B1 
                CUS=GRDR(2,1,1)+GRDR(2,2,1)*B1 
                A3=DDP*CUP 
                A4=DDS*CUS 
                VR20=ER(2,1)+ER(2,2)*B1+ER(2,3)*A3+ER(2,4)*A4 
                CDP=1.D0/VR20 
                CDS=B1*CDP 
              ELSE 
                CDPU=CDP 
                CDP=GTDR(1,1,LAY-1)*CDPU+GTDR(1,2,LAY-1)*CDS 
                CDS=GTDR(2,1,LAY-1)*CDPU+GTDR(2,2,LAY-1)*CDS 
              END IF 
              CUP=GRDR(1,1,LAY)*CDP+GRDR(1,2,LAY)*CDS 
              CUS=GRDR(2,1,LAY)*CDP+GRDR(2,2,LAY)*CDS 
            END IF 
C 
            IF(LAY.EQ.NL)THEN 
              A1=0.5D0/CMP 
              A2=0.5D0/CMS 
              A3=1.D0/(CMP+CMS) 
C 
              B1=ER(1,1)*ER(1,1)*CDP*CDP 
              B2=ER(1,2)*ER(1,2)*CDS*CDS 
              B3=2.D0*ER(1,1)*ER(1,2)*CDP*CDS 
              DR1=B1*A1+B2*A2+B3*A3 
              B1=ER(2,1)*ER(2,1)*CDP*CDP 
              B2=ER(2,2)*ER(2,2)*CDS*CDS 
              B3=2.D0*ER(2,1)*ER(2,2)*CDP*CDS 
              DR2=B1*A1+B2*A2+B3*A3 
              B1=ER(1,1)*ER(4,1)*CDP*CDP 
              B2=ER(1,2)*ER(4,2)*CDS*CDS 
              B3=(ER(1,1)*ER(4,2)+ER(1,2)*ER(4,1))*CDP*CDS 
              DR14=B1*A1+B2*A2+B3*A3 
              B1=ER(2,1)*ER(3,1)*CDP*CDP 
              B2=ER(2,2)*ER(3,2)*CDS*CDS 
              B3=(ER(2,1)*ER(3,2)+ER(2,2)*ER(3,1))*CDP*CDS 
              DR23=B1*A1+B2*A2+B3*A3 
            ELSE 
              DDP=CDEXP(-2.D0*CMP*THK(LAY)) 
              A1=0.5D0*(1.D0-DDP)/CMP               
              DDS=CDEXP(-2.D0*CMS*THK(LAY)) 
              A2=0.5D0*(1.D0-DDS)/CMS               
              DDS=CDEXP(-(CMP+CMS)*THK(LAY)) 
              A3=(1.D0-DDS)/(CMP+CMS)               
              DDP=CDEXP(-CMP*THK(LAY)) 
              A4=THK(LAY)*DDP               
              DDS=CDEXP(-CMS*THK(LAY)) 
              A5=(DDS-DDP)/(CMP-CMS) 
              A6=THK(LAY)*DDS               
C 
              B1=ER(1,1)*ER(1,1)*CDP*CDP+ER(1,3)*ER(1,3)*CUP*CUP 
              B2=ER(1,2)*ER(1,2)*CDS*CDS+ER(1,4)*ER(1,4)*CUS*CUS 
              B3=ER(1,1)*ER(1,2)*CDP*CDS+ER(1,3)*ER(1,4)*CUP*CUS 
              B4=ER(1,1)*ER(1,3)*CDP*CUP 
              B5=ER(1,1)*ER(1,4)*CDP*CUS+ER(1,2)*ER(1,3)*CUP*CDS 
              B6=ER(1,2)*ER(1,4)*CDS*CUS 
              DR1=B1*A1+B2*A2+2.D0*(B3*A3+B4*A4+B5*A5+B6*A6) 
C 
              B1=ER(2,1)*ER(2,1)*CDP*CDP+ER(2,3)*ER(2,3)*CUP*CUP 
              B2=ER(2,2)*ER(2,2)*CDS*CDS+ER(2,4)*ER(2,4)*CUS*CUS 
              B3=ER(2,1)*ER(2,2)*CDP*CDS+ER(2,3)*ER(2,4)*CUP*CUS 
              B4=ER(2,1)*ER(2,3)*CDP*CUP 
              B5=ER(2,1)*ER(2,4)*CDP*CUS+ER(2,2)*ER(2,3)*CUP*CDS 
              B6=ER(2,2)*ER(2,4)*CDS*CUS 
              DR2=B1*A1+B2*A2+2.D0*(B3*A3+B4*A4+B5*A5+B6*A6) 
C 
              B1=ER(1,1)*ER(4,1)*CDP*CDP+ER(1,3)*ER(4,3)*CUP*CUP 
              B2=ER(1,2)*ER(4,2)*CDS*CDS+ER(1,4)*ER(4,4)*CUS*CUS 
              B3=(ER(1,1)*ER(4,2)+ER(1,2)*ER(4,1))*CDP*CDS 
     *          +(ER(1,3)*ER(4,4)+ER(1,4)*ER(4,3))*CUP*CUS 
              B4=(ER(1,1)*ER(4,3)+ER(1,3)*ER(4,1))*CDP*CUP 
              B5=(ER(1,1)*ER(4,4)+ER(1,4)*ER(4,1))*CDP*CUS 
     *          +(ER(1,2)*ER(4,3)+ER(1,3)*ER(4,2))*CUP*CDS 
              B6=(ER(1,2)*ER(4,4)+ER(1,4)*ER(4,2))*CDS*CUS 
              DR14=B1*A1+B2*A2+B3*A3+B4*A4+B5*A5+B6*A6 
C 
              B1=ER(2,1)*ER(3,1)*CDP*CDP+ER(2,3)*ER(3,3)*CUP*CUP 
              B2=ER(2,2)*ER(3,2)*CDS*CDS+ER(2,4)*ER(3,4)*CUS*CUS 
              B3=(ER(2,1)*ER(3,2)+ER(2,2)*ER(3,1))*CDP*CDS 
     *          +(ER(2,3)*ER(3,4)+ER(2,4)*ER(3,3))*CUP*CUS 
              B4=(ER(2,1)*ER(3,3)+ER(2,3)*ER(3,1))*CDP*CUP 
              B5=(ER(2,1)*ER(3,4)+ER(2,4)*ER(3,1))*CDP*CUS 
     *          +(ER(2,2)*ER(3,3)+ER(2,3)*ER(3,2))*CUP*CDS 
              B6=(ER(2,2)*ER(3,4)+ER(2,4)*ER(3,2))*CDS*CUS 
              DR23=B1*A1+B2*A2+B3*A3+B4*A4+B5*A5+B6*A6 
            END IF 
C 
            SM=DNS(LAY)*VS(LAY)*VS(LAY) 
            RAM=DNS(LAY)*VP(LAY)*VP(LAY)-2.D0*SM 
            A=RAM+2.D0*SM 
C      
            E1R(IO,MD)=dreal(E1R(IO,MD)+DNS(LAY)*(DR1+DR2))
            E2R =dreal(E2R +A*DR1+SM*DR2 )
            E3R =dreal(E3R +RAM*(DR14-WN*RAM*DR1)/A-SM*WN*DR2-DR23)
C 
  260     CONTINUE 
C 
          E1R(IO,MD)=E1R(IO,MD)/2.D0 
          E2R       =E2R       /2.D0 
          UUR(IO,MD)=(E2R+E3R/(2.D0*WN))/(E1R(IO,MD)*CCR(IO,MD)) 
C          
  210   CONTINUE 
c 
        do MD=0,MDR(IO) 
          VDIS=0.001d0
          VSTR=0.001d0
          do IZ=1,NZINT+1 
            if(dabs(VR1(IZ,MD)).gt.VDIS)VDIS=dabs(VR1(IZ,MD)) 
            if(dabs(VR2(IZ,MD)).gt.VDIS)VDIS=dabs(VR2(IZ,MD)) 
            if(dabs(VR3(IZ,MD)).gt.VSTR)VSTR=dabs(VR3(IZ,MD)) 
            if(dabs(VR4(IZ,MD)).gt.VSTR)VSTR=dabs(VR4(IZ,MD)) 
          end do 
          do IZ=1,NZINT+1 
            VR1(IZ,MD)=VR1(IZ,MD)/VDIS 
            VR2(IZ,MD)=VR2(IZ,MD)/VDIS 
            VR3(IZ,MD)=VR3(IZ,MD)/VSTR 
            VR4(IZ,MD)=VR4(IZ,MD)/VSTR 
          end do 
        end do 
c 
        IF(IVCT.EQ.1)THEN 
        do IZ=1,NZINT+1 
          Z=DZ*(IZ-1) 
          write(26,2500)z/1000.d0,(VR1(IZ,MD)/2.d0+MD
     *                           ,VR2(IZ,MD)/2.d0+MD,MD=MDR(IO),0,-1)
          write(27,2500)z/1000.d0,(VR3(IZ,MD)/2.d0+MD
     *                           ,VR4(IZ,MD)/2.d0+MD,MD=MDR(IO),0,-1)
         end do
        END IF 
C 
  100 CONTINUE 
C 
      IF(IVCT.EQ.1)THEN 
      CLOSE(10) 
      CLOSE(11) 
      CLOSE(20) 
      CLOSE(21) 
      CLOSE(22) 
      CLOSE(23) 
      CLOSE(24) 
      CLOSE(25) 
      CLOSE(26) 
      CLOSE(27) 
      END IF 
 2500 FORMAT(f8.4,500(',',e12.4))
c
c *** Output of Dispersion Curves and Medium Response from Aki and Richards,
c   * MedRes-L (A & K, p302, (7.106)) -> |Gyy|=1/(4*k*c*U*I1)=1/(4*om*U*I1)
c   * MedRes-R (A & K, p303, (7.110)) -> |Gzz|=1/(4*k*c*U*I1)=1/(4*om*U*I1)
c
      OPEN(44,FILE='Ldisper.csv',STATUS='unknown')
      OPEN(45,FILE='Rdisper.csv',STATUS='unknown')
      OPEN(46,FILE='Lmedres.csv',STATUS='unknown')
      OPEN(47,FILE='Rmedres.csv',STATUS='unknown')
      WRITE(44,2600)(MD,MD,MD=0,MDL(NOM))
      WRITE(45,2600)(MD,MD,MD=0,MDR(NOM))
      WRITE(46,2650)(MD,MD=0,MDL(NOM))
      WRITE(47,2650)(MD,MD=0,MDR(NOM))
      do IO=IOM,NOM
	om=DOM*DFLOAT(IO)
        Hz=DOM*DFLOAT(IO)/6.283185307d0
        WRITE(44,2620)Hz,(CCL(IO,MD),UUL(IO,MD),MD=0,MDL(NOM))
        WRITE(45,2620)Hz,(CCR(IO,MD),UUR(IO,MD),MD=0,MDR(NOM))
	do MD=0,MDL(NOM)
	    res(MD)=0.d0
	    if (IO.ne.0.and.UUL(IO,MD).gt.0.00001d0)then
	      res(MD)=0.25d0/(om*UUL(IO,MD)*E1L(IO,MD))
	    end if
	end do
	WRITE(46,2670)Hz,(res(MD),MD=0,MDL(NOM))
	do MD=0,MDR(NOM)
	    res(MD)=0.d0
	    if (IO.ne.0.and.UUR(IO,MD).gt.0.00001d0)then
	      res(MD)=0.25d0/(om*UUR(IO,MD)*E1R(IO,MD))
	    end if
	end do
	WRITE(47,2670)Hz,(res(MD),MD=0,MDR(NOM))
      end do
 2600 FORMAT('Freq (Hz)',500(',C',I3,',U',I3))
 2620 FORMAT(f8.4,500(',',e15.7))
 2650 FORMAT('Freq (Hz)',500(',M',I3))
 2670 FORMAT(f8.4,500(',',e15.7))
      CLOSE(44)
      CLOSE(45)
      CLOSE(46)
      CLOSE(47)
C
      RETURN
      END 
C 
C 
C *** SH GENERALIZED R/T MATIX *** 
C 
      SUBROUTINE GNRTLV(OM,WN) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 GRU0,GTD,GRD,TD,TU,RD,RU 
      COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
      COMMON /GRTLV/GRU0,GTD(20),GRD(20) 
     *    /MRTLV/TD(20),TU(20),RD(20),RU(0:20) 
C 
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM GRTLV: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** MODIFIED R/T MATRIX *** 
C 
      CALL MDRTLV(OM,WN) 
C 
C *** GENERALIZED R/T MATRIX *** 
C 
      GRU0=RU(0) 
      GTD(NL-1)=TD(NL-1) 
      GRD(NL-1)=RD(NL-1) 
C 
      DO 200 I=NL-2,1,-1 
        if(cdabs(1.D0-RU(I)*GRD(I+1)).lt.1.d-15)then
          write(10,*)'from GNRTLV:om,wn',om,wn
          write(10,*)' 1.D0-RU(I)*GRD(I+1) < 1.d-15 !!'
          write(10,*)1.D0-RU(I)*GRD(I+1)
          write(10,*)' Try a different number of Partitions !!'
        end if
        GTD(I)=TD(I)/(1.D0-RU(I)*GRD(I+1))
        GRD(I)=RD(I)+TU(I)*GRD(I+1)*GTD(I) 
  200 CONTINUE 
C 
      RETURN 
      END 
C 
C
C
C *** SH MODIFIED R/T MATIX *** 
C 
      SUBROUTINE MDRTLV(OM,WN) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 CMS,EA(2,2),EB(2,2),TD,TU,RD,RU,DDA,DUB,EAB 
      COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
      COMMON /MRTLV/TD(20),TU(20),RD(20),RU(0:20) 
C 
C *** CHECK FOR THE DIMENSION OF TD, TU, RD, AND RU *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM MDRTSH: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** MODIFIED R/T MATRIX *** 
C 
      WN2=WN*WN 
      OM2=OM*OM 
C 
      VS2=VS(1)*VS(1) 
      WS2=OM2/VS2 
      CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
C 
      SM=DNS(1)*VS2 
      EA(1,1)=(1.D0,0.D0) 
      EA(1,2)=(1.D0,0.D0) 
      EA(2,1)=-SM*CMS 
      EA(2,2)= SM*CMS 
C 
      DUB=CDEXP(-CMS*THK(1)) 
      DDA=DUB 
C 
      RU(0)=DUB 
C 
      DO 100 I=1, NL-2 
C 
        VS2=VS(I+1)*VS(I+1) 
        WS2=OM2/VS2 
        CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
C 
        SM=DNS(I+1)*VS2 
        EB(1,1)=(1.D0,0.D0) 
        EB(1,2)=(1.D0,0.D0) 
        EB(2,1)=-SM*CMS 
        EB(2,2)= SM*CMS 
C 
        DUB=CDEXP(-CMS*THK(I+1)) 
C 
        EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2)  
        TD(I)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB 
        RU(I)=(EB(1,2)*EA(2,2)-EB(2,2)*EA(1,2))*DUB/EAB 
        RD(I)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB 
        TU(I)=(EB(1,2)*EB(2,1)-EB(2,2)*EB(1,1))*DUB/EAB 
C 
        EA(1,1)=EB(1,1) 
        EA(1,2)=EB(1,2) 
        EA(2,1)=EB(2,1) 
        EA(2,2)=EB(2,2) 
C 
        DDA=DUB 
C 
  100 CONTINUE 
C         
      VS2=VS(NL)*VS(NL) 
      WS2=OM2/VS2 
      CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
      IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
C 
      SM=DNS(NL)*VS2 
      EB(1,1)=(1.D0,0.D0) 
      EB(1,2)=(1.D0,0.D0) 
      EB(2,1)=-SM*CMS 
      EB(2,2)= SM*CMS 
C 
      EAB=EB(2,1)*EA(1,2)-EB(1,1)*EA(2,2)  
C 
      TD(NL-1)=(EA(1,2)*EA(2,1)-EA(1,1)*EA(2,2))*DDA/EAB 
      RD(NL-1)=(EB(1,1)*EA(2,1)-EA(1,1)*EB(2,1))*DDA/EAB 
C 
      RETURN 
      END 
C 
C 
C 
C *** P-SV GENERALIZED R/T MATIX *** 
C 
      SUBROUTINE GNRTRY(OM,WN) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 GTD,GRD,TD,TU,RD,RU,E1(2,2),E2(2,2),DT 
      COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
      COMMON /MRTRY/TD(2,2,20),TU(2,2,20),RD(2,2,20),RU(2,2,20) 
     *  /GRTRY/GTD(2,2,20),GRD(2,2,20) 
C 
C *** CHECK FOR THE DIMENSION OF GTD, GTU, GRD, AND GRU *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM GRTPRY: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** MODIFIED R/T MATRIX *** 
C 
      IF(NL.EQ.1)RETURN 
      CALL MDRTRY(OM,WN) 
C 
C *** GENERALIZED R/T MATRIX *** 
C 
      GTD(1,1,NL-1)=TD(1,1,NL-1) 
      GTD(1,2,NL-1)=TD(1,2,NL-1) 
      GTD(2,1,NL-1)=TD(2,1,NL-1) 
      GTD(2,2,NL-1)=TD(2,2,NL-1) 
      GRD(1,1,NL-1)=RD(1,1,NL-1) 
      GRD(1,2,NL-1)=RD(1,2,NL-1) 
      GRD(2,1,NL-1)=RD(2,1,NL-1) 
      GRD(2,2,NL-1)=RD(2,2,NL-1) 
C 
      DO 200 I=NL-2,1,-1 
C 
        E1(1,1)=1.D0-RU(1,1,I)*GRD(1,1,I+1)-RU(1,2,I)*GRD(2,1,I+1) 
        E1(1,2)=    -RU(1,1,I)*GRD(1,2,I+1)-RU(1,2,I)*GRD(2,2,I+1) 
        E1(2,1)=    -RU(2,1,I)*GRD(1,1,I+1)-RU(2,2,I)*GRD(2,1,I+1) 
        E1(2,2)=1.D0-RU(2,1,I)*GRD(1,2,I+1)-RU(2,2,I)*GRD(2,2,I+1) 
        DT=E1(1,1)*E1(2,2)-E1(1,2)*E1(2,1) 
        E2(1,1)= E1(2,2)/DT 
        E2(1,2)=-E1(1,2)/DT 
        E2(2,1)=-E1(2,1)/DT 
        E2(2,2)= E1(1,1)/DT 
        GTD(1,1,I)=E2(1,1)*TD(1,1,I)+E2(1,2)*TD(2,1,I) 
        GTD(1,2,I)=E2(1,1)*TD(1,2,I)+E2(1,2)*TD(2,2,I) 
        GTD(2,1,I)=E2(2,1)*TD(1,1,I)+E2(2,2)*TD(2,1,I) 
        GTD(2,2,I)=E2(2,1)*TD(1,2,I)+E2(2,2)*TD(2,2,I) 
C 
        E1(1,1)=TU(1,1,I)*GRD(1,1,I+1)+TU(1,2,I)*GRD(2,1,I+1) 
        E1(1,2)=TU(1,1,I)*GRD(1,2,I+1)+TU(1,2,I)*GRD(2,2,I+1) 
        E1(2,1)=TU(2,1,I)*GRD(1,1,I+1)+TU(2,2,I)*GRD(2,1,I+1) 
        E1(2,2)=TU(2,1,I)*GRD(1,2,I+1)+TU(2,2,I)*GRD(2,2,I+1) 
        E2(1,1)=E1(1,1)*GTD(1,1,I)+E1(1,2)*GTD(2,1,I) 
        E2(1,2)=E1(1,1)*GTD(1,2,I)+E1(1,2)*GTD(2,2,I) 
        E2(2,1)=E1(2,1)*GTD(1,1,I)+E1(2,2)*GTD(2,1,I) 
        E2(2,2)=E1(2,1)*GTD(1,2,I)+E1(2,2)*GTD(2,2,I) 
        GRD(1,1,I)=RD(1,1,I)+E2(1,1) 
        GRD(1,2,I)=RD(1,2,I)+E2(1,2) 
        GRD(2,1,I)=RD(2,1,I)+E2(2,1) 
        GRD(2,2,I)=RD(2,2,I)+E2(2,2) 
C 
  200 CONTINUE 
C 
      RETURN 
      END 
C 
C 
C 
C *** P-SV MODIFIED R/T MATIX *** 
C 
      SUBROUTINE MDRTRY(OM,WN) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 CMS,CMP,EA(4,4),EB(4,4),DDA(2,2),DUB(2,2)
      COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
C
C *** CHECK FOR THE DIMENSION OF COMMON *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM MDRTPSV: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** MODIFIED R/T MATRIX *** 
C 
      CALL MTXRY(OM,WN,1,CMS,CMP,EA) 
      CALL DDURY(THK(1),CMS,CMP,DDA)
C 
      DO 100 I=1, NL-1 
C 
        CALL MTXRY(OM,WN,I+1,CMS,CMP,EB) 
        CALL DDURY(THK(I+1),CMS,CMP,DUB)
C 
        CALL MTRRY(EA,EB,DDA,DUB,I)  
C 
        DO 110 J=1,4 
          DO 120 K=1,4 
            EA(J,K)=EB(J,K) 
  120     CONTINUE 
  110   CONTINUE 
C 
        DDA(1,1)=DUB(1,1) 
        DDA(1,2)=DUB(1,2) 
        DDA(2,1)=DUB(2,1) 
        DDA(2,2)=DUB(2,2) 
C 
  100 CONTINUE 
C 
      RETURN 
      END 
C 
C 
C 
C *** P-SV MODIFIED R/T MATRIX FOR Ith LAYER *** 
C 
        SUBROUTINE MTRRY(EA,EB,DDA,DUB,L)  
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 EA(4,4),EB(4,4),DDA(2,2),DUB(2,2),TD,TU,RD,RU 
     *          ,E1(4,4),E2(4,4),E3(4,4) 
      COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
     *      /MRTRY/TD(2,2,20),TU(2,2,20),RD(2,2,20),RU(2,2,20) 
C 
C *** CHECK FOR THE DIMENSION OF COMMON *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM MTRRY: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** MATRIX *** 
C 
      DO 100 I=1,2 
        DO 110 J=1,2 
          E1(I,J)= EB(I,J) 
          E1(I,J+2)=-EA(I,J+2) 
          E1(I+2,J)= EB(I+2,J) 
          E1(I+2,J+2)=-EA(I+2,J+2) 
          E2(I,J)= EA(I,J) 
          E2(I,J+2)=-EB(I,J+2) 
          E2(I+2,J)= EA(I+2,J) 
          E2(I+2,J+2)=-EB(I+2,J+2) 
  110   CONTINUE 
  100 CONTINUE 
C 
      CALL DCINVA(E1) 
C 
      DO 200 I=1,4 
        DO 210 J=1,4 
          E3(I,J)=(0.D0,0.D0) 
          DO 220 K=1,4 
            E3(I,J)=E3(I,J)+E1(I,K)*E2(K,J)  
  220     CONTINUE 
  210   CONTINUE 
  200 CONTINUE 
C 
      TD(1,1,L)=E3(1,1)*DDA(1,1)+E3(1,2)*DDA(2,1)   
      TD(2,1,L)=E3(2,1)*DDA(1,1)+E3(2,2)*DDA(2,1)   
      TD(1,2,L)=E3(1,1)*DDA(1,2)+E3(1,2)*DDA(2,2)   
      TD(2,2,L)=E3(2,1)*DDA(1,2)+E3(2,2)*DDA(2,2)   
      RD(1,1,L)=E3(3,1)*DDA(1,1)+E3(3,2)*DDA(2,1)   
      RD(2,1,L)=E3(4,1)*DDA(1,1)+E3(4,2)*DDA(2,1)   
      RD(1,2,L)=E3(3,1)*DDA(1,2)+E3(3,2)*DDA(2,2)   
      RD(2,2,L)=E3(4,1)*DDA(1,2)+E3(4,2)*DDA(2,2)   
C 
      IF(L.EQ.NL)RETURN 
C 
      RU(1,1,L)=E3(1,3)*DUB(1,1)+E3(1,4)*DUB(2,1)   
      RU(2,1,L)=E3(2,3)*DUB(1,1)+E3(2,4)*DUB(2,1)   
      RU(1,2,L)=E3(1,3)*DUB(1,2)+E3(1,4)*DUB(2,2)   
      RU(2,2,L)=E3(2,3)*DUB(1,2)+E3(2,4)*DUB(2,2)   
      TU(1,1,L)=E3(3,3)*DUB(1,1)+E3(3,4)*DUB(2,1)   
      TU(2,1,L)=E3(4,3)*DUB(1,1)+E3(4,4)*DUB(2,1)   
      TU(1,2,L)=E3(3,3)*DUB(1,2)+E3(3,4)*DUB(2,2)   
      TU(2,2,L)=E3(4,3)*DUB(1,2)+E3(4,4)*DUB(2,2)   
C 
      RETURN 
      END 
C 
C 
C *** Inverse of 4x4 Matrix *** 
C 
        SUBROUTINE DCINVA(E)  
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 E(4,4),ee(4,4),aa 
     *          ,a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3   
C 
      do 100 i=1,4 
      do 100 j=1,4 
      ee(i,j)=e(i,j) 
 100  continue 
C 
      DO 200 I=1,4 
C 
        IF(I.EQ.1)THEN 
          J=2 
          K=3 
          M=4 
        ELSE 
          IF(I.EQ.2)THEN 
            J=3 
            K=4 
            M=1 
          ELSE 
            IF(I.EQ.3)THEN 
              J=4 
              K=1 
              M=2 
            ELSE 
              J=1 
              K=2 
              M=3 
            END IF 
          END IF 
        END IF 
C 
        A1=EE(I,4)*EE(J,1)-EE(I,1)*EE(J,4)   
        A2=EE(I,4)*EE(J,2)-EE(I,2)*EE(J,4)   
        A3=EE(I,4)*EE(J,3)-EE(I,3)*EE(J,4)   
        B1=EE(I,3)*EE(J,2)-EE(I,2)*EE(J,3)   
        B2=EE(I,1)*EE(J,3)-EE(I,3)*EE(J,1)   
        B3=EE(I,2)*EE(J,1)-EE(I,1)*EE(J,2)   
        C1=EE(K,3)*EE(M,2)-EE(K,2)*EE(M,3)   
        C2=EE(K,1)*EE(M,3)-EE(K,3)*EE(M,1)   
        C3=EE(K,2)*EE(M,1)-EE(K,1)*EE(M,2)   
        D1=EE(K,4)*EE(M,1)-EE(K,1)*EE(M,4)   
        D2=EE(K,4)*EE(M,2)-EE(K,2)*EE(M,4)   
        D3=EE(K,4)*EE(M,3)-EE(K,3)*EE(M,4)   
C 
        AA=A1*C1+A2*C2+A3*C3+B1*D1+B2*D2+B3*D3 
C 
        E(1,I)=(-EE(J,2)*D3+EE(J,3)*D2-EE(J,4)*C1)/AA 
        E(2,I)=( EE(J,1)*D3-EE(J,3)*D1-EE(J,4)*C2)/AA 
        E(3,I)=(-EE(J,1)*D2+EE(J,2)*D1-EE(J,4)*C3)/AA 
        E(4,I)=( EE(J,1)*C1+EE(J,2)*C2+EE(J,3)*C3)/AA 
C 
 200  continue 
C 
      return 
      end 
c 
C 
C *** P-SV ELEMENT MATIX *** 
C 
      SUBROUTINE MTXRY(OM,WN,L,CMS,CMP,EE) 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 CMS,CMP,EE(4,4),CMSG,CMPG 
      COMMON /NSTRC/DNS(20),VP(20),VS(20),THK(20),NL 
C 
C *** CHECK FOR THE DIMENSION OF COMMOM *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM MTXRY: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** ELEMENT MATRIX *** 
C 
      VS2=VS(L)*VS(L) 
      VP2=VP(L)*VP(L) 
      SM=DNS(L)*VS2 
C 
      IF(OM.EQ.0.D0)THEN 
C 
        X0=(1.D0+VS2/VP2)/(1.D0-VS2/VP2) 
        AA=WN*SM 
C 
        EE(1,1)=(1.D0,0.D0) 
        EE(1,2)=(1.D0,0.D0) 
        EE(1,3)=(1.D0,0.D0) 
        EE(1,4)=(1.D0,0.D0) 
        EE(2,1)=-(X0-1.D0) 
        EE(2,2)=(1.D0,0.D0) 
        EE(2,3)= X0-1.D0 
        EE(2,4)=(-1.D0,0.D0) 
        EE(3,1)= (X0-3.D0)*AA 
        EE(3,2)= -2.D0*AA 
        EE(3,3)=-(X0-3.D0)*AA 
        EE(3,4)= 2.D0*AA 
        EE(4,1)= (X0-1.D0)*AA 
        EE(4,2)= -2.D0*AA 
        EE(4,3)= (X0-1.D0)*AA 
        EE(4,4)=-2.D0*AA 
C 
      ELSE 
C 
        WN2=WN*WN 
        OM2=OM*OM 
C 
        WS2=OM2/VS2 
        CMS=CDSQRT(DCMPLX(WN2-WS2,0.D0)) 
        IF(DREAL(CMS).LT.0.D0)CMS=-CMS 
        CMSG=CMS/WN 
C 
        WP2=OM2/VP2 
        CMP=CDSQRT(DCMPLX(WN2-WP2,0.D0)) 
        IF(DREAL(CMP).LT.0.D0)CMP=-CMP 
        CMPG=CMP/WN 
C 
        X0=2.D0*WN-WS2/WN 
C 
        EE(1,1)=(-1.D0,0.D0) 
        EE(1,2)= CMSG 
        EE(1,3)=(-1.D0,0.D0) 
        EE(1,4)= CMSG 
        EE(2,1)=-CMPG 
        EE(2,2)=(1.D0,0.D0) 
        EE(2,3)= CMPG 
        EE(2,4)=(-1.D0,0.D0) 
        EE(3,1)= 2.D0*SM*CMP  
        EE(3,2)=-SM*X0  
        EE(3,3)=-2.D0*SM*CMP  
        EE(3,4)= SM*X0  
        EE(4,1)= SM*X0  
        EE(4,2)=-2.D0*SM*CMS  
        EE(4,3)= SM*X0  
        EE(4,4)=-2.D0*SM*CMS  
C 
      END IF 
C 
      RETURN 
      END 
C 
C 
C *** P-SV UP/DOWN-GOING MATRIX *** 
C 
      SUBROUTINE DDURY(THK,CMS,CMP,DU)
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      COMPLEX*16 DU(2,2),CMS,CMP 
C 
        DU(1,1)= CDEXP(-CMP*THK) 
        DU(1,2)= (0.D0,0.D0) 
        DU(2,1)= (0.D0,0.D0) 
        DU(2,2)= CDEXP(-CMS*THK) 
C 
      RETURN 
      END 

