C
C *** RICKER WAVES IN THE 3D LAYERED HALF-SPACE USING
C     THE EXACT GREEN FUNCTION BASED ON THE RECIPROCITY RELATION ***
C
      PROGRAM MRICK
C
      PARAMETER(MOM=100,ML=30,MOB=6)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML),AA,CC,CPX
     *          ,GQ(MOM,MOB),GV(MOM,MOB),GW(MOM,MOB)
     *          ,QX1,QX2,QZ1,QZ2,VX1,VX2,WX1,WX2,WZ1,WZ2
      DIMENSION DNS(ML),THK(ML),ROB(MOB)
      COMMON / CN / C1,N1,C2,N2 /FDAT/ OM,NL,RR
      CHARACTER FNAME*15
C
      WRITE(*,1000)
 1000 FORMAT(' MRICKN: RICKER WAVES IN THE 3D LAYERED HALF-SPACE'/
     *       '   USING THE EXACT GREEN FUNCTION BASED ON'/
     *       '   THE RECIPLOCITY RELATION '
     *       '  (SP AND OB ARE IN ANY LAYER),')
C
      CALL DATIN2(DOM,NOM,ML,NL1,DNS,VP,VS,THK,LS,SDEP,LO,ODEP)
      NL=NL1
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(*,*)ROB(I)
  100 CONTINUE
      WRITE(*,1430)
      READ(*,*)RAD
      RAD=RAD*0.01745329252D0
C
 1410 FORMAT(/' Enter Number of OB Point               : ',\)
 1420 FORMAT( ' Enter R-axis of OB Point of Num. =',I3,'  : ',\)
 1430 FORMAT(/' Enter Angle of OB from X-axis (degree) : ',\)
C
      WRITE(*,1521)
      READ(*,*)TP,TS
C
 1521 FORMAT(/' Tp(Dominant Period) and Ts(Delay Time) '
     *         'of Ricker Wavelet(sec)?: ',\)
C
      WRITE(*,1550)
      READ(*,*)QX,QY,QZ
C
 1550 FORMAT(/' Enter QX, QY and QZ of SP: ',\)
C
      AA=DNS(1)*VS(1)*VS(1)
      WRITE(*,1551)AA
      WRITE(*,1552)  
      READ(*,*)POW
      AA=AA**POW
      WRITE(*,1553)AA
C
 1551 FORMAT(/' Amplification Factor (AA) = ',2D12.4)
 1552 FORMAT( ' AA**N; N=?: ',\)
 1553 FORMAT(/' New Amplification Factor  = ',2D12.4)
C
C *** GREEN FUNC ***
C
      CPX=(0.D0,1.D0)
      BB=TP*TP*0.0253302D0
      CR=DCOS(RAD)
      SR=DSIN(RAD)
      WRITE(*,*)
C
      DO 200 IOB=1,NOB
C
        RR=ROB(IOB)
C
        DO 210 IOM=1,NOM
C
          OM=DOM*DFLOAT(IOM)
          OM2=OM*OM
          WRITE(*,*)' OM =',OM
C
          WNB=OM/C2
          WNA=WNB/(1000.D0*N1)
          WNC=OM/C1
C
          IF((LO.LT.LS).OR.((LO.EQ.LS).AND.(ODEP.LE.SDEP)))THEN
C
            CALL QVW(WNA,WNB,N1,DNS,VP,VS,THK,ML,IER,LS,SDEP,LO,ODEP
     *              ,1,QX1,QZ1,VX1,WX1,WZ1)
            IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
            CALL QVW(WNB,WNC,N2,DNS,VP,VS,THK,ML,IER,LS,SDEP,LO,ODEP
     *              ,1,QX2,QZ2,VX2,WX2,WZ2)
            IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          ELSE
C
            CALL QVW(WNA,WNB,N1,DNS,VP,VS,THK,ML,IER,LO,ODEP,LS,SDEP
     *              ,2,QX1,QZ1,VX1,WX1,WZ1)
            IF(IER.NE.0)WRITE(*,*)' 1:IOB,IOM,IER',IOB,IOM,IER
            CALL QVW(WNB,WNC,N2,DNS,VP,VS,THK,ML,IER,LO,ODEP,LS,SDEP
     *              ,2,QX2,QZ2,VX2,WX2,WZ2)
            IF(IER.NE.0)WRITE(*,*)' 2:IOB,IOM,IER',IOB,IOM,IER
C
          END IF
C
          CC=AA*OM2*DEXP(-OM2*BB)*CDEXP(DCMPLX(0.D0,-OM*TS))
          GQ(IOM,IOB)=CC*((QX1+QX2)*(CR*QX+SR*QY)+(QZ1+QZ2)*QZ)
          GV(IOM,IOB)=CC*((VX1+VX2)*(SR*QX-CR*QY))
          GW(IOM,IOB)=CC*((WX1+WX2)*(CR*QX+SR*QY)+(WZ1+WZ2)*QZ)
C
  210   CONTINUE
C
  200 CONTINUE
