C
      PROGRAM MLVRICK
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 ML    : NL IN MAIN PROGRAM DIMENSION (LAYER NUMBER)
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 
C VS (I): VS OF Ith LAYER      
C DEP(I): DEPTH ( THICKNESS ) OF Ith LAYER  
C NL    : TOTAL NUMBER OF LAYERS            
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 ** 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 ALM(I,J): AMPLITUDE SPEC (i*ALM/K) OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : TO OBTAIN AN ACTUAL ALM, MULTIPLY (-i) BY ALM
C
      PARAMETER(MOM=150,MMD=6,ML=20,MOB=6)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),ALM(MOM,0:MMD)
     *         ,DNS(ML),VS(ML),DEP(ML),XOB(MOB)
      COMPLEX*16 VDO(MOM,MOB),EX,BBB
      CHARACTER FNAME*15
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
 1000 FORMAT(/' From MLVRICK :' 
     *       /' Cal. of LOVE WAVES IN A 2D MEDIUM'
     *       /' SOURCE: P=RICKER*SM1')
C
C *** DATA INPUT ***
C
      CALL LVINN(DOM,NOM,ML,NL,DNS,VS,DEP,NDOM,NSMD,NEMD,NDC,NADD)
C
C *** ENTER OF SOURCE DEPTH IN ANY LAYER ***
C
   10 WRITE(*,1010)NL
      READ(*,*)MLS
      IF(MLS.EQ.0)GOTO 10
      IF(MLS.GT.NL)GOTO 10
      WRITE(*,1020)DEP(MLS)
      READ(*,*)SDEP
 1010 FORMAT(/' From MLVRICK : layers number =',I5,
     *       /' Enter source layer number (from top layer) : ',\)
 1020 FORMAT( ' Source layer depth =',F10.4,
     *       /' Enter source depth from upper boundary : ',\)
C
C *** ENTER OF OBS. DEPTH IN ANY LAYER ***
C
   20 WRITE(*,2010)NL
      READ(*,*)MLO
      IF(MLO.EQ.0)GOTO 20
      IF(MLO.GT.NL)GOTO 20
      WRITE(*,2020)DEP(MLO)
      READ(*,*)ODEP
 2010 FORMAT(/' From MLVRICK : layers number =',I5,
     *       /' Enter obs. layer number (from top layer) : ',\)
 2020 FORMAT( ' obs. layer depth =',F10.4,
     *       /' Enter obs. depth from upper boundary : ',\)
C
C *** INPUT OF OBSERVATION POINTS ON BASIN ***
C
   21 WRITE(*,1011)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 21
      END IF
      IF(NOB.EQ.0)GOTO 21
C
      DO 101 I=1,NOB
        WRITE(*,1021)I
        READ(*,*)XOB(I)
  101 CONTINUE
C
 1011 FORMAT(/' Enter Number of Observation Point : ',\)
 1021 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
C
C *** INPUT OF Q ***
C
      WRITE(*,1040)
      READ(*,*)IQF
      IF(IQF.EQ.0)THEN
        WRITE(*,1050)
        READ(*,*)QFL
      END IF
 1040 FORMAT(/' Include Quality Factor (yes=0)?   : ',\)
 1050 FORMAT( ' Enter Quality Factor for Love Wave: ',\)
C
C *** TP OF RICKER WAVELET ***
C
      WRITE(*,1221)
      READ(*,*)TP,TS
C
 1221 FORMAT(/' Enter Tp and Ts (sec) of Ricker wavelet : ',\)
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 LVUAL4(DOM,NOM,MOM,ML,NL,DNS,VS,DEP
     *           ,CC,NMODE,CFOM,NSMD,MMD,MLS,SDEP,MLO,ODEP,UU,ALM)
C
C *** LOVE WAVE AMPLITUDE ***
C
      IF(IQF.NE.0)QFL=9999.D0
      AA=DNS(1)*VS(1)*VS(1)
      BB=TP*TP*0.0253302D0
C
      DO 200 IOM=1,NOM
C
        OM=DFLOAT(IOM)*DOM
        BBB=CDEXP(DCMPLX(0.D0,-OM*TS))
C
        DO 210 IOB=1,NOB
C
          VDO(IOM,IOB)=(0.D0,0.D0)
C
          DO 220 MODE=NSMD,NCFOM
            IF(CC(IOM,MODE).EQ.0.D0)THEN
              EX=(0.D0,0.D0)
              EXQ=1.D0
            ELSE
              CCC=CC(IOM,MODE)
              EX=CDEXP(DCMPLX(0.D0,-OM*XOB(IOB)/CCC))
              EXQ=DEXP(-OM*XOB(IOB)/(2.D0*QFL*UU(IOM,MODE)))
            END IF
            IF(IQF.NE.0)EXQ=1.D0
            WNI=CC(IOM,MODE)/OM
            VDO(IOM,IOB)=VDO(IOM,IOB)-DCMPLX(0.D0,1.D0)
     *                   *ALM(IOM,MODE)*WNI*EX*EXQ
  220     CONTINUE
C
          OM2=OM*OM
          VDO(IOM,IOB)=VDO(IOM,IOB)*AA*BBB*OM2*DEXP(-OM2*BB)
C
  210   CONTINUE
C
  200 CONTINUE
C
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,1400)DOM,NOM
        DO 300 IOB=1,NOB
        WRITE(*,1410)XOB(IOB)
          DO 310 K=1,NOM,4
            WRITE(*,1420)(VDO(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
C
 1400 FORMAT(/' DOM,NOM=',F6.2,I5)
 1410 FORMAT(/' X=',F10.2)
 1420 FORMAT(8D10.3)
C
C *** OUTPUT OF DATA FILE ***
C
      WRITE(*,*)
   30 WRITE(*,1310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 30
      WRITE(*,1320)
 1310 FORMAT( ' Enter file name : ',\)
 1320 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VV roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1500)DOM,NOM
        WRITE(8,1510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,5
            WRITE(8,1520)(VDO(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 1500 FORMAT(F15.8,I15)
 1510 FORMAT(///)
 1520 FORMAT(10F8.5)
C
C *** DATA OUTPUT ON SCREEN ***
C
      CALL LVOTN(DOM,NOM,MOM,ML,NL,VS,CFOM,NCFOM,CC
     *          ,NSMD,MMD,UU,ALM,NADD)
C
      STOP
      END
