CC *** PHASE & GROUP VELOCITIES 
C     AND AMPLITUDE SPECTRA OF RAYLEIGH WAVES 
C     DUE TO VERTICAL & HORIZONTAL POINT-FORCES 
C     IN A MULTI-LAYERD MEDIUM ***
C     
      PROGRAM MRCUA
C
C ** MARAMETERS **
C MOM   : NUMBER OF DELTA OMEGAS (FOR NOM)
C MMD   : NUMBER OF RAYLEIGH MODES (FOR NEMD)
C ML    : NUMBER OF SOIL LAYERS
C ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY (DELTA OMEGA: RAD/SEC)
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCIES
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 BETWEEN MIN VS AND MAX VS FOR PHASE VELOCITY
C NSMD  : FIRST RAYLEIGH MODE TO OUTPUT
C NEMD  : FINAL RAYLEIGH MODE TO OUTPUT
C NDOM  : DIVISION NUMBER BETWEEN 0 AND MAX OMEGA FOR CUT-OFF OMEGAS
C NLS   : SOURCE LOCATED LAYER
C SDEP  : SOURCE DEPTH FROM THE UPPER BOUNDARY OF THE SOURCE LAYER
C ** CALCULATED DATA **
C CFOM(I) : CUT-OFF OMEGA OF Ith RAYLEIGH MODE
C NCFOM   : NUMBER OF CUT-OFF OMEGAS
C CC(I,J) : PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): NUMBER OF RAYLEIGH MODES AT Ith OMEGA
C UU(I,J) : GROUP VELOCITY OF Jth MODE AT Ith OMEGA
C UW0(I,J): AMPLITUDE RATIO (IMAG(RK/RL)) OF Jth MODE AT Ith OMEGA
C ARV(I,J): AMPLITUDE SPECTRUM (DISPLACEMENT) DUE TO VERTICAL UNIT POINT-FOURCE
c           OF Jth MODE AT Ith OMEGA (REAL)
C ARH(I,J): AMPLITUDE SPECTRUM (DISPLACEMENT) DUE TO HORIZONTAL 
C           UNIT POINT-FOURCE OF Jth MODE AT Ith OMEGA (REAL)
C
      PARAMETER(MOM=200,MMD=5,ML=35)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),UW0(MOM,0:MMD)
     *         ,ARH(MOM,0:MMD),ARV(MOM,0:MMD)
     *         ,DNS(ML),VP(ML),VS(ML),DEP(ML)
C
C *** EXPLANATION OF THIS PROGRAM ***
C
      WRITE(*,1000)
C
 1000 FORMAT(/' From MRCUA (Main Program) :'
     */'   Phase & Group Velocities and Amplitude Spectra'
     */'   of Rayleigh Waves Due to Vertical & Horizontal Point Forces'
     */'   in a multi-layerd Medium (Programed by Y.HISADA)')
C
C *** INPUT OF DATA ***
C
      CALL RINPUT(DOM,NOM,ML,NL,DNS,VP,VS,DEP,NDOM,NSMD,NEMD,NDC)
C
C *** DATA CHECK ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MRCUA :'
        WRITE(*,*)'MOM < NOM (OMEGGA NUMBER), THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(ML.LT.NL)THEN
        WRITE(*,*)'FROM MAIN MRCUA :'
        WRITE(*,*)'ML < NL (SOIL LAYER NUMBER), THUS STOP !!!'
        WRITE(*,*)'ML,NL =',ML,NL
        STOP
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MRCUA :'
        WRITE(*,*)'MMD < NEMD (MODE NUMBER), THUS STOP !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        STOP
      END IF
C
C *** ENTER OF SOURCE LAYER NUMBER AND SOURCE DEPTH ***
C
   20 WRITE(*,1010)NL
      READ(*,*,ERR=20)NLS
      IF(NLS.EQ.0)GOTO 20
      IF(NLS.GT.NL)GOTO 20
   30 WRITE(*,1020)DEP(NLS)
      READ(*,*,ERR=30)SDEP
      IF((NLS.NE.NL).AND.(SDEP.GT.DEP(NLS)))GOTO 30
C
 1010 FORMAT(/' From MRCUA : layers number =',I5,
     *       /' Enter source layer number (from top layer) : ',\)
 1020 FORMAT( ' Source layer thickness     =',F10.4,
     *       /' Enter source depth from upper boundary     : ',\)
C
C *** CUT-OFF OMEGAS ***
C
      CALL RCUTOM(DOM,NOM,ML,NL,DNS,VP,VS,DEP,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITIES ***
C
      CALL RPHASE(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *           ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRA ***
C
      CALL RGRAMP(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *           ,CC,NMODE,CFOM,NSMD,MMD
     *           ,NLS,SDEP,UU,UW0,ARV,ARH)
C
C *** DATA OUTPUT ON SCREEN ***
C
      CALL ROUT(DOM,NOM,MOM,ML,NL,VS
     *         ,CFOM,NCFOM,CC,NSMD,MMD,UU,ARV,ARH,UW0)
C
C *** DATA OUTPUT ON PLOTTER ***
C
      CALL RPLOTV(DOM,NOM,MOM,ML,NL,VS
     *           ,CFOM,NCFOM,CC,NSMD,MMD,UU,ARV,UW0)
      CALL RPLOTH(DOM,NOM,MOM,CFOM,NCFOM,NSMD,MMD,ARH,UW0)
C
      STOP
      END
