C 
C Program GPMIG for 2-D Gaussian packet prestack depth migration
C
C Version: 6.00
C Date: 2006, May 30
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: zacek@karel.troja.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings, read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'SEP'/
C     'SEP'...Name of the file with input parameters.
C             Description of file SEP
C             If blank, default values of the corresponding data are
C             considered.
C Default: 'SEP'=' '.
C
C                                                     
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file
C     'sep.for'.
C Names of input files with the results of ray tracing from a receiver:
C     RAYALL='string'... Name of the input file CRT-R with the
C             quantities stored along rays (see C.R.T.5.5.1).
C             Default: RAYALL='r01.out' (usually sufficient)
C     INIALL='string'... Name of the input file CRT-I with the
C             quantities at the initial points of rays, corresponding
C             to file RAYALL (see C.R.T.6.1).
C             Default: INIALL='r01i.out' (usually sufficient)
C Names of input formatted files with quantities corresponding to the
C             common source, interpolated at the nodes of input grid:
C     NUM='string'...Name of the input file with N1*N2*N3 array
C             of integer values, containing the number of arrivals at
C             each gridpoint of the target grid.
C             Default: NUM='mtt-num.out'
C     MTT='string'...Name of the input file with the array of real
C             values, containing for each gridpoint the travel times
C             of all determined arrivals.
C             Default: MTT='mtt-tt.out'
C     AMP='string'...Name of the input file with the array of real
C             values, containing for each gridpoint the amplitudes.
C             Default: AMP='mtt-ap.out'
C Names of input formatted files with amplitudes of Gaussian packets:
C     GPAMPR='string'...Name of the input file with ...
C             Default: GPAMPR='gpampr.out'
C     GPAMPI='string'...Name of the input file with ...
C             Default: GPAMPR='gpampi.out'
C Names of output files:
C     GPMSEC='string'... Name of the output file containing
C             migrated section.
C             Default: GPMSEC='gpmsec.out'
C Data specifying the dimensions of the input grid:
C     O1=real, O2=real ... Coordinates of the origin of the grid.
C             Default: O1=0., O2=0.
C     D1=real... Grid interval along the X1 axis.
C             Default: D1=1.
C     D2=real... Grid interval along the X2 axis.
C             Default: D2=1.
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C
C Input unformatted file RAYALL:
C     See the description within source code file 'writ.for'.
C     Description of files CRT-R
C
C Input unformatted file INIALL:
C     See the description within source code file 'writ.for'.
C     Description of files CRT-I
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL AP00
C     AP00... File 'ap.for'.
      EXTERNAL RARRAY
C     RARRAY...File 'forms.for'.
      EXTERNAL RARRAI
C     RARRAI...File 'forms.for'.
      EXTERNAL ERROR
C     ERROR...File 'error.for'.
      EXTERNAL RSEP1
C     RSEP1...File 'sep.for'.
      EXTERNAL RSEP3T
C     RSEP3T...File 'sep.for'.
      EXTERNAL RSEP3I
C     RSEP3I...File 'sep.for'.
      EXTERNAL RSEP3R
C     RSEP3R...File 'sep.for'.

C-----------------------------------------------------------------------

C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc

C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.

      INCLUDE 'sep.inc'
C     sep.inc

C-----------------------------------------------------------------------

C     Auxiliary storage locations:
      REAL O1,O2,D1,D2,PI,KAPPA,X1,X2,R1,R2,R4,R4A,W1,W2,W4
      REAL HI(18),H(18),HUI(9),HH(3,3)
      REAL VI1,VI2,VI3,VI,V1,V2,V3,V
      REAL GPVAL,GPMIN,GPMAX,RFAR,RFAI,AK0,DETT
      REAL TR,GPAR,GPAI,DXR,DTR,DOMEGA,DPR,OXR,OTR,OOMEGA,OPR,DT,D
      COMPLEX M0(4,4),N(4,4),M(4,4),C,Q1,Q2,Q3,Q4,P1,P2,P3,P4,DET,N44
      COMPLEX GPIMAGE,GPEXP,GP1,GP2,GPA
      INTEGER I,J,K,L,IGRD
      INTEGER LU1,LU2,LU3,LU4,LU5,LU6,LU7,LU8,LU9,LU10,LU11
      INTEGER IRAM(MRAM)
      INTEGER ITR,NTR,IOMEGA,NOMEGA,N1,N2,N1N2,ISH
      INTEGER IXR,NXR,IPR,NPR,NXPT,NXPTO,IGRDA,IGRDB,IGRDC,IGRDG
      INTEGER IHISTORY,NHISTORY,IXTZL,IXTZR,IXTZU,IXTZD,IPROLD,IXROLD
      INTEGER IX1,IX2,IX1U,IX1D,IX2L,IX2R
      INTEGER INUM,IMTT,IAPR,IAPI,IGPMSEC,IGPAPTZ
      INTEGER IGPAR,IGPAI,IFGBR,IFGBI,IX1NEW,IX2NEW,IGRDNEW,ID
      LOGICAL TZONE
      PARAMETER (LU1=21,LU2=22,LU3=23,LU4=24,LU5=25,LU6=26,LU7=27)
      PARAMETER (LU8=28,LU9=29,LU10=30,LU11=31,LU12=32,LU13=33)
      PARAMETER (LU14=34,LU15=35)
      CHARACTER*80 FILSEP,RAYALL,INIALL,NUM,MTT,APR,API,GPK0
      CHARACTER*80 GPAMPR,GPAMPI,FGBY22,FGBR22,GPGRID,GPMSEC,GPAPTZ
      EQUIVALENCE (IRAM,RAM)
