C
C *** CUT-OFF OMEGA OF LOVE WAVE (C=VS(NL))
C
      SUBROUTINE LOVECT(DOM,NOM,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 : HIGHEST OUTPUT MODE OF LOVE WAVE
C MMD  : NMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CFOM(I): CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM  : NUMBER OF CUT-OFF OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD)
      EXTERNAL FNCA21
C
C *** 0 th MODE CUT-OFF OMEGA ***
C
      ICO=0
      CFOM(ICO)=0.0D0
      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 D2ROOT ***
C
      DXEPS=100.0D0
      XEPSO=DDOM/DXEPS
      DFEPS=1000.0D0
      ITMAX=100
      WRITE(*,*)
      WRITE(*,*)'From SUBROUTINE LOVECT :'
      WRITE(*,*)
      WRITE(*,*)'Convergence Criteria for Cut-Off Omega :'
      WRITE(*,1000)XEPSO
      WRITE(*,1010)DFEPS
      WRITE(*,1020)ITMAX
      WRITE(*,1030)
      READ(*,*)IYES
      IF(IYES.NE.0)THEN
        WRITE(*,*)
        WRITE(*,1040)
        READ(*,*)XEPSO,DFEPS,ITMAX
      END IF
      WRITE(*,*)
 1000 FORMAT(/' XEPS = DDOM/DXEPS           : XEPS =',F15.10)
 1010 FORMAT( ' FEPS = DABS(FL1-FL2)/DFEPS : DFEPS =',F15.3)
 1020 FORMAT( ' Maximum Iteration Number   : ITMAX =',I15/)
 1030 FORMAT( ' Are These Criteria OK ? ( yes = 0 ) :',\)
 1040 FORMAT( ' Enter New DDOM, DXEPS, DFEPS, ITMAX :',\)
C
C *** LOWER TRIAL OMEGA AND EIGEN VALUE EQUATION ***
C
      OM1=DDOM
      FL1=FNCA21(OM1)
C
C *** LOWER TRIAL OMEGA CHECK ***
C
      ICO=1
      IFOM=2
      IF(FL1.EQ.0.0D0)THEN
        CFOM(1)=OM1
        IF(NEMD.EQ.1)GOTO 40
        OM1=2.0D0*DDOM
        FL1=FNCA21(OM1)
        ICO=2
        IFOM=3
      ENDIF
C
C *** SEARCH OF CUT-OFF OMEGA OF HIGHER LOVE MODE ***
C
      DO 100 IOM=IFOM,NDOM
C
C  ** HIGHER TRIAL OMEGA AND EIGEN VALUE EQUATION **
C
        OM2=DFLOAT(IOM)*DDOM
        FL2=FNCA21(OM2)
C
C  ** LOWER AND HIGER OMEGA CHECK **
C
        IF(FL1*FL2)10,20,30
C
C   * CASE OF LOWER OM < CUT-OFF OM < HIGHER OM *
C
C     * SEARCH OF CUT-OFF OM BY D2ROOT *
C
   10     XEPS=XEPSO
          FEPS=DABS(FL1-FL2)/DFEPS
          ITR=ITMAX
          OML=OM1
          OMR=OM2
          CALL D2ROOT(FNCA21,OML,OMR,XEPS,FEPS,ITR,OMLR,IER)
C
C     * D2ROOT ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE LOVECT :'
            WRITE(*,*)' D2ROOT ERROR, IER =',IER
            WRITE(*,*)' XEPS,FEPS,ITMAX =',XEPS,FEPS,ITR
            WRITE(*,*)' OM1,OM2 =',OM1,OM2
            WRITE(*,*)' FL1,FL2 =',FL1,FL2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' OML,OMR =',OML,OMR
              FLL=FNCA21(OML)
              FLR=FNCA21(OMR)
              WRITE(*,*)' FLL,FLR =',FLL,FLR
              WRITE(*,*)' OMEGA IS CHANGED BY MEAN OMEGA !!!'
              OMLR=(OMR+OML)/2.0D0
              FLLR=FNCA21(OMLR)
              WRITE(*,*)' NEW OM, FL =',OMLR,FLLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FL1*FL2 > 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
          FL2=FNCA21(OM2)
