C
C *** APROXIMATE LOVE WAVE IN SEDIMENT BASIN (ASUUMED AT (0,0))
C     BY 2-D SH HASKELL MODEL ***
C
      PROGRAM MLVKINEM
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 VDO(I,J): LOVE WAVE GENERATED ON SEDIMENT EDGE OF Ith OMEGA AT Jth 
C         : POINT ON BASIN SURFACE
C
      PARAMETER(MOM=100,MMD=6,MOB=8,ML=35)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM),XOB(MOB)
     *         ,UU(MOM,0:MMD),DNS(ML),VS(ML),DEP(ML)
      COMPLEX*16 VDO(MOM,MOB),VD(MOM,0:MMD),EX
      CHARACTER FNAME*15
C
C *** INPUT OF SOIL DATA ***
C
   10 CALL LVIN(DOM,NOM,ML,NL,DNS,VS,DEP,NDOM,NSMD,NEMD,NDC)
C
C *** CHECK OF INPUT DATA ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MLVKINEM'
        WRITE(*,*)'MOM < NOM RE-ENTRER !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        GOTO 10
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MLVKINEM :'
        WRITE(*,*)'MMD < NEMD RE-ENTER !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        GOTO 10
      END IF
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
 1000 FORMAT(/' From MLVKINEM : Basin edge is located at (0,0)'
     * /' Approximate Lave Wave in multi-layered basin by 2-D Haskell'
     * /' Fundamental solusion in semi-infinite medium is used'/)
 1010 FORMAT( ' Enter Number of Observation Point : ',\)
 1020 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
C
C *** ENTER OF LAYERS OF OUTPUT DALM ***
C
   15 WRITE(*,1110)NL
      READ(*,*)MLS
      IF(MLS.EQ.0)GOTO 15
      IF(MLS.GT.NL)GOTO 15
   25 WRITE(*,1120)
      READ(*,*)MLE
      IF(MLE.EQ.0)GOTO 25
      IF(MLE.GT.NL)GOTO 25
 1110 FORMAT(/' For output of VD : layers number =',I5,
     *       /' Enter start layer number (MLS>= 1 ): ',\)
 1120 FORMAT( ' Enter   end layer number (MLE>=MLS): ',\)
C
C *** CUT-OFF OMEGA ***
C
      CALL LVCT(DOM,NOM,ML,NL,DNS,VS,DEP,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITY ***
C
      CALL LVCC(DOM,NOM,MOM,ML,NL,DNS,VS,DEP
     *         ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      CALL LVUVD(DOM,NOM,MOM,ML,NL,DNS,VS,DEP,MLS,MLE
     *          ,CC,NMODE,CFOM,NSMD,MMD,UU,VD)
C
C *** CALCULATION OF LOVE WAVES ON BASIN ***
C
      DO 200 K=1,NOM
        OM=DOM*DFLOAT(K)
        DO 210 IOB=1,NOB
          VDO(K,IOB)=(0.D0,0.D0)
          DO 220 MODE=NSMD,NCFOM
            IF(CC(K,MODE).EQ.0.D0)THEN
              EX=(0.D0,0.D0)
            ELSE
              CCC=CC(K,MODE)
              EX=CDEXP(DCMPLX(0.D0,-OM*XOB(IOB)/CCC))
            END IF
            VDO(K,IOB)=VDO(K,IOB)+VD(K,MODE)*EX
  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:VDO 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)(VDO(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
      STOP
      END