C      PARAMETER (GPMIN=1.0E-10)
      PARAMETER (GPMIN=0.0)
      PARAMETER (PI=3.141592654)

C.......................................................................

C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPMIG: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPMIG-01
        CALL ERROR('GPMIG-01: No input file specified')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
      WRITE(*,'(A)')'+GPMIG: Reading, calculating... '

C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)

C     Reading filenames of the files with computed rays and
C     multivalued travel times
      CALL RSEP3T('RAYALL',RAYALL,'r01.out')
      CALL RSEP3T('INIALL',INIALL,'r01i.out')
      CALL RSEP3T('NUM',NUM,'mtt-num.out')
      CALL RSEP3T('MTT',MTT,'mtt-tt.out')
      CALL RSEP3T('APR',APR,'mtt-apr.out')
C      CALL RSEP3T('API',API,'mtt-api.out')
      CALL RSEP3T('GPAMPR',GPAMPR,'gpampr.out')
      CALL RSEP3T('GPAMPI',GPAMPI,'gpampi.out')
C      CALL RSEP3T('FGBR22',FGBR22,'fgbr22.out')
C      CALL RSEP3T('FGBY22',FGBY22,'fgby22.out')
      CALL RSEP3T('GPGRID',GPGRID,'gpgrid.out')
      CALL RSEP3T('GPK0',GPK0,'gpm-par.out')

C     Reading filenames of the output files
      CALL RSEP3T('GPMSEC',GPMSEC,'gpmsec.out')
C      CALL RSEP3T('GPAPTZ',GPAPTZ,'gpaptz.out')

      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3R('D1',D1,4.)
      CALL RSEP3R('D2',D2,4.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)

      CALL RSEP3R('DT',DT,0.1)

      CALL RSEP3I('IXTZU',IXTZU,1)
      CALL RSEP3I('IXTZD',IXTZD,1)
      CALL RSEP3I('IXTZL',IXTZL,1)
      CALL RSEP3I('IXTZR',IXTZR,1)

      CALL RSEP3R('KAPPA',KAPPA,1.253314137)

      OPEN(LU2,FILE=RAYALL,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=INIALL,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU4,FILE=NUM,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU5,FILE=MTT,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU6,FILE=APR,FORM='FORMATTED',STATUS='OLD')
C      OPEN(LU7,FILE=API,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU8,FILE=GPAMPR,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU9,FILE=GPAMPI,FORM='FORMATTED',STATUS='OLD')
C      OPEN(LU10,FILE=FGBR22,FORM='FORMATTED',STATUS='OLD')
C      OPEN(LU11,FILE=FGBY22,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU12,FILE=GPGRID,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU13,FILE=GPK0,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU14,FILE=GPMSEC,FORM='FORMATTED')
C      OPEN(LU15,FILE=GPAPTZ,FORM='FORMATTED')

      N1N2=N1*N2
      ISH=2*N1N2+1

      READ(LU12,*) OXR,OPR,OTR,OOMEGA,DXR,DPR,DTR,DOMEGA,
     *             NXR,NPR,NTR,NOMEGA
      NXPT =NXR*NPR*NTR
      NXPTO=NXPT*NOMEGA

      READ(LU13,*) AK0

      CALL RARRAI(LU4,' ','FORMATTED',.FALSE.,0,N1N2,IRAM)
      IRAM(N1N2+1)=0
      DO 10 IGRD=1,N1N2
        IRAM(N1N2+IGRD+1)=IRAM(N1N2+IGRD)+IRAM(IGRD)
   10 CONTINUE

      IF(3*N1N2+2*IRAM(ISH)+1.GT.MRAM) THEN
C       GPMIG-02
        CALL ERROR('GPMIG-02: Too small array RAM(MRAM)')
C       Too small array RAM(MRAM).  If possible, increase dimension MRAM
C       in include file ram.inc.
      END IF

      CALL RARRAY(LU5,' ','FORMATTED',.FALSE.,0.,IRAM(ISH),
     *            RAM(ISH+1))
      CALL RARRAY(LU6,' ','FORMATTED',.FALSE.,1.,IRAM(ISH),
     *            RAM(ISH+IRAM(ISH)+1))
