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 phs3 
C 
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER(MOM=700) 
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 1995)' 
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"' 
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,CTOL,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 "grflt01s.in"' 
c      read(*,'(a20)')infile 
 
      infile='grflt01s.in' 
      OPEN (UNIT=7,FILE=infile,STATUS='OLD') 
C 
      read(7,1100)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,1110)period 
 1100 format(/f8.0,i8) 
 1110 format(f8.0) 
      duration=dt*nt 
      pi2=datan(1.d0)*8.d0 
      dom=pi2/duration 
      iom=0 
      nom=duration/period 
c 
      write(*,*) 
      write(*,*)'Check Secular Functions and Eigen Vectors? (yes=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,1200)NL 
      DO 100 IL=1,NL 
        READ(7,1210)DNS(IL),VP(IL),QP,VS(IL),QS,THK(IL) 
  100 CONTINUE 
 1200 FORMAT(/I10) 
 1210 FORMAT(6E8.0) 
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., 100):') 
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 
      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
   26 WRITE(*,1650) 
 1650 FORMAT
     *( ' Enter Number of Partitions from Cmin to Cmax (ex.,200):') 
      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) 
 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 
 1800 FORMAT(/' Output Eigen Vectors? (yes=1)              :') 
 1810 FORMAT( ' Enter Depth to Output the Eigen Vectors (m):') 
 1820 FORMAT( ' Enter Number of Partitions along the Depth :') 
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,CTOL,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:701,0:230),CCR(0:701,0:230),MDL(0:701),MDR(0:701) 
     *   /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=701 
      IF(MOM.GT.NC)THEN 
         WRITE(*,*)' FROM PHASE: MOM.GT.701 !!!!!!!!!!!!!!!' 
         WRITE(10,*)' FROM PHASE: MOM.GT.701 !!!!!!!!!!!!!!!' 
         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) 
        if(io.eq.0)om=0.001d0
        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 
c        IF(IO.EQ.0)THEN 
c          MDL(IO)=0 
c          CCL(IO,0)=VSMAX 
c          GOTO 20 
c        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') 
        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=CMIN 
        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) 
        end if 
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) 
          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) 
          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=230 
        IF(MDL(IO).GT.NC)THEN 
           WRITE(*,*)' FROM PHASE: LOVE MODE.GT.230 !!!!!!!!!!!!!!' 
           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) 
        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) 
          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) 
          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=230 
        IF(MDR(IO).GT.NC)THEN 
           WRITE(*,*)' FROM PHASE: RAYLEIGH MODE.GT.230 !!!!!!!!!!!!!!!' 
           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='phsl_dat',STATUS='unknown') 
      OPEN(23,FILE='phsr_dat',STATUS='unknown') 
        DO 200 MD=0,MDL(NOM) 
        DO 200 IO=IOM,NOM 
          WRITE(22,1100)CCL(IO,MD) 
  200   CONTINUE 
        DO 210 MD=0,MDR(NOM) 
        DO 210 IO=IOM,NOM 
          WRITE(23,1100)CCR(IO,MD) 
  210   CONTINUE 
      CLOSE(22) 
      CLOSE(23) 
c 
      OPEN(25,FILE='mdl_dat',STATUS='unknown') 
      OPEN(26,FILE='mdr_dat',STATUS='unknown') 
        DO 300 IO=IOM,NOM 
          WRITE(25,1000)MDL(IO) 
          WRITE(26,1000)MDR(IO) 
  300   CONTINUE 
 1000 FORMAT(I5) 
      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 
      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,CTOL,NITER,NCC 
     *   /GRTLV/GRU0L,GTDL(20),GRDL(20) 
     *   /GRTRY/GTDR(2,2,20),GRDR(2,2,20) 
     *   /CCMD/CCL(0:701,0:230),CCR(0:701,0:230),MDL(0:701),MDR(0:701) 
     *   /VCTL/VL1(702,0:230),VL2(702,0:230) 
     * /VCLR/VR1(702,0:230),VR2(702,0:230),VR3(702,0:230),VR4(702,0:230) 
     *   /UEGL/UUL(0:701,0:230),E1L(0:701,0:230) 
     *   /UEGR/UUR(0:701,0:230),E1R(0:701,0:230) 
