C
      PROGRAM MRICKR
C
C ** MARAMETER **
C MOM   : NOM IN MAIN PROGRAM DIMENSION
C       : IF MOM > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C MMD   : NEMD IN MAIN PROGRAM DIMENSION
C ML    : NL (LAYER NUMBER)
C ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCY
C DNS(I): DENSITY OF Ith LAYER 
C VS (I): VS OF Ith LAYER      
C DEP(I): DEPTH ( THICKNESS ) OF Ith LAYER   
C NL    : TOTAL NUMBER OF LAYERS             
C NDC   : DIVISION NUMBER FROM MIN VS TO MAX VS TO PHASE VELOCITY
C NSMD  : START OUTPUT MODE OF RAYLEIGH WAVE
C NEMD  : END OUTPUT MODE OF RAYLEIGH WAVE
C NDOM  : DIVISION NUMBER FROM 0 TO MAX OMEGA TO CUT-OFF OMEGA
C ** CALCULATE DATA **
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE RAYLEIGH WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : RAYLEIGH WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): OUTPUT HIGHEST RAYLEIGH MODE NUMBER AT Ith OMEGA
C UU(I,J) : GROUP VELOCITY OF Jth MODE RAYLEIGH WAVE AT Ith OMEGA
C
      PARAMETER(MOM=50,MMD=7,ML=10,MOB=5)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 FW0(MOM,MOB),FU0(MOM,MOB),CPX,EKX,FWW,FUU,CCC
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),XOB(MOB)
     *         ,ARWZ(MOM,0:MMD),ARWX(MOM,0:MMD)
     *         ,ARUZ(MOM,0:MMD),ARUX(MOM,0:MMD)
     *         ,DNS(ML),VP(ML),VS(ML),THK(ML) 
      COMMON /SODAT/ SDEP,LS,ODEP,LO
      CHARACTER FNAME*15
C
      WRITE(*,1000)
 1000 FORMAT(' MRICKR: RICKER WAVE IN 2-D IN-PLANE LAYERED'/
     *       '         HALF-SPACE USING NORMAL MODE SOLUTION'/
     *       '        (SP AND OB ARE IN ANY LAYER),')
C
C *** DATA INPUT ***
C
      CALL RAYINN(DOM,NOM,ML,NL,DNS,VP,VS,THK,NDOM,NSMD,NEMD,NDC)
C
C *** PARAMETER CHECK ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MRCUAV :'
        WRITE(*,*)'MOM < NOM THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(ML.LT.NL)THEN
        WRITE(*,*)'FROM MAIN MRCUAV :'
        WRITE(*,*)'ML < NL THUS STOP !!!'
        WRITE(*,*)'ML,NL =',ML,NL
        STOP
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MRCUAV :'
        WRITE(*,*)'MMD < NEMD THUS STOP !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        STOP
      END IF
C
C
C *** INPUT OF SP AND OP DEPTH IN ANY LAYER ***
C
   10 WRITE(*,1005)NL
      READ(*,*)LS
      IF(LS.EQ.0)GOTO 10
      IF(LS.GT.NL)GOTO 10
      WRITE(*,1010)THK(LS)
      READ(*,*)SDEP
 1005 FORMAT(/' Total layers number =',I5,
     *       /' Enter source layer number (from top layer) : ',$)
 1010 FORMAT( ' Source layer thickness =',F10.4,
     *       /' Enter source depth from upper boundary : ',$)
C
   20 WRITE(*,1100)NL
      READ(*,*)LO
      IF(LO.EQ.0)GOTO 20
      IF(LO.GT.NL)GOTO 20
      WRITE(*,1110)THK(LO)
      READ(*,*)ODEP
 1100 FORMAT(/' Total layers number =',I5,
     *       /' Enter OP layer number (from top layer) : ',$)
 1110 FORMAT( ' Source layer thickness =',F10.4,
     *       /' Enter OP depth from upper boundary : ',$)
C
   30 WRITE(*,1410)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 30
      END IF
      IF(NOB.EQ.0)GOTO 30
C
      DO 100 I=1,NOB
        WRITE(*,1420)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 1410 FORMAT(/' Enter Number of Observation Point : ',$)
 1420 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',$)
C
      WRITE(*,1521)
      READ(*,*)TP,TS
C
 1521 FORMAT(/' Enter Tp and Ts (sec) of Ricker wavelet : ',$)
C
      WRITE(*,1550)
      READ(*,*)QX,QZ
C
 1550 FORMAT(/' Enter QX and QZ : ',$)
C
      WRITE(*,1640)
      READ(*,*)IQF
      IF(IQF.EQ.0)THEN
        WRITE(*,1650)
        READ(*,*)QFL
      END IF
 1640 FORMAT(/' Include Quality Factor (yes=0)?   : ',$)
 1650 FORMAT( ' Enter Quality Factor for Rayleigh Wave: ',$)
