C
C *** PHASE AND GROUP VELOCITY AND 
C     AMPLITUDE SPECTRUM DUE TO VERTICAL POINT-FOORCE ON FREE SURFACE 
C     OF RAYLEIGH WAVE
C
      PROGRAM MRKIN
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 (LAYERS NUMBER)
C MOB   : NUMBER OF OBSERVATION POINTS
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 RAYLEIGH WAVE
C NEMD  : END OUTPUT MODE OF RAYLEIGH 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 RAYLEIGH WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : RAYLEIGH WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): OUTPUT HIGHEST RAYLEIGH MODE NUMBER AT Ith OMEGA
C XOB(I)  : X-AXIS OF OBSERBATION POINTS ON THE BASIN
C UU(I,J) : GROUP VELOCITY OF Jth MODE RAYLEIGH WAVE AT Ith OMEGA
C ARK(I,J): AMPLITUDE SPECTRUM (AL/K) OF Jth MODE RAYLEIGH WAVE
C           AT Ith OMEGA
C UW0(I,J): AMPLITUDE RATIO ( IMAG(RK/RL) ) OF Jth MODE RAYLEIGH WAVE
C           AT Ith OMEGA
C U0(I)   : APPROXIMATE RAYLEIGH WAVE(HORIZONTAL DISPLACEMENT
C           ON FREE-SURFACE) IN A SEDIMENTARY BASIN 
C           DUE TO 2-D SV HASKELL'S FAULT MODEL
C W0(I)   : APPROXIMATE RAYLEIGH WAVE(VERTICAL DISPLACEMENT
C           ON FREE-SURFACE) IN A SEDIMENTARY BASIN 
C           DUE TO 2-D SV HASKELL'S FAULT MODEL
C
      PARAMETER(MOM=10,MMD=2,ML=2,MOB=2)
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*15
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM),XOB(MOB)
     *         ,UU(MOM,0:MMD),ARK(MOM,0:MMD),UW0(MOM,0:MMD)
     *         ,DNS(ML),VP(ML),VS(ML),DEP(ML) 
     *         ,DGS(ML),DHS(ML),DRS(ML),DSS(ML)
     *         ,DGSZ(ML),DHSZ(ML),DRSZ(ML),DSSZ(ML)
      COMPLEX*16 U0(MOM),W0(MOM),UO(MOM,MOB),WO(MOM,MOB)
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
      WRITE(*,1010)
 1000 FORMAT(/' From MRKIN :' 
     *   /' Cal. of C,U and W0 & U0 of Rayleigh wave in multi-layers'
     *   /'  : W0 and U0 are displacements on the free-surface' 
     *   /'    in a sedimentary basin') 
 1010 FORMAT('    They are generated by 2-D Haskell'' fault model'
     *   /'    in the bedrock'
     *   /'    Basin egde is located at (0,0)')
C
C *** DATA INPUT ***
C
      CALL RINPUT(DOM,NOM,ML,NL,DNS,VP,VS,DEP,NDOM,NSMD,NEMD,NDC)
C
C *** PARAMETER CHECK ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MRKIN :'
        WRITE(*,*)'MOM < NOM THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(ML.LT.NL)THEN
        WRITE(*,*)'FROM MAIN MRKIN :'
        WRITE(*,*)'ML < NL THUS STOP !!!'
        WRITE(*,*)'ML,NL =',ML,NL
        STOP
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MRKIN :'
        WRITE(*,*)'MMD < NEMD THUS STOP !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        STOP
      END IF
C
C *** INPUT OF OBSERVATION POINTS ON THE BASIN ***
C
      WRITE(*,*)
    1 WRITE(*,2010)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-enter NOB : MOB=',MOB
        GOTO  1
      END IF
      IF(NOB.EQ.0)GOTO  1
C
      DO 100 I=1,NOB
        WRITE(*,2020)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 2010 FORMAT( ' Enter Number of Observation Points                : ',\)
 2020 FORMAT( ' Enter X-axis of Obsevation Points of Num. =',I3,' : ',\)