C 
C *** CHECK OF DIMENSION *** 
C 
      NC=20 
      IF(NL.GT.NC)THEN 
         WRITE(*,*)' FROM ENGUU: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
      NC=702 
      IF((NZINT+2).GT.NC)THEN 
         WRITE(*,*)' FROM ENGUU: NZINT+2.GT.702 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
      NC=701 
      IF(MOM.GT.NC)THEN 
         WRITE(*,*)' FROM ENGUU: MOM.GT.701 !!!!!!!!!!!!!!!' 
         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',STATUS='unknown') 
      OPEN(25,FILE='Lstr',STATUS='unknown') 
      OPEN(26,FILE='Rdis',STATUS='unknown') 
      OPEN(27,FILE='Rstr',STATUS='unknown') 
      END IF 
 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)=DDS*CDS 
              VL2(IZ,MD)=EL(2,1)*DDS*CDS 
            ELSE 
              DUS=CDEXP(-CMS*(BDEP-Z)) 
              VL1(IZ,MD)=DDS*CDS+DUS*CUS 
              VL2(IZ,MD)=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)=E1L(IO,MD)+DNS(LAY)*DR1 
            E2L       =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.001 
          VSTR=0.001 
          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 
        do IZ=1,NZINT+1 
          Z=DZ*(IZ-1) 
          write(24,2500)z/1000.0,(VL1(IZ,MD)/2.0+MD,MD=0,MDL(IO)) 
          write(25,2500)z/1000.0,(VL2(IZ,MD)/2.0+MD,MD=0,MDL(IO)) 
c          write(24,2500)z/1000.0,(VL1(IZ,MD),MD=0,MDL(IO)) 
c          write(25,2500)z/1000.0,(VL2(IZ,MD),MD=0,MDR(IO)) 
        end do 
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)=ER(1,1)*A1+ER(1,2)*A2 
              VR2(IZ,MD)=ER(2,1)*A1+ER(2,2)*A2 
              VR3(IZ,MD)=ER(3,1)*A1+ER(3,2)*A2 
              VR4(IZ,MD)=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)=ER(1,1)*A1+ER(1,2)*A2 
                VR2(IZ,MD)=ER(2,1)*A1+ER(2,2)*A2 
                VR3(IZ,MD)=ER(3,1)*A1+ER(3,2)*A2 
                VR4(IZ,MD)=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)=ER(1,1)*A1+ER(1,2)*A2+ER(1,3)*A3+ER(1,4)*A4 
                VR2(IZ,MD)=ER(2,1)*A1+ER(2,2)*A2+ER(2,3)*A3+ER(2,4)*A4 
                VR3(IZ,MD)=ER(3,1)*A1+ER(3,2)*A2+ER(3,3)*A3+ER(3,4)*A4 
                VR4(IZ,MD)=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)=E1R(IO,MD)+DNS(LAY)*(DR1+DR2) 
            E2R =E2R +A*DR1+SM*DR2 
            E3R =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.001 
          VSTR=0.001
          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.0,(VR1(IZ,MD)/2.0+MD 
     *                           ,VR2(IZ,MD)/2.0+MD,MD=0,MDR(IO)) 
          write(27,2500)z/1000.0,(VR3(IZ,MD)/2.0+MD 
     *                           ,VR4(IZ,MD)/2.0+MD,MD=0,MDR(IO)) 
c          write(26,2500)z/1000.0,(VR1(IZ,MD),VR2(IZ,MD),MD=0,MDL(IO)) 
c          write(27,2500)z/1000.0,(VR3(IZ,MD),VR4(IZ,MD),MD=0,MDR(IO)) 
         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 
c 
      OPEN(32,FILE='grpl_dat',STATUS='unknown') 
      OPEN(33,FILE='grpr_dat',STATUS='unknown') 
C 
      DO 200 MD=0,MDL(NOM) 
      DO 200 IO=IOM,NOM 
        WRITE(32,2100)UUL(IO,MD) 
  200 CONTINUE 
C 
      DO 300 MD=0,MDR(NOM) 
      DO 300 IO=IOM,NOM 
        WRITE(33,2100)UUR(IO,MD) 
  300 CONTINUE 
 2100 format(e15.7) 
