C
C *** CUT-OFF OMEGAS OF RAYLEIGH WAVE (C=VS(NL)) ***
C
      SUBROUTINE RCUTOM(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 ML    : NL IN MAIN DIMENSION
C NL    : TOTAL NUMBER OF LAYERS
C DNS(I): DENSITY OF Ith LAYER
C VS (I): VS OF Ith LAYER
C DEP(I): DEPTH (THICKNESS) OF Ith LAYER
C NDOM  : DIVISION NUMBER from 0 to MAX OMEGA TO GET CUT-OFF OMEGAS
C NEMD  : FINAL RAYLEIGH MODE
C MMD   : NEMD IN MAIN DIMENSION
C ** OUTPUT DATA **
C CFOM(I): CUT-OFF OMEGAS OF Ith RAYLEIGH MODE
C NCFOM  : NUMBER OF CUT-OFF OMEGAS
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 *** CUT-OFF OMEGA OF 0th MODE ***
C
      ICO=0
      CFOM(ICO)=0.D0
      IF(NEMD.EQ.0)GOTO 40
C
C *** OMEGA INCREMENT FOR ITERATION ***
C
      DDOM=DFLOAT(NOM)*DOM/DFLOAT(NDOM)
C
C *** CRITERIA NUMBERS FOR ITERATION ***
C
      DXEPS=100.0D0
      XEPSO=DDOM/DXEPS
      DFEPS=1000.0D0
      ITMAX=100
C
C *** LOWER TRIAL OMEGA AND ITS EIGEN VALUE EQUATION ***
C
      OM1=DDOM
      FR1=FNCCUT(OM1,ML,NL,DNS,VP,VS,DEP)
C
C *** CHECK OF LOWER TRIAL OMEGA ***
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=FNCCUT(OM1,ML,NL,DNS,VP,VS,DEP)
        ICO=2
        IFOM=3
      ENDIF
C
C *** CUT-OFF OMEGAS OF HIGHER RAYLEIGH MODES ***
C
      DO 100 IOM=IFOM,NDOM
C
C    ** HIGHER TRIAL OMEGA AND ITS EIGEN VALUE EQUATION **
C
        OM2=DFLOAT(IOM)*DDOM
        FR2=FNCCUT(OM2,ML,NL,DNS,VP,VS,DEP)
C
C    ** CHECK OF LOWER AND HIGER OMEGAS **
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(FNCCUT,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=FNCCUT(OML,ML,NL,DNS,VS,VP,DEP)
              FRR=FNCCUT(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=FNCCUT(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       * 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=FNCCUT(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 OMEGAS ***
C
      NCFOM=ICO-1
      RETURN
C
   40 NCFOM=ICO
C
      RETURN
      END
C
C
C
C *** EIGEN VALUE EQUATION OF RAYLEIGH WAVE IN THE CASE OF CC=VS(NL)
C
      REAL FUNCTION FNCCUT*8 (OM,ML,NL,DNS,VP,VS,DEP)
C
C INPUT DATA
C OM     : CIRCULAR FREQUENCY (OMEGA)
C OUTPUT DATA
C FNCCUT : 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
      CV=VS(NL)
      WN=OM/CV
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=CV*CV/(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=CV*CV/(AM*AM)
        CALL RYSC(WN,AM,CV,CA2,DEP(IL),CP,STP,SDP)
C
C    ** C11-C16, D11-D15 AND SB MATRIX **
C
        GM=2.D0/CB2
        CALL RYSB(GM,DNS(IL),CV,CP,CQ,SDP,SDQ,STP,STQ
     *           ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15,SB)
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)/(CV*CV)
      RAN=-DCMPLX(0.D0,DSQRT(1.0D0-CV*CV/(VP(NL)*VP(NL))))
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*CV*CV)
      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 *** FNCCUT (EIGEN VALURE EQUATION) **
C
      FNCCUT=REAL(RN*RK-RL*RM)
C                     
      RETURN
      END
C
C
C
      SUBROUTINE RYSC(WN,BM,CV,CB2,DEP,CQ,STQ,SDQ)
      IMPLICIT REAL*8(A-H,O-Z)
C
      IF(BM-CV)10,20,30
C
   10 RBM=DSQRT(CB2-1.D0)
      QM=WN*RBM*DEP
      CQ=DCOS(QM)
      SQ=DSIN(QM)
      STQ=RBM*SQ
      SDQ=SQ/RBM
      GOTO 40
C
   20 CQ=1.D0
      STQ=0.D0
      SDQ=WN*DEP
      GOTO 40
C
   30 RBM=DSQRT(1.D0-CB2)
      QM=WN*RBM*DEP
      EX1=DEXP(QM)
      EX2=DEXP(-QM)
      CQ=(EX1+EX2)/2.0D0
      EX12=(EX1-EX2)/2.0D0
      STQ=-RBM*EX12
      SDQ=EX12/RBM
C
   40 RETURN
      END
C
C
C
      SUBROUTINE RYSB(GM,DNS,CV,CP,CQ,SDP,SDQ,STP,STQ
     *               ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15,SB)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SB(4,4)
C
      C11=GM
      C12=GM-1.D0
      C13=GM*GM
      C14=C12*C12
      C16=DNS*CV*CV
      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
      RETURN
      END
C
C
C
C
C *** ROOT OF FUNC(X)= 0 USING TWO DIVISION AND LINEAR
C     INTERPOLATION COMBINED METHOD
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 TO JUGDE CONVERGENCE OF X (XEPS>0)
C FEPS  : CRITERION TO JUGDE CONVERGENCE OF FUNC(X) (FEPS>0)
C ITMAX : MAXIMUM ITERATION NUMBER
C ** OUTPUT DATA **
C ITMAX : FINAL ITERATION NUMBER IN THIS 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 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 *** ERROR 1 (ITERATION OVER) ***
C
      IER=3000
      N = ITMAX
      GO TO 555
C
C *** ERROR 2 (INPUT DATA) ***
C
  260 IER=2000
      N = ITMAX
C
C *** ITERATION NUMBER ***
C
  555 CONTINUE
      ITMAX = N
C
      RETURN
      END
C
C
C
C *** PHASE VELOCITY (CC) OF RAYLEIGH WAVE
C
      SUBROUTINE RPHASE(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 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 FNCFRN 
C
      WRITE(*,*)
C
C *** INCREMENT OF CC (DC) BETWEEN MAX VS AND MIN CC FOR ITERATION ***
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
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=FNCFRN(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=FNCFRN(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=FNCFRN(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(FNCFRN,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=FNCFRN(CL,ML,NL,DNS,VP,VS,DEP)
              FRR=FNCFRN(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=FNCFRN(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=FNCFRN(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
C
C *** EIGEN VALUE EQUATION OF RAYLEIGH WAVE IN A MULTI-LAYERED MEDIUM ***
C
      REAL FUNCTION FNCFRN*8 (CV,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/CV
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 & SDQ AND CP, STP & SDP **
C
        BM=VS(IL)
        CB2=CV*CV/(BM*BM)
        CALL RYSC(WN,BM,CV,CB2,DEP(IL),CQ,STQ,SDQ)
C
        AM=VP(IL)
        CA2=CV*CV/(AM*AM)
        CALL RYSC(WN,AM,CV,CA2,DEP(IL),CP,STP,SDP)
C
C    ** C11-C16, D11-D15 AND SB MATRIX **
C
        GM=2.D0/CB2
        CALL RYSB(GM,DNS(IL),CV,CP,CQ,SDP,SDQ,STP,STQ
     *           ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15,SB)
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)/(CV*CV)
      RAN=-DCMPLX(0.D0,DSQRT(1.0D0-CV*CV/(VP(NL)*VP(NL))))
      RBN=-DCMPLX(0.D0,DSQRT(1.0D0-CV*CV/(VS(NL)*VS(NL))))
C
      B11=GN*RAN
      B12=GN-1.D0
      B14=1.D0/(DNS(NL)*CV*CV)
      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 *** FNCFRN ***
C
      FNCFRN=REAL(RN*RK-RL*RM)
C                     
      RETURN
      END