C
C *** LAYERS AND INTERGRATION DATA ***
C
   10 WRITE(*,3010)NL
      READ(*,*)NFSL
      IF(NFSL.EQ.0)GOTO 10
      IF(NFSL.GT.NL)GOTO 10
   20 WRITE(*,3020)
      READ(*,*)NFNL
      IF(NFNL.EQ.0)GOTO 20
      IF(NFNL.GT.NL)GOTO 20
   30 WRITE(*,3030)
      READ(*,*)DDEP
      IF(DDEP.EQ.0.D0)GOTO 30
   40 WRITE(*,3040)
      READ(*,*)NNL
      IF(NNL.EQ.0)GOTO 40
C
 3010 FORMAT(/' Layers number =',I5,
     *       /' Enter first (top) layer number    : ',\)
 3020 FORMAT( ' Enter final (bottom) layer number : ',\)
 3030 FORMAT( ' Enter delta thickness in bottom layer'
     *       /'       for numerical integration   : ',\)
 3040 FORMAT( ' Enter number of delta depth       : ',\)
C
C *** CUT-OFF OMEGA ***
C
      CALL RAYCTN(DOM,NOM,ML,NL,DNS,VP,VS,DEP,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITY ***
C
      CALL RAYCCN(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *           ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      CALL RAYKIN(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *           ,CC,NMODE,CFOM,NSMD,MMD,UU,ARK,UW0
     *           ,NFSL,NFNL,DDEP,NNL,U0,W0
     *           ,DGS,DHS,DRS,DSS,DGSZ,DHSZ,DRSZ,DSSZ)

C
C *** CAL.OF UO AND WO ***
C
      DO 200 K=1,NOM
C
        OM=DOM*DFLOAT(K)
C
        DO 210 IOB=1,NOB
C
          WO(K,IOB)=(0.D0,0.D0)
          UO(K,IOB)=(0.D0,0.D0)
C
          DO 220 MODE=NSMD,NCFOM
C
            IF(CC(K,MODE).EQ.0.D0)THEN
              EX=(0.D0,0.D0)
            ELSE
              EX=CDEXP(DCMPLX(0.D0,-OM*XOB(IOB)/CC(K,MODE)))
            END IF
C
            UO(K,IOB)=UO(K,IOB)+U0(K)*EX
            WO(K,IOB)=WO(K,IOB)+W0(K)*EX
C
  220     CONTINUE
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** WAVE DATA OUTPUT TO FILE ***
C
      WRITE(*,*)
   50 WRITE(*,1200)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 50
      WRITE(*,1210)
 1200 FORMAT(/' Enter name of wave data output file for PLOTMS : ',\)
 1210 FORMAT( '  : Data format for PLOTMS 1:OB => 2:UO => 3:WO roop')
C      
      OPEN(7,FILE=FNAME)
C
        WRITE(7,1300)DOM,NOM
        WRITE(7,1310)
C
        DO 300 IOB=1,NOB
C
          DO 310 K=1,NOM,5
            WRITE(7,1320)(UO(K+J,IOB),J=0,4)
  310     CONTINUE
C
          DO 320 K=1,NOM,5
            WRITE(7,1320)(WO(K+J,IOB),J=0,4)
  320     CONTINUE
C
  300   CONTINUE
C
      CLOSE(7)
C
 1300 FORMAT(F15.8,I15)
 1310 FORMAT(///)
 1320 FORMAT(10F8.5)
C
C *** DATA OUTPUT ON SCREEN ***
C
      WRITE(*,1400)
      READ(*,*)NYES
      IF(NYES.NE.0)GOTO 60
 1400 FORMAT(/' Output Rayleigh wave data on screen ? (yes=0) :',\) 
C
      CALL RAYOTN(DOM,NOM,MOM,ML,NL,VS
     *           ,CFOM,NCFOM,CC,NSMD,MMD,UU,ARK,UW0)
C
C *** DATA OUTPUT ON PLOTTER ***
C
   60 WRITE(*,1500)
      READ(*,*)NYES
      IF(NYES.NE.0)GOTO 70
 1500 FORMAT(/' Output Rayleigh wave data to plotter ? (yes=0) :',\) 
C
      CALL RAYPTN(DOM,NOM,MOM,ML,NL,VS
     *           ,CFOM,NCFOM,CC,NSMD,MMD,UU,ARK,UW0)
C
   70 STOP
      END