C 
      CLOSE(32) 
      CLOSE(33) 
c 
      OPEN(24,FILE='Ldisper',STATUS='unknown') 
      OPEN(25,FILE='Rdisper',STATUS='unknown') 
      do IO=IOM,NOM 
        Hz=DOM*DFLOAT(IO)/6.283185307d0 
        WRITE(24,2500)Hz,(CCL(IO,MD),UUL(IO,MD),MD=0,MDL(NOM)) 
        WRITE(25,2500)Hz,(CCR(IO,MD),UUR(IO,MD),MD=0,MDR(NOM)) 
      end do 
 2500 FORMAT(f8.4,200e12.4) 
      CLOSE(24) 
      CLOSE(25) 
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),TD,TU,RD,RU 
     *       ,DDA(2,2),DUB(2,2) 
      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 MDRTPSV: NL.GT.20 !!!!!!!!!!!!!!!' 
         STOP 
      END IF 
C 
C *** MODIFIED R/T MATRIX *** 
C 
      CALL MTXRY(OM,WN,1,CMS,CMP,EA) 
      CALL DDURY(OM,WN,THK(1),CMS,CMP,DDA) 
C 
      DO 100 I=1, NL-1 
C 
        CALL MTXRY(OM,WN,I+1,CMS,CMP,EB) 
        CALL DDURY(OM,WN,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(OM,WN,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 
C 
C 
C *** BESSEL FUNCTION J0(X) AND N0(X) *** 
C 
      SUBROUTINE BESJN0(X,BJ0,BN0) 
C 
      IMPLICIT REAL*8 (A-H,O-Z) 
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10 
     * /  -2.2270614         D-9,  7.85940046        D-8, 
     *    -1.9891965786      D-6,  4.03109603859     D-5, 
     *    -6.449978225950    D-4,  7.9012341828039   D-3, 
     *    -7.11111110046304  D-2,  4.444444444257695 D-1, 
     *    -1.7777777777758446D 0,  3.9999999999998966D 0, 
     *    -3.9999999999999978D 0                         / 
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13 
     * /   3.315018          D-11, 4.697063          D-11, 
     *    -5.66336210        D-9, -4.88429867        D-9, 
     *     6.7810966225      D-7,  2.1332909858      D-7, 
     *    -5.379500918201    D-5,  1.120221502384    D-5, 
     *     2.52021197018629  D-3, -1.70738759690737  D-3, 
     *    -5.614869347449754 D-2,  5.604047189802616 D-2, 
     *     3.2757913759146509D-1, -1.7759677131433834D-1 / 
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14 
     * /  -1.98650           D-12, -1.517802          D-11, 
     *     3.9649790         D-10,  2.40512241        D-9, 
     *    -5.707143476       D-8,  -2.5797893105      D-7, 
     *     5.66074614693     D-6,   1.722612595802    D-5, 
     *    -3.5284910023231   D-4,  -5.9314897388202   D-4, 
     *     1.179012935710246 D-2,   6.39612989784324  D-3, 
     *    -1.5037412265137390D-1,   4.68282348234590  D-3, 
     *     3.0007927051955560D-1                          / 
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10 
     * /   2.8228605          D-10, -1.94155134      D-9, 
     *     6.51602613         D-9,  -1.533982923     D-8, 
     *     3.267679187        D-8,  -8.099284883     D-8, 
     *     2.8872150072       D-7,  -1.74249896802   D-6, 
     *     2.184678334973     D-5,  -8.7658215906110 D-4, 
     *     0.79788456080286526D 0                         / 
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12 
     * /  -8.9331052         D-10,  6.65327442      D-9, 
     *    -2.322479945       D- 8,  5.187839807     D-8, 
     *    -8.765163373       D- 8,  1.3089752200    D-7, 
     *    -2.0632292049      D- 7,  4.0856346318    D-7, 
     *    -1.15933062842     D- 6,  5.25865147508   D-6, 
     *    -4.423973781374    D- 5,  9.1310641570497 D-4, 
     *    -9.973557010035815 D- 2                         / 
      DATA H0,H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11 
     * /   3.4432403         D-9,  -1.175592825       D-7, 
     *     2.8512737759      D-6,  -5.49344073691     D-5, 
     *     8.276591200190    D-4,  -9.4202584419331   D-3, 
     *     7.72372081288392  D-2,  -4.261441276304766 D-1, 
     *     1.4216343895861494D 0,  -2.3498510134188587D 0, 
     *     1.0766114686838581D 0,   3.6746690519661597D-1  / 
      DATA G0,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
     * /   1.50410           D-12, -5.72117           D-12, 
     *     4.050146          D-11, -6.4960948         D-10, 
     *     5.3845950         D-10,  5.922131699       D- 8, 
     *     6.874770130       D- 8, -6.56197858169     D- 6, 
     *     5.9100539363      D- 7,  3.9226240684696   D- 4, 
     *    -1.5291708227790   D- 4, -1.348097236217831 D- 2, 
     *     1.238842781168767 D- 2,  1.6904512696363842D- 1, 
     *    -1.4786314339122686D- 1, -3.0851762524903378D- 1  / 
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14 
     * /  -1.14460           D-12,  2.930515          D-11, 
     *     1.9780007         D-10, -4.93772861        D- 9, 
     *    -2.641221656       D- 8,  5.9964479378      D- 7, 
     *     2.24627737107     D- 6, -4.765550355443    D- 5, 
     *    -1.1239446152872   D- 4,  2.24134018235420  D- 3, 
     *     2.36753744518725  D- 3, -4.900342975478215 D- 2, 
     *    -8.64421637526570  D- 3,  3.0266723702418476D- 1, 
     *    -2.5949743967209265D- 2                           / 
      DATA PI2,PI4 
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 / 
C 
      AX=DABS(X) 
      IF (AX .EQ. 0.0D0 ) GO TO 50 
      IF (AX .GE. 1.0D12) GO TO 40 
      IF (AX .GE. 8.0D0 ) GO TO 30 
      IF (AX .GE. 6.0D0 ) GO TO 20 
      IF (AX .GE. 4.0D0 ) GO TO 10 
C 
C     !X! < 4 
C 
      Z=X*X*0.0625D0 
      BJ0=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z 
     *    +C7)*Z +C8)*Z +C9)*Z +C10)*Z + 1.0D0 
      BN0=((((((((((H0*Z +H1)*Z +H2)*Z +H3)*Z +H4)*Z +H5)*Z +H6)*Z 
     *    +H7)*Z +H8)*Z +H9)*Z +H10)*Z +H11 
      BN0=PI2*BJ0*DLOG(AX*0.5D0) +BN0 
      GO TO 9999 