C
C   * CASE OF LOWER OM < HIGHER OM < CUT-OFF OM *
C
C     * HIGHER OM ==> LOWER OM AND NEXT STEP *
C
   30     OM1=OM2
          FL1=FL2
C
  100 CONTINUE
C
C *** NUMBER OF CUT-OFF OMEGA ***
C
      NCFOM=ICO-1
      RETURN
C
   40 NCFOM=ICO
C
      RETURN
      END
C
C
C
C *** EIGEN VALUE EQUATION OF LOVE WAVE, IN CASE OF C=VS(NL)
C
      REAL FUNCTION FNCA21*8 (OM)
C
C INPUT DATA
C OM     : CIRCULAR FREQUENCY ( OMEGA )
C OUTPUT DATA
C FNCA21 : EIGEN VALUE EQUATION OF LOVE WAVE, IN CASE OF C=VS(NL)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON / LAYER / NL,VS(35),DEP(35),SM(35)
      COMPLEX*16 R11,R12,R21,R22,SA11,SA12,SA21,SA22
     *         ,SB11,SB12,SB21,SB22
C
      BN=VS(NL)
      WN=OM/BN
C
      R11=(1.0D0,0.0D0)
      R12=(0.0D0,0.0D0)
      R21=(0.0D0,0.0D0)
      R22=(1.0D0,0.0D0)
C
      DO 100 ML=1,NL-1
C
        SA11=R11
        SA12=R12
        SA21=R21
        SA22=R22
C
        BM=VS(ML)
        RBM=DSQRT(BN*BN/(BM*BM)-1.0D0)
        QM=WN*RBM*DEP(ML)
        CQ=DCOS(QM)
        SQ=DSIN(QM)
C
        SB11=DCMPLX(CQ,0.0D0)
        SB12=DCMPLX(0.0D0,SQ/(SM(ML)*RBM))
        SB21=DCMPLX(0.0D0,SQ*SM(ML)*RBM)
        SB22=DCMPLX(CQ,0.0D0)
C
        R11=SB11*SA11+SB12*SA21
        R12=SB11*SA12+SB12*SA22
        R21=SB21*SA11+SB22*SA21
        R22=SB21*SA12+SB22*SA22
C
  100 CONTINUE
C
      FNCA21=IMAG(R21)
C
      RETURN
      END
C
C
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 D2ROOT(FUNC,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)
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)
      IF( FL ) 100,555,100
  100 X   = XR
      FR  =  FUNC(X)
      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)= 0 ***
C
  130 DO  250 N=1,ITMAX
C
C   * TRIAL X AND FUNC(X) BY TWO DIVISION METHOD *
C
        X   = (XL+XR)*.5D0
        F   = FUNC(X)
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) BY LINEAR INTERPOLATION METHOD *
C
  170   DF  = FL-FR
        DX  = (XR-XL)*FL
        DX  = DX/DF
        X   = XL+DX
        F   = FUNC(X)
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) *
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
C
C *** PHASE VELOCITY (CC) OF LOVE WAVE
C
      SUBROUTINE LOVECC(DOM,NOM,MOM,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 VS TO MAX VS TO GET CC
C NSMD   : START OUTPUT MODE OF LOVE WAVE
C NEMD   : END OUTPUT MODE OF LOVE WAVE
C CFOM(I): CUT-OFF OMEGA OF Ith MODE
C MMD    : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA
C CC(I,J): LOVE 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)
      COMMON / LAYER / NL,VS(35),DEP(35),SM(35) / OMEGA / OM
     *       / MINMAX / VSMIN, VSMAX
      DIMENSION CC(MOM,0:MMD),CFOM(0:MMD),NMODE(MOM)
      EXTERNAL FUNCFL
C
C *** INCREMENT OF CC (DC) BETWEEN VS(NL)(= MAX VS) AND MINIMUM VS ***
C
      BN=VS(NL)
      DC=(BN-VSMIN)/DFLOAT(NDC)
