C
C *** LOVE WAVE GROUP VELOCITY AND MEDIUM RESPONSE ***
C
      SUBROUTINE LVUALN(DOM,NOM,MOM,ML,NL,DNS,VS,DEP,CC,NMODE,CFOM
     *                 ,NSMD,MMD,MLS,SDEP,MLO,ODEP,UU,ALM)
C
C ** INPUT DATA **
C DOM     : DELTA OMEGA
C NOM     : TOTAL NUMBER OF OMEGA
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C ML    : NL IN MAIN 
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 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 MLS     : NUMBER OF LAYER WHERE A SOURCE IS LOCATED (FROM TOP LAYER)
C SDEP    : SOURCE DEPTH IN THE SOURCE LAYER (FROM UPPER BOUNDARY)
C MLO     : NUMBER OF LAYER WHERE A OBS. IS LOCATED (FROM TOP LAYER)
C ODEP    : SOURCE DEPTH IN THE OBS. LAYER (FROM UPPER BOUNDARY)
C ** OUTPUT DATA **
C UU(I,J) : GROUP VELOCITY OF Jth MODE LOVE WAVE AT Ith OMEGA
C ALM(I,J): NEW AMPLITUDE SPEC OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : BY REGAN AND HARKRIDER (1989)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CC(MOM,0:MMD),NMODE(MOM),UU(MOM,0:MMD),ALM(MOM,0:MMD)
     *         ,CFOM(0:MMD),DNS(ML),VS(ML),DEP(ML)
      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 IL=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(IL)
            SM=DNS(IL)*BM*BM
            CB2=CV*CV/(BM*BM)
            CALL LVSC(WN,BM,CV,CB2,DEP(IL),CQ,STQ,SDQ,RBM2)
C
            DSG=DEP(IL)*SDQ
            DC=DEP(IL)*CQ
            SGK=SDQ/WN
            SG2=SM*RBM2
C
            SB11=DCMPLX(CQ,0.0D0)
            SB12=DCMPLX(0.0D0,SDQ/SM)
            SB21=DCMPLX(0.0D0,STQ*SM)
            SB22=DCMPLX(CQ,0.0D0)
C
            SBK11=DCMPLX(DSG,0.0D0)
            SBK12=DCMPLX(0.0D0,(CB2*SGK-DC)/SG2)
            SBK21=DCMPLX(0.0D0,-SM*(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*(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
            IF(IL.EQ.MLS)THEN
              CALL LVSC(WN,BM,CV,CB2,SDEP,CQ,STQ,SDQ,RBM2)
              SB11=DCMPLX(CQ,0.0D0)
              SB12=DCMPLX(0.0D0,SDQ/SM)
              AS11=REAL(SB11*SA11+SB12*SA21)
            END IF
C
            IF(IL.EQ.MLO)THEN
              CALL LVSC(WN,BM,CV,CB2,ODEP,CQ,STQ,SDQ,RBM2)
              SB11=DCMPLX(CQ,0.0D0)
              SB12=DCMPLX(0.0D0,SDQ/SM)
              AO11=REAL(SB11*SA11+SB12*SA21)
            END IF
C
  120     CONTINUE
C
          BM=VS(NL)
          CB2=CV*CV/(BM*BM)
          RRBN=DSQRT(1.D0-CB2)
          RBN=DCMPLX(0.D0,-RRBN)
          SM=DNS(NL)*BM*BM
C
          IF(NL.EQ.MLS)THEN
            CALL LVSC(WN,BM,CV,CB2,SDEP,CQ,STQ,SDQ,RBM2)
            SB11=DCMPLX(CQ,0.0D0)
            SB12=DCMPLX(0.0D0,SDQ/SM)
            AS11=REAL(SB11*R11+SB12*R21)
          END IF
C
          IF(NL.EQ.MLO)THEN
            CALL LVSC(WN,BM,CV,CB2,ODEP,CQ,STQ,SDQ,RBM2)
            SB11=DCMPLX(CQ,0.0D0)
            SB12=DCMPLX(0.0D0,SDQ/SM)
            AO11=REAL(SB11*R11+SB12*R21)
          END IF
C
          FLK= -IMAG(RK21+SM*(RBN*RK11-CB2*R11/(RBN*WN)))
          FLW= -IMAG(RW21+SM*(RBN*RW11+CB2*R11/(CV*RBN*WN)))
          UU(IOM,MODE)=-FLK/FLW
          AL11=REAL(R11)
          ALM(IOM,MODE)=AS11*AO11/(AL11*FLK)
C
  110   CONTINUE
C
  100 CONTINUE
C
      RETURN
      END
C
C
C
      SUBROUTINE LVSC(WN,BM,CV,CB2,DEP,CQ,STQ,SDQ,RBM2)
      IMPLICIT REAL*8(A-H,O-Z)
C
      IF(BM-CV)10,20,30
C
   10   RBM=DSQRT(CB2-1.0D0)
        QM=WN*RBM*DEP
        CQ=DCOS(QM)
        SQ=DSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
        RBM2=RBM*RBM
        GOTO 40
C
   20   WRITE(*,*)' FROM SUBROUTINE LVSC, LVUAM :'
        WRITE(*,*)' RBM=0 THEN STOP !!!!'
        WRITE(*,*)' DATA : OM,IL,CV=',OM,IL,CV
        STOP
C
   30   RBM=DSQRT(1.0D0-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
        RBM2=-RBM*RBM
C
   40 RETURN
C
      END