C 
C     4 <= !X! < 6 
   10 Z=AX-5.0D0 
      BJ0=((((((((((((D0*Z+ D1)*Z+ D2)*Z+ D3)*Z+ D4)*Z+ D5)*Z 
     *       +D6)*Z+ D7)*Z+ D8)*Z+ D9)*Z+ D10)*Z+ D11)*Z+ D12)*Z+ D13 
      BN0=((((((((((((((G0*Z +G1)*Z +G2)*Z +G3)*Z +G4)*Z +G5)*Z 
     *       +G6)*Z +G7)*Z +G8)*Z +G9)*Z +G10)*Z +G11)*Z +G12)*Z 
     *       +G13)*Z +G14)*Z +G15 
      GO TO 9999 
C 
C     6 <= !X! < 8 
C 
   20 Z=AX-7.0D0 
      BJ0=(((((((((((((E0*Z+ E1)*Z+ E2)*Z+ E3)*Z+ E4)*Z+ E5)*Z 
     *       + E6)*Z+ E7)*Z+ E8)*Z+ E9)*Z+ E10)*Z+ E11)*Z+ E12)*Z 
     *       + E13)*Z+ E14 
      BN0=(((((((((((((F0*Z +F1)*Z +F2)*Z +F3)*Z +F4)*Z +F5)*Z 
     *       +F6)*Z +F7)*Z +F8)*Z +F9)*Z +F10)*Z +F11)*Z +F12)*Z 
     *       +F13)*Z +F14 
      GO TO 9999 
