C
C *** WAVE READ, FOURIER TRANS., AMPLITUDE OF LOVE WAVE
C     SPECTRUM PLOT, INVERSE FOURIER TRANS., AND WAVE PLOT
C
      PROGRAM MLVWV
C
C ** PARAMETER **
C MWV : WAVE NUMBER
C MOM : OMEGA NUMBER (> NFAST > MWV/2+1)
C ML  : LAYER NUMBER
C MMD : MODE NUMBERT
C ** INPUT DATA **
C TITL  : TITLE OF EACH CHANNEL
C WV    : WAVE AMPLITUDE DATA
C NWV   : WAVE DATA NUMBER
C NOM   : OMEGA NUMBER (> NFAST > NWV/2+1)
C DT    : DELTA TIME ( SEC )
C DOM   : INCREMENT OF CIRCULAR FREQUENCY (DELTA OMEGA)
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  : FIRST MODE OF LOVE WAVE
C NEMD  : FINAL MODE OF LOVE WAVE
C NDOM  : DIVISION NUMBER FROM 0 TO MAX OMEGA TO CUT-OFF OMEGA
C ** CALCURATED 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 DAL(I,J) : AMPLITUDE SPEC OF Jth MODE LOVE WAVE AT Ith OMEGA
C FOUR(I)  : FOURIER COEFFICIENTS OR SPECTRUM
C          : FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0
C
      PARAMETER(MWV=1025,MOM=513,MMD=5,ML=10)
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER TITL*50,FNAME*30
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),DAL(MOM,0:MMD)
      DIMENSION DNS(ML),DEP(ML),VS(ML),QF(ML),WV(MWV)
      COMPLEX*16 FOUR(MWV),FF,EXC
C
C *** HEAD ***
C
      WRITE(*,1000)MWV,MOM,MMD,ML
 1000 FORMAT(/' From MLVWV (main):' 
     *  /' Cal. of analytical Love waves in multi-layers'
     *  /' MWV =',I5,' MOM =',I5,' ,MMD =',I5,' ,ML =',I5)  
C
C *** WAVE DATA READ AND FOURIER COEFF. TRANS. ***
C
   10 CALL WVRD(TITL,WV,NWV,MWV,DT)
      CALL WVFR(WV,NWV,MWV,FOUR,NFAST,DT)
C
      DUR=DFLOAT(NFAST)*DT
      FOM=0.0D0
      DOM=6.2831853070D0/DUR
      NOM=NFAST/2+1
C
C *** INPUT SOIL DATA AND LOVE WAVE DATA INPUT ***
C
   15 CALL LVWVIN(ML,NL,DNS,VS,DEP,QF,NDOM,NSMD,NEMD,NDC)
C
C *** INPUT DATA CHECK ***
C
      IF(NFAST.GT.MWV)THEN
        WRITE(*,*)'NFAST (FFT wave numbers) > MWV !!!, MWV =',MWV
        GOTO 10
      ENDIF
      IF(NOM.GT.MOM)THEN
        WRITE(*,*)'NOM (omega numbers) > MOM !!!, MOM =',MOM
        GOTO 10
      ENDIF
      IF(NL.GT.ML)THEN
        WRITE(*,*)'NL (layer numbers) > ML !!!, ML =',ML
        GOTO 15
      ENDIF
      IF(NEND.GT.MMD)THEN
        WRITE(*,*)'NEMD (highest mode) > MMD !!!, MMD =',MMD
        GOTO 15
      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
C
 1010 FORMAT(/' From MLVWV : 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 *** INPUT OF X-AXIS OF A OBSERBATION  POINT ***
C
      WRITE(*,1030)
      READ(*,*)XOB
 1030 FORMAT(/' Enter X-Axis of OB Point from basin edge: ',\)
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 *** CUT-OFF OMEGA ***
C
      CALL LVCT(DOM,NOM-1,ML,NL,DNS,VS,DEP,NDOM,CFOM,NEMD,MMD,NCFOM)
      WRITE(*,*)'Higest mode of Love waves =',NCFOM
C
C *** PHASE VELOCITY ***
C
      CALL LVCC(DOM,NOM-1,MOM,ML,NL,DNS,VS,DEP
     *         ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      CALL LVDAL(DOM,NOM-1,MOM,ML,NL,DNS,VS,DEP,MLS,MLE
     *          ,CC,NMODE,CFOM,NSMD,MMD,UU,DAL)
      CALL LVOTN(DOM,NOM-1,MOM,ML,NL,VS
     *          ,CFOM,NCFOM,CC,NSMD,MMD,UU,DAL,NADD)
C
C *** CALCULATION OF LOVE WAVES ON BASIN ***
C
      FOUR(1)=FOUR(1)
      IF(MLE.NE.NL)FOUR(1)=(0.D0,0.D0)
C
      IF(IQF.NE.0)QFL=9999.D0
      DO 100 K=2,NOM
        OM=DOM*DFLOAT(K-1)
        FF=(0.D0,0.D0)
        DO 110 MODE=NSMD,NCFOM
          IF(CC(K-1,MODE).EQ.0.D0)THEN
            EXC=(0.D0,0.D0)
            EXQ=1.D0
          ELSE
            EXC=CDEXP(DCMPLX(0.D0,-OM*XOB/CC(K-1,MODE)))
            EXQ=DEXP(-OM*XOB/(2.D0*QFL*UU(K-1,MODE)))
          END IF
          IF(IQF.NE.0)EXQ=1.D0
          FF=FF+DAL(K-1,MODE)*EXC*EXQ/2.D0
  110   CONTINUE
        FOUR(K)=FOUR(K)*FF
  100 CONTINUE
C
C *** OUTPUT OF SPEC DATA TO FILE FOR LAZER SHOT ***
C
      WRITE(*,*)
   11 WRITE(*,1001)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 11
 1001 FORMAT( ' From MLVWV : '
     *       /' Enter spec file name for Lazer Shot : ',\)
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1101)DOM,NOM-1
        WRITE(8,1111)
        DO 101 K=2,NOM,5
          WRITE(8,1121)(FOUR(K+J)*DUR,J=0,4)
  101   CONTINUE
      CLOSE(8)
C
 1101 FORMAT(F15.8,I15)
 1111 FORMAT(///)
 1121 FORMAT(10F8.2)
C
C *** TRANS. FROM FOURIER COEFF. TO AMPLITUDE SPECTRUM ***
C
      DO 150 I=1,NOM
        WV(I)=CDABS(FOUR(I))*DUR
  150 CONTINUE
C
C *** SPECTRUM PLOT ***
C
      CALL SPPT(FOM,DOM,NOM,MWV,WV)
C
C *** INVERSE FOURIER TRANSFORM ***
C
      DO 200 I=2,NFAST/2
        FOUR(NFAST+2-I)=DCMPLX(DREAL(FOUR(I)),-DIMAG(FOUR(I)))
  200 CONTINUE
C           
      CALL WVFT(NFAST,FOUR,MWV,1)
C
      DO 300 I=1,NFAST
        WV(I)=DREAL(FOUR(I))
  300 CONTINUE
C
      CALL WVOT(WV,NFAST,MWV,DT)
      CALL WVPT(TITL,WV,NFAST,MWV,DT)
C
      STOP
      END
C
C *** OUTPUT OF LOVE WAVE PHASE AND GROUP VELOCITY AND 
C     AMPLITUDE SPECTRUM ON SCREEN
C
      SUBROUTINE LVOTN(DOM,NOM,MOM,ML,NL,VS
     *                ,CFOM,NCFOM,CC,NSMD,MMD,UU,ALK,NADD)
C
C ** INPUT DATA **
C DOM     : DELTA OMEGA
C NOM     : TOTAL NUMBER OF CIRCULAR FREQUENCY
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C         : IF MOM > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C ML      : NL IN MAIN 
C NL      : TOTAL NUMBER OF LAYERS           
C VS (I)  : VS OF Ith LAYER      
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM   : FINAL MODE
C CC(I,J) : LOVE WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NSMD    : OUTPUT START LOVE MODE OF LOVE WAVE
C MMD     : NMD IN MAIN PROGRAM DIMENSION
C UU(I,J) : GROUP VELOCITY OF Jth MODE LOVE WAVE AT Ith OMEGA
C ALK(I,J): AMPLITUDE SPECTRUM (AL/K) OF Jth MODE LOVE WAVE AT Ith OMEGA
C         : TO OBTAIN AN ACTUAL AL/K, MULTIPLY AL/K BY (-i)
C NADD    : =1: ADD FROM FIRST MODE TO FINAL MODE
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),UU(MOM,0:MMD),ALK(MOM,0:MMD)
     *         ,VS(ML)
C
      WRITE(*,1000)
      WRITE(*,1010)NCFOM
      IF(NADD.EQ.1)THEN
        WRITE(*,1020)NSMD
      ENDIF
 1000 FORMAT(/' From SUBROUTINE LVOT :')
 1010 FORMAT(/' Highest Mode of Cut-Off Omega =',I5)
 1020 FORMAT( ' DAMP of First Mode is DAMP adding'
     *        '   from First to Final Mode: First Mode =',I5)
C
      DO 100 MODE=NSMD,NCFOM
C
        WRITE(*,1100)MODE
        WRITE(*,1110)CFOM(MODE),VS(NL)
        WRITE(*,1120)
C
        IOMMAX=0
        ALKMAX=0.D0
        DO 110 IOM=1,NOM
          ALKI=DABS(ALK(IOM,MODE))
          IF(ALKI.GT.ALKMAX)THEN
          IOMMAX=IOM
          ALKMAX=ALKI
          END IF
  110   CONTINUE
C
        ICN=1
        DO 120 IOM=1,NOM
C
         OM=DOM*FLOAT(IOM)
         IF(OM.LT.CFOM(MODE))GOTO 120
         PER=6.2831853070D0/OM
         WRITE(*,1130)ICN,OM,PER,CC(IOM,MODE),UU(IOM,MODE)
     *               ,ALK(IOM,MODE),ALK(IOM,MODE)/ALK(IOMMAX,MODE)
         ICN=ICN+1
C
  120   CONTINUE
C
  100 CONTINUE
C
 1100 FORMAT(//' Love Mode =',I5)
 1110 FORMAT(' Cut-Off Omega =',F10.5,',  Phase and Group Vel. =',F9.3/)
 1120 FORMAT(' Num., Omega, Period, Phase Vel., Group Vel., Amp.spc*i'
     *       ,', Amp(Max=1)'/)
 1130 FORMAT(I4,2F8.3,F11.3,F12.3,1X,D12.4,F8.3)  
C
      WRITE(*,*)
C
      RETURN
      END
