C
C *** APROXIMATE LOVE WAVE IN SEDIMENT BASIN (ASUUMED AT (0,0))
C     BY 2-D SH HASKELL MODEL ***
C
      PROGRAM MLVKIN2S
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 ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCY
C DNS(I): DENSITY OF Ith LAYER ==> SM(I):SHEAR MODULUS => COMMOM SM(35)
C VS (I): VS OF Ith LAYER                              => COMMON VS(35)
C DEP(I): DEPTH ( THICKNESS ) OF Ith LAYER             => COMMON DEP(35)
C NL    : TOTAL NUMBER OF LAYERS                       => COMMON NL
C NDC   : DIVISION NUMBER FROM MIN VS TO MAX VS TO PHASE VELOCITY
C NSMD  : START OUTPUT MODE OF LOVE WAVE
C NEMD  : END OUTPUT MODE OF LOVE WAVE
C NDOM  : DIVISION NUMBER FROM 0 TO MAX OMEGA TO CUT-OFF OMEGA
C XF    : FAULT LOCATION ON X-AXIS 
C ZF    : FAULT LOCATION ON Z-AXIS 
C XO    : LOCATION OF OBSERVATION POINT ON X-AXIS 
C ZO    : LOCATION OF OBSERVATION POINT ON Z-AXIS 
C W     : WIDTH OF FAULT
C DA    : DIP ANGLE OF FAULT (DEGREE--->RADIAN IN THIS SUBROUTINE)
C D     : DISPLACEMENT DIFFELENCE OF FAULT
C VRR   : RUPTURE VELOCITY OF FAULT
C RT    : RISE TIME
C ** CALCULATE DATA **
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : LOVE WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): OUTPUT HIGHEST LOVE MODE NUMBER AT Ith OMEGA
C UU(I,J) : GROUP VELOCITY OF Jth MODE LOVE WAVE AT Ith OMEGA
C AL(I,J) : AMPLITUDE SPEC (i*AL/K) OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : TO OBTAIN AN ACTUAL AL/K, MULTIPLE AL/K BY (-i)
C VF      : Y DISPLACEMENT AT (XO,ZO) BY FAULT
C VLS(I,J): LOVE WAVE GENERATED ON SEDIMENT EDGE OF Ith OMEGA AT Jth 
C           POINT ON BASIN SURFACE
C VLB(I,J): LOVE WAVE GENERATED IN BEDROCK OF Ith OMEGA AT Jth POINT
C           ON BASIN SURFACE
C VSH(I,J): SH WAVE OF Ith OMEGA AT Jth POINT ON BAISN (SH)
C AMP(I)  : AMPLITUDE FACTOR ON SURFACE OF Ith OMEGA
C
      PARAMETER(MOM=50,MMD=10,MOB=10)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM),XOB(MOB)
     *         ,UU(MOM,0:MMD),AL(MOM,0:MMD)
      COMPLEX*16 VF,VLS(MOM,MOB),VLB(MOM,MOB),VSH(MOM,MOB)
     *          ,EX,VLSA,AMP(MOM)
      CHARACTER FNAME*15
      COMMON /LAYER/NL,VS(35),DEP(35),SM(35) /DNSITY/DNS(35)
C
C *** INPUT OF FAULT RARAMETERS ***
C
      CALL H2DIN(XF,ZF,W,DA,D,VRR,RT)
C
C *** INPUT OF SOIL DATA ***
C
   10 CALL LOVEIN(DOM,NOM,NDOM,NSMD,NEMD,NDC)
C
C *** INPUT OF OBSERVATION POINTS ON BASIN ***
C
      WRITE(*,*)
      WRITE(*,1000)
   20 WRITE(*,1010)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 20
      END IF
      IF(NOB.EQ.0)GOTO 20
C
      DO 100 I=1,NOB
        WRITE(*,1020)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
      WRITE(*,1030)
      READ(*,*)SDP,NH
C
 1000 FORMAT(/' From MLVKIN2S : Basin edge is located at (0,0)'
     *       /' Approximate Lave Wave on Basin by 2-D Haskell'/)
 1010 FORMAT( ' Enter Number of Observation Point : ',\)
 1020 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
 1030 FORMAT(/' Enter depth and number of depth' 
     *       ,' for integration in lower layer : ',\)
C
C *** CHECK OF INPUT DATA ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MLVKINE'
        WRITE(*,*)'MOM < NOM THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MLVKINE :'
        WRITE(*,*)'MMD < NEMD THUS STOP !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        STOP
      END IF