C 
C     8 <= !X! < 1.0D12 
C 
   30 Z=64.0D0/(AX*AX) 
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z 
     *  +A7)*Z +A8)*Z +A9)*Z +A10 
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z 
     *  +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX 
      Q=AX-PI4 
      CC= DCOS(Q) 
      CS= DSIN(Q) 
      CQ= DSQRT(AX) 
      BJ0=(A*CC - B*CS) /CQ 
      BN0=(A*CS + B*CC) /CQ 
      GO TO 9999 
C 
C     PARAMETER ERROR 
C 
   40 BJ0= 0.0D0 
      BN0= 0.0D0 
      GO TO 9999 
C 
   50 BJ0=  1.0D 0 
C      BN0= -1.0D75 
      BN0= -1.0D30 
C 
 9999 RETURN 
      END 
C 
C 
C 
C 
C *** BESSEL FUNCTION J1(X), N1(X) 
C 
      SUBROUTINE BESJN1(X,BJ1,BN1) 
C 
      IMPLICIT REAL*8 (A-H,O-Z) 
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9 
     * /  1.22494472        D-8, -3.929965711       D-7, 
     *    8.9514968720      D-6, -1.612440490740    D-4, 
     *    2.2574925903842   D-3, -2.37037026763209  D-2, 
     *    1.777777775582676 D-1, -8.888888888615124 D-1, 
     *    2.6666666666649150D 0, -3.9999999999999563D 0  / 
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11 
     * /  2.3353113         D-9, -8.71498135        D-8, 
     *    2.3316556575      D-6, -5.00978325295     D-5, 
     *    8.533226991090    D-4, -1.11766285197065  D-2, 
     *    1.080130262837708 D-1, -7.271013435638566 D-1, 
     *    3.1262109002217823D 0, -7.3980378533256488D 0, 
     *    6.8529249642057360D 0,  3.9325615210260295D-1  / 
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14,E15 
     * /   4.84682           D-12, -2.280090          D-11, 
     *     6.313256          D-11, -5.2574110         D-10, 
     *     7.81894200        D- 9, -5.92404185        D- 9, 
     *    -5.9222983079      D- 7, -6.1872869462      D- 7, 
     *     5.249583490136    D- 5, -4.13703795176     D- 6, 
     *    -2.35357444227453  D- 3,  7.6458541141893   D- 4, 
     *     5.392388944881264 D- 2, -3.716528343506466 D- 2, 
     *    -3.3809025392727920D- 1,  1.4786314339122685D- 1  / 
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14 
     * /   1.78737           D-12,  1.601135          D-11, 
     *    -3.8677585         D-10, -2.37356163        D- 9, 
     *     5.432238759       D- 8,  2.6412212075      D- 7, 
     *    -5.39680775206     D- 6, -1.797021894410    D- 5, 
     *     3.3358852634720   D- 4,  6.7436676916589   D- 4, 
     *    -1.120670091199092 D- 2, -9.47014978074827  D- 3, 
     *     1.4701028926435867D- 1,  1.728843275053138 D- 2, 
     *    -3.0266723702418487D- 1                           / 
      DATA G0,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14 
     * /  -1.07463           D-12,  2.778003          D-11, 
     *     2.0080 677        D-10, -4.75788198        D- 9, 
     *    -2.6460 77937      D- 8,  5.7071424129      D- 7, 
     *     2.32181314997     D- 6, -4.528596911743    D- 5, 
     *    -1.2058288258765   D- 4,  2.11709460137864  D- 3, 
     *     2.96574486954234  D- 3, -4.716051742840817 D- 2, 
     *    -1.918838969353708 D- 2,  3.0074824530274775D- 1, 
     *    -4.68282348234583  D- 3                           / 
      DATA H0,H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12,H13,H14 
     * /   2.21931           D-12,  4.11634           D-12, 
     *    -4.3816515         D-10, -5.7599652         D-10, 
     *     6.230613774       D- 8,  4.885713666       D- 8, 
     *    -6.10299268189     D- 6, -1.70664050678     D- 6, 
     *     3.7656506609458   D- 4, -6.721328811703    D- 5, 
     *    -1.260105985120454 D- 2,  6.82955038740438  D- 3, 
     *     1.6844608042350780D- 1, -1.1208094379604528D- 1, 
     *    -3.2757913759146522D- 1                           / 
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10 
     * /  -3.0207732         D-10,  2.08062912       D-9, 
     *    -7.00035699        D- 9,  1.655768548      D-8, 
     *    -3.558232648       D- 8,  8.955040706      D-8, 
     *    -3.2722783265      D- 7,  2.05931893593    D-6, 
     *    -2.808872167333    D- 5,  1.46097026511633 D-3, 
     *     7.9788456080286546D- 1                           / 
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12 
     * /   9.4650569         D-10, -7.05369184       D-9, 
     *     2.464612368       D- 8, -5.514232274      D-8, 
     *     9.343285420       D- 8, -1.4024484043     D-7, 
     *     2.2297669562      D- 7, -4.4752033584     D-7, 
     *     1.29573186799     D- 6, -6.06767597238    D-6, 
     *     5.407079074732    D- 5, -1.27834898198997 D-3, 
     *     2.9920671030107448D- 1                           / 
      DATA PI2,PI4 
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 / 
C 
      AX=DABS(X) 
