C
C *** RICKER WAVE IN 2-D IN-PLANE LAYERED HALF-SPACE
C     USING EXACT GREEN FUNCTION
C     (SOURCE IS (0,SDEP IN ANY LAYER), AND OB IS (XOB,0) ***
C
      PROGRAM MRICK1
C
      PARAMETER(MOM=100,ML=30,MOB=6)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML),FW0(MOM,MOB),FU0(MOM,MOB)
     *          ,FW1,FU1,FW2,FU2,AA,CC
      DIMENSION DNS(ML),THK(ML),XOB(MOB)
      COMMON / CN / C1,N1,C2,N2 /FDAT/ OM,NL,SDEP,LS,QX,QZ
      CHARACTER FNAME*15
C
      WRITE(*,1000)
 1000 FORMAT(' MRICK1: RICKER WAVE IN 2-D IN-PLANE LAYERED'/
     *       '         HALF-SPACE USING EXACT GREEN FUNCTION'/
     *       '        (SOURCE IS (0,SDEP IN ANY LAYER),'/
     *       '         AND OB IS (XOB,0)')
C
      CALL DATAIN(DOM,NOM,ML,NL1,DNS,VP,VS,THK,LS1,SDEP1)
      NL=NL1
      LS=LS1
      SDEP=SDEP1
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
C *** GREEN FUNC ***
C
      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
C
          WNB=OM/C2
          WNA=WNB/(1000.D0*N1)
          WNC=OM/C1
C
          CALL NEWTN1(WNA,WNB,N1,FW1,FU1,XX,DNS,VP,VS,THK,ML,IER)
          IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
          CALL NEWTN1(WNB,WNC,N2,FW2,FU2,XX,DNS,VP,VS,THK,ML,IER)
          IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          CC=AA*OM2*DEXP(-OM2*BB)*CDEXP(DCMPLX(0.D0,-OM*TS))
          FW0(IOM,IOB)=CC*(FW1+FW2)*DCMPLX(0.D0,-0.3183198D0)
          FU0(IOM,IOB)=CC*(FU1+FU2)*DCMPLX(0.D0,+0.3183198D0)
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
