C
      PROGRAM MRICKR
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 (LAYER NUMBER)
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 UU(I,J) : GROUP VELOCITY OF Jth MODE RAYLEIGH WAVE AT Ith OMEGA
C
      PARAMETER(MOM=50,MMD=7,ML=10,MOB=5)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 FW0(MOM,MOB),FU0(MOM,MOB),CPX,EKX,FWW,FUU,CCC
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD),NMODE(MOM)
     *         ,UU(MOM,0:MMD),XOB(MOB)
     *         ,ARWZ(MOM,0:MMD),ARWX(MOM,0:MMD)
     *         ,ARUZ(MOM,0:MMD),ARUX(MOM,0:MMD)
     *         ,DNS(ML),VP(ML),VS(ML),THK(ML) 
      COMMON /SODAT/ SDEP,LS,ODEP,LO
      CHARACTER FNAME*15
C
      WRITE(*,1000)
 1000 FORMAT(' MRICKR: RICKER WAVE IN 2-D IN-PLANE LAYERED'/
     *       '         HALF-SPACE USING NORMAL MODE SOLUTION'/
     *       '        (SP AND OB ARE IN ANY LAYER),')
C
C *** DATA INPUT ***
C
      CALL RAYINN(DOM,NOM,ML,NL,DNS,VP,VS,THK,NDOM,NSMD,NEMD,NDC)
C
C *** PARAMETER CHECK ***
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MRCUAV :'
        WRITE(*,*)'MOM < NOM THUS STOP !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        STOP
      END IF
C
      IF(ML.LT.NL)THEN
        WRITE(*,*)'FROM MAIN MRCUAV :'
        WRITE(*,*)'ML < NL THUS STOP !!!'
        WRITE(*,*)'ML,NL =',ML,NL
        STOP
      END IF
C
      IF(MMD.LT.NEMD)THEN
        WRITE(*,*)'FROM MAIN MRCUAV :'
        WRITE(*,*)'MMD < NEMD THUS STOP !!!'
        WRITE(*,*)'MMD,NEMD =',MMD,NEMD
        STOP
      END IF
C
C
C *** INPUT OF SP AND OP DEPTH IN ANY LAYER ***
C
   10 WRITE(*,1005)NL
      READ(*,*)LS
      IF(LS.EQ.0)GOTO 10
      IF(LS.GT.NL)GOTO 10
      WRITE(*,1010)THK(LS)
      READ(*,*)SDEP
 1005 FORMAT(/' Total layers number =',I5,
     *       /' Enter source layer number (from top layer) : ',\)
 1010 FORMAT( ' Source layer thickness =',F10.4,
     *       /' Enter source depth from upper boundary : ',\)
C
   20 WRITE(*,1100)NL
      READ(*,*)LO
      IF(LO.EQ.0)GOTO 20
      IF(LO.GT.NL)GOTO 20
      WRITE(*,1110)THK(LO)
      READ(*,*)ODEP
 1100 FORMAT(/' Total layers number =',I5,
     *       /' Enter OP layer number (from top layer) : ',\)
 1110 FORMAT( ' Source layer thickness =',F10.4,
     *       /' Enter OP depth from upper boundary : ',\)
C
   30 WRITE(*,1410)
      READ(*,*)NOB
      IF(NOB.GT.MOB)THEN
        WRITE(*,*)'NOB > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 30
      END IF
      IF(NOB.EQ.0)GOTO 30
C
      DO 100 I=1,NOB
        WRITE(*,1420)I
        READ(*,*)XOB(I)
  100 CONTINUE
C
 1410 FORMAT(/' Enter Number of Observation Point : ',\)
 1420 FORMAT( ' Enter X-axis of Observation Point of Num. =',I3,' : ',\)
C
      WRITE(*,1521)
      READ(*,*)TP,TS
C
 1521 FORMAT(/' Enter Tp and Ts (sec) of Ricker wavelet : ',\)
C
      WRITE(*,1550)
      READ(*,*)QX,QZ
C
 1550 FORMAT(/' Enter QX and QZ : ',\)