C      IF (AX .LT. 8.797D-77) GO TO 60 
      IF (AX .LT. 8.797D-30) GO TO 60 
      IF (AX .GE. 1.0D12) GO TO 50 
      IF (AX .GE. 8.0D0) GO TO 30 
      IF (AX .GE. 6.0D0) GO TO 20 
      IF (AX .GE. 4.0D0) GO TO 10 
C 
C     0 < !X! < 4 
C 
      Y=X*0.25D0 
      Z=Y*Y 
      BJ1=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z 
     *       +C7)*Z +C8)*Z +C9)*Z + 2.0D0) *Y 
      BN1=((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z +D6)*Z 
     *       +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z -PI2) /X 
      BN1=PI2*BJ1*DLOG(AX*0.5D0) +BN1 
      GO TO 9999 
C 
C     4 <= !X! < 6 
C 
   10 Z=X-5.0D0 
      BJ1=(((((((((((((H0*Z +H1)*Z +H2)*Z +H3)*Z +H4)*Z +H5)*Z 
     *       +H6)*Z +H7)*Z +H8)*Z +H9)*Z +H10)*Z +H11)*Z +H12)*Z 
     *       +H13)*Z +H14 
      BN1=((((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z 
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z 
     *       +E13)*Z +E14)*Z +E15 
      GO TO 40 
C 
C     6 <= !X! < 8 
C 
   20 Z=X-7.0D0 
      BJ1=(((((((((((((G0*Z +G1)*Z +G2)*Z +G3)*Z +G4)*Z +G5)*Z 
     *       +G6)*Z +G7)*Z +G8)*Z +G9)*Z +G10)*Z +G11)*Z +G12)*Z 
     *       +G13)*Z +G14 
      BN1=(((((((((((((F0*Z +F1)*Z +F2)*Z +F3)*Z +F4)*Z +F5)*Z 
     *       +F6)*Z +F7)*Z +F8)*Z +F9)*Z +F10)*Z +F11)*Z +F12)*Z 
     *       +F13)*Z +F14 
      GO TO 40 
C 
C     8 <= !X! < 1.0D12 
C 
   30 Z=64.0D0/(X*X) 
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z +A7)*Z 
     *       +A8)*Z +A9)*Z +A10 
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z 
     *       +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX 
      Q=X-PI4 
      CC= DCOS(Q) 
      CS= DSIN(Q) 
      CQ= DSQRT(AX) 
      BJ1=(A*CS + B*CC) /CQ 
      BN1=(B*CS - A*CC) /CQ 
C 
   40 IF (AX .LT. 0.0D0) THEN 
      BJ1=-BJ1 
      BN1=-BN1 
      END IF 
      GO TO 9999 
C 
C     PARAMETER ERROR 
C 
   50 BJ1= 0.D0 
      BN1= 0.D0 
      GO TO 9999 
C 
   60 BJ1= 0.0D0 
C      BN1=-1.0D75 
      BN1=-1.0D30 
C 
 9999 RETURN 
      END 

