C
C *** PHASE AND GROUP VELOCITY AND MEDIUM RESPONSE SPECTRUM OF LOVE WAVE
C     NEW VERSION OF MLOVE :
C     WE CAN OBTAIN MEDIUM RESPONSE (ALM) OF DISP, VEL OR ACC
C
      PROGRAM MLVCUAM
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 NADD  : =1: ADD FROM MSMD TO NEMD FOR DALM
C ** CALCULATE DATA **
C CFOM(I)    : CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM      : HIGHEST MODE CALCURATED IN LVCT
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 DALM(I,J,M): AMPLITUDE SPEC (i*ALM) OF Jth MODE LOVE WAVE AT Ith OMEGA
C            : OF Mth LAYERS
C            : TO OBTAIN AN ACTUAL DALM, MULTIPLY (-i) BY DALM
C DAL(I,J)   : AMPLITUDE SPEC (i*ALM) OF Jth MODE LOVE WAVE AT Ith OMEGA
C            : TO OBTAIN AN ACTUAL DAL, MULTIPLY (-i) BY DAL
C
      PARAMETER(MOM=300,MMD=10,ML=35)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),DAL(MOM,0:MMD)
     *         ,DNS(ML),VS(ML),DEP(ML)
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)MOM,MMD,ML
 1000 FORMAT(/' From MLVCUAM :' 
     *  /' Cal. of C,U and DALM*i of Love wave in multi-layers'
     *  /'   C    : phase velocity'
     *  /'   U    : group velocity'
     *  /'   DALM : amplification factor'
     *  /'          due to unit horizontal displacement'
     *  /' MOM =',I5,' ,MMD =',I5,' ,ML =',I5)  
C
C *** DATA INPUT ***
C
   10 CALL LVINN(DOM,NOM,ML,NL,DNS,VS,DEP,NDOM,NSMD,NEMD,NDC,NADD)
C
C *** INPUT DATA CHECK ***
C
      IF(NOM.GT.MOM)THEN
        WRITE(*,*)'NOM (omega number) > MOM !!!, MOM =',MOM
        GOTO 10
      ENDIF
      IF(NL.GT.ML)THEN
        WRITE(*,*)'NL (layer number) > ML !!!, ML =',ML
        GOTO 10
      ENDIF
      IF(NEND.GT.MMD)THEN
        WRITE(*,*)'NEMD (output highest mode) > MMD !!!, MMD =',MMD
        GOTO 10
      ENDIF
C
C *** ENTER OF LAYERS OF OUTPUT DALM ***
C
   20 WRITE(*,1010)NL
      READ(*,*)MLS
      IF(MLS.EQ.0)GOTO 20
      IF(MLS.GT.NL)GOTO 20
   30 WRITE(*,1020)
      READ(*,*)MLE
      IF(MLE.EQ.0)GOTO 30
      IF(MLE.GT.NL)GOTO 30
 1010 FORMAT(/' From MLVCUAM : total layers number =',I5,
     *       /' To get DALM by adding from first to final layer''DALM'
     *       /'   Enter first layer number (MLS>= 1 ): ',\)
 1020 FORMAT( '   Enter final 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 LVUDAM(DOM,NOM,MOM,ML,NL,DNS,VS,DEP,MLS,MLE
     *           ,CC,NMODE,CFOM,NSMD,MMD,UU,DAL)
C
C *** SELECTION OF MEDIUM RESPONSE (ALM) : DISP, VEL OR ACC ? ***
C
      WRITE(*,*)
      WRITE(*,1110)
      READ(*,*)ISEL
 1110 FORMAT(' Enter type of DALM (disp=0,vel=1,acc=2) : ',\)
C
C *** LOVE WAVE AMPLITUDE ***
C
      DO 100 IOM=1,NOM
        OM=DFLOAT(IOM)*DOM
        IF(ISEL.EQ.0)OMFAC=1.D0
        IF(ISEL.EQ.1)OMFAC=OM
        IF(ISEL.EQ.2)OMFAC=OM*OM
        IF(NADD.EQ.1)THEN
          DALC=0.D0
          DO 110 MODE=NSMD,NCFOM
            DALC=DALC+DAL(IOM,MODE)
  110     CONTINUE
          DAL(IOM,NSMD)=DALC*OMFAC
        ELSE
          DO 120 MODE=NSMD,NCFOM
            DAL(IOM,MODE)=DAL(IOM,MODE)*OMFAC
  120     CONTINUE
        ENDIF
  100 CONTINUE
C
C *** DATA OUTPUT ON SCREEN ***
C
      CALL LVOTN(DOM,NOM,MOM,ML,NL,VS,CFOM
     *          ,NCFOM,CC,NSMD,MMD,UU,DAL,NADD)
C
C *** DATA OUTPUT ON PLOTTER ***
C
      WRITE(*,*)
      WRITE(*,2000)
      READ(*,*)IHZ 
 2000 FORMAT(' Plot out (HZ=0 or rad/sec=1): ',\)
C
      IF(IHZ.EQ.0)THEN
      CALL LVPTHZ(DOM,NOM,MOM,ML,NL,VS,MLE,CFOM,NCFOM
     *           ,CC,NSMD,MMD,UU,DAL,NADD)
      ELSE
      CALL LVPTN(DOM,NOM,MOM,ML,NL,VS,MLE,CFOM,NCFOM
     *          ,CC,NSMD,MMD,UU,DAL,NADD)
      END IF
C
      STOP
      END