C      CALL RARRAY(LU7,' ','FORMATTED',.FALSE.,1.,IRAM(ISH),
C     *            RAM(ISH+2*IRAM(ISH)+1))
      CALL RARRAY(LU8,' ','FORMATTED',.FALSE.,0.,NXPTO,
     *            RAM(ISH+3*IRAM(ISH)+1))
      CALL RARRAY(LU9,' ','FORMATTED',.FALSE.,0.,NXPTO,
     *            RAM(ISH+3*IRAM(ISH)+NXPTO+1))
C      CALL RARRAY(LU10,' ','FORMATTED',.FALSE.,0.,NXPT,
C     *            RAM(ISH+3*IRAM(ISH)+2*NXPTO+1))
C      CALL RARRAY(LU11,' ','FORMATTED',.FALSE.,0.,NXPT,
C     *            RAM(ISH+3*IRAM(ISH)+2*NXPTO+NXPT+1))

      INUM   =N1N2
      IMTT   =ISH
      IAPR   =ISH+IRAM(ISH)
      IAPI   =ISH+2*IRAM(ISH)
      IGPAR  =ISH+3*IRAM(ISH)
      IGPAI  =ISH+3*IRAM(ISH)+NXPTO
      IFGBR  =ISH+3*IRAM(ISH)+2*NXPTO
      IFGBI  =ISH+3*IRAM(ISH)+2*NXPTO+NXPT
      IGPMSEC=ISH+3*IRAM(ISH)+2*NXPTO+2*NXPT
C      IGPAPTZ=IGPMSEC+N1N2
      IGPAPTZ=IGPMSEC+N1N2+2301*751

      IF((IXTZU.LT.1).OR.(IXTZU.GT.N1)) THEN
        IXTZU=1
      END IF
      IF((IXTZD.LT.1).OR.(IXTZD.GT.N1)) THEN
        IXTZD=N1
      END IF
      IF((IXTZL.LT.1).OR.(IXTZL.GT.N2)) THEN
        IXTZL=1
      END IF
      IF((IXTZR.LT.1).OR.(IXTZR.GT.N2)) THEN
        IXTZR=N2
      END IF

      IXROLD=0
      IPROLD=0

      DO 12 I=1,N1N2
        RAM(IGPMSEC+I)=0.0
C        RAM(IGPMSEC+I)=1.0E20
   12 CONTINUE

      DO 15 I=IXTZU,IXTZD
        DO 14 J=IXTZL,IXTZR
          RAM(IGPMSEC+N2*(N1-I)+J)=0.0
   14   CONTINUE
   15 CONTINUE

      DO 18 I=1,NXTPO
        RAM(IGPAPTZ+I)=0.0
   18 CONTINUE

C     Loop for the points of rays
   20 CONTINUE
C       Reading the results of the complete ray tracing
        CALL AP00(0,LU2,LU3)
        CALL AP03(0,HI,H,HUI)

        IF(IWAVE.LT.1)THEN
C         End of rays
          GO TO 100
        END IF

        IF (IPT.LE.1) THEN
C         New ray
          IXR=NINT((YI(4)-OXR)/DXR)+1
          IPR=NPR-NINT((YI(7)-OPR)/DPR)

C          IF ((YLI(1)*YI(4)).LE.-1.000001) THEN
C            IPR=IPROLD-1
C          ELSE IF ((YLI(1)*YI(4)).GE.0.999999) THEN
C            IPR=IPROLD-1
C          END IF
C          IPROLD=IPR
C          IF (IPR.EQ.1) THEN
C            IPROLD=NPR+1
C          END IF

          IF ((IPR.EQ.IPROLD).AND.(IXR.EQ.IXROLD)) THEN
            GO TO 20
          END IF
          IXROLD=IXR
          IPROLD=IPR

C          WRITE(*,'(A,1(G16.6,1X))') '(YI(4)-OXR)/DXR=',(YI(4)-OXR)/DXR
C          WRITE(*,'(A,1(G16.6,1X))') '(YI(7)-OPR)/DPR=',(YI(7)-OPR)/DPR
C          WRITE(*,'(A,1(G16.6,1X))') 'YI(4) =',YI(4)
C          WRITE(*,'(A,1(G16.6,1X))') 'YI(7) =',YI(7)
C          WRITE(*,'(A,I6)')          'IXR   =',IXR
C          WRITE(*,'(A,I6)')          'IPR   =',IPR
          WRITE(*,'(A,I6,A,I6)')      'IXR   =',IXR,'      IPR   =',IPR

          HH(1,1)= HI(1)
          HH(2,1)= HI(2)
          HH(3,1)= HI(3)
          HH(1,2)= HI(4)
          HH(2,2)= HI(5)
          HH(3,2)= HI(6)
          HH(1,3)=-HI(7)
          HH(2,3)=-HI(8)
          HH(3,3)=-HI(9)
          VI     =YLI(1)
          VI1    =YLI(4)*HH(1,1)+YLI(5)*HH(2,1)+YLI(6)*HH(3,1)
          VI2    =YLI(4)*HH(1,2)+YLI(5)*HH(2,2)+YLI(6)*HH(3,2)
          VI3    =YLI(4)*HH(1,3)+YLI(5)*HH(2,3)+YLI(6)*HH(3,3)
        END IF

        IX1=NINT((Y(5)-O1)/D1)+1
        IX2=NINT((Y(4)-O2)/D2)+1