C
C *** OUTPUT OF DATA (GQ, GV AND GW)***
C
        WRITE(*,2400)DOM,NOM
        DO 300 IOB=1,NOB
        WRITE(*,2410)ROB(IOB)
          DO 310 K=1,NOM,2
            AA=GQ(K  ,IOB)
            CC=GQ(K+1,IOB)
            DAA=CDABS(AA)
            DCC=CDABS(CC)
            WRITE(*,2420)AA,DAA,CC,DCC
  310     CONTINUE
  300   CONTINUE
        DO 320 IOB=1,NOB
          WRITE(*,2415)ROB(IOB)
          DO 330 K=1,NOM,2
            AA=GV(K  ,IOB)
            CC=GV(K+1,IOB)
            DAA=CDABS(AA)
            DCC=CDABS(CC)
            WRITE(*,2420)AA,DAA,CC,DCC
  330     CONTINUE
  320   CONTINUE
        DO 340 IOB=1,NOB
          WRITE(*,2417)ROB(IOB)
          DO 350 K=1,NOM,2
            AA=GW(K  ,IOB)
            CC=GW(K+1,IOB)
            DAA=CDABS(AA)
            DCC=CDABS(CC)
            WRITE(*,2420)AA,DAA,CC,DCC
  350     CONTINUE
  340   CONTINUE
C
 2400 FORMAT(/' DOM,NOM=',F6.2,I5)
 2410 FORMAT(/' X=',F10.2,'; GQ*2 (Real, Imag, Abs)')
 2415 FORMAT(/' X=',F10.2,'; GV*2 (Real, Imag, Abs)')
 2417 FORMAT(/' X=',F10.2,'; GW*2 (Real, Imag, Abs)')
 2420 FORMAT(2(3D10.3,3X))
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: GQ, GV and GW')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,2500)DOM,NOM
        WRITE(8,2510)
        DO 500 IOB=1,NOB
          DO 510 K=1,NOM,4
            WRITE(8,2520)(GQ(K+J,IOB),J=0,3)
  510     CONTINUE
  500   CONTINUE
        DO 520 IOB=1,NOB
          DO 530 K=1,NOM,4
            WRITE(8,2520)(GV(K+J,IOB),J=0,3)
  530     CONTINUE
  520   CONTINUE
        DO 540 IOB=1,NOB
          DO 550 K=1,NOM,4
            WRITE(8,2520)(GW(K+J,IOB),J=0,3)
  550     CONTINUE
  540   CONTINUE
      CLOSE(8)
C
 2500 FORMAT(F15.8,I15)
 2510 FORMAT(///)
 2520 FORMAT(8D10.3)
C
C *** OUTPUT OF DATA (GX, GY AND GZ)***
C
      WRITE(*,3000)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 9999
C
 3000 FORMAT(/' Output of GX, GY AND GZ? (YES=0): ',\)
C
        DO 600 IOB=1,NOB
          DO 610 K=1,NOM
            AA=GQ(K,IOB)
            CC=GV(K,IOB)
            GQ(K,IOB)=AA*CR-CC*SR
            GV(K,IOB)=AA*SR+CC*CR
  610     CONTINUE
  600   CONTINUE
C
        WRITE(*,3400)DOM,NOM
        DO 700 IOB=1,NOB
        WRITE(*,3410)ROB(IOB)
          DO 710 K=1,NOM,2
            AA=GQ(K  ,IOB)
            CC=GQ(K+1,IOB)
            DAA=CDABS(AA)
            DCC=CDABS(CC)
            WRITE(*,3420)AA,DAA,CC,DCC
  710     CONTINUE
  700   CONTINUE
        DO 720 IOB=1,NOB
          WRITE(*,3415)ROB(IOB)
          DO 730 K=1,NOM,2
            AA=GV(K  ,IOB)
            CC=GV(K+1,IOB)
            DAA=CDABS(AA)
            DCC=CDABS(CC)
            WRITE(*,3420)AA,DAA,CC,DCC
  730     CONTINUE
  720   CONTINUE
        DO 740 IOB=1,NOB
          WRITE(*,3417)ROB(IOB)
          DO 750 K=1,NOM,2
            AA=GW(K  ,IOB)
            CC=GW(K+1,IOB)
            DAA=CDABS(AA)
            DCC=CDABS(CC)
            WRITE(*,3420)AA,DAA,CC,DCC
  750     CONTINUE
  740   CONTINUE
C
 3400 FORMAT(/' DOM,NOM=',F6.2,I5)
 3410 FORMAT(/' X=',F10.2,'; GX*2 (Real, Imag, Abs)')
 3415 FORMAT(/' X=',F10.2,'; GY*2 (Real, Imag, Abs)')
 3417 FORMAT(/' X=',F10.2,'; GZ*2 (Real, Imag, Abs)')
 3420 FORMAT(2(3D10.3,3X))
C
      WRITE(*,*)
   55 WRITE(*,3310)
      READ(*,'(A)')FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 55
      WRITE(*,3320)
 3310 FORMAT( ' Enter file name : ',\)
 3320 FORMAT( '   Data format for PLOTMS 1: OB => 2: GX, GY and GZ')
C      
      OPEN(8,FILE=FNAME)
        WRITE(8,3500)DOM,NOM
        WRITE(8,3510)
        DO 800 IOB=1,NOB
          DO 810 K=1,NOM,4
            WRITE(8,3520)(GQ(K+J,IOB),J=0,3)
  810     CONTINUE
  800   CONTINUE
        DO 820 IOB=1,NOB
          DO 830 K=1,NOM,4
            WRITE(8,3520)(GV(K+J,IOB),J=0,3)
  830     CONTINUE
  820   CONTINUE
        DO 840 IOB=1,NOB
          DO 850 K=1,NOM,4
            WRITE(8,3520)(GW(K+J,IOB),J=0,3)
  850     CONTINUE
  840   CONTINUE
      CLOSE(8)
C
 3500 FORMAT(F15.8,I15)
 3510 FORMAT(///)
 3520 FORMAT(8D10.3)
C
 9999 STOP
      END