C
      WRITE(*,1640)
      READ(*,*)IQF
      IF(IQF.EQ.0)THEN
        WRITE(*,1650)
        READ(*,*)QFL
      END IF
 1640 FORMAT(/' Include Quality Factor (yes=0)?   : ',\)
 1650 FORMAT( ' Enter Quality Factor for Rayleigh Wave: ',\)
C
C *** CUT-OFF OMEGA ***
C
       CALL RAYCTN(DOM,NOM,ML,NL,DNS,VP,VS,THK
     *            ,NDOM,CFOM,NEMD,MMD,NCFOM)
C
C *** PHASE VELOCITY ***
C
      CALL RAYCCN(DOM,NOM,MOM,ML,NL,DNS,VP,VS,THK
     *           ,NDC,CFOM,NSMD,NEMD,MMD,CC,NMODE)
C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      CALL RAYUAR(DOM,NOM,MOM,ML,NL,DNS,VP,VS,THK
     *           ,CC,NMODE,CFOM,NSMD,MMD,UU,ARWZ,ARWX,ARUZ,ARUX)
C
C *** RICKER WAVE ***
C
      IF(IQF.NE.0)QFL=9999.D0
      CPX=DCMPLX(0.D0,1.D0)
      AA=DNS(1)*VS(1)*VS(1)
      BB=TP*TP*0.0253302D0
C
      DO 200 IOB=1,NOB
C
        XX=XOB(IOB)
C
        DO 210 IOM=1,NOM
C
          OM=DOM*DFLOAT(IOM)
          OM2=OM*OM
          CCC=AA*OM2*DEXP(-OM2*BB)*CDEXP(DCMPLX(0.D0,-OM*TS))
C
          FW0(IOM,IOB)=(0.D0,0.D0)
          FU0(IOM,IOB)=(0.D0,0.D0)
C
          DO 220 MODE=0,NMODE(IOM)
C
            IF(CC(IOM,MODE).EQ.0.D0)THEN
              EKX=DCMPLX(1.D0,0.D0)
              EXQ=1.D0
            ELSE
              EKX=CDEXP(-CPX*OM/CC(IOM,MODE)*XX)
              EXQ=DEXP(-OM*XX/(2.D0*QFL*UU(IOM,MODE)))
            END IF
C
            IF(IQF.NE.0)EXQ=1.D0
C
            FWW=(CPX*ARWZ(IOM,MODE)*QZ+ARWX(IOM,MODE)*QX)*EKX*EXQ
            FUU=(ARUZ(IOM,MODE)*QZ+CPX*ARUX(IOM,MODE)*QX)*EKX*EXQ
C
            FW0(IOM,IOB)=FW0(IOM,IOB)+FWW
            FU0(IOM,IOB)=FU0(IOM,IOB)+FUU
C
  220     CONTINUE
C
          FW0(IOM,IOB)=FW0(IOM,IOB)*CCC
          FU0(IOM,IOB)=FU0(IOM,IOB)*CCC
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,2400)DOM,NOM
        DO 300 IOB=1,NOB
        WRITE(*,2410)XOB(IOB)
          DO 310 K=1,NOM,4
            WRITE(*,2420)(FW0(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
        DO 320 IOB=1,NOB
          WRITE(*,2415)XOB(IOB)
          DO 330 K=1,NOM,4
            WRITE(*,2420)(FU0(K+J,IOB),J=0,3)
  330     CONTINUE
  320   CONTINUE
C
 2400 FORMAT(/' DOM,NOM=',F6.2,I5)
 2410 FORMAT(/' X=',F10.2,'; FW0*4')
 2415 FORMAT(/' X=',F10.2,'; FU0*4')
 2420 FORMAT(8D10.3)
C
C *** OUTPUT OF DATA FILE ***
C
      WRITE(*,*)
   50 WRITE(*,2310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 50
      WRITE(*,2320)
 2310 FORMAT( ' Enter file name : ',\)
 2320 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VV roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,2500)DOM,NOM
        WRITE(8,2510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,5
            WRITE(8,2520)(FW0(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
        DO 520 IOB=1,NOB
          DO 530 K=1,NOM,5
            WRITE(8,2520)(FU0(K+J,IOB),J=0,4)
  530     CONTINUE
  520   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(F15.8,I15)
 2510 FORMAT(///)
 2520 FORMAT(10F8.5)
C
      STOP
      END