C        IF(IX1.GT.IXTZD) THEN
C          GO TO 20
C        ELSE IF(IX1.LT.IXTZU) THEN
C          GO TO 20
C        ELSE IF(IX2.GT.IXTZR) THEN
C          GO TO 20
C        ELSE IF(IX2.LT.IXTZL) THEN
C          GO TO 20
C        END IF

        IF(IX1.GT.(IXTZD+20)) THEN
          GO TO 20
        ELSE IF(IX1.LT.(IXTZU-20)) THEN
          GO TO 20
        ELSE IF(IX2.GT.(IXTZR+20)) THEN
          GO TO 20
        ELSE IF(IX2.LT.(IXTZL-20)) THEN
          GO TO 20
        END IF

C        IX1U=IX1-9
C        IX1D=IX1+10
C        IX2L=IX2-9
C        IX2R=IX2+10

        IX1U=IX1-14
        IX1D=IX1+15
        IX2L=IX2-14
        IX2R=IX2+15

C        IX1U=IX1-49
C        IX1D=IX1+50
C        IX2L=IX2-49
C        IX2R=IX2+50

        IF(IX1D.GT.IXTZD) THEN
          IX1D=IXTZD
        END IF
        IF(IX1U.LT.IXTZU) THEN
          IX1U=IXTZU
        END IF
        IF(IX2R.GT.IXTZR) THEN
          IX2R=IXTZR
        END IF
        IF(IX2L.LT.IXTZL) THEN
          IX2L=IXTZL
        END IF

C        IX1U=IXTZU
C        IX1D=IXTZD
C        IX2L=IXTZL
C        IX2R=IXTZR

C        WRITE(*,'(A,I6)')          'IX1    =',IX1
C        WRITE(*,'(A,I6)')          'IX2    =',IX2
C        WRITE(*,'(A,I6)')          'IX1U   =',IX1U
C        WRITE(*,'(A,I6)')          'IX1D   =',IX1D
C        WRITE(*,'(A,I6)')          'IX2L   =',IX2L
C        WRITE(*,'(A,I6)')          'IX2R   =',IX2R

        W1=-Y(8)
        W2=-Y(7)
        W4=-1.0

C       Loop over time in the common-shot gather
        DO 90 ITR=1,NTR
C          WRITE(*,'(A,I6)')          'ITR   =',ITR
C          WRITE(*,'(A,I6)')          'NTR   =',NTR
          TR=OTR+DTR*(ITR-1)

          R4A=Y(1)-TR

C         Loop over frequency
          DO 80 IOMEGA=1,NOMEGA

            OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)

C            IF((1.5*D1*OMEGA).GE.(PI*YL(1))) THEN
C            IF(ABS(D1*OMEGA*YI(7)).GE.PI) THEN
C              GO TO 79
C            END IF

C            IF(((IXR.NE.20).OR.(IPR.NE.12)).OR.
C     *         ((ITR.NE.8).OR.(IOMEGA.NE.40))) THEN
C              GO TO 79
C            END IF
C            AK0=30.0

C            IF(((IXR.NE.10).OR.(IPR.NE.14)).OR.
C     *          (IOMEGA.NE.10)) THEN
C              GO TO 79
C            END IF

C            IF(((IXR.NE.10).OR.(IPR.NE.8)).OR.
C     *          (ITR.NE.20)) THEN
C              GO TO 79
C            END IF

C            WRITE(*,'(A,1(G16.6,1X))') 'Y(1)    =',Y(1)
C            IF(NINT(Y(1)/DT).NE.15) THEN
C              GO TO 79
C            END IF

C            IF((IXR.NE.10).OR.(IPR.NE.8)) THEN
C              GO TO 79
C            END IF

C            IF(IXR.NE.25) THEN
C              GO TO 79
C            END IF

C           Initial shape of GPs
            DO 25 I=1,4
              DO 24 J=1,4
                N(I,J) =CMPLX(0.0,0.0)
                M0(I,J)=CMPLX(0.0,0.0)
   24         CONTINUE
   25       CONTINUE

            IGRD=((IXR-1)*NPR+IPR-1)*NTR+ITR