C
C *** INPUT OF CONVERGENCE JUDGMENT CRITERIA FOR D2ROOT ***
C
      DCEPS=1000.0D0
      CEPSO=DC/DCEPS
      DFEPS=100000.0D0
      ITMAX=100
      WRITE(*,*)
      WRITE(*,*)'From SUBROUTINE LOVECC :'
      WRITE(*,*)
      WRITE(*,*)'Convergence Criteria for Phase Velocity :'
      WRITE(*,1000)CEPSO
      WRITE(*,1010)DFEPS
      WRITE(*,1020)ITMAX
      WRITE(*,1030)
      READ(*,*)IYES
      IF(IYES.NE.0)THEN
        WRITE(*,*)
        WRITE(*,1040)
        READ(*,*)CEPSO,DFEPS,ITMAX
      END IF
      WRITE(*,*)
      WRITE(*,*)
 1000 FORMAT(/' CEPS = DC/DCEPS            : CEPS  =',F15.10)
 1010 FORMAT( ' FEPS = DABS(FL1-FL2)/DFEPS : DFEPS =',F15.3)
 1020 FORMAT( ' Maximum Iteration Number   : ITMAX =',I15/)
 1030 FORMAT( ' Are These Criteria OK ? ( yes = 0 ) :',\)
 1040 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=VSMIN-DC
        FL1=FUNCFL(CC1)
C
C  ** CHECK OF LOWER TRIAL CC **
C
        IFC=2
        ICC=0
        IF(FL1.EQ.0.0D0)THEN
          CC(IOM,0)=CC1
          IF(NEMD.EQ.0)GOTO 40
          CC1=VSMIN
          FL1=FUNCFL(CC1)
          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=VSMIN+DFLOAT(IC-2)*DC
          FL2=FUNCFL(CC2)
C
C   * CHECK OF TRIAL CC ( LOWER AND HIGER CC ) *
C
          IF(FL1*FL2)10,20,30
C
C     * CASE OF LOWER CC < CORRECT CC < HIGHER CC *
C
C       * SEARCH OF CORRECT CC BY D2ROOT *
C
   10     CEPS=CEPSO
          FEPS=DABS(FL1-FL2)/DFEPS
          ITR=ITMAX
          CL=CC1
          CR=CC2
          CALL D2ROOT(FUNCFL,CL,CR,CEPS,FEPS,ITR,CCLR,IER)
C
C       * D2ROOT ERROR OR NOT CONVERGENCE WORK *
C
          IF(IER.NE.0)THEN
            WRITE(*,*)
            WRITE(*,*)' FROM SUBROUTINE LOVECC :'
            WRITE(*,*)' D2ROOT ERROR, IER =',IER
            WRITE(*,*)' OMEGA =',OM
            WRITE(*,2000)CEPS,FEPS,ITR
            WRITE(*,*)' CC1,CC2 =',CC1,CC2
            WRITE(*,*)' FL1,FL2 =',FL1,FL2
            IF(IER.EQ.3000)THEN
              WRITE(*,*)' NOT CONVERGENCE (IER=3000) !!!'
              WRITE(*,*)' CL,CR =',CL,CR
              FLL=FUNCFL(CL)
              FLR=FUNCFL(CR)
              WRITE(*,*)' FLL,FLR =',FLL,FLR
              WRITE(*,*)' CC IS CHANGED BY MEAN CC !!!'
              CCLR=(CR+CL)/2.0D0
              FLLR=FUNCFL(CCLR)
              WRITE(*,*)' NEW CC, FL =',CCLR,FLLR
            ELSE
              WRITE(*,*)' IER=2000 INPUT DATA ERROR  !!!'
              WRITE(*,*)' IER=2100 FL1*FL2 > 0       !!!'
              WRITE(*,*)' STOP  !!!'
              STOP
            ENDIF
          ENDIF
 2000     FORMAT(' CEPS,FEPS,ITRATION =',2D15.7,I5)
C
C       * CORRECT CC BY D2ROOT *
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
          FL2=FUNCFL(CC2)
C
C    * CASE OF LOWER CC < HIGER CC < CORRECT CC *
C
C      * LOWER CC ==> HIGHER CC AND NEXT STEP *
C
   30     CC1=CC2
          FL1=FL2
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 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 LOVE WAVE
C
      REAL FUNCTION FUNCFL*8 (CC)
