C
      PROGRAM MRICKM
C
      PARAMETER(MOM=200,MOB=21)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VS1,VS2,SM1,SM2,WB1,WB2,S1,S2,VV(MOM,MOB)
     *          ,FNC1,FNC2,FNC3,FNC4,FNC5,FNC6,CC
      DIMENSION XOB(MOB)
      COMMON /DATA/ X,WB1,WB2,SM1,SM2,THK,SDEP,ODEP
      CHARACTER FNAME*15
      EXTERNAL FNC1,FNC2,FNC3,FNC4,FNC5,FNC6
C
C *** HEAD ***
C
      WRITE(*,*)
      WRITE(*,1000)
 1000 FORMAT(/' From MRICKM: '
     * /' EXACT GREEN FUNTION IN 2-LAYER SH MEDIUM BY RICKER WAVE'
     * /' 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 M'
        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(*,1110)
      READ(*,*)C1
      WRITE(*,1115)
      READ(*,*)C2
      WRITE(*,1120)
      READ(*,*)NS1
      WRITE(*,1125)
      READ(*,*)NS2
C
 1110 FORMAT(/' Enter the lowest phase vel. C1        : ',\)
 1115 FORMAT( ' Enter the second phase vel. C2(>C1)   : ',\)
 1120 FORMAT( ' Enter div. num. for int. (0-OM/C2)    : ',\)
 1125 FORMAT( ' Enter div. num. for int. (OM/C2-OM/C1): ',\)
C
C *** ENTER OF SOURCE DEPTH IN ANY LAYER ***
C
   20 WRITE(*,1210)2
      READ(*,*)MLS
      IF(MLS.EQ.0)GOTO 20
      IF(MLS.GT.2)GOTO 20
      WRITE(*,1220)THK
      READ(*,*)SDEP
 1210 FORMAT(/' From MRICKM : layers number =',I5,
     *       /' Enter source layer number (from top layer) : ',\)
 1220 FORMAT( ' Surface layer depth =',F10.4,
     *       /' Enter source depth from upper boundary : ',\)
C
C *** ENTER OF OBS. DEPTH IN ANY LAYER ***
C
   25 WRITE(*,1310)2
      READ(*,*)MLO
      IF(MLO.EQ.0)GOTO 25
      IF(MLO.GT.2)GOTO 25
      WRITE(*,1320)THK
      READ(*,*)ODEP
 1310 FORMAT(/' From MRICKM : layers number =',I5,
     *       /' Enter obs. layer number (from top layer) : ',\)
 1320 FORMAT( ' Surface layer depth =',F10.4,
     *       /' Enter obs. 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
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 200 IOB=1,NOB
C
        X=XOB(IOB)
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
          WNA=0.D0
          WNB=OM/C2
          WNC=OM/C1
          IF(OM.EQ.0.D0)WRITE(*,*)'OM=0 !!!!!!!!!!!!!!!!!!!!'
C
          IF(MLS.EQ.1)THEN
            IF(MLO.EQ.1)THEN
              IF(ODEP.LT.SDEP)THEN
                CALL NEWTON(FNC1,WNA,WNB,NS1,S1,IER)
                IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
                CALL NEWTON(FNC1,WNB,WNC,NS2,S2,IER)
                IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
              ELSE
                CALL NEWTON(FNC2,WNA,WNB,NS1,S1,IER)
                IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
                CALL NEWTON(FNC2,WNB,WNC,NS2,S2,IER)
                IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
              END IF
            ELSE
              CALL NEWTON(FNC3,WNA,WNB,NS1,S1,IER)
              IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
              CALL NEWTON(FNC3,WNB,WNC,NS2,S2,IER)
              IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
            END IF
          ELSE
            IF(MLO.EQ.1)THEN
              CALL NEWTON(FNC4,WNA,WNB,NS1,S1,IER)
              IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
              CALL NEWTON(FNC4,WNB,WNC,NS2,S2,IER)
              IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
            ELSE
              IF(ODEP.LT.SDEP)THEN
                CALL NEWTON(FNC5,WNA,WNB,NS1,S1,IER)
                IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
                CALL NEWTON(FNC5,WNB,WNC,NS2,S2,IER)
                IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
              ELSE
                CALL NEWTON(FNC6,WNA,WNB,NS1,S1,IER)
                IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
                CALL NEWTON(FNC6,WNB,WNC,NS2,S2,IER)
                IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
              END IF
            END IF
          END IF
C
          CC=AA*OM2*DEXP(-OM2*BB)*CDEXP(DCMPLX(0.D0,-OM*TS))
          VV(IOM,IOB)=CC*(S1+S2)/3.1415927D0
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)(VV(K+J,IOB),J=0,3)
  310     CONTINUE
  300   CONTINUE
C
 2400 FORMAT(/' DOM,NOM=',F6.2,I5)
 2410 FORMAT(/' X=',F10.2)
 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)(VV(K+J,IOB),J=0,4)
  510     CONTINUE
  500   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(F15.8,I15)
 2510 FORMAT(///)
 2520 FORMAT(10F8.5)
C
      STOP
      END