C            N44=CMPLX(0.0,0.25E-6)*CMPLX(0.0,AK0)
C     *         /CMPLX(0.0,0.25E-6*OMEGA-AK0*YI(7)**2)
            N44=CMPLX(-0.25E-6,0.25E-6)*CMPLX(0.0,AK0)
     *         /CMPLX(-0.25E-6*OMEGA,0.25E-6*OMEGA-AK0*YI(7)**2)

C            N(2,2)=CMPLX(RAM(IFGBR+IGRD),RAM(IFGBI+IGRD))+N44*YI(7)**2
C            N(2,2)=CMPLX(-0.25E-6,0.25E-6)+N44*YI(7)**2
C            N(2,4)=-N44*YI(7)
C            N(4,2)=N(2,4)
C            N(4,4)=N44

            M0(4,4)= N44
            M0(3,3)=(M0(4,4)-VI3)/VI**2
            M0(3,2)=-VI2/VI**2
            M0(2,3)= M0(3,2)
            M0(4,3)=-M0(4,4)/VI
            M0(3,4)= M0(4,3)
C            M0(2,2)=(CMPLX(RAM(IFGBR+IGRD),RAM(IFGBI+IGRD))
C     *             +2.0*VI2*HH(2,2)*HH(2,3)/VI**2
C     *                 +VI3*HH(3,2)*HH(3,2)/VI**2)/HH(2,2)**2
C            M0(2,2)=(CMPLX(-0.25E-6,0.25E-6)
C     *             +2.0*VI2*HH(2,2)*HH(2,3)/VI**2
C     *                 +VI3*HH(3,2)*HH(3,2)/VI**2)/HH(2,2)**2

C            M0(2,2)=(CMPLX(0.0,0.25E-6)
C     *             -2.0*VI2*HI(5)*HI(8)/VI**2
C     *                 +VI3*HI(8)*HI(8)/VI**2)/HI(5)**2
            M0(2,2)=(CMPLX(-0.25E-6,0.25E-6)
     *             -2.0*VI2*HI(5)*HI(8)/VI**2
     *                 +VI3*HI(8)*HI(8)/VI**2)/HI(5)**2

C            M0(1,1)=CMPLX(-REAL(M0(1,1)), AIMAG(M0(1,1)))
C            M0(1,2)=CMPLX(-REAL(M0(1,2)), AIMAG(M0(1,2)))
C            M0(1,3)=CMPLX(-REAL(M0(1,3)), AIMAG(M0(1,3)))
C            M0(1,4)=CMPLX( REAL(M0(1,4)),-AIMAG(M0(1,4)))
C            M0(2,1)=CMPLX(-REAL(M0(2,1)), AIMAG(M0(2,1)))
C            M0(2,2)=CMPLX(-REAL(M0(2,2)), AIMAG(M0(2,2)))
C            M0(2,3)=CMPLX(-REAL(M0(2,3)), AIMAG(M0(2,3)))
C            M0(2,4)=CMPLX( REAL(M0(2,4)),-AIMAG(M0(2,4)))
C            M0(3,1)=CMPLX(-REAL(M0(3,1)), AIMAG(M0(3,1)))
C            M0(3,2)=CMPLX(-REAL(M0(3,2)), AIMAG(M0(3,2)))
C            M0(3,3)=CMPLX(-REAL(M0(3,3)), AIMAG(M0(3,3)))
C            M0(3,4)=CMPLX( REAL(M0(3,4)),-AIMAG(M0(3,4)))
C            M0(4,1)=CMPLX( REAL(M0(4,1)),-AIMAG(M0(4,1)))
C            M0(4,2)=CMPLX( REAL(M0(4,2)),-AIMAG(M0(4,2)))
C            M0(4,3)=CMPLX( REAL(M0(4,3)),-AIMAG(M0(4,3)))
C            M0(4,4)=CMPLX(-REAL(M0(4,4)), AIMAG(M0(4,4)))

            HH(1,1)= H(1)
            HH(2,1)= H(2)
            HH(3,1)= H(3)
            HH(1,2)= H(4)
            HH(2,2)= H(5)
            HH(3,2)= H(6)
            HH(1,3)=-H(7)
            HH(2,3)=-H(8)
            HH(3,3)=-H(9)

            V =YL(1)
            V1=YL(4)*HH(1,1)+YL(5)*HH(2,1)+YL(6)*HH(3,1)
            V2=YL(4)*HH(1,2)+YL(5)*HH(2,2)+YL(6)*HH(3,2)
            V3=YL(4)*HH(1,3)+YL(5)*HH(2,3)+YL(6)*HH(3,3)

