C
      SUBROUTINE DATIN(DOM,MOM,NOM,ML,NL,DNS,VP,VS,THK,NDOM,NDC
     *                ,MMD,NSR,NER,NSL,NEL
     *                ,LS,SDEP,LO,ODEP,MOB,NOB,ROB,RAD
     *                ,TP,TS,QX,QY,QZ,AMP,IQF,QFL
     *                ,VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL)
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30
      DIMENSION DNS(ML),VP(ML),VS(ML),THK(ML),ROB(MOB)
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 DATIN :')
 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(IL),VS(IL),THK(IL),QF
  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),VP(IL),VS(IL),THK(IL)
  200 CONTINUE
 2000 FORMAT(/' Total Layers Number            =',I5)
 2010 FORMAT( ' Layer Num., Density, Vp, Vs, Depth =',I3,4F10.3)
C
C *** MIN AND MAX VS AND SOIL DATA CHECK ***
C
      VSMIN=VS(NL)
      VSMAX=0.0D0
      VPMAX=0.0D0 
      MINIL=1
      MAXIL=1
      DO 300 IL=1,NL
        IF(VS(IL).LE.VSMIN)THEN
          VSMIN=VS(IL)
          MINIL=IL
        END IF
        IF(VS(IL).GE.VSMAX)THEN
          VSMAX=VS(IL)
          MAXIL=IL
        END IF
        VPMAX=MAX(VPMAX,VP(IL))
  300 CONTINUE
