C
C *** LOVE WAVE GROUP VELOCITY AND MEDIUM RESPONSE ***
C
      SUBROUTINE LVUAN(DOM,NOM,MOM,ML,NL,DNS,VS,DEP,MLS,MLE
     *                ,DH,NH,CC,NMODE,CFOM,NSMD,MMD,UU,VD)
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 MLS     : OUTPUT START LAYER FOR DAL  
C MLE     : OUTPUT END LAYER FOR DAL  
C DH      : THICKNESS FOR INTEGRATION IN THE LOWEST LAYER
C NH      : NUMBER OF DH FOR INTEGRATION
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 VD(I,J) : AMPLITUDE SPEC (i*AL) OF LOVE WAVE GENERATED ON X=0 
C           OF Jth MODE AT Ith OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CC(MOM,0:MMD),NMODE(MOM),UU(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,DS1,DS2,FL1,FL2,VF,VD(MOM,0:MMD)
C
      
      SMN=DNS(NL)*VS(NL)*VS(NL)
      CALL H2DIN(XF,ZF,W,DA,D,VRR,RT)
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
          CB2N=CV*CV/(VS(NL)*VS(NL))
          RRBN=DSQRT(1.0D0-CB2N)
          RBN=DCMPLX(0.0D0,-RRBN)
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
          VD(IOM,MODE)=(0.D0,0.D0)
C
          DO 120 IL=NL-1,1,-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
            SM=DNS(IL)*VS(IL)*VS(IL)
            BM=VS(IL)
            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=SA11*SB11+SA12*SB21
            R12=SA11*SB12+SA12*SB22
            R21=SA21*SB11+SA22*SB21
            R22=SA21*SB12+SA22*SB22
C
            RK11=SAK11*SB11+SAK12*SB21+SA11*SBK11+SA12*SBK21
            RK12=SAK11*SB12+SAK12*SB22+SA11*SBK12+SA12*SBK22
            RK21=SAK21*SB11+SAK22*SB21+SA21*SBK11+SA22*SBK21
            RK22=SAK21*SB12+SAK22*SB22+SA21*SBK12+SA22*SBK22
C
            RW11=SAW11*SB11+SAW12*SB21+SA11*SBW11+SA12*SBW21
            RW12=SAW11*SB12+SAW12*SB22+SA11*SBW12+SA12*SBW22
            RW21=SAW21*SB11+SAW22*SB21+SA21*SBW11+SA22*SBW21
            RW22=SAW21*SB12+SAW22*SB22+SA21*SBW12+SA22*SBW22
C
            IF((MLS.LE.IL).AND.(IL.LE.MLE))THEN
              FL1=SMN*RBN*R12+R22
              FL2=SMN*RBN*R11+R21
              GB2=CB2-1.D0
              DS1=DCMPLX(SDQ/WN,0.D0)
              DS2=DCMPLX(0.D0,(1.D0-CQ)/(WN*SM*GB2))
              TM=2.D0/(1.D0+VS(NL)*SM/(CC(IOM,MODE)*SMN))
C              TM=2.D0/(1.D0+DNS(IL)*VS(IL)/(DNS(NL)*VS(NL)))
              XO=0.D0
              ZO=0.D0
              DO 130 IIL=1,IL-1
                ZO=ZO+DEP(IIL)
  130         CONTINUE
              VSN=VS(NL)
              CALL HKSHS1(OM,XF,ZF,XO,ZO,VSN,W,DA,D,VRR,RT,VF)
              VD(IOM,MODE)=VD(IOM,MODE)-TM*SM*(FL1*DS1-FL2*DS2)*VF
            END IF
C
  120     CONTINUE
C
          FLKIM=IMAG(RK21+SMN*(RBN*RK11-CB2N*R11/(RBN*WN)))
          FLWIM=IMAG(RW21+SMN*(RBN*RW11+CB2N*R11/(CV*RBN*WN)))
          UU(IOM,MODE)=-FLKIM/FLWIM
C
          WRBN=WN*RRBN
          IF(MLE.EQ.NL)THEN
            XO=0.D0
            ZO=0.D0
            DO 140 IIL=1,NL-1
              ZO=ZO+DEP(IIL)
  140       CONTINUE
            HN1=ZO
            DO 150 IH=1,NH
              ZO=HN1+DH*DFLOAT(IH-1)
              ZO1=ZO+DH
              CALL HKSHS1(OM,XF,ZF,XO,ZO,VS(NL),W,DA,D,VRR,RT,VF)
              EX1=DEXP(-WRBN*(ZO-HN1))
              EX2=DEXP(-WRBN*(ZO1-HN1))
              VD(IOM,MODE)=VD(IOM,MODE)-SMN*VF*(EX1-EX2)/WRBN
  150       CONTINUE
          END IF
C
          VD(IOM,MODE)=2.D0*VD(IOM,MODE)/FLKIM
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
      RBM2=CB2-1.0D0
      IF(BM-CV)10,20,30
C
   10   RBM=DSQRT(RBM2)
        QM=WN*RBM*DEP
        CQ=DCOS(QM)
        SQ=DSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/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(-RBM2)
        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
C
      END