C            C=M0(4,4)+CMPLX(0.0,0.5)*(M0(1,4)**2+M0(2,4)**2)
C     *                                       /AIMAG(M0(2,2))
            C=M0(4,4)

            Q1= Y(12)-Y(20)*M0(1,1)-Y(24)*M0(2,1)
            Q2= Y(13)-Y(21)*M0(1,1)-Y(25)*M0(2,1)
            Q3= Y(16)-Y(20)*M0(1,2)-Y(24)*M0(2,2)
            Q4= Y(17)-Y(21)*M0(1,2)-Y(25)*M0(2,2)
            P1=-Y(14)+Y(22)*M0(1,1)+Y(26)*M0(2,1)
            P2=-Y(15)+Y(23)*M0(1,1)+Y(27)*M0(2,1)
            P3=-Y(18)+Y(22)*M0(1,2)+Y(26)*M0(2,2)
            P4=-Y(19)+Y(23)*M0(1,2)+Y(27)*M0(2,2)

            DET   = Q1*Q4-Q2*Q3
            M(1,1)=          (P1*Q4-P3*Q2)/DET
            M(2,1)=          (P2*Q4-P4*Q2)/DET
            M(1,2)= M(2,1)
            M(2,2)=          (P4*Q1-P2*Q3)/DET
            M(1,4)=(Q4*M0(1,4)-Q2*M0(2,4))/DET
            M(2,4)=(Q1*M0(2,4)-Q3*M0(1,4))/DET
            M(4,1)= M(1,4)
            M(4,2)= M(2,4)
            M(1,3)=-M(1,4)/V-V1/V**2
            M(2,3)=-M(2,4)/V-V2/V**2
            M(3,1)= M(1,3)
            M(3,2)= M(2,3)
            M(4,4)= C-CMPLX(0.0,0.5)*(M0(1,4)**2*(Q3**2+Q4**2)
     *                      -2.0*M0(1,4)*M0(2,4)*(Q1*Q3+Q2*Q4)
     *                               +M0(2,4)**2*(Q1**2+Q2**2))
     *                                  /(AIMAG(M(2,2))*DET**2)
            M(3,4)=-M(4,4)/V
            M(4,3)= M(3,4)
            M(3,3)=(M(4,4)-V3)/V**2

C            M(1,1)=CMPLX(-REAL(M(1,1)), AIMAG(M(1,1)))
C            M(1,2)=CMPLX(-REAL(M(1,2)), AIMAG(M(1,2)))
C            M(1,4)=CMPLX( REAL(M(1,4)),-AIMAG(M(1,4)))
C            M(2,1)=CMPLX(-REAL(M(2,1)), AIMAG(M(2,1)))
C            M(2,2)=CMPLX(-REAL(M(2,2)), AIMAG(M(2,2)))
C            M(2,4)=CMPLX( REAL(M(2,4)),-AIMAG(M(2,4)))
C            M(4,1)=CMPLX( REAL(M(4,1)),-AIMAG(M(4,1)))
C            M(4,2)=CMPLX( REAL(M(4,2)),-AIMAG(M(4,2)))
C            M(4,4)=CMPLX(-REAL(M(4,4)), AIMAG(M(4,4)))

C            M(1,3)=-M(1,4)/V-V1/V**2
C            M(2,3)=-M(2,4)/V-V2/V**2
C            M(3,1)= M(1,3)
C            M(3,2)= M(2,3)
C            M(3,3)=(M(4,4)-V3)/V**2
C            M(3,4)=-M(4,4)/V
C            M(4,3)= M(3,4)

C            M(1,3)=CMPLX(-REAL(M(1,3)), AIMAG(M(1,3)))
C            M(2,3)=CMPLX(-REAL(M(2,3)), AIMAG(M(2,3)))
C            M(3,1)=CMPLX(-REAL(M(3,1)), AIMAG(M(3,1)))
C            M(3,2)=CMPLX(-REAL(M(3,2)), AIMAG(M(3,2)))
C            M(3,3)=CMPLX(-REAL(M(3,3)), AIMAG(M(3,3)))
C            M(3,4)=CMPLX( REAL(M(3,4)),-AIMAG(M(3,4)))
C            M(4,3)=CMPLX( REAL(M(4,3)),-AIMAG(M(4,3)))

            N(4,4)=M(4,4)
            DO 43 I=1,3
              N(I,4)=CMPLX(0.0,0.0)
              DO 42 J=1,3
                N(I,4)=N(I,4)+M(J,4)*HH(I,J)
                N(4,I)=N(I,4)
                N(I,J)=CMPLX(0.0,0.0)
                DO 41 K=1,3
                  DO 40 L=1,3
                    N(I,J)=N(I,J)+M(K,L)*HH(I,K)*HH(J,L)
   40             CONTINUE
   41           CONTINUE
   42         CONTINUE
   43       CONTINUE