C
C INPUT DATA
C OM     : CIRCULAR FREQUENCY ( OMEGA ) FROM COMMON
C CC     : PHASE VELOCITY OF LOVE WAVE
C OUTPUT DATA
C FUNCFL : EIGEN VALUE EQUATION OF LOVE WAVE
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON / LAYER / NL,VS(35),DEP(35),SM(35) / OMEGA /OM
      COMPLEX*16 R11,R12,R21,R22,SA11,SA12,SA21,SA22
     *          ,SB11,SB12,SB21,SB22
C
      WN=OM/CC
C
      R11=(1.0D0,0.0D0)
      R12=(0.0D0,0.0D0)
      R21=(0.0D0,0.0D0)
      R22=(1.0D0,0.0D0)
C
      DO 100 ML=1,NL-1
C
        SA11=R11
        SA12=R12
        SA21=R21
        SA22=R22
C
        BM=VS(ML)
        CB2=CC*CC/(BM*BM)
        IF(BM-CC)10,20,30
C
   10     RBM=DSQRT(CB2-1.0D0)
          QM=WN*RBM*DEP(ML)
          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(ML)
          GOTO 40
C
   30     RBM=DSQRT(1.0D0-CB2)
          QM=WN*RBM*DEP(ML)
          EX1=DEXP(QM)
          EX2=DEXP(-QM)
          CQ=(EX1+EX2)/2.0D0
          EX12=(EX1-EX2)/2.0D0
          STQ=-RBM*EX12
          SDQ=EX12/RBM
C
   40   SB11=DCMPLX(CQ,0.0D0)
        SB12=DCMPLX(0.0D0,SDQ/SM(ML))
        SB21=DCMPLX(0.0D0,STQ*SM(ML))
        SB22=DCMPLX(CQ,0.0D0)
C
        R11=SB11*SA11+SB12*SA21
        R12=SB11*SA12+SB12*SA22
        R21=SB21*SA11+SB22*SA21
        R22=SB21*SA12+SB22*SA22
C
  100 CONTINUE
C
      RBN=-DSQRT(1.0D0-CC*CC/(VS(NL)*VS(NL)))
      FUNCFL=IMAG(R21)+SM(NL)*RBN*REAL(R11)
C
      RETURN
      END
C
C
C
C *** LOVE WAVE GROUP VELOCITY AND MEDIUM RESPONSE ***
C
      SUBROUTINE LVUA(DOM,NOM,MOM,CC,NMODE,CFOM,NSMD,MMD,UU,AL)
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 LOVE WAVE 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 LOVE WAVE AT Ith OMEGA
C AL(I,J) : AMPLITUDE SPEC (i*AL) OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : TO OBTAIN AN ACTUAL AL, MULTIPLY (-i) BY AL
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON / LAYER / NL,VS(35),DEP(35),SM(35)
      DIMENSION CC(MOM,0:MMD),NMODE(MOM),UU(MOM,0:MMD),AL(MOM,0:MMD)
     *         ,CFOM(0:MOM)
      COMPLEX*16 R11,R12,R21,R22,SA11,SA12,SA21,SA22
     *         ,SB11,SB12,SB21,SB22
     *         ,RK11,RK12,RK21,RK22,SAK11,SAK12,SAK21,SAK22
     *         ,SBK11,SBK12,SBK21,SBK22
     *         ,RW11,RW12,RW21,RW22,SAW11,SAW12,SAW21,SAW22
     *         ,SBW11,SBW12,SBW21,SBW22
     *         ,RBN
C
      DO 100 IOM=1,NOM
C
        OM=DOM*DFLOAT(IOM)
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
        DO 110 MODE=0,NMODE(IOM)
C
          CV=CC(IOM,MODE)
          WN=OM/CV
C
          R11=(1.0D0,0.0D0)
          R12=(0.0D0,0.0D0)
          R21=(0.0D0,0.0D0)
          R22=(1.0D0,0.0D0)
C
          RK11=(0.0D0,0.0D0)
          RK12=(0.0D0,0.0D0)
          RK21=(0.0D0,0.0D0)
          RK22=(0.0D0,0.0D0)
C
          RW11=(0.0D0,0.0D0)
          RW12=(0.0D0,0.0D0)
          RW21=(0.0D0,0.0D0)
          RW22=(0.0D0,0.0D0)
