C
      PROGRAM MRICKF
C
      PARAMETER(MOM=100,MWN=10000,MOB=6)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VS1,VS2,SM1,SM2,WB1,WB2,SS,FNCF,FOUR(MWN),VV(MOM,MOB)
      DIMENSION NOB(MOB)
      COMMON /DATA/ OM,WB1,WB2,SM1,SM2,THK
      CHARACTER FNAME*15
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
 1000 FORMAT(/' From MRICKF : SOURCE POINT is located at (0,0)'
     * /' EXACT RICKER WAVE IN 2-LAYER SH MEDIUM BY FFT'
     * /' SOURCE: P=RICKER*SM1')
C
      Q1=100.D0
      Q2=100.D0
      DN1=2.0D0
      DN2=2.5D0
      V1=1000.D0
      V2=3000.D0
      THK=2500.D0
      WRITE(*,1010)Q1,DN1,V1,THK
      WRITE(*,1020)Q2,DN2,V2
 1010 FORMAT(/' Q1,DN1,V1,THK=',4F8.1)
 1020 FORMAT( ' Q2,DN2,V2    =',4F8.1)
C
C *** DATA INPUT ***
C
   10 WRITE(*,*)'OMEGA NUMBER?'
      READ(*,*)NOM
C
      IF(MOM.LT.NOM)THEN
        WRITE(*,*)'FROM MAIN MRICKF'
        WRITE(*,*)'MOM < NOM RE-ENTRER !!!'
        WRITE(*,*)'MOM,NOM =',MOM,NOM
        GOTO 10
      END IF
C
      WRITE(*,*)'FIRST OM?'
      READ(*,*)FOM
      WRITE(*,*)'DELTA OM?'
      READ(*,*)DOM
C
      WRITE(*,1210)
      READ(*,*)SLL
      DWN=6.2831853D0/SLL
   19 WRITE(*,1220)
      READ(*,*)NFT
      NFFT=2**NFT
      IF(NFFT.GT.MWN)GOTO 19
      DX=SLL/DFLOAT(NFFT)
      WRITE(*,1230)DWN,DWN*(NFFT/2+1),NFFT,DX
C
 1210 FORMAT(/' Enter X-length for FFT    : ',\)
 1220 FORMAT( ' Enter NFT for FFT (2**NFT): ',\)
 1230 FORMAT(/' Delta Wave Number =',D12.4,', Final Wave Number =',D12.4,
     *       /' Num. of FFT Data  =',I12,', Delta X           =',D12.4)
C
   20 WRITE(*,1110)
      READ(*,*)NOBT
      IF(NOBT.GT.MOB)THEN
        WRITE(*,*)'NOBT > MOB !!! : Re-Enter NOB : MOB=',MOB
        GOTO 20
      END IF
      IF(NOBT.EQ.0)GOTO 20
C
      DO 100 I=1,NOBT
   21   WRITE(*,1120)I
        READ(*,*)NOB(I)
        OB=NOB(I)*DX
        WRITE(*,1130)OB
        READ(*,*)IYES
        IF(IYES.NE.0)GOTO 21
  100 CONTINUE
C
 1110 FORMAT(/' Enter Number of Observation Point: ',\)
 1120 FORMAT( ' Enter NX (Dx*NX= X-axis of OP)=',I3,': ',\)
 1130 FORMAT( ' X =',F10.3,' Is this X OK ? (yes=0)',\)
C
      WRITE(*,1221)
      READ(*,*)TP
C
 1221 FORMAT(/' Enter Tp(sec) of Ricker wavelet : ',\)
C
C *** DATA CAL ***
C
      VS1=V1*DCMPLX(1.D0,0.5D0/Q1)
      VS2=V2*DCMPLX(1.D0,0.5D0/Q2)
      SM1=DN1*VS1*VS1
      SM2=DN2*VS2*VS2
C
C *** GREEN FUNC ***
C
      AA=DN1*V1*V1
      BB=TP*TP*0.0253302D0
C
      DO 210 IOM=1,NOM
C
        OM=FOM+DOM*DFLOAT(IOM-1)
        OM2=OM*OM
        WB1=OM2/(VS1*VS1)
        WB2=OM2/(VS2*VS2)
C
        DO 220 K=0,NFFT/2+1
          WN=DWN*DFLOAT(K)
          FOUR(K+1)=FNCF(WN)
  220   CONTINUE
C
        DO 230 I=2,NFFT/2
          FOUR(NFFT+2-I)=DCMPLX(DREAL(FOUR(I)),-DIMAG(FOUR(I)))
  230   CONTINUE
C
        CALL FFT(NFFT,FOUR,MWN,-1)
C
        CC=AA*OM2*DEXP(-OM2*BB)
        DO 240 IOB=1,NOBT
          VV(IOM,IOB)=CC*FOUR(NOB(IOB)+1)/SLL
  240   CONTINUE
C
  210 CONTINUE
C
C *** OUTPUT OF DATA ***
C
        WRITE(*,1400)DOM,NOM
        DO 300 IOB=1,NOBT
        WRITE(*,1410)NOB(IOB)*DX
          DO 310 K=1,NOM,4
            WRITE(*,1420)(VV(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
C
 1400 FORMAT(/' DOM,NOM=',F6.2,I5)
 1410 FORMAT(/' X=',F10.2)
 1420 FORMAT(8D10.3)
C
C *** OUTPUT OF DATA FILE ***
C
      WRITE(*,*)
   30 WRITE(*,1310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 30
      WRITE(*,1320)
 1310 FORMAT( ' Enter file name : ',\)
 1320 FORMAT( '  : Data format for PLOTMS 1:OB => 2:VV roop')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,1500)DOM,NOM
        WRITE(8,1510)
        DO 500 IOB=1,NOBT
          DO 510 K=1,NOM,5
            WRITE(8,1520)(VV(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 1500 FORMAT(F15.8,I15)
 1510 FORMAT(///)
 1520 FORMAT(10F8.5)
C
      STOP
      END