C
C *** CUT-OFF OMEGA ***
C
       CALL RAYCTN(DOM,NOM,ML,NL,DNS,VP,VS,THK
     *            ,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITY ***
C
      CALL RAYCCN(DOM,NOM,MOM,ML,NL,DNS,VP,VS,THK
     *           ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      CALL RAYUAR(DOM,NOM,MOM,ML,NL,DNS,VP,VS,THK
     *           ,CC,NMODE,CFOM,NSMD,MMD,UU,ARWZ,ARWX,ARUZ,ARUX)
C
C *** RICKER WAVE ***
C
      IF(IQF.NE.0)QFL=9999.D0
      CPX=DCMPLX(0.D0,1.D0)
      AA=DNS(1)*VS(1)*VS(1)
      BB=TP*TP*0.0253302D0
C
      DO 200 IOB=1,NOB
C
        XX=DABS(XOB(IOB))
C
        DO 210 IOM=1,NOM
C
          OM=DOM*DFLOAT(IOM)
          OM2=OM*OM
          CCC=AA*OM2*DEXP(-OM2*BB)*CDEXP(DCMPLX(0.D0,-OM*TS))
C
          FW0(IOM,IOB)=(0.D0,0.D0)
          FU0(IOM,IOB)=(0.D0,0.D0)
C
          DO 220 MODE=0,NMODE(IOM)
C
            IF(CC(IOM,MODE).EQ.0.D0)THEN
              EKX=DCMPLX(1.D0,0.D0)
              EXQ=1.D0
            ELSE
              EKX=CDEXP(-CPX*OM/CC(IOM,MODE)*XX)
              EXQ=DEXP(-OM*XX/(2.D0*QFL*UU(IOM,MODE)))
            END IF
C
            IF(IQF.NE.0)EXQ=1.D0
C
            IF(XOB(IOB).LT.0.D0)THEN
              ARWX(IOM,MODE)=-ARWX(IOM,MODE)
              ARUZ(IOM,MODE)=-ARUZ(IOM,MODE)
            END IF
C
            FWW=(CPX*ARWZ(IOM,MODE)*QZ+ARWX(IOM,MODE)*QX)*EKX*EXQ
            FUU=(ARUZ(IOM,MODE)*QZ+CPX*ARUX(IOM,MODE)*QX)*EKX*EXQ
C
            FW0(IOM,IOB)=FW0(IOM,IOB)+FWW
            FU0(IOM,IOB)=FU0(IOM,IOB)+FUU
C
  220     CONTINUE
C
          FW0(IOM,IOB)=FW0(IOM,IOB)*CCC
          FU0(IOM,IOB)=FU0(IOM,IOB)*CCC
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,2400)DOM,NOM
        DO 300 IOB=1,NOB
        WRITE(*,2410)XOB(IOB)
          DO 310 K=1,NOM,4
            WRITE(*,2420)(FW0(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
        DO 320 IOB=1,NOB
          WRITE(*,2415)XOB(IOB)
          DO 330 K=1,NOM,4
            WRITE(*,2420)(FU0(K+J,IOB),J=0,3)
  330     CONTINUE
  320   CONTINUE
C
 2400 FORMAT(/' DOM,NOM=',F6.2,I5)
 2410 FORMAT(/' X=',F10.2,'; FW0*4')
 2415 FORMAT(/' X=',F10.2,'; FU0*4')
 2420 FORMAT(1X,8D10.3)
C
C *** OUTPUT OF DATA FILE ***
C
      WRITE(*,*)
   50 WRITE(*,2310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 50
      WRITE(*,2320)
 2310 FORMAT( ' Enter file name : ',$)
 2320 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VV roop')
C      
      OPEN(8,FILE=FNAME,STATUS='NEW')
        WRITE(8,2500)DOM,NOM
        WRITE(8,2510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,5
            WRITE(8,2520)(FW0(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
        DO 520 IOB=1,NOB
          DO 530 K=1,NOM,5
            WRITE(8,2520)(FU0(K+J,IOB),J=0,4)
  530     CONTINUE
  520   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(1X,F15.8,I15)
 2510 FORMAT(///)
 2520 FORMAT(1X,10F8.5)
C
      STOP
      END
C
C *** READ OF SOIL DATA FROM DATA FILE
C
      SUBROUTINE RAYINN(DOM,NOM,ML,NL,DNS,VP,VS,DEP,NDOM,NSMD,NEMD,NDC)
C
C ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCY
C NDOM  : DIVISION NUMBER FROM 0 TO MAX OMEGA TO GET CUT-OFF OMEGA
C NSMD  : START OUTPUT MODE OF RAYLEIGH WAVE
C NEMD  : END OUTPUT MODE OF RAYLEIGH WAVE
C NDC   : DIVISION NUMBER FROM MIN TO MAX VS TO GET PHASE VELOCITY
C DNS(I): DENSITY OF Ith LAYER 
C VP (I): VP OF Ith LAYER                              
C VS (I): VS OF Ith LAYER                              
C DEP(I): DEPTH ( THICKNESS ) OF Ith LAYER             
C ML    : NL IN MAIN PARAMETER
C NL    : TOTAL NUMBER OF LAYERS                       
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30
      DIMENSION DNS(ML),VP(ML),VS(ML),DEP(ML)
      COMMON / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
C
C *** SOIL DATA INPUT FROM FILE DATA ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1000 FORMAT( ' From SUBROUTINE RAYINN :')
 1010 FORMAT(/' Enter Soil Data File Name ',$)

      OPEN (7,FILE=FNAME,STATUS='OLD')
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VP(IL),VS(IL),DEP(IL),QF
  100 CONTINUE
      CLOSE(7)
 1200 FORMAT(I10)
 1210 FORMAT(5F10.0)
C
C *** DATA ECHO ***
C
      WRITE(*,2000)NL
      DO 200 IL=1,NL
        WRITE(*,2010)IL,DNS(IL),VP(IL),VS(IL),DEP(IL)
  200 CONTINUE
 2000 FORMAT(/' Total Layers Number            =',I5)
 2010 FORMAT( ' Layer Num., Density, Vp, Vs, Depth =',I3,4F10.3)
C
C *** MIN AND MAX VS AND SOIL DATA CHECK ***
C
      VSMIN=VS(NL)
      VSMAX=0.0D0
      VPMAX=0.0D0 
      MINIL=1
      MAXIL=1
      DO 300 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))
  300 CONTINUE
C
      IF(VSMAX.NE.VS(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE RAYINN :'
        WRITE(*,*)' VS(NL).NE.VSMAX, THUS STOP !!!'
        STOP
      END IF
      IF(VPMAX.NE.VP(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE RAYINN :'
        WRITE(*,*)' VP(NL).NE.VPMAX, THUS STOP !!!'
        STOP
      END IF
C
C *** MINUMUM AND MAXIMUM PHASE VELOCITY ***
C
      CCMIN=FNCCR(VP(MINIL),VS(MINIL))
      CCMAX=FNCCR(VP(MAXIL),VS(MAXIL))
C
C *** DATA INPUT FROM SCREEN ***
C
      WRITE(*,3000)
      READ(*,*)DOM,NOM
      OMMAX=DOM*DFLOAT(NOM)
      WRITE(*,3010)0.D0,OMMAX
      WRITE(*,3020)
      READ(*,*)NDOM
      WRITE(*,3030)CCMIN,VSMAX
      WRITE(*,3040)
      READ(*,*)NDC
      WRITE(*,3050)
      READ(*,*)NSMD
      WRITE(*,3060)
      READ(*,*)NEMD
      WRITE(*,3070)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 10
 3000 FORMAT(/' Enter Delta Omega and Omega Number        ? :',$)
 3010 FORMAT(/' Minimum Omega and Maximum Omega = ',2F8.3)
 3020 FORMAT( ' To Get Cut-Off Omega,  Enter',
     *       /' Division Number between Min and Max Omega ? :',$)
 3030 FORMAT(/' Minimum C and Maximum Vs = ',2F10.3)
 3040 FORMAT( ' To Get General Phase Velocity (C), Enter',
     *       /' Division Number between Min and Max C     ? :',$)
 3050 FORMAT(/' Enter Output First Rayleigh Mode ( >= 0 ) ? :',$)
 3060 FORMAT( ' Enter Output Final Rayleigh Mode          ? :',$)
 3070 FORMAT(/' Input Data OK ? ( yes = Enter 0 )         ? :',$)
      WRITE(*,*)
C
      RETURN
      END
C
C *** CUT-OFF OMEGA OF RAYLEIGH WAVE (C=VS(NL))
C
      SUBROUTINE RAYCTN(DOM,NOM,ML,NL,DNS,VP,VS,DEP
     *                 ,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C ** INPUT DATA
C DOM  : INCREMENT OF CIRCULAR FREQUENCY (DELTA OMEGA)
C NOM  : TOTAL NUMBER OF CIRCULAR FREQUENCY
C NDOM : DIVISION NUMBER from 0 to MAX OMEGA TO GET CUT-OFF OMEGA
C NEMD : END OUTPUT MODE OF RAYLEIGH WAVE
C MMD  : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CFOM(I): CUT-OFF OMEGA OF Ith MODE RAYLEIGH WAVE
C NCFOM  : NUMBER OF CUT-OFF OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD)
     *         ,DNS(ML),VP(ML),VS(ML),DEP(ML) 
      COMMON / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
      EXTERNAL FCFRN
C
C *** 0 th MODE CUT-OFF OMEGA ***
C
      ICO=0
      CFOM(ICO)=0.D0
      IF(NEMD.EQ.0)GOTO 40
C
C *** INCREMENT OF OMEGA ***
C
      DDOM=DFLOAT(NOM)*DOM/DFLOAT(NDOM)
C
C *** INPUT OF CONVERGENCE JUDGMENT CRITERIA FOR ROOTN ***
C
      DXEPS=100.0D0
      XEPSO=DDOM/DXEPS
      DFEPS=1000.0D0
      ITMAX=100
C     WRITE(*,*)
C     WRITE(*,*)'From SUBROUTINE RAYCTN :'
C     WRITE(*,*)
C     WRITE(*,*)'Convergence Criteria for Cut-Off Omega :'
C     WRITE(*,1000)XEPSO
C     WRITE(*,1010)DFEPS
C     WRITE(*,1020)ITMAX
C     WRITE(*,1030)
C     READ(*,*)IYES
C     IF(IYES.NE.0)THEN
C       WRITE(*,*)
C       WRITE(*,1040)
C       READ(*,*)XEPSO,DFEPS,ITMAX
C     END IF
C     WRITE(*,*)
C1000 FORMAT(/' XEPS = DDOM/DXEPS          : XEPS =',F15.10)
C1010 FORMAT( ' FEPS = DABS(FR1-FR2)/DFEPS : DFEPS =',F15.3)
C1020 FORMAT( ' Maximum Iteration Number   : ITMAX =',I15/)
C1030 FORMAT( ' Are These Criteria OK ? ( yes = 0 ) :',\)
C1040 FORMAT( ' Enter New DDOM, DXEPS, DFEPS, ITMAX :',\)
C
C *** LOWER TRIAL OMEGA AND EIGEN VALUE EQUATION ***
C
      OM1=DDOM
      FR1=FCFRN(OM1,ML,NL,DNS,VP,VS,DEP)
C
C *** LOWER TRIAL OMEGA CHECK ***
C
      ICO=1
      IFOM=2
      IF(FR1.EQ.0.0D0)THEN
        CFOM(1)=OM1
        IF(NEMD.EQ.1)GOTO 40
        OM1=2.0D0*DDOM
        FR1=FCFRN(OM1,ML,NL,DNS,VP,VS,DEP)
        ICO=2
        IFOM=3
      ENDIF
C
C *** SEARCH OF CUT-OFF OMEGA OF HIGHER RAYLEIGH MODE ***
C
      DO 100 IOM=IFOM,NDOM
C
C  ** HIGHER TRIAL OMEGA AND EIGEN VALUE EQUATION **
C
        OM2=DFLOAT(IOM)*DDOM
        FR2=FCFRN(OM2,ML,NL,DNS,VP,VS,DEP)
C
C  ** LOWER AND HIGER OMEGA CHECK **
C
        IF(FR1*FR2)10,20,30
C
C   * CASE OF LOWER OM < CUT-OFF OM < HIGHER OM *
C
C     * SEARCH OF CUT-OFF OM BY ROOTN *
C
   10     XEPS=XEPSO
          FEPS=DABS(FR1-FR2)/DFEPS
          ITR=ITMAX
          OML=OM1
          OMR=OM2
          CALL ROOTN(FCFRN,ML,NL,DNS,VP,VS,DEP
     *              ,OML,OMR,XEPS,FEPS,ITR,OMLR,IER)
C
C     * ROOTN ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE RAYCTN :'
            WRITE(*,*)' ROOTN ERROR, IER =',IER
            WRITE(*,*)' XEPS,FEPS,ITMAX =',XEPS,FEPS,ITR
            WRITE(*,*)' OM1,OM2 =',OM1,OM2
            WRITE(*,*)' FR1,FR2 =',FR1,FR2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' OML,OMR =',OML,OMR
              FRL=FCFRN(OML,ML,NL,DNS,VS,VP,DEP)
              FRR=FCFRN(OMR,ML,NL,DNS,VP,VS,DEP)
              WRITE(*,*)' FRL,FRR =',FRL,FRR
              WRITE(*,*)' OMEGA IS CHANGED BY MEAN OMEGA !!!'
              OMLR=(OMR+OML)/2.0D0
              FRLR=FCFRN(OMLR,ML,NL,DNS,VS,VP,DEP)
              WRITE(*,*)' NEW OM, FR =',OMLR,FRLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FR1*FR2 > 0       !!!'
              WRITE(*,*)' STOP  !!!'
              STOP
            ENDIF
          ENDIF
C
C     * ENTER OF CUT-OFF OMEGA *
C
          CFOM(ICO)=OMLR
          IF(NEMD.EQ.ICO)GOTO 40
          ICO=ICO+1
          GOTO 30
C
C   * CASE OF HIGER TRIAL OM = CUT-OFF OM *
C
   20     CFOM(ICO)=OM2
          IF(NEMD.EQ.ICO)GOTO 40
          ICO=ICO+1
          OM2=OM1+DDOM/2.0D0
          FR2=FCFRN(OM2,ML,NL,DNS,VP,VS,DEP)
C
C   * CASE OF LOWER OM < HIGHER OM < CUT-OFF OM *
C
C     * HIGHER OM ==> LOWER OM AND NEXT STEP *
C
   30     OM1=OM2
          FR1=FR2
C
  100 CONTINUE
C
C *** NUMBER OF CUT-OFF OMEGA ***
C
      NCFOM=ICO-1
      RETURN
C
   40 NCFOM=ICO
C
      RETURN
      END
C
C *** ROOT FINDING OF FUNC(X)= 0 BY TWO DIVISION AND LINEAR
C     INTERPOLATION COMBINED METHOD
C *** SOURCE PROGRAM IS $DNTBM (MSL2) AND ARRENGED IN 1988
C
      SUBROUTINE ROOTN(FUNC,ML,NL,DNS,VP,VS,DEP
     *                ,XL,XR,XEPS,FEPS,ITMAX,X,IER)
C
C ** INPUT DATA **
C FUNC  : FUNCTION OF X ( NEED EXTERNAL IN MAIN PROGRAM !!! )
C XL    : LEFT(LOWER LIMIT) X (XL<XS)
C XR    : RIGHT(HIGHER LIMIT) X (XS<XR)
C XEPS  : CRITERION VALUE TO JUGDE CONVERGENCE OF X (XEPS>0)
C FEPS  : CRITERION VALUE TO JUGDE CONVERGENCE OF FUNC(X) (FEPS>0)
C ITMAX : MAXIMUM ITERATION NUMBER
C ** OUTPUT DATA **
C ITMAX : FINAL ITERATION NUMBER IN THE SUBROUTINE
C X     : ROOT (SOLUSION) OF FUNC(X)=0
C IER   : RETURN CODE (IF ERROR, IER.NE.0.  SEE $DNTBM,MSL2)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),VP(ML),VS(ML),DEP(ML)
C
C *** INPUT DATA CHECK ***
C
      IER = 0
      IF(XEPS.LT.0.0D0 .OR. FEPS.LT.0.0D0 .OR. ITMAX.LT.1) GO TO 260
C
      N   = 1
      X   = XL
      FL  = FUNC(X,ML,NL,DNS,VP,VS,DEP)
      IF( FL ) 100,555,100
  100 X   = XR
      FR  =  FUNC(X,ML,NL,DNS,VP,VS,DEP)
      IF( FR  ) 110,555,110
  110 IF( FR*FL ) 130,555,120
  120 X   = (XL+XR)*.5D0
  125 CONTINUE
      IER=2100
      N   = ITMAX
      GO TO 555
C
C *** ROOT FINDING OF FUNC(X,ML,NL,DNS,VP,VS,DEP)= 0 ***
C
  130 DO  250 N=1,ITMAX
C
C   * TRIAL X AND FUNC(X,ML,NL,DNS,VP,VS,DEP) BY TWO DIVISION METHOD *
C
        X   = (XL+XR)*.5D0
        F   = FUNC(X,ML,NL,DNS,VP,VS,DEP)
C
C   * REPLACE OF XL OR XR BY TRIAL X *
C
        IF( F ) 140,555,140
C
  140     IF( F*FL ) 150,555,160
C
  150       XR  = X
            FR  = F
            GO TO 170
  160       XL  = X
            FL  = F
C
C   * TRIAL X AND FUNC(X,ML,NL,DNS,VP,VS,DEP) BY LINEAR INTERPOLATION METHOD *
C
  170   DF  = FL-FR
        DX  = (XR-XL)*FL
        DX  = DX/DF
        X   = XL+DX
        F   = FUNC(X,ML,NL,DNS,VP,VS,DEP)
C
C   * REPLACE OF XL OR XR BY TRIAL X *
C
        IF( F ) 175,555,175
C
  175     IF( FL*F ) 180,190,190
C
  180       XR  = X
            FR  = F
            GO TO 200
  190       XL  = X
            FL  = F
C
C   * CONVERGENCE JUDGMENT OF TRIAL X AND FUNC(X,ML,NL,DNS,VP,VS,DEP) *
C
  200   SX  = XR-XL
        IF( ABS(X).LE.1.0D0 ) GO TO 210
          SX  = SX/X
  210   IF( ABS(SX)-XEPS ) 220,220,250
  220     X   = XR
        IF( ABS(FR)-FEPS ) 555,555,230
  230     X   = XL
        IF( ABS(FL)-FEPS ) 555,555,250
C
  250 CONTINUE
C
C *** ITERATION OVER ERROR ***
C
      IER=3000
      N = ITMAX
      GO TO 555
C
C *** INPUT DATA ERROR ***
C
  260 IER=2000
      N = ITMAX
C
C *** FINAL WORK ***
C
  555 CONTINUE
      ITMAX = N
C
      RETURN
      END
C
C *** CAL. OF RAYLEIGH PHASE WAVE VELOCITY IN ONE LAYER MEDIUM ***
C
      REAL FUNCTION FNCCR*8(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=FNCFR1(SP,Z)
  100 Z=Z-DZ
      F1=FNCFR1(SP,Z)
      AF=F0/F1
      F00=F0
      F0 =F1
      IF(AF.GE.0.D0)GOTO 100
      FZ=Z+DZ*F0/(F0-F00)
C
      FNCCR=VS*DSQRT(FZ)
C
      RETURN
      END
C
C
      REAL FUNCTION FNCFR1*8(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
      FNCFR1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0
C
      RETURN
      END
C
C *** PHASE VELOCITY (CC) OF RAYLEIGH WAVE
C
      SUBROUTINE RAYCCN(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *                 ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C ** INPUT DATA
C DOM    : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM    : TOTAL NUMBER OF CIRCULAR FREQUENCY
C MOM    : NOM IN MAIN PROGRAM DIMENSION
C        : IF MOM > 2**13(=8192), CAN NOT USE !!!
C NDC    : DIVISION NUMVER FROM MIN CC TO MAX CC TO GET ANY CC
C NSMD   : START OUTPUT MODE OF RAYLEIGH WAVE
C NEMD   : END OUTPUT MODE OF RAYLEIGH WAVE
C CFOM(I): CUT-OFF OMEGA OF Ith MODE
C MMD    : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CC(I,J) : RAYLEIGH WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): TOTAL NUMBER OF PHASE VELOCITY AT Ith OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),DEP(ML),DNS(ML)
      COMMON / OMEGA / OM / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
      DIMENSION CC(MOM,0:MMD),CFOM(0:MMD),NMODE(MOM)
      EXTERNAL FNCFR 
C
C *** INCREMENT OF CC (DC) BETWEEN MAX VS AND MIN CC ***
C
      DC=(VSMAX-CCMIN)/DFLOAT(NDC)
C
C *** INPUT OF CONVERGENCE JUDGMENT CRITERIA FOR ROOTN ***
C
      DCEPS=1000.0D0
      CEPSO=DC/DCEPS
      DFEPS=100000.0D0
      ITMAX=100
C     WRITE(*,*)
C     WRITE(*,*)'From SUBROUTINE RAYCCN :'
C     WRITE(*,*)
C     WRITE(*,*)'Convergence Criteria for Phase Velocity :'
C     WRITE(*,1000)CEPSO
C     WRITE(*,1010)DFEPS
C     WRITE(*,1020)ITMAX
C     WRITE(*,1030)
C     READ(*,*)IYES
C     IF(IYES.NE.0)THEN
C       WRITE(*,*)
C       WRITE(*,1040)
C       READ(*,*)CEPSO,DFEPS,ITMAX
C     END IF
C     WRITE(*,*)
C     WRITE(*,*)
C1000 FORMAT(/' CEPS = DC/DCEPS            : CEPS  =',F15.10)
C1010 FORMAT( ' FEPS = DABS(FR1-FR2)/DFEPS : DFEPS =',F15.3)
C1020 FORMAT( ' Maximum Iteration Number   : ITMAX =',I15/)
C1030 FORMAT( ' Are These Criteria OK ? ( yes = 0 ) :',\)
C1040 FORMAT( ' Enter New CEPS, DFEPS, ITMAX :',\)
C
C *** PHASE VELOCITY FOR EACH OMEGA ***
C
      DO 100 IOM=1,NOM
C
C  ** OMEGA AND THE LOWEST USED OMEGA **
C
        OM=DFLOAT(IOM)*DOM
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
C  ** LOWER TRIAL CC AND EIGEN VALUE EQUATION **
C
        CC1=CCMIN-DC
        FR1=FNCFR(CC1,ML,NL,DNS,VP,VS,DEP)
C
C  ** CHECK OF LOWER TRIAL CC **
C
        IFC=2
        ICC=0
        IF(FR1.EQ.0.0D0)THEN
          CC(IOM,0)=CC1
          IF(NEMD.EQ.0)GOTO 40
          CC1=CCMIN
          FR1=FNCFR(CC1,ML,NL,DNS,VP,VS,DEP)
          IFC=3
          ICC=1
        ENDIF
C
C  ** SEARCH OF CORRECT CC FROM MINIMUM VS TO MAXIMUM VS(=VS(NL)) **
C
        DO 110 IC=IFC,NDC+2
C
C   * HIGHER TRIAL CC AND EIGEN VALUE EQUARION *
C
          CC2=CCMIN+DFLOAT(IC-2)*DC
          FR2=FNCFR(CC2,ML,NL,DNS,VP,VS,DEP)
C
C   * CHECK OF TRIAL CC ( LOWER AND HIGER CC ) *
C
          IF(FR1*FR2)10,20,30
C
C     * CASE OF LOWER CC < CORRECT CC < HIGHER CC *
C
C       * SEARCH OF CORRECT CC BY ROOTN *
C
   10     CEPS=CEPSO
          FEPS=DABS(FR1-FR2)/DFEPS
          ITR=ITMAX
          CL=CC1
          CR=CC2
          CALL ROOTN(FNCFR,ML,NL,DNS,VP,VS,DEP
     *              ,CL,CR,CEPS,FEPS,ITR,CCLR,IER)
C
C       * ROOTN ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE RAYCCN :'
            WRITE(*,*)' ROOTN ERROR, IER =',IER
            WRITE(*,*)' OMEGA =',OM
            WRITE(*,2000)CEPS,FEPS,ITR
            WRITE(*,*)' CC1,CC2 =',CC1,CC2
            WRITE(*,*)' FR1,FR2 =',FR1,FR2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' CL,CR =',CL,CR
              FRL=FNCFR(CL,ML,NL,DNS,VP,VS,DEP)
              FRR=FNCFR(CR,ML,NL,DNS,VP,VS,DEP)
              WRITE(*,*)' FRL,FRR =',FRL,FRR
              WRITE(*,*)' CC IS CHANGED BY MEAN CC !!!'
              CCLR=(CR+CL)/2.0D0
              FRLR=FNCFR(CCLR,ML,NL,DNS,VP,VS,DEP)
              WRITE(*,*)' NEW CC, FR =',CCLR,FRLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FR1*FR2 > 0       !!!'
              WRITE(*,*)' STOP  !!!'
              STOP
            ENDIF
          ENDIF
 2000     FORMAT(' CEPS,FEPS,ITRATION =',2D15.7,I5)
C
C       * CORRECT CC BY ROOTN *
C
          CC(IOM,ICC)=CCLR
          IF(NEMD.EQ.ICC)GOTO 40
          ICC=ICC+1
          GOTO 30
C
C    * CASE OF CORRECT CC = HIGHER CC *
C
   20     CC(IOM,ICC)=CC2
          IF(NEMD.EQ.ICC)GOTO 40
          ICC=ICC+1
          CC2=CC1+DC/2.0D0
          FR2=FNCFR(CC2,ML,NL,DNS,VP,VS,DEP)
C
C    * CASE OF LOWER CC < HIGER CC < CORRECT CC *
C
C      * LOWER CC ==> HIGHER CC AND NEXT STEP *
C
   30     CC1=CC2
          FR1=FR2
C
  110   CONTINUE
C
C  ** MODE NUMBER AT EACH OMEGA **
C
        NMODE(IOM)=ICC-1
        GOTO 50
   40   NMODE(IOM)=NEMD
C
C *** IMFORMATION OUTPUT ON SCREEN ***
C
   50   WRITE(*,3000)OM,DOM*DFLOAT(NOM)
 3000   FORMAT('+Now Caluculating C at OMEGA / Final OMEGA =',2F8.3)
C
C *** NEXT OMEGA ***
C
  100 CONTINUE
C
      WRITE(*,*)
C
      RETURN
      END
C
C *** EIGEN VALUE EQUATION OF RAYLEIGH WAVE IN CASE OF CC=VS(NL)
C
      REAL FUNCTION FCFRN*8 (OM,ML,NL,DNS,VP,VS,DEP)
C
C INPUT DATA
C OM     : CIRCULAR FREQUENCY ( OMEGA ) FROM COMMON
C OUTPUT DATA
C FNCFRN : EIGEN VALUE EQUATION OF RAYLEIGH WAVE
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),DEP(ML),DNS(ML)
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4)
     *          ,RAN,B11,B13,RL,RK,RN,RM
C
C *** WAVE NUMBER ***
C
      BN=VS(NL)
      WN=OM/BN
C
C *** INITIALIZATION OF AA MATRIX ***
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J)=(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
C *** AA MATRIX ***
C
      DO 200 IL=NL-1,1,-1
C
        DO 210 I=1,4
          DO 220 J=1,4
            SA(I,J)=AA(I,J)
  220     CONTINUE
  210   CONTINUE
C
C    ** CQ, STQ AND SDQ **
C
        BM=VS(IL)
        CB2=BN*BN/(BM*BM)
C
        RBM=DSQRT(CB2-1.0D0)
        QM=WN*RBM*DEP(IL)
        CQ=DCOS(QM)
        SQ=DSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
C    ** CP, STP AND SDP **
C
        AM=VP(IL)
        CA2=BN*BN/(AM*AM)
        IF(AM-BN)10,20,30
C
   10     RAM=DSQRT(CA2-1.0D0)
          PM=WN*RAM*DEP(IL)
          CP=DCOS(PM)
          SP=DSIN(PM)
          STP=RAM*SP
          SDP=SP/RAM
          GOTO 40
C
   20     CP=1.0D0
          STP=0.0D0
          SDP=WN*DEP(IL)
          GOTO 40
C
   30     RAM=DSQRT(1.0D0-CA2)
          PM=WN*RAM*DEP(IL)
          EX1=DEXP(PM)
          EX2=DEXP(-PM)
          CP=(EX1+EX2)/2.0D0
          EX12=(EX1-EX2)/2.0D0
          STP=-RAM*EX12
          SDP=EX12/RAM
C
C    ** C11-C16, D11-D15 AND SB MATRIX **
C
   40   GM=2.D0/CB2
C
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*BN*BN
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= DCMPLX(C11*CP-C12*CQ,0.0D0)
        SB(1,2)= DCMPLX(0.0D0,C12*SDP+C11*STQ)
        SB(1,3)=-DCMPLX(C15*D11,0.D0)
        SB(1,4)= DCMPLX(0.D0,C15*D12)
        SB(2,1)=-DCMPLX(0.0D0,C11*STP+C12*SDQ)
        SB(2,2)= DCMPLX(-C12*CP+C11*CQ,0.D0)
        SB(2,3)= DCMPLX(0.D0,C15*D13)
        SB(2,4)= SB(1,3)
        SB(3,1)= DCMPLX(C16*C11*C12*D11,0.0D0)
        SB(3,2)= DCMPLX(0.D0,C16*D14)
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= DCMPLX(0.D0,C16*D15)
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
C    ** AA MATRIX **
C
        DO 250 I=1,4
          DO 260 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 270 K=1,4
              AA(I,J)=AA(I,J)+SA(I,K)*SB(K,J)
  270       CONTINUE
  260     CONTINUE
  250   CONTINUE
C
  200 CONTINUE
C
C *** GN,RAN,B11-B24 AND RL-RM ***
C
      GN=2.D0*VS(NL)*VS(NL)/(BN*BN)
      RAN=-DCMPLX(0.D0,DSQRT(1.0D0-BN*BN/(VP(NL)*VP(NL))))
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*BN*BN)
      B13=-RAN*B14
      B21=-B12
      B23=B14
C
      RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
      RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
      RN=B21*AA(1,1)+B23*AA(3,1)
      RM=B21*AA(1,2)+B23*AA(3,2)
C
C *** FNCFRN ***
C
      FCFRN=REAL(RN*RK-RL*RM)
C                     
      RETURN
      END
C
C *** EIGEN VALUE EQUATION OF RAYLEIGH WAVE
C
      REAL FUNCTION FNCFR*8 (CC,ML,NL,DNS,VP,VS,DEP)
C
C INPUT DATA
C OM     : CIRCULAR FREQUENCY ( OMEGA ) FROM COMMON
C OUTPUT DATA
C FUNCFR : EIGEN VALUE EQUATION OF RAYLEIGH WAVE
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),DEP(ML),DNS(ML)
      COMMON / OMEGA /OM
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4)
     *          ,RAN,RBN,B11,B22,B13,B24,RL,RK,RN,RM
C
C *** WAVE NUMBER ***
C
      WN=OM/CC
C
C *** INITIALIZATION OF AA MATRIX ***
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J)=(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
C *** AA MATRIX ***
C
      DO 200 IL=NL-1,1,-1
C
        DO 210 I=1,4
          DO 220 J=1,4
            SA(I,J)=AA(I,J)
  220     CONTINUE
  210   CONTINUE
C
C    ** CQ, STQ AND SDQ **
C
        BM=VS(IL)
        CB2=CC*CC/(BM*BM)
        IF(BM-CC)10,20,30
C
   10     RBM=DSQRT(CB2-1.0D0)
          QM=WN*RBM*DEP(IL)
          CQ=DCOS(QM)
          SQ=DSIN(QM)
          STQ=RBM*SQ
          SDQ=SQ/RBM
          GOTO 40
C
   20     CQ=1.0D0
          STQ=0.0D0
          SDQ=WN*DEP(IL)
          GOTO 40
C
   30     RBM=DSQRT(1.0D0-CB2)
          QM=WN*RBM*DEP(IL)
          EX1=DEXP(QM)
          EX2=DEXP(-QM)
          CQ=(EX1+EX2)/2.0D0
          EX12=(EX1-EX2)/2.0D0
          STQ=-RBM*EX12
          SDQ=EX12/RBM
C
C    ** CP, STP AND SDP **
C
   40   AM=VP(IL)
        CA2=CC*CC/(AM*AM)
        IF(AM-CC)50,60,70
C
   50     RAM=DSQRT(CA2-1.0D0)
          PM=WN*RAM*DEP(IL)
          CP=DCOS(PM)
          SP=DSIN(PM)
          STP=RAM*SP
          SDP=SP/RAM
          GOTO 80
C
   60     CP=1.0D0
          STP=0.0D0
          SDP=WN*DEP(IL)
          GOTO 80
C
   70     RAM=DSQRT(1.0D0-CA2)
          PM=WN*RAM*DEP(IL)
          EX1=DEXP(PM)
          EX2=DEXP(-PM)
          CP=(EX1+EX2)/2.0D0
          EX12=(EX1-EX2)/2.0D0
          STP=-RAM*EX12
          SDP=EX12/RAM
C
C    ** C11-C16, D11-D15 AND SB MATRIX **
C
   80   GM=2.D0/CB2
C
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CC*CC
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= DCMPLX(C11*CP-C12*CQ,0.0D0)
        SB(1,2)= DCMPLX(0.0D0,C12*SDP+C11*STQ)
        SB(1,3)=-DCMPLX(C15*D11,0.D0)
        SB(1,4)= DCMPLX(0.D0,C15*D12)
        SB(2,1)=-DCMPLX(0.0D0,C11*STP+C12*SDQ)
        SB(2,2)= DCMPLX(-C12*CP+C11*CQ,0.D0)
        SB(2,3)= DCMPLX(0.D0,C15*D13)
        SB(2,4)= SB(1,3)
        SB(3,1)= DCMPLX(C16*C11*C12*D11,0.0D0)
        SB(3,2)= DCMPLX(0.D0,C16*D14)
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= DCMPLX(0.D0,C16*D15)
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
C    ** AA MATRIX **
C
        DO 250 I=1,4
          DO 260 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 270 K=1,4
              AA(I,J)=AA(I,J)+SA(I,K)*SB(K,J)
  270       CONTINUE
  260     CONTINUE
  250   CONTINUE
C
  200 CONTINUE
C
C *** GN,RAN,RBN,B11-B24 AND RL-RM ***
C
      GN=2.D0*VS(NL)*VS(NL)/(CC*CC)
      RAN=-DCMPLX(0.D0,DSQRT(1.0D0-CC*CC/(VP(NL)*VP(NL))))
      RBN=-DCMPLX(0.D0,DSQRT(1.0D0-CC*CC/(VS(NL)*VS(NL))))
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*CC*CC)
      B22=GN*RBN
      B13=-RAN*B14
      B21=-B12
      B23=B14
      B24=RBN*B14
C
      RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
      RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
      RN=B21*AA(1,1)+B22*AA(2,1)+B23*AA(3,1)+B24*AA(4,1)
      RM=B21*AA(1,2)+B22*AA(2,2)+B23*AA(3,2)+B24*AA(4,2)
C
C *** FUNCFR ***
C
      FNCFR=REAL(RN*RK-RL*RM)
C                     
      RETURN
      END
C
C *** RAYLEIGH WAVE GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      SUBROUTINE RAYUAR(DOM,NOM,MOM,ML,NL,DNS,VP,VS,THK
     *          ,CC,NMODE,CFOM,NSMD,MMD,UU,ARWZ,ARWX,ARUZ,ARUX)
C
C ** INPUT DATA **
C DOM     : DELTA OMEGA
C NOM     : TOTAL NUMBER OF OMEGA
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C CC(I,J) : PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): OUTPUT HIGHEST MODE AT Ith OMEGA
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE
C NSMD    : START MODE
C MMD     : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA **
C UU(I,J) : GROUP VELOCITY OF Jth MODE AT Ith OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),THK(ML),DNS(ML)
      DIMENSION CC(MOM,0:MMD),NMODE(MOM),UU(MOM,0:MMD),CFOM(0:MOM)
     *         ,ARWZ(MOM,0:MMD),ARWX(MOM,0:MMD)
     *         ,ARUZ(MOM,0:MMD),ARUX(MOM,0:MMD)
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4),AAK(4,4),SAK(4,4),SBK(4,4)
     *          ,AAW(4,4),SAW(4,4),SBW(4,4)
     *          ,RAN,RBN,B11,B22,B13,B24,RL,RK,RN,RM
     *          ,RANK,RBNK,B11K,B22K,B13K,B24K,RLK,RKK,RNK,RMK
     *          ,RANW,RBNW,B11W,B22W,B13W,B24W,RLW,RKW,RNW,RMW
     *          ,AR1(4,4),AS2(4,4),RS2,GS2,HS2,SS2
     *          ,FWZ0,FWX0,FUZ0,FUX0,FWZR,FWXR,FUZR,FUXR
     *          ,CJ11,CJ12,CJ31,CJ32,DNZ,DNX,ONZ,ONX
     *          ,FUZN,FUXN,FNZN,FNXN
      COMMON /SODAT/ SDEP,LS,ODEP,LO
C
C *** OMEGA ***
C
      DO 100 IOM=1,NOM
C
        OM=DOM*DFLOAT(IOM)
        OM2=OM*OM
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
C *** MODE ***
C
        DO 110 MODE=0,NMODE(IOM)

C      ** WAVE NUMBER **
C
          CV=CC(IOM,MODE)
          CV2=CV*CV
          WN=OM/CV
C
C      ** CK AND CW **
C
          CK=-CV/WN
          CW=1.D0/WN
C
C      ** INITIALIZATION OF AA,AAK AND AAW MATRIX **
C
          DO 120 I=1,4
            DO 130 J=1,4
              AA(I,J) =(0.0D0,0.0D0)
              AAK(I,J)=(0.0D0,0.0D0)
              AAW(I,J)=(0.0D0,0.0D0)
  130       CONTINUE
            AA(I,I)=(1.0D0,0.0D0)
  120     CONTINUE
C
C      ** AA MATRIX **
C
          DO 200 IL=NL-1,1,-1
C
            DO 210 I=1,4
              DO 220 J=1,4
                SA(I,J) =AA(I,J)
                SAK(I,J)=AAK(I,J)
                SAW(I,J)=AAW(I,J)
  220         CONTINUE
  210       CONTINUE
C
C         * CQ, STQ AND SDQ *
C
            BM=VS(IL)
            BM2=BM*BM
            CB2=CV2/BM2
            IF(BM-CV)10,20,30
C
   10         RBM=DSQRT(CB2-1.D0)
              QM=WN*RBM*THK(IL)
              CQ=DCOS(QM)
              SQ=DSIN(QM)
              STQ=RBM*SQ
              SDQ=SQ/RBM
              GOTO 40
C
   20         CQ=1.D0
              STQ=0.D0
              SDQ=WN*THK(IL)
C
   30         RBM=DSQRT(1.D0-CB2)
              QM=WN*RBM*THK(IL)
              EX1=DEXP(QM)
              EX2=DEXP(-QM)
              CQ=(EX1+EX2)/2.0D0
              EX12=(EX1-EX2)/2.0D0
              STQ=-RBM*EX12
              SDQ=EX12/RBM
C
C         * CP, STP AND SDP *
C
   40       AM=VP(IL)
            AM2=AM*AM
            CA2=CV2/AM2
            IF(AM-CV)50,60,70
C
   50         RAM=DSQRT(CA2-1.D0)
              PM=WN*RAM*THK(IL)
              CP=DCOS(PM)
              SP=DSIN(PM)
              STP=RAM*SP
              SDP=SP/RAM
              GOTO 80
C
   60         CP=1.D0
              STP=0.D0
              SDP=WN*THK(IL)
C
   70         RAM=DSQRT(1.D0-CA2)
              PM=WN*RAM*THK(IL)
              EX1=DEXP(PM)
              EX2=DEXP(-PM)
              CP=(EX1+EX2)/2.0D0
              EX12=(EX1-EX2)/2.0D0
              STP=-RAM*EX12
              SDP=EX12/RAM
C
C         * CPK,CPW,-,SDQK,SDQW *
C
   80       CVAK=CV*CK/AM2
            CVAW=CV*CW/AM2
            CVBK=CV*CK/BM2
            CVBW=CV*CW/BM2
C
            RB2=CB2-1.D0
            RA2=CA2-1.D0
            CPK=-THK(IL)*(STP+WN*CVAK*SDP)
            CPW=-THK(IL)*WN*CVAW*SDP
            CQK=-THK(IL)*(STQ+WN*CVBK*SDQ)
            CQW=-THK(IL)*WN*CVBW*SDQ
            STPK=CVAK*SDP+THK(IL)*CP*(RA2+WN*CVAK)
            STPW=CVAW*(SDP+THK(IL)*CP*WN)
            STQK=CVBK*SDQ+THK(IL)*CQ*(RB2+WN*CVBK)
            STQW=CVBW*(SDQ+THK(IL)*CQ*WN)
            SDPK=-CVAK*SDP/RA2+THK(IL)*CP*(1.D0+WN*CVAK/RA2)
            SDPW=CVAW*(-SDP+THK(IL)*CP*WN)/RA2
            SDQK=-CVBK*SDQ/RB2+THK(IL)*CQ*(1.D0+WN*CVBK/RB2)
            SDQW=CVBW*(-SDQ+THK(IL)*CQ*WN)/RB2
C
C         * C11-C16, D11-D15 AND SB MATRIX *
C
            GM=2.D0/CB2
C
            C11=GM
            C12=GM-1.D0
            C13=GM*GM
            C14=C12*C12
            C16=DNS(IL)*CV2
            C15=1.D0/C16
C
            D11=CP-CQ
            D12=SDP+STQ
            D13=STP+SDQ
            D14=C14*SDP+C13*STQ
            D15=C13*STP+C14*SDQ
C
            SB(1,1)= DCMPLX(C11*CP-C12*CQ,0.0D0)
            SB(1,2)= DCMPLX(0.0D0,C12*SDP+C11*STQ)
            SB(1,3)=-DCMPLX(C15*D11,0.D0)
            SB(1,4)= DCMPLX(0.D0,C15*D12)
            SB(2,1)=-DCMPLX(0.0D0,C11*STP+C12*SDQ)
            SB(2,2)= DCMPLX(-C12*CP+C11*CQ,0.D0)
            SB(2,3)= DCMPLX(0.D0,C15*D13)
            SB(2,4)= SB(1,3)
            SB(3,1)= DCMPLX(C16*C11*C12*D11,0.0D0)
            SB(3,2)= DCMPLX(0.D0,C16*D14)
            SB(3,3)= SB(2,2)
            SB(3,4)= SB(1,2)
            SB(4,1)= DCMPLX(0.D0,C16*D15)
            SB(4,2)= SB(3,1)
            SB(4,3)= SB(2,1)
            SB(4,4)= SB(1,1)
C
C         * C11K-C16W, D11K-D15W, SBK AND SBW MATRIX *
C
            GMK=-2.D0*GM*CK/CV
            GMW=-2.D0*GM*CW/CV
C
            C11K=GMK
            C11W=GMW
            C12K=GMK
            C12W=GMW
            C13K=2.D0*GM*GMK
            C13W=2.D0*GM*GMW
            C14K=2.D0*C12*C12K
            C14W=2.D0*C12*C12W
            C16K=2.D0*DNS(IL)*CV*CK
            C16W=2.D0*DNS(IL)*CV*CW
            C15K=-C16K/(C16*C16)
            C15W=-C16W/(C16*C16)
C
            D11K=CPK-CQK
            D11W=CPW-CQW
            D12K=SDPK+STQK
            D12W=SDPW+STQW
            D13K=STPK+SDQK
            D13W=STPW+SDQW
            D14K=C14K*SDP+C14*SDPK+C13K*STQ+C13*STQK
            D14W=C14W*SDP+C14*SDPW+C13W*STQ+C13*STQW
            D15K=C13K*STP+C13*STPK+C14K*SDQ+C14*SDQK
            D15W=C13W*STP+C13*STPW+C14W*SDQ+C14*SDQW
C
            SBK(1,1)= DCMPLX(C11K*CP+C11*CPK-C12K*CQ-C12*CQK,0.0D0)
            SBK(1,2)= DCMPLX(0.0D0,C12K*SDP+C12*SDPK+C11K*STQ+C11*STQK)
            SBK(1,3)=-DCMPLX(C15K*D11+C15*D11K,0.D0)
            SBK(1,4)= DCMPLX(0.D0,C15K*D12+C15*D12K)
            SBK(2,1)=-DCMPLX(0.0D0,C11K*STP+C11*STPK+C12K*SDQ+C12*SDQK)
            SBK(2,2)= DCMPLX(-C12K*CP-C12*CPK+C11K*CQ+C11*CQK,0.D0)
            SBK(2,3)= DCMPLX(0.D0,C15K*D13+C15*D13K)
            SBK(2,4)= SBK(1,3)
            SBK(3,1)= DCMPLX((C16K*C11+C16*C11K)*C12*D11,0.0D0)
            SBK(3,1)= SBK(3,1)+DCMPLX(C16*C11*(C12K*D11+C12*D11K),0.0D0)
            SBK(3,2)= DCMPLX(0.D0,C16K*D14+C16*D14K)
            SBK(3,3)= SBK(2,2)
            SBK(3,4)= SBK(1,2)
            SBK(4,1)= DCMPLX(0.D0,C16K*D15+C16*D15K)
            SBK(4,2)= SBK(3,1)
            SBK(4,3)= SBK(2,1)
            SBK(4,4)= SBK(1,1)
C
            SBW(1,1)= DCMPLX(C11W*CP+C11*CPW-C12W*CQ-C12*CQW,0.0D0)
            SBW(1,2)= DCMPLX(0.0D0,C12W*SDP+C12*SDPW+C11W*STQ+C11*STQW)
            SBW(1,3)=-DCMPLX(C15W*D11+C15*D11W,0.D0)
            SBW(1,4)= DCMPLX(0.D0,C15W*D12+C15*D12W)
            SBW(2,1)=-DCMPLX(0.0D0,C11W*STP+C11*STPW+C12W*SDQ+C12*SDQW)
            SBW(2,2)= DCMPLX(-C12W*CP-C12*CPW+C11W*CQ+C11*CQW,0.D0)
            SBW(2,3)= DCMPLX(0.D0,C15W*D13+C15*D13W)
            SBW(2,4)= SBW(1,3)
            SBW(3,1)= DCMPLX((C16W*C11+C16*C11W)*C12*D11,0.0D0)
            SBW(3,1)= SBW(3,1)+DCMPLX(C16*C11*(C12W*D11+C12*D11W),0.0D0)
            SBW(3,2)= DCMPLX(0.D0,C16W*D14+C16*D14W)
            SBW(3,3)= SBW(2,2)
            SBW(3,4)= SBW(1,2)
            SBW(4,1)= DCMPLX(0.D0,C16W*D15+C16*D15W)
            SBW(4,2)= SBW(3,1)
            SBW(4,3)= SBW(2,1)
            SBW(4,4)= SBW(1,1)
            
C
C         * AA MATRIX *
C
            DO 250 I=1,4
              DO 260 J=1,4
                AA(I,J)=(0.D0,0.D0)
                AAK(I,J)=(0.D0,0.D0)
                AAW(I,J)=(0.D0,0.D0)
                DO 270 K=1,4
                  AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
                  AAK(I,J)=AAK(I,J)+SAK(I,K)*SB(K,J)+SA(I,K)*SBK(K,J)
                  AAW(I,J)=AAW(I,J)+SAW(I,K)*SB(K,J)+SA(I,K)*SBW(K,J)
  270           CONTINUE
  260         CONTINUE
  250       CONTINUE
C
  200     CONTINUE
C
C       * GN,RAN,RBN,B11-B24 AND RL-RM ***
C
          GN=2.D0*VS(NL)*VS(NL)/CV2
          RANI=DSQRT(1.0D0-CV2/(VP(NL)*VP(NL)))
          RBNI=DSQRT(1.0D0-CV2/(VS(NL)*VS(NL)))
          RAN=-DCMPLX(0.D0,RANI)
          RBN=-DCMPLX(0.D0,RBNI)
C
          B11=GN*RAN
          B12=GN-1.D0
          B14=1.D0/(DNS(NL)*CV2)
          B22=GN*RBN
          B13=-RAN*B14
          B21=-B12
          B23=B14
          B24=RBN*B14
C
          RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
          RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
          RN=B21*AA(1,1)+B22*AA(2,1)+B23*AA(3,1)+B24*AA(4,1)
          RM=B21*AA(1,2)+B22*AA(2,2)+B23*AA(3,2)+B24*AA(4,2)
C
C       * GNK,RANK,RBNK,B11K-B24K AND RLK-RMK ***
C
          GNK=-2.D0*GN*CK/CV
          RANK=CV*CK/(VP(NL)*VP(NL)*RAN)
          RBNK=CV*CK/(VS(NL)*VS(NL)*RBN)
C
          B11K=GNK*RAN+GN*RANK
          B12K=GNK
          B14K=-2.D0*CK*B14/CV
          B22K=GNK*RBN+GN*RBNK
          B13K=-RANK*B14-RAN*B14K
          B21K=-B12K
          B23K=B14K
          B24K=RBNK*B14+RBN*B14K
C
          RLK=B11K*AA(1,1)+B11*AAK(1,1)+B12K*AA(2,1)+B12*AAK(2,1)
          RLK=RLK+B13K*AA(3,1)+B13*AAK(3,1)+B14K*AA(4,1)+B14*AAK(4,1)
          RKK=B11K*AA(1,2)+B11*AAK(1,2)+B12K*AA(2,2)+B12*AAK(2,2)
          RKK=RKK+B13K*AA(3,2)+B13*AAK(3,2)+B14K*AA(4,2)+B14*AAK(4,2)
          RNK=B21K*AA(1,1)+B21*AAK(1,1)+B22K*AA(2,1)+B22*AAK(2,1)
          RNK=RNK+B23K*AA(3,1)+B23*AAK(3,1)+B24K*AA(4,1)+B24*AAK(4,1)
          RMK=B21K*AA(1,2)+B21*AAK(1,2)+B22K*AA(2,2)+B22*AAK(2,2)
          RMK=RMK+B23K*AA(3,2)+B23*AAK(3,2)+B24K*AA(4,2)+B24*AAK(4,2)
C
C       * GNW,RANW,RBNW,B11W-B24W AND RLW-RMW ***
C
          GNW=-2.D0*GN*CW/CV
          RANW=CV*CW/(VP(NL)*VP(NL)*RAN)
          RBNW=CV*CW/(VS(NL)*VS(NL)*RBN)
C
          B11W=GNW*RAN+GN*RANW
          B12W=GNW
          B14W=-2.D0*CW*B14/CV
          B22W=GNW*RBN+GN*RBNW
          B13W=-RANW*B14-RAN*B14W
          B21W=-B12W
          B23W=B14W
          B24W=RBNW*B14+RBN*B14W
C
          RLW=B11W*AA(1,1)+B11*AAW(1,1)+B12W*AA(2,1)+B12*AAW(2,1)
          RLW=RLW+B13W*AA(3,1)+B13*AAW(3,1)+B14W*AA(4,1)+B14*AAW(4,1)
          RKW=B11W*AA(1,2)+B11*AAW(1,2)+B12W*AA(2,2)+B12*AAW(2,2)
          RKW=RKW+B13W*AA(3,2)+B13*AAW(3,2)+B14W*AA(4,2)+B14*AAW(4,2)
          RNW=B21W*AA(1,1)+B21*AAW(1,1)+B22W*AA(2,1)+B22*AAW(2,1)
          RNW=RNW+B23W*AA(3,1)+B23*AAW(3,1)+B24W*AA(4,1)+B24*AAW(4,1)
          RMW=B21W*AA(1,2)+B21*AAW(1,2)+B22W*AA(2,2)+B22*AAW(2,2)
          RMW=RMW+B23W*AA(3,2)+B23*AAW(3,2)+B24W*AA(4,2)+B24*AAW(4,2)
C
C      ** FRK,FRW, RG AND RH **
C
          FRK=REAL(RNK*RK+RN*RKK-RLK*RM-RL*RMK)
          FRW=REAL(RNW*RK+RN*RKW-RLW*RM-RL*RMW)
C
          RGI=DIMAG(B11*AA(1,3)+B12*AA(2,3)+B13*AA(3,3)+B14*AA(4,3))
          RNR=REAL(RN)
          RLI=DIMAG(RL)
          RHR=REAL(B21*AA(1,3)+B22*AA(2,3)+B23*AA(3,3)+B24*AA(4,3))
C
C      ** UU **
C
          UU(IOM,MODE) =-FRK/FRW
C
C      ** ARWZ - ARUX **
C
          IF(LS.EQ.NL)THEN
C
            EPD=DEXP(-WN*RANI*SDEP)
            ESD=DEXP(-WN*RBNI*SDEP)
            RS2= B14*EPD
            GS2=-RAN*RS2
            HS2= B14*ESD
            SS2= RBN*HS2
C
            FWZ0=GS2*RN-HS2*RL
            FWX0=RS2*RN-SS2*RL
            FUZ0=GS2*RM-HS2*RK
            FUX0=RS2*RM-SS2*RK
C
            IF((LO.LT.LS).OR.((LO.EQ.LS).AND.(ODEP.LE.SDEP)))THEN
C
              ODEP1=ODEP
              LO1=LO
              CALL FNC21(AR1,LO1,ODEP1,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
              FWZR=AR1(2,2)*FWZ0-AR1(2,1)*FUZ0
              FWXR=AR1(2,2)*FWX0-AR1(2,1)*FUX0
              FUZR=AR1(1,1)*FUZ0-AR1(1,2)*FWZ0
              FUXR=AR1(1,1)*FUX0-AR1(1,2)*FWX0
C
            ELSE
C
              CJ11=B11*AA(1,1)+B13*AA(3,1)
              CJ12=B11*AA(1,2)+B13*AA(3,2)
              CJ31=B21*AA(1,1)+B23*AA(3,1)
              CJ32=B21*AA(1,2)+B23*AA(3,2)
C
              DNZ=-(CJ11*FUZ0-CJ12*FWZ0)/RAN
              DNX=-(CJ11*FUX0-CJ12*FWX0)/RAN
              ONZ=+(CJ31*FUZ0-CJ32*FWZ0)/RBN
              ONX=+(CJ31*FUX0-CJ32*FWX0)/RBN
C
              EPZ=DEXP(-WN*RANI*ODEP)
              ESZ=DEXP(-WN*RBNI*ODEP)
C
              FWZR= RAN*DNZ*EPZ+ONZ*ESZ
              FWXR= RAN*DNX*EPZ+ONX*ESZ
              FUZR=-DNZ*EPZ+RBN*ONZ*ESZ
              FUXR=-DNX*EPZ+RBN*ONX*ESZ
C
            END IF
C
          ELSE
C
            SDEP1=SDEP
            LS1=LS
            CALL FNC22(AS2,LS1,SDEP1,OM2,CV2,WN,DNS,
     *                 VP,VS,THK,ML,NL)
C
            GS2=B11*AS2(1,3)+B12*AS2(2,3)+B13*AS2(3,3)+B14*AS2(4,3)
            RS2=B11*AS2(1,4)+B12*AS2(2,4)+B13*AS2(3,4)+B14*AS2(4,4)
            HS2=B21*AS2(1,3)+B22*AS2(2,3)+B23*AS2(3,3)+B24*AS2(4,3)
            SS2=B21*AS2(1,4)+B22*AS2(2,4)+B23*AS2(3,4)+B24*AS2(4,4)
C
            FWZ0=GS2*RN-HS2*RL
            FWX0=RS2*RN-SS2*RL
            FUZ0=GS2*RM-HS2*RK
            FUX0=RS2*RM-SS2*RK
C
            IF(LO.NE.NL)THEN
C
              ODEP1=ODEP
              LO1=LO
              CALL FNC21(AR1,LO1,ODEP1,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
              FWZR=AR1(2,2)*FWZ0-AR1(2,1)*FUZ0
              FWXR=AR1(2,2)*FWX0-AR1(2,1)*FUX0
              FUZR=AR1(1,1)*FUZ0-AR1(1,2)*FWZ0
              FUXR=AR1(1,1)*FUX0-AR1(1,2)*FWX0
C
            ELSE
C
              FUZN=AA(1,1)*FUZ0-AA(1,2)*FWZ0
              FUXN=AA(1,1)*FUX0-AA(1,2)*FWX0
              FNZN=AA(3,1)*FUZ0-AA(3,2)*FWZ0
              FNXN=AA(3,1)*FUX0-AA(3,2)*FWX0
C
              DNZ=B14*FNZN-GN*FUZN
              DNX=B14*FNXN-GN*FUXN
              ONZ=(B21*FUZN+B14*FNZN)/RBN
              ONX=(B21*FUXN+B14*FNXN)/RBN
C
              EPZ=DEXP(-WN*RANI*ODEP)
              ESZ=DEXP(-WN*RBNI*ODEP)
C
              FWZR= RAN*DNZ*EPZ+ONZ*ESZ
              FWXR= RAN*DNX*EPZ+ONX*ESZ
              FUZR=-DNZ*EPZ+RBN*ONZ*ESZ
              FUXR=-DNX*EPZ+RBN*ONX*ESZ
C
            END IF
C
          END IF
C
          ARWZ(IOM,MODE)=-DIMAG(FWZR)/(WN*FRK)
          ARWX(IOM,MODE)=-REAL(FWXR)/(WN*FRK)
          ARUZ(IOM,MODE)= REAL(FUZR)/(WN*FRK)
          ARUX(IOM,MODE)= DIMAG(FUXR)/(WN*FRK)
C
  110   CONTINUE
C
  100  CONTINUE
C                     
      RETURN
      END
C
      SUBROUTINE FNC21(AA,LL,DEP,OM2,CV2,WN,DNS,VP,VS,THK,ML)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML),VP(ML),VS(ML)
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4),CPX
     *          ,CMP,CMS,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
C
      CPX=DCMPLX(0.D0,1.D0)
      WN2=WN*WN
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=LL,1,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-DCMPLX(WN2,0.D0))
        CMS=CDSQRT(WS2-DCMPLX(WN2,0.D0))
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        IF(IL.EQ.LL)PM=CMP*DEP
        CP=CDCOS(PM)
        SP=CDSIN(PM)
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        IF(IL.EQ.LL)QM=CMS*DEP
        CQ=CDCOS(QM)
        SQ=CDSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
      RETURN
      END
