C
C *** PHASE AND GROUP VELOCITY AND AMPLITUDE SPECTRUM OF LOVE WAVE
C
      PROGRAM MLOVE
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 ** 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 ALK(I,J): AMPLITUDE SPEC (i*AL/K) OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : TO OBTAIN AN ACTUAL AL/K, MULTIPLY (-i) BY AL/K 
C
      PARAMETER(MOM=500,MMD=10)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),ALK(MOM,0:MMD)
C
C *** DATA INPUT ***
C
      CALL LOVEIN(DOM,NOM,NDOM,NSMD,NEMD,NDC)
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 LOVEUA(DOM,NOM,MOM,CC,NMODE,CFOM,NSMD,MMD,UU,ALK)
C
C *** DATA OUTPUT ON SCREEN ***
C
      CALL LOVEOT(DOM,NOM,MOM,CFOM,NCFOM,CC,NSMD,MMD,UU,ALK)
C
C *** DATA OUTPUT ON PLOTTER ***
C
      CALL LOVEPT(DOM,NOM,MOM,CFOM,NCFOM,CC,NSMD,MMD,UU,ALK)
C
      STOP
      END