C
          DO 120 ML=1,NL-1
C
            SA11=R11
            SA12=R12
            SA21=R21
            SA22=R22
C
            SAK11=RK11
            SAK12=RK12
            SAK21=RK21
            SAK22=RK22
C
            SAW11=RW11
            SAW12=RW12
            SAW21=RW21
            SAW22=RW22
C
            BM=VS(ML)
            CB2=CV*CV/(BM*BM)
            IF(BM-CV)10,20,30
C
   10         RBM=DSQRT(CB2-1.0D0)
              QM=WN*RBM*DEP(ML)
              CQ=DCOS(QM)
              SQ=DSIN(QM)
              STQ=RBM*SQ
              SDQ=SQ/RBM
              RBM2=RBM*RBM
              GOTO 40
C
   20         WRITE(*,*)' FROM SUBROUTINE LOVEUA :'
              WRITE(*,*)' RBM=0 THEN STOP !!!!'
              WRITE(*,*)' DATA : OM,ML,CV=',OM,ML,CV
              STOP
C
   30         RBM=DSQRT(1.0D0-CB2)
              QM=WN*RBM*DEP(ML)
              EX1=DEXP(QM)
              EX2=DEXP(-QM)
              CQ=(EX1+EX2)/2.0D0
              EX12=(EX1-EX2)/2.0D0
              STQ=-RBM*EX12
              SDQ=EX12/RBM
              RBM2=-RBM*RBM
C
   40       DSG=DEP(ML)*SDQ
            DC=DEP(ML)*CQ
            SGK=SDQ/WN
            SG2=SM(ML)*RBM2
C
            SB11=DCMPLX(CQ,0.0D0)
            SB12=DCMPLX(0.0D0,SDQ/SM(ML))
            SB21=DCMPLX(0.0D0,STQ*SM(ML))
            SB22=DCMPLX(CQ,0.0D0)
C
            SBK11=DCMPLX(DSG,0.0D0)
            SBK12=DCMPLX(0.0D0,(CB2*SGK-DC)/SG2)
            SBK21=DCMPLX(0.0D0,-SM(ML)*(CB2*SGK+DC))
            SBK22=DCMPLX(DSG,0.0D0)
C
            CVM2=CV/(BM*BM)
            SBW11=CVM2*DCMPLX(-DSG,0.0D0)
            SBW12=CVM2*DCMPLX(0.0D0,(-SGK+DC)/SG2)
            SBW21=CVM2*DCMPLX(0.0D0,SM(ML)*(SGK+DC))
            SBW22=CVM2*DCMPLX(-DSG,0.0D0)
C
            R11=SB11*SA11+SB12*SA21
            R12=SB11*SA12+SB12*SA22
            R21=SB21*SA11+SB22*SA21
            R22=SB21*SA12+SB22*SA22
C
            RK11=SBK11*SA11+SBK12*SA21+SB11*SAK11+SB12*SAK21
            RK12=SBK11*SA12+SBK12*SA22+SB11*SAK12+SB12*SAK22
            RK21=SBK21*SA11+SBK22*SA21+SB21*SAK11+SB22*SAK21
            RK22=SBK21*SA12+SBK22*SA22+SB21*SAK12+SB22*SAK22
C
            RW11=SBW11*SA11+SBW12*SA21+SB11*SAW11+SB12*SAW21
            RW12=SBW11*SA12+SBW12*SA22+SB11*SAW12+SB12*SAW22
            RW21=SBW21*SA11+SBW22*SA21+SB21*SAW11+SB22*SAW21
            RW22=SBW21*SA12+SBW22*SA22+SB21*SAW12+SB22*SAW22
C
  120     CONTINUE
C
          CB2=CV*CV/(VS(NL)*VS(NL))
          RBN=DCMPLX(0.0D0,-DSQRT(1.0D0-CB2))
C
          FLK=IMAG(RK21+SM(NL)*(RBN*RK11-CB2*R11/(RBN*WN)))
          FLW=IMAG(RW21+SM(NL)*(RBN*RW11+CB2*R11/(CV*RBN*WN)))
          UU(IOM,MODE)=-FLK/FLW
          AL(IOM,MODE)=1.0D0/(R11*FLK)
C
  110   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