C
      SUBROUTINE FNC22(AA,LL,DEP,OM2,CV2,WN,DNS,VP,VS,THK,ML,NL)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(ML),THK(ML),VP(ML),VS(ML)
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4),CPX
     *          ,CMP,CMS,PM,QM,RBM,RAM
     *          ,CQ,SQ,STQ,SDQ,CP,SP,STP,SDP,GM
     *          ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15
C
      CPX=DCMPLX(0.D0,1.D0)
      WN2=WN*WN
C
      DO 100 I=1,4
        DO 110 J=1,4
          AA(I,J) =(0.0D0,0.0D0)
  110   CONTINUE
        AA(I,I)=(1.0D0,0.0D0)
  100 CONTINUE
C
      DO 150 IL=NL-1,LL,-1
C
        DO 160 I=1,4
          DO 170 J=1,4
            SA(I,J) =AA(I,J)
  170     CONTINUE
  160   CONTINUE
C
        BM2=VS(IL)*VS(IL)
        WP2=OM2/(VP(IL)*VP(IL))
        WS2=OM2/BM2
C
        CMP=CDSQRT(WP2-DCMPLX(WN2,0.D0))
        CMS=CDSQRT(WS2-DCMPLX(WN2,0.D0))
        IF(DIMAG(CMP).GT.0.D0)CMP=-CMP
        IF(DIMAG(CMS).GT.0.D0)CMS=-CMS