C           Loop over ray history
            IX1=NINT((Y(5)-O1)/D1)+1
            IX2=NINT((Y(4)-O2)/D2)+1
            NHISTORY=IRAM(N1*(IX2-1)+IX1)
            DO 45 IHISTORY=1,NHISTORY
              R4=R4A+RAM(IMTT+IRAM(INUM+N1*(IX2-1)+IX1)+IHISTORY)
              IF(ABS(R4).LT.2.0*SQRT(1.0/(OMEGA*AIMAG(N(4,4))))) THEN
                GO TO 46
              END IF
   45       CONTINUE
C            WRITE(*,'(A,1(G16.6,1X))') 'huh   =',R4
            GO TO 79

   46       CONTINUE

C           Initial amplitudes of GPs
            IGRD=(((IXR-1)*NTR+ITR-1)*NPR+IPR-1)*NOMEGA+IOMEGA
            GPAR=RAM(IGPAR+IGRD)
            GPAI=RAM(IGPAI+IGRD)
C            GPAR=500.0
C            WRITE(*,'(A,1(G16.6,1X))') 'GPAR   =',GPAR
C            WRITE(*,'(A,1(G16.6,1X))') 'GPAI   =',GPAI
            IF(((GPAR.GT.-100.0).AND.(GPAR.LT.100.0)).AND.
     *         ((GPAI.GT.-100.0).AND.(GPAI.LT.100.0))) THEN
              GO TO 79
            END IF

C           Amplitudes of GPs
            GPA=SQRT(CMPLX(0.0,1.0)*(Y(20)*Y(25)-Y(21)*Y(24))
     *                    /(ABS(Y(20)*Y(25)-Y(21)*Y(24))*DET))
            GPA=CMPLX(ABS(REAL(GPA)),AIMAG(GPA))
            GPA=GPA*CMPLX(GPAR,GPAI)*SQRT(VI/V)
     *         *CEXP(CMPLX(0.0,1.0)*PI/2.0*(KMAH-0.5))
            GPAR= REAL(GPA)
            GPAI=AIMAG(GPA)

C            WRITE(*,'(A,I6,A,I6,A,I6,A,I6)') 'IXR=',IXR,' ITR=',ITR,
C     *                                ' IPR=',IPR,' IOMEGA=',IOMEGA
C            WRITE(*,'(A,1(G16.6,1X),A,1(G16.6,1X))') 'GPAR=',GPAR,
C     *                                              ' GPAI=',GPAI

            IGRDG=IGPAPTZ+IGRD

C            IX1=IXTZU
            IX1=IX1U

C           Loops over the target zone
   50       CONTINUE

              IX1=IX1+1

C              IF(IX1.GT.IXTZD) THEN
C                GO TO 79
C              END IF

              IF(IX1.GT.IX1D) THEN
                GO TO 79
              END IF

              IGRDA=N2*(IX1-1)

              X1=O1+D1*(IX1-1)
              R1=X1-Y(5)

C              IX2=IXTZL
              IX2=IX2L

   60         CONTINUE

                IX2=IX2+1

C                IF(IX2.GT.IXTZR) THEN
C                  GO TO 50
C                END IF

                IF(IX2.GT.IX2R) THEN
                  GO TO 50
                END IF

C                IGRDB=IGRDA+IX2
                IGRDB=N1*(IX2-1)+IX1
C                IGRDC=N2*(IX1-1)+IX2
                IGRDC=N2*(N1-IX1)+IX2

                X2=O2+D2*(IX2-1)
                R2=X2-Y(4)

C                GPEXP=CMPLX(0.0,OMEGA)*(R1*W1+R2*W2
C     *               +0.5*R1*R1*N(3,3)+R1*R2*N(3,2)+0.5*R2*R2*N(2,2))
C     *                                   -(KAPPA*(R1*W1+R2*W2)/DT)**2
                GPEXP=CMPLX(0.0,OMEGA)*(R1*W1+R2*W2
     *               +0.5*R1*R1*N(3,3)+R1*R2*N(3,2)+0.5*R2*R2*N(2,2))
     *                           -(KAPPA*(R1*Y(8)+R2*Y(7))/DT)**2


                GP1=CMPLX(0.0,0.0)
                GP2=CMPLX(0.0,0.0)
                NHISTORY=IRAM(IGRDB)
C                WRITE(*,'(A,I6)')          'NHISTORY    =',NHISTORY

C               Loop over ray history
                DO 70 IHISTORY=1,NHISTORY

                  J=IRAM(INUM+IGRDB)+IHISTORY

                  R4=R4A+RAM(IMTT+J)
C                  R4=R4A+1.0
C                  R4=-R4
C                  WRITE(*,'(A,1(G16.6,1X))') 'R4   =',R4

                  GPIMAGE=CEXP(GPEXP+CMPLX(0.0,OMEGA)*(R4*W4
     *                              +R1*R4*N(3,4)+R2*R4*N(2,4)
     *                                       +0.5*R4*R4*N(4,4)))

                  RFAR=RAM(IAPR+J)*1.0E10