C
      IF(VSMAX.NE.VS(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE DATIN :'
        WRITE(*,*)' VS(NL).NE.VSMAX, THUS STOP !!!'
        STOP
      END IF
      IF(VPMAX.NE.VP(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE DATIN :'
        WRITE(*,*)' VP(NL).NE.VPMAX, THUS STOP !!!'
        STOP
      END IF
C
C *** MINUMUM AND MAXIMUM PHASE VELOCITY FOR RAYLEIGH WAVE ***
C
      CCMIN=FRHF(VP(MINIL),VS(MINIL))
      CCMAX=FRHF(VP(MAXIL),VS(MAXIL))
C
C *** DATA INPUT FROM SCREEN ***
C
   15 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(*,*)NSR
      WRITE(*,3060)
      READ(*,*)NER
      WRITE(*,3055)
      READ(*,*)NSL
      WRITE(*,3065)
      READ(*,*)NEL
 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',
     *       /' Number of Partition between Min and Max Omega ? :',\)
 3030 FORMAT(/' Minimum C and Maximum Vs = ',2F10.3)
 3040 FORMAT( ' To Get General Phase Velocity (C), Enter',
     *       /' Number of Partition between Min and Max C     ? :',\)
 3050 FORMAT(/' Enter Output Starting Rayleigh Mode ( >= 0 )  ? :',\)
 3060 FORMAT( ' Enter Output Ending Rayleigh Mode             ? :',\)
 3055 FORMAT(/' Enter Output Starting Love Mode     ( >= 0 )  ? :',\)
 3065 FORMAT( ' Enter Output Ending Love Mode                 ? :',\)
C
C *** PARAMETER CHECK ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM DATIN  :'
        WRITE(*,*)'MOM < NOM THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(ML.LT.NL)THEN
        WRITE(*,*)'FROM DATIN  :'
        WRITE(*,*)'ML < NL THUS STOP !!!'
        WRITE(*,*)'ML,NL =',ML,NL
        STOP
      END IF
C
      IF(MMD.LT.NER)THEN
        WRITE(*,*)'FROM DATIN  :'
        WRITE(*,*)'MMD < NER THUS STOP !!!'
        WRITE(*,*)'MMD,NER =',MMD,NER
        STOP
      END IF
C
      IF(MMD.LT.NEL)THEN
        WRITE(*,*)'FROM DATIN  :'
        WRITE(*,*)'MMD < NEL THUS STOP !!!'
        WRITE(*,*)'MMD,NEL =',MMD,NEL
        STOP
      END IF
C
C *** INPUT OF DEPTH OF SP AND OP, AND SO ON ***
C
   20 WRITE(*,4000)NL
      READ(*,*)LS
      IF(LS.EQ.0)GOTO 20
      IF(LS.GT.NL)GOTO 20
      WRITE(*,4010)THK(LS)
      READ(*,*)SDEP
 4000 FORMAT(/' Number of total layers    =',I5,
     *       /' Enter number of the SP layer (from top layer): ',\)
 4010 FORMAT( ' Thickness of the SP layer =',F10.4,
     *       /' Enter SP depth from the upper boundary       : ',\)
C
   30 WRITE(*,4100)NL
      READ(*,*)LO
      IF(LO.EQ.0)GOTO 30
      IF(LO.GT.NL)GOTO 30
      WRITE(*,4110)THK(LO)
      READ(*,*)ODEP
 4100 FORMAT(/' Number of total layers    =',I5,
     *       /' Enter number of the OP layer (from top layer): ',\)
 4110 FORMAT( ' Thickness of the OP layer =',F10.4,
     *       /' Enter OP depth from the upper boundary       : ',\)
C
   40 WRITE(*,4410)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 40
      END IF
      IF(NOB.EQ.0)GOTO 40
C
      DO 400 I=1,NOB
        WRITE(*,4420)I
        READ(*,*)ROB(I)
  400 CONTINUE
      WRITE(*,4430)
      READ(*,*)RAD
      RAD=RAD*0.01745329252D0
C
 4410 FORMAT(/' Enter Number of OB Point               : ',\)
 4420 FORMAT( ' Enter R-axis of OB Point of Num. =',I3,'  : ',\)
 4430 FORMAT(/' Enter Angle of OB from X-axis (degree) : ',\)
C
      WRITE(*,4521)
      READ(*,*)TP,TS
C
 4521 FORMAT(/' Peak Period(Tp) and Delay(Ts) Time of'
     *        ' Ricker wavelet (sec) ?: ',\)
C
      WRITE(*,4550)
      READ(*,*)QX,QY,QZ
C
 4550 FORMAT(/' Enter QX, QY and QZ : ',\)
C
      AMP=DNS(1)*VS(1)*VS(1)
      WRITE(*,4551)AMP
      WRITE(*,4552)  
      READ(*,*)POW
      AMP=AMP**POW
      WRITE(*,4553)AMP
C
 4551 FORMAT(/' Amplification Factor (AMP) = ',2D12.4)
 4552 FORMAT( ' AMP**N; N=?: ',\)
 4553 FORMAT(/' New Amplification Factor   = ',2D12.4)
C
      WRITE(*,4640)
      READ(*,*)IQF
      IF(IQF.EQ.0)THEN
        WRITE(*,4650)
        READ(*,*)QFL
      END IF
 4640 FORMAT(/' Include Quality Factor (yes=0)?        : ',\)
 4650 FORMAT( ' Enter Quality Factor for Surface Wave  : ',\)
C
      WRITE(*,4700)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 15
 4700 FORMAT(/' All Input Data OK ? ( yes = Enter 0 ) ?:',\)
      WRITE(*,*)
C
      RETURN
      END
C
C
C
C
C *** PHASE VELOCITY OF THE RAYLEIGHWAVE IN THE HOMOGENEOUS HALF-SPACE ***
C
      REAL FUNCTION FRHF*8(VP,VS)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DZ=0.0001D0
C
      SP=VS/VP
      IF(SP.GE.0.0D0) Z= 0.91263D0
      IF(SP.GE.0.2D0) Z= 0.90767D0
      IF(SP.GE.0.4D0) Z= 0.88899D0
      IF(SP.GE.0.6D0) Z= 0.83575D0
      IF(SP.GE.0.7D0) Z= 0.77062D0
      IF(SP.GE.0.8D0) Z= 0.63338D0
C
      F0=FRHF1(SP,Z)
  100 Z=Z-DZ
      F1=FRHF1(SP,Z)
      AF=F0/F1
      F00=F0
      F0 =F1
      IF(AF.GE.0.D0)GOTO 100
      FZ=Z+DZ*F0/(F0-F00)
C
      FRHF=VS*DSQRT(FZ)
C
      RETURN
      END
C
C
      REAL FUNCTION FRHF1*8(SP,Z)
      IMPLICIT REAL*8(A-H,O-Z)
C
      Z2=Z*Z
      Z3=Z2*Z
C
      F1=16.D0*(Z-1.D0)*SP*SP
      FRHF1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0
C
      RETURN
      END