C
        RAM=CMP/WN
        PM=CMP*THK(IL)
        IF(IL.EQ.LL)PM=CMP*(THK(IL)-DEP)
        CP=CDCOS(PM)
        SP=CDSIN(PM)
        STP=RAM*SP
        SDP=SP/RAM
C
        RBM=CMS/WN
        QM=CMS*THK(IL)
        IF(IL.EQ.LL)QM=CMS*(THK(IL)-DEP)
        CQ=CDCOS(QM)
        SQ=CDSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
C
        GM=2.D0*BM2/CV2
        C11=GM
        C12=GM-1.D0
        C13=GM*GM
        C14=C12*C12
        C16=DNS(IL)*CV2
        C15=1.D0/C16
C
        D11=CP-CQ
        D12=SDP+STQ
        D13=STP+SDQ
        D14=C14*SDP+C13*STQ
        D15=C13*STP+C14*SDQ
C
        SB(1,1)= C11*CP-C12*CQ
        SB(1,2)= (C12*SDP+C11*STQ)*CPX
        SB(1,3)=-C15*D11
        SB(1,4)= C15*D12*CPX
        SB(2,1)=-(C11*STP+C12*SDQ)*CPX
        SB(2,2)=-C12*CP+C11*CQ
        SB(2,3)= C15*D13*CPX
        SB(2,4)= SB(1,3)
        SB(3,1)= C16*C11*C12*D11
        SB(3,2)= C16*D14*CPX
        SB(3,3)= SB(2,2)
        SB(3,4)= SB(1,2)
        SB(4,1)= C16*D15*CPX
        SB(4,2)= SB(3,1)
        SB(4,3)= SB(2,1)
        SB(4,4)= SB(1,1)
C
        DO 180 I=1,4
          DO 190 J=1,4
            AA(I,J)=(0.D0,0.D0)
            DO 195 K=1,4
              AA(I,J) =AA(I,J) +SA(I,K)*SB(K,J)
  195       CONTINUE
  190     CONTINUE
  180   CONTINUE
C
  150 CONTINUE
C
      RETURN
      END