C
C *** CUT-OFF OMEGA ***
C
       CALL LOVECT(DOM,NOM,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITY ***
C
      CALL LOVECC(DOM,NOM,MOM,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
       CALL LVUA(DOM,NOM,MOM,CC,NMODE,CFOM,NSMD,MMD,UU,AL)
C
C *** CALCULATION OF LOVE WAVES GENERATED AT BASIN EDGE (0,0)***
C
      XO=0.D0
      ZO=0.D0
      VSN=VS(NL)
      TOKA=2.D0/(1.D0+DNS(1)*VS(1)/(DNS(NL)*VS(NL)))
      DO 200 K=1,NOM
        OM=DOM*DFLOAT(K)
        CALL HKSHS1(OM,XF,ZF,XO,ZO,VSN,W,DA,D,VRR,RT,VF)
        DO 210 IOB=1,NOB
          VLS(K,IOB)=(0.D0,0.D0)
          DO 220 MODE=NSMD,NCFOM
            IF(CC(K,MODE).EQ.0.D0)THEN
              SR1=0.D0
              EX=(0.D0,0.D0)
            ELSE
              CCC=CC(K,MODE)
              R1=OM*DSQRT(CCC*CCC-VS(1)*VS(1))/(VS(1)*CCC)
              SR1=DSIN(R1*DEP(1))/R1
              EX=CDEXP(DCMPLX(0.D0,-OM*XOB(IOB)/CCC))
            END IF
            VLSA=-2.D0*TOKA*SM(1)*SR1*VF*AL(K,MODE)*EX
            VLS(K,IOB)=VLS(K,IOB)+VLSA
  220     CONTINUE
  210   CONTINUE
  200 CONTINUE
C
C *** OUTPUT OF LOVE DATA TO FILE ***
C
      WRITE(*,*)
   30 WRITE(*,1210)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 30
      WRITE(*,1220)
 1210 FORMAT( ' Enter file name of sediment Love wave (PLOTMS) : ',\)
 1220 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VLS roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1300)DOM,NOM
        WRITE(8,1310)
        DO 300 IOB=1,NOB
          DO 310 K=1,NOM,5
            WRITE(8,1320)(VLS(K+J,IOB),J=0,4)
  310     CONTINUE
  300   CONTINUE
      CLOSE(8)
C
 1300 FORMAT(F15.8,I15)
 1310 FORMAT(///)
 1320 FORMAT(10F8.5)
C
C *** CALCULATION OF LOVE WAVES GENERATED ON BEDROCK (X=0) ***
C
      DO 330 K=1,NOM
        DO 340 IOB=1,NOB
          VLB(K,IOB)=(0.D0,0.D0)
  340   CONTINUE
  330 CONTINUE
C
      XO=0.D0
      VSN=VS(NL)
      DO 350 IH=1,NH
        AA=SDP*DFLOAT(IH-1)
        BB=AA+SDP
        ZO=DEP(1)+AA
        DO 360 K=1,NOM
          OM=DOM*DFLOAT(K)
          CALL HKSHS1(OM,XF,ZF,XO,ZO,VSN,W,DA,D,VRR,RT,VF)
          DO 370 IOB=1,NOB
            DO 380 MODE=NSMD,NCFOM
              IF(CC(K,MODE).EQ.0.D0)THEN
                CR1=0.D0
                ER2=0.D0
                EX=(0.D0,0.D0)
              ELSE
                CCC=CC(K,MODE)
                R1=OM*DSQRT(CCC*CCC-VS(1)*VS(1))/(VS(1)*CCC)
                R2=OM*DSQRT(VS(NL)*VS(NL)-CCC*CCC)/(VS(NL)*CCC)
                CR1=DCOS(R1*DEP(1))
                ER2=(DEXP(-R2*AA)-DEXP(-R2*BB))/R2
                EX=CDEXP(DCMPLX(0.D0,-OM*XOB(IOB)/CCC))
              END IF
              VLSA=-2.D0*SM(NL)*VF*AL(K,MODE)*CR1*ER2*EX
              VLB(K,IOB)=VLB(K,IOB)+VLSA
  380       CONTINUE
  370     CONTINUE
  360   CONTINUE
  350 CONTINUE
C
C *** OUTPUT OF LOVE DATA TO FILE ***
C
      WRITE(*,*)
      WRITE(*,1350)
   35 READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 35
      WRITE(*,1360)
 1350 FORMAT( ' Enter file name of bedrock Love wave (PLOTMS) : ',\)
 1360 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VLB roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1370)DOM,NOM
        WRITE(8,1380)
        DO 305 IOB=1,NOB
          DO 315 K=1,NOM,5
            WRITE(8,1390)(VLB(K+J,IOB),J=0,4)
  315     CONTINUE
  305   CONTINUE
      CLOSE(8)
C
 1370 FORMAT(F15.8,I15)
 1380 FORMAT(///)
 1390 FORMAT(10F8.5)
C
C *** multiple reflection of SH wave ***
C
      THETA=0.D0
      CALL SBMLTPL(0.D0,DOM,NOM,MOM,THETA,1,AMP)
C
C *** CALCULATION OF SH WAVES ***
C
      ZO=DEP(NL)
      VSN=VS(NL)
      DO 400 IOB=1,NOB
        XO=XOB(IOB)
        DO 410 K=1,NOM
          OM=DOM*DFLOAT(K)
          CALL HKSHS1(OM,XF,ZF,XO,ZO,VSN,W,DA,D,VRR,RT,VF)
          VSH(K,IOB)=VF*AMP(K)
  410   CONTINUE
  400 CONTINUE
C
C *** OUTPUT OF SH DATA TO FILE ***
C
      WRITE(*,*)
   40 WRITE(*,1410)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 40
      WRITE(*,1420)
 1410 FORMAT( ' Enter file name of SH wave data (PLOTMS) : ',\)
 1420 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VSH roop')
C      
      OPEN(7,FILE=FNAME)
        WRITE(7,1500)DOM,NOM
        WRITE(7,1510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,5
            WRITE(7,1520)(VSH(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
      CLOSE(7)
C
 1500 FORMAT(F15.8,I15)
 1510 FORMAT(///)
 1520 FORMAT(10F8.5)
C
C *** LOVES + LOVEB + SH ***
C
      WRITE(*,1550)
      READ(*,*)NY
      IF (NY.EQ.0)THEN
        WRITE(*,1560)
        READ(*,*)NAY
        IF (NAY.EQ.0)THEN
          AFC=1.D0
        ELSE
          WRITE(*,1570)
          READ(*,*)AFC
        END IF
      END IF
C
 1550 FORMAT(/' Love + SH wave ? (enter yes = 0) : ',\)
 1560 FORMAT(/' Amp factor of SH wave = 1.0 ? (enter yes = 0) : ',\)
 1570 FORMAT(/' Enter amp factor of SH wave : ',\)
C
   45 DO 600 K=1,NOM
        OM=DOM*DFLOAT(K)
        DO 610 IOB=1,NOB
          IF (NY.EQ.0)THEN
            VLS(K,IOB)=VLS(K,IOB)+VLB(K,IOB)+AFC*VSH(K,IOB)
          ELSE
            VLS(K,IOB)=VLS(K,IOB)+VLB(K,IOB)
          END IF
  610   CONTINUE
  600 CONTINUE
C
C *** OUTPUT OF SH + LOVE DATA TO FILE ***
C
      WRITE(*,*)
   50 WRITE(*,1600)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 50
      WRITE(*,1610)
 1600 FORMAT(/' Enter file name of Love s+b+SH Love data (PLOTMS) : ',\)
 1610 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VLS roop')
C      
      OPEN(9,FILE=FNAME)
        WRITE(9,1700)DOM,NOM
        WRITE(9,1710)
        DO 700 IOB=1,NOB
          DO 710 K=1,NOM,5
            WRITE(9,1720)(VLS(K+J,IOB),J=0,4)
  710     CONTINUE
  700   CONTINUE
      CLOSE(9)
C
 1700 FORMAT(F15.8,I15)
 1710 FORMAT(///)
 1720 FORMAT(10F8.5)
C
C *** OUTPUT OF DATA OF LOVE WAVE ON SCREEN ***
C
      WRITE(*,1800)
      READ(*,*)NYES
      IF(NYES.NE.0)GOTO 60
 1800 FORMAT(/' Output Love wave data on screen ? (yes=0) :',\) 
      CALL LOVEOT(DOM,NOM,MOM,CFOM,NCFOM,CC,NSMD,MMD,UU,AL)
C
   60 STOP
      END