C                  RFAI=RAM(IAPI+J)*1.0E10
C                  RFAR=1.0
                  RFAI=0.0

                  GP1=GP1+CMPLX(RFAR,-RFAI)*CMPLX(GPAR,GPAI)*GPIMAGE
                  GP2=GP2+CMPLX(RFAR,-RFAI)*CMPLX(RFAR,RFAI)

   69             CONTINUE
   70           CONTINUE

                IF((REAL(GP2).NE.0.0).OR.(AIMAG(GP2).NE.0.0)) THEN
                  I         =IGPMSEC+IGRDC
                  RAM(I)    =RAM(I)+REAL(GP1/GP2)
C                  RAM(I)    =RAM(I)+AIMAG(GP1/GP2)
C                  RAM(IGRDG)=RAM(IGRDG)+RAM(I)
                END IF

              GO TO 60

            GO TO 50

   79       CONTINUE
   80     CONTINUE

   90   CONTINUE

      GO TO 20

  100 CONTINUE

C     Interpolating output grid values
      WRITE(*,'(A)')'+GPMIG: Interpolating... '

      DO 105 I=1,1728051
        RAM(IGPMSEC+N1N2+I)=0.0
  105 CONTINUE

      D= D1/4.0
      ID=NINT(D)
      IF(ID.EQ.1) THEN
        DO 107 I=1,1728051
          RAM(IGPMSEC+N1N2+I)=RAM(IGPMSEC+I)
  107   CONTINUE
        GO TO 160
      END IF

      DO 115 IX2=1,N2
        IX2NEW=(IX2-1)*ID+1
        DO 110 IX1=1,N1
          IX1NEW=(IX1-1)*ID+1
          IGRD   =  N2*(IX1   -1)+IX2
          IGRDNEW=2301*(IX1NEW-1)+IX2NEW
          RAM(IGPMSEC+N1N2+IGRDNEW)=RAM(IGPMSEC+IGRD)
  110   CONTINUE
  115 CONTINUE

C      GO TO 160

      DO 130 IX2=1,N2
        IX2NEW=(IX2-1)*ID+1
        DO 125 IX1=1,N1
          IX1NEW=(IX1-1)*ID+1
C          IGRDNEW=751*(IX2NEW-1)+IX1NEW
          IGRDNEW=2301*(IX1NEW-1)+IX2NEW
          DO 120 I=1,751
            IF(I.NE.IX1NEW) THEN
C              IGRD=751*(IX2NEW-1)+I
              IGRD=2301*(I-1)+IX2NEW
              RAM(IGPMSEC+N1N2+IGRD)=RAM(IGPMSEC+N1N2+IGRD)
     *                           +RAM(IGPMSEC+N1N2+IGRDNEW)
     *              *SIN(PI*(I-IX1NEW)/D)/(PI*(I-IX1NEW)/D)
            END IF
  120     CONTINUE
  125   CONTINUE
  130 CONTINUE

      DO 150 I=1,751
        DO 145 IX2=1,N2
          IX2NEW=(IX2-1)*ID+1
C          IGRDNEW=751*(IX2NEW-1)+I
          IGRDNEW=2301*(I-1)+IX2NEW
          DO 140 J=1,2301
            IF(J.NE.IX2NEW) THEN
C              IGRD=751*(J-1)+I
              IGRD=2301*(I-1)+J
              RAM(IGPMSEC+N1N2+IGRD)=RAM(IGPMSEC+N1N2+IGRD)
     *                           +RAM(IGPMSEC+N1N2+IGRDNEW)
     *              *SIN(PI*(J-IX2NEW)/D)/(PI*(J-IX2NEW)/D)
            END IF
  140     CONTINUE
  145   CONTINUE
  150 CONTINUE

  160 CONTINUE

C     Writing output grid values
C      IF(GPMSEC.NE.' ') THEN
C        CALL WARRAY(LU14,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
C     *                                     N1N2,RAM(IGPMSEC+1))
C      END IF

      IF(GPMSEC.NE.' ') THEN
        CALL WARRAY(LU14,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                             1728051,RAM(IGPMSEC+N1N2+1))
      END IF

C     IF(GPAPTZ.NE.' ') THEN
C       CALL WARRAY(LU15,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
C    *                                    NXTPO,RAM(IGPAPTZ+1))
C     END IF

      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
C      CLOSE(LU7)
      CLOSE(LU8)
      CLOSE(LU9)
C      CLOSE(LU10)
C      CLOSE(LU11)
      CLOSE(LU12)
      CLOSE(LU13)
      CLOSE(LU14)
C      CLOSE(LU15)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'forms.for'
C     forms.forC
C
C=======================================================================
C