C
C *** READ OF SOIL DATA FROM DATA FILE
C
      SUBROUTINE LOVEIN(DOM,NOM,NDOM,NSMD,NEMD,NDC)
C
C ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCY
C NDOM  : DIVISION NUMBER FROM 0 TO MAX OMEGA TO GET CUT-OFF OMEGA
C NSMD  : FIRST OUTPUT MODE OF LOVE WAVE
C NEMD  : END OUTPUT MODE OF LOVE WAVE
C NDC   : DIVISION NUMBER FROM MIN TO MAX VS TO GET PHASE VELOCITY
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
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30
      COMMON / LAYER / NL,VS(35),DEP(35),SM(35) / DNSITY / DNS(35)
     *       / MINMAX / VSMIN, VSMAX
C
C *** SOIL DATA INPUT FROM FILE DATA ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1000 FORMAT( ' From SUBROUTINE LOVEIN :')
 1010 FORMAT(/' Enter Soil Data File Name ',\)

      OPEN (7,FILE=FNAME)
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VP,VS(IL),DEP(IL),QF
        SM(IL)=DNS(IL)*VS(IL)*VS(IL)
  100 CONTINUE
      CLOSE(7)
 1200 FORMAT(I10)
 1210 FORMAT(5F10.0)
C
C *** DATA ECHO ***
C
      WRITE(*,2000)NL
      DO 200 IL=1,NL
        WRITE(*,2010)IL,DNS(IL),VS(IL),DEP(IL)
  200 CONTINUE
 2000 FORMAT(/' Total Layers Number            =',I5)
 2010 FORMAT( ' Layer Num., Density, Vs, Depth =',I3,3F10.3)
C
C *** MIN AND MAX VS AND SOIL DATA CHECK ***
C
      VSMIN=VS(NL)
      VSMAX=0.0D0
      DO 300 IL=1,NL
        VSMIN=MIN(VSMIN,VS(IL))
        VSMAX=MAX(VSMAX,VS(IL))
  300 CONTINUE
      IF(VSMAX.NE.VS(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE LOVEIN'
        WRITE(*,*)' VS(NS).NE.VSMAX, THUS STOP !!!'
        STOP
      END IF
C
      IF(NL.GT.50)THEN
        WRITE(*,*)' FROM SUBROUTINE LOVEIN'
        WRITE(*,*)' SOIL LAYERS NUMBER > 50, THUS STOP !!!'
        STOP
      END IF
C
C *** DATA INPUT FROM SCREEN ***
C
      WRITE(*,3000)
      READ(*,*)DOM,NOM
      OMMAX=DOM*DFLOAT(NOM)
      WRITE(*,3010)0.D0,OMMAX
      WRITE(*,3020)
      READ(*,*)NDOM
      WRITE(*,3030)VSMIN,VSMAX
      WRITE(*,3040)
      READ(*,*)NDC
      WRITE(*,3050)
      READ(*,*)NSMD
      WRITE(*,3060)
      READ(*,*)NEMD
      WRITE(*,3070)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 10
 3000 FORMAT(/' Enter Delta Omega and Omega Number        ?    :',\)
 3010 FORMAT(/' Minimum Omega and Maximum Omega = ',2F8.3)
 3020 FORMAT( ' To Get Cut-Off Omega,  Enter',
     *       /' Division Number between Min and Max Omega ?    :',\)
 3030 FORMAT(/' Minimum Vs and Maximum Vs = ',2F10.3)
 3040 FORMAT( ' To Get Phase Velocity, Enter',
     *       /' Division Number between Min and Max Vs       ? :',\)
 3050 FORMAT(/' Enter Output Start Love Mode (Fundamental=0) ? :',\)
 3060 FORMAT( ' Enter Output End Love Mode ( >= Start Mode ) ? :',\)
 3070 FORMAT(/' Input Data OK ? ( yes = Enter 0 )            ? :',\)
      WRITE(*,*)
C
      RETURN
      END
C
C
C
C *** OUTPUT OF LOVE WAVE PHASE AND GROUP VELOCITY AND 
C     AMPLITUDE SPECTRUM ON SCREEN
C
      SUBROUTINE LOVEOT(DOM,NOM,MOM,CFOM,NCFOM,CC,NSMD,MMD,UU,ALK)
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 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 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
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD)
     *         ,UU(MOM,0:MMD),ALK(MOM,0:MMD)
      COMMON / LAYER / NL,VS(35),DEP(35),SM(35)
C
      WRITE(*,1000)
      WRITE(*,1010)NCFOM
 1000 FORMAT(/' From SUBROUTINE LOVEOT :')
 1010 FORMAT(/' Highest Mode of Cut-Off Omega =',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
