C
C *** PHASE AND GROUP VELOCITY AND AMP OF GREEN FUNC ***
C     NEW VERSION OF MLVCUA4 :
C     WE CAN AMP FUNC (ALM/K) OR ALM OF DISP, VEL OR ACC
C     WE CAN OBTAIN ARBITRARY SOURCE AND OBSERVATION POINTS 
C
      PROGRAM MLVAMP
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=300,MMD=10,ML=40)
      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)
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
 1000 FORMAT(/' From MLVAMP :' 
     *       /' Cal. of C,U and ALM*i/K or AMP*I of Love wave'
     *       /' Source and obs. points can be located in any layer'
     *       /' If source and obs. is located on free surface,'
     *       /' ALM is equal to Harkrider''s AL')
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 MLVAMP : 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 MLVAMP : 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 *** ENTER OF AMP OR AMP/K ***
C
      WRITE(*,2050)
      READ(*,*)IAMP
 2050 FORMAT(/' From MLVAMP : ALM(=0) OR ALM/K(=1) : ',\)
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 *** SELECTION OF MEDIUM RESPONSE (ALM) : DISP, VEL OR ACC ? ***
C
      WRITE(*,*)
      WRITE(*,3110)
      READ(*,*)ISEL
 3110 FORMAT(' Enter type of ALM (disp=0,vel=1,acc=2) : ',\)
C
C *** LOVE WAVE AMPLITUDE ***
C
      DO 110 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
          ALC=0.D0
          DO 120 MODE=NSMD,NCFOM
            AAA=1.0 
            IF(IAMP.EQ.1)AAA=CC(IOM,MODE)/OM
            ALC=ALC+ALM(IOM,MODE)*AAA
  120     CONTINUE
          ALM(IOM,NSMD)=ALC*OMFAC
        ELSE    
          DO 100 MODE=NSMD,NCFOM
            AAA=1.0 
            IF(IAMP.EQ.1)AAA=CC(IOM,MODE)/OM
            ALM(IOM,MODE)=ALM(IOM,MODE)*OMFAC*AAA
  100     CONTINUE
        END IF
  110 CONTINUE
C
C *** DATA OUTPUT ON SCREEN ***
C
      CALL LVOTN(DOM,NOM,MOM,ML,NL,VS,CFOM,NCFOM,CC
     *          ,NSMD,MMD,UU,ALM,NADD)
C
C *** DATA OUTPUT ON PLOTTER ***
C
      WRITE(*,*)
      WRITE(*,3500)
      READ(*,*)IHZ
 3500 FORMAT(' Enter type of PLOT (rad/sec=0 or HZ=1): ',\)
C
      MLE=0
      IF(IHZ.EQ.0)CALL LVPTN(DOM,NOM,MOM,ML,NL,VS,MLE,CFOM
     *                      ,NCFOM,CC,NSMD,MMD,UU,ALM,NADD)
      IF(IHZ.EQ.1)CALL LVPTHZ(DOM,NOM,MOM,ML,NL,VS,MLE,CFOM
     *                       ,NCFOM,CC,NSMD,MMD,UU,ALM,NADD)
C
      STOP
      END
