C 
C Program GPMIG for 2-D Gaussian packet prestack depth migration
C
C Version: 8.00
C Date: 2022, April 24
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
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       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       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     APR='string'... Name of the input file with the array of real
C             values, containing for each gridpoint the amplitudes.
C             Default: APR='mtt-apr.out'
C Names of other input files:
C     GPAMPR='string'... Name of the output file with
C             amplitudes of Gaussian packets - real part.
C             Default: GPAMPR='gpampr.out'
C     GPAMPI='string'... Name of the output file with
C             amplitudes of Gaussian packets - imaginary part.
C             Default: GPAMPR='gpampi.out'
C     FGBR22='string'... Name of the input file with the resampled
C             smoothed optimized parameters of Gaussian beams Re(N0).
C             Default: FGBR22=' '
C     FGBY22='string'... Name of the input file with the resampled
C             smoothed optimized parameters of Gaussian beams Im(N0).
C             Default: FGBY22=' '
C     GPSTEP='string'... Name of the input file containing
C             discretization steps NXR,DXR,OXR,NTR,DTR,OTR,NP,DP,OP,
C             NOMEGA,DOMEGA,OOMEGA, complex-valued constants
C             N0, K0, ~K0, Im(N)min (same as Y0min).
C             Default: GPSTEP=' '
C Data specifying the dimensions of the input grid:
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     N1MAX=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1MAX=1
C     N2MAX=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2MAX=1
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     O1=real, O2=real ... Coordinates of the origin of the grid.
C             Default: O1=0., O2=0.
C Input SEP parameters:
C     NX=positive integer... Number of gridpoints along the profile
C             axis X of the decomposed wavefield.
C             Default: NX=1
C     NT=positive integer... Number of gridpoints along the time axis T
C             of the decomposed wavefield.
C             Default: NT=1
C     DX=real... Grid interval along the profile axis X of
C             the decomposed wavefield.
C             Default: DX=1.
C     DTRAY=real... Time step along the central ray.
C             Default: DTRAY=1.
C     OT=real, OX=real ... Coordinates of the origin of the grid.
C             Default: OT=0., OX=0.
C     IXTZU=positive integer... Limits of target zone (up).
C             Default: IXTZU=1
C     IXTZD=positive integer... Limits of target zone (down)
C             Default: IXTZD=1
C     IXTZL=positive integer... Limits of target zone (left).
C             Default: IXTZL=1
C     IXTZR=positive integer... Limits of target zone (right).
C             Default: IXTZR=1
C     KAPPA=real... Parameter kappa has the meaning of the ratio of
C             the step between two subsequent beams to their
C             characteristic half-width.
C             Default: KAPPA=1.253314137
C     RELAMP=real... Relative decay of the Gaussian envelope at which
C             the loop over the points of the seismograms is terminated.
C             The relative error due to this economizing roughly
C             corresponds to the value of RELAMP.
C             Default: RELAMP=0.1
C Name of output file:
C     GPMSEC='string'... Name of the output file containing
C             migrated section.
C             Default: GPMSEC='gpmsec.out'
C Form of the input/output files:
C     FORM='string'... Form of the input/output files: either
C             'FORMATTED' or 'UNFORMATTED' (case insensitive).
C             Default: FORM='FORMATTED'
C     FORMW='string'... Form of the output file. If both FORM and FORMW
C             are specified, FORMW is used for output files.
C             Default: FORMW=FORM
C     FORMR='string'... Form of the input file. If both FORM and FORMR
C             are specified, FORMR is used for output files.
C             Default: FORMR=FORM
C
C Input unformatted file RAYALL:
C     See the description within source code file 'writ.for'.
C     Description of files
C     CRT-R
C
C Input unformatted file INIALL:
C     See the description within source code file 'writ.for'.
C     Description of files
C     CRT-I
C
C Optional parameters specifying the form of the quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             
C             forms.for.
C     NUMLIN=positive integer ... Number of reals or integers in one
C             line of the output file. See the description in file
C             
C             forms.for.
C Value of undefined quantities:
C     UNDEF=real... The value to be used for undefined real quantities.
C             Default: UNDEF=undefined value used in forms.for
C
C=======================================================================
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.
C
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL AP00,OARRAY,RARRAY,RARRAI,ERROR,RSEP1,RSEP3T,RSEP3I,
C    *         RSEP3R
C
C     Filenames and parameters:
      CHARACTER*80 FILSEP,RAYALL,INIALL,NUM,MTT,APR
      CHARACTER*80 GPAMPR,GPAMPI,FGBY22,FGBR22,GPSTEP,GPMSEC
      INTEGER LU1,LU2
      REAL PI
      PARAMETER (LU1=1,LU2=2,PI=3.141592654)
C
C
C     Auxiliary storage locations:
      REAL O1,O2,D1,D2,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 RFAR,OMEGA
      REAL TR,CGPAR,CGPAI,DXR,DTR,DOMEGA,DP,OXR,OTR,OOMEGA,OP,D
      REAL RELAMP,Y22,Y23,Y33,Y24,Y34,Y44,YDET,X22,X33
      REAL AK0,AK0T,AN0R,AN0I,ANIMIN
      REAL OX,OT,DX,DTRAY
      COMPLEX CM0(4,4),CN(4,4),CM(4,4),CC,CQ1,CQ2,CQ3,CQ4
      COMPLEX CP1,CP2,CP3,CP4,CDET,CN44,CGPA
      REAL GPIMGR,GPIMGI,CGP1R,CGP2R
      REAL AUX,CAUX2R,CAUX2I
      REAL GPEXPR,GPEXPI,CN33R,CN33I,CN32R,CN32I,CN22R,CN22I
      REAL CN24R,CN24I,CN34R,CN34I,CN44R,CN44I
      REAL OW1,OW2,OW4,OCN22R,OCN22I,OCN32R,OCN32I,OCN33R,OCN33I
      REAL OCN24R,OCN24I,OCN34R,OCN34I,OCN44R,OCN44I,CAUXR1,CAUXI1
      REAL YINVY3,AUX1,AUX2,EXPR0,EXPR1,EXPR2,EXPI0,EXPI1,EXPI2
      REAL EXPR00,EXPR10,PID
      INTEGER I,J,K,L,IGRD
      INTEGER IRAM(MRAM)
      INTEGER ITR,NTR,IOMEGA,NOMEGA,N1,N2,N1N2,ISH,N2MAX,N1MAX,N1N2MA
      INTEGER IXR,NXR,IP,NP,NXPT,NXPTO,IGRDB,IGRDC
      INTEGER IHIST,NHIST,IXTZL,IXTZR,IXTZU,IXTZD,IPOLD,IXROLD
      INTEGER IX1,IX2,IX1U,IX1D,IX2L,IX2R,ID1,ID2
      INTEGER INUM,IMTT,IAPR,IAPI,IGPSEC
      INTEGER ICGPAR,ICGPAI,IFGBR,IFGBI,IX1NEW,IX2NEW,IGRDNE,ID
      INTEGER NX,NT
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
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
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading filenames of the input files
      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')
      CALL RSEP3T('GPAMPR',GPAMPR,'gpampr.out')
      CALL RSEP3T('GPAMPI',GPAMPI,'gpampi.out')
      CALL RSEP3T('FGBR22',FGBR22,'fgbr22.out')
      CALL RSEP3T('FGBY22',FGBY22,'fgby22.out')
      CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.out')
C
C     Reading input parameters
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N2MAX',N2MAX,1)
      CALL RSEP3I('N1MAX',N1MAX,1)
      CALL RSEP3R('D1',D1,4.)
      CALL RSEP3R('D2',D2,4.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
C
      CALL RSEP3I('NX',NX,116)
      CALL RSEP3R('DX',DX,25.)
      CALL RSEP3R('OX',OX,0.)
      CALL RSEP3I('NT',NT,750)
      CALL RSEP3R('DTRAY',DTRAY,0.1)
      CALL RSEP3R('OT',OT,0.)
C
      CALL RSEP3I('IXTZU',IXTZU,1)
      CALL RSEP3I('IXTZD',IXTZD,1)
      CALL RSEP3I('IXTZL',IXTZL,1)
      CALL RSEP3I('IXTZR',IXTZR,1)
C
      CALL RSEP3R('KAPPA',KAPPA,1.253314137)
      CALL RSEP3R('RELAMP',RELAMP,0.1)
C
C     Reading filenames of the output files
      CALL RSEP3T('GPMSEC',GPMSEC,'gpmsec.out')
      CLOSE(LU1)
C
      RELAMP=-ALOG(RELAMP)
C
C
      N1N2=N1*N2
      N1N2MA=N2MAX*N1MAX
      ISH=2*N1N2+1
C
      OPEN(LU1,FILE=GPSTEP,FORM='FORMATTED',STATUS='OLD')
      READ(LU1,'(A)') TEXT
      READ(LU1,*) NXR,DXR,OXR
      READ(LU1,*) NTR,DTR,OTR
      READ(LU1,*) NP, DP, OP
      READ(LU1,*) NOMEGA,DOMEGA,OOMEGA
C     Reading complex-valued constants
      READ(LU1,*) AN0R,AN0I
      READ(LU1,*) AK0,AK0T,ANIMIN
      CLOSE(LU1)
C
C     Writing discretization steps to the screen
      WRITE(*,'(A,I6)')          ' NXR   =',NXR
      WRITE(*,'(A,1(G16.6,1X))') ' DXR   =',DXR
      WRITE(*,'(A,1(G16.6,1X))') ' OXR   =',OXR
      WRITE(*,'(A,I6)')          ' NTR   =',NTR
      WRITE(*,'(A,1(G16.6,1X))') ' DTR   =',DTR
      WRITE(*,'(A,1(G16.6,1X))') ' OTR   =',OTR
      WRITE(*,'(A,I6)')          ' NP    =',NP
      WRITE(*,'(A,1(G16.6,1X))') ' DP    =',DP
      WRITE(*,'(A,1(G16.6,1X))') ' OP    =',OP
      WRITE(*,'(A,I6)')          ' NOMEGA=',NOMEGA
      WRITE(*,'(A,1(G16.6,1X))') ' DOMEGA=',DOMEGA
      WRITE(*,'(A,1(G16.6,1X))') ' OOMEGA=',OOMEGA
C
C     Writing complex-valued constants to the screen
      WRITE(*,'(A,1(G16.6,1X))') ' N0R   =',AN0R
      WRITE(*,'(A,1(G16.6,1X))') ' N0I   =',AN0I
      WRITE(*,'(A,1(G16.6,1X))') ' K0    =',AK0
      WRITE(*,'(A,1(G16.6,1X))') ' ~K0   =',AK0T
      WRITE(*,'(A,1(G16.6,1X))') ' ANIMIN=',ANIMIN
C
C     Writing constants and grid values to the screen
      WRITE(*,'(A,1(G16.6,1X))') ' KAPPA =',KAPPA
      WRITE(*,'(A,I6)')          ' NX    =',NX
      WRITE(*,'(A,1(G16.6,1X))') ' DX    =',DX
      WRITE(*,'(A,1(G16.6,1X))') ' OX    =',OX
      WRITE(*,'(A,I6)')          ' NT    =',NT
      WRITE(*,'(A,1(G16.6,1X))') ' DTRAY =',DTRAY
      WRITE(*,'(A,1(G16.6,1X))') ' OT    =',OT
C
      NXPT =NXR*NP*NTR
      NXPTO=NXPT*NOMEGA
C
      CALL OARRAY(LU1,NUM,'R')
      CALL RARRAI(LU1,' ',.FALSE.,0,N1N2,IRAM)
      CLOSE(LU1)
      IRAM(N1N2+1)=0
      DO 10 IGRD=1,N1N2
        IRAM(N1N2+IGRD+1)=IRAM(N1N2+IGRD)+IRAM(IGRD)
   10 CONTINUE
C
      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
C       ram.inc.
      END IF
C
      CALL OARRAY(LU1,MTT,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,IRAM(ISH),
     *            RAM(ISH+1))
      CLOSE(LU1)
      CALL OARRAY(LU1,APR,'R')
      CALL RARRAY(LU1,' ',.FALSE.,1.,IRAM(ISH),
     *            RAM(ISH+IRAM(ISH)+1))
      CLOSE(LU1)
C?    CALL RARRAY(LU7,' ',.FALSE.,1.,IRAM(ISH),
C?   *            RAM(ISH+2*IRAM(ISH)+1))
      CALL OARRAY(LU1,GPAMPR,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXPTO,
     *            RAM(ISH+3*IRAM(ISH)+1))
      CLOSE(LU1)
      CALL OARRAY(LU1,GPAMPI,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXPTO,
     *            RAM(ISH+3*IRAM(ISH)+NXPTO+1))
      CLOSE(LU1)
      CALL OARRAY(LU1,FGBR22,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXPT,
     *            RAM(ISH+3*IRAM(ISH)+2*NXPTO+1))
      CLOSE(LU1)
      CALL OARRAY(LU1,FGBY22,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXPT,
     *            RAM(ISH+3*IRAM(ISH)+2*NXPTO+NXPT+1))
      CLOSE(LU1)
C
      INUM  =N1N2
      IMTT  =ISH
      IAPR  =ISH+IRAM(ISH)
C     IAPI  =ISH+2*IRAM(ISH)
      ICGPAR=ISH+3*IRAM(ISH)
      ICGPAI=ISH+3*IRAM(ISH)+NXPTO
      IFGBR =ISH+3*IRAM(ISH)+2*NXPTO
      IFGBI =ISH+3*IRAM(ISH)+2*NXPTO+NXPT
      IGPSEC=ISH+3*IRAM(ISH)+2*NXPTO+2*NXPT
C
      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
C
      IXROLD=0
      IPOLD=0
C
      DO 12 I=1,N1N2
        RAM(IGPSEC+I)=0.
   12 CONTINUE
      DO 15 I=IXTZU,IXTZD
        DO 14 J=IXTZL,IXTZR
          RAM(IGPSEC+N2*(N1-I)+J)=0.
   14   CONTINUE
   15 CONTINUE
C
C     Loop for the points of rays
      OPEN(LU1,FILE=RAYALL,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=INIALL,FORM='UNFORMATTED',STATUS='OLD')
   20 CONTINUE
C       Reading the results of the complete ray tracing
        CALL AP00(0,LU1,LU2)
C       Change of model coordinates X1->X3
C       Array YL
        AUX=YL(6)
        YL(6)=YL(4)
        YL(4)=AUX
C       Array Y
        AUX=Y(5)
        Y(5)=Y(3)
        Y(3)=AUX
        AUX=Y(8)
        Y(8)=Y(6)
        Y(6)=AUX
        AUX=Y(11)
        Y(11)=Y(9)
        Y(9)=AUX
C       Array YLI
        AUX=YLI(6)
        YLI(6)=YLI(4)
        YLI(4)=AUX
C       Array YI
        AUX=YI(5)
        YI(5)=YI(3)
        YI(3)=AUX
        AUX=YI(8)
        YI(8)=YI(6)
        YI(6)=AUX
        AUX=YI(11)
        YI(11)=YI(9)
        YI(9)=AUX
        CALL AP03(0,HI,H,HUI)
C
        IF(IWAVE.LT.1)THEN
C         End of rays
          WRITE(*,*)' End of rays'
          CLOSE(LU1)
          CLOSE(LU2)
          GO TO 100
        END IF
C
        IF (IPT.LE.1) THEN
C         New ray
          IXR=NINT((YI(4)-OXR)/DXR)+1
          IP=NP-NINT((YI(7)-OP)/DP)
          IF ((IP.EQ.IPOLD).AND.(IXR.EQ.IXROLD)) THEN
            GO TO 20
          END IF
          IXROLD=IXR
          IPOLD=IP
          WRITE(*,'(4(A,I3))') '+GPMIG: Calculating NXR=',IXR,'/',NXR,
     *' NP=',IP,'/',NP
          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
        W1=-Y(8)
        W2=-Y(7)
        W4=-1.
C
C       Loop over corresponding arrival time of Gaussian packet
        DO 90 ITR=1,NTR
          TR=OTR+DTR*(ITR-1)
          R4A=Y(1)-TR
C
C         Loop over positive circular frequency of Gaussian packet
          DO 80 IOMEGA=1,NOMEGA
            OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
C
C           Initial shape of Gausian packet
            DO 25 I=1,4
              DO 24 J=1,4
                CN(I,J) =CMPLX(0.,0.)
                CM0(I,J)=CMPLX(0.,0.)
   24         CONTINUE
   25       CONTINUE
C
            IGRD=((IXR-1)*NP+IP-1)*NTR+ITR
            CN44=CMPLX(RAM(IFGBR+IGRD),RAM(IFGBI+IGRD))*CMPLX(0.,AK0)
     *         /CMPLX(RAM(IFGBR+IGRD)*OMEGA,RAM(IFGBI+IGRD)*OMEGA
     *         -AK0*YI(7)**2)
            CM0(4,4)= CN44
            CM0(3,3)=(CM0(4,4)-VI3)/VI**2
            CM0(3,2)=-VI2/VI**2
            CM0(2,3)= CM0(3,2)
            CM0(4,3)=-CM0(4,4)/VI
            CM0(3,4)= CM0(4,3)
            CM0(2,2)=(CMPLX(RAM(IFGBR+IGRD),RAM(IFGBI+IGRD))
     *             -2.*VI2*HI(5)*HI(8)/VI**2
     *                 +VI3*HI(8)*HI(8)/VI**2)/HI(5)**2
            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)
            CC=CM0(4,4)
C-C         CQ1= Y(12)-Y(20)*CM0(1,1)-Y(24)*CM0(2,1)
            CQ1= Y(12)
C-C         CQ2= Y(13)-Y(21)*CM0(1,1)-Y(25)*CM0(2,1)
            CQ2= Y(13)
C-C         CQ3= Y(16)-Y(20)*CM0(1,2)-Y(24)*CM0(2,2)
            CQ3= Y(16)-Y(24)*CM0(2,2)
C-C         CQ4= Y(17)-Y(21)*CM0(1,2)-Y(25)*CM0(2,2)
            CQ4= Y(17)-Y(25)*CM0(2,2)
C-C         CP1=-Y(14)+Y(22)*CM0(1,1)+Y(26)*CM0(2,1)
            CP1=-Y(14)
C-C         CP2=-Y(15)+Y(23)*CM0(1,1)+Y(27)*CM0(2,1)
            CP2=-Y(15)
C-C         CP3=-Y(18)+Y(22)*CM0(1,2)+Y(26)*CM0(2,2)
            CP3=-Y(18)+Y(26)*CM0(2,2)
C-C         CP4=-Y(19)+Y(23)*CM0(1,2)+Y(27)*CM0(2,2)
            CP4=-Y(19)+Y(27)*CM0(2,2)
            CDET   = CQ1*CQ4-CQ2*CQ3
            CM(1,1)=          (CP1*CQ4-CP3*CQ2)/CDET
            CM(2,1)=          (CP2*CQ4-CP4*CQ2)/CDET
            CM(1,2)= CM(2,1)
            CM(2,2)=          (CP4*CQ1-CP2*CQ3)/CDET
            CM(1,4)=(CQ4*CM0(1,4)-CQ2*CM0(2,4))/CDET
            CM(2,4)=(CQ1*CM0(2,4)-CQ3*CM0(1,4))/CDET
            CM(4,1)= CM(1,4)
            CM(4,2)= CM(2,4)
            CM(1,3)=-CM(1,4)/V-V1/V**2
            CM(2,3)=-CM(2,4)/V-V2/V**2
            CM(3,1)= CM(1,3)
            CM(3,2)= CM(2,3)
            CM(4,4)= CC-CMPLX(0.,0.5)*(CM0(1,4)**2*(CQ3**2+CQ4**2)
     *                      -2.*CM0(1,4)*CM0(2,4)*(CQ1*CQ3+CQ2*CQ4)
     *                               +CM0(2,4)**2*(CQ1**2+CQ2**2))
     *                                  /(AIMAG(CM(2,2))*CDET**2)
            CM(3,4)=-CM(4,4)/V
            CM(4,3)= CM(3,4)
            CM(3,3)=(CM(4,4)-V3)/V**2
            CN(4,4)=CM(4,4)
            DO 43 I=1,3
              CN(I,4)=CMPLX(0.,0.)
              DO 42 J=1,3
                CN(I,4)=CN(I,4)+CM(J,4)*HH(I,J)
                CN(4,I)=CN(I,4)
                CN(I,J)=CMPLX(0.,0.)
                DO 41 K=1,3
                  DO 40 L=1,3
                    CN(I,J)=CN(I,J)+CM(K,L)*HH(I,K)*HH(J,L)
   40             CONTINUE
   41           CONTINUE
   42         CONTINUE
   43       CONTINUE
            CN22R=REAL(CN(2,2))
            CN22I=AIMAG(CN(2,2))
            CN32R=REAL(CN(3,2))
            CN32I=AIMAG(CN(3,2))
            CN33R=REAL(CN(3,3))
            CN33I=AIMAG(CN(3,3))
            CN24R=REAL(CN(2,4))
            CN24I=AIMAG(CN(2,4))
            CN34R=REAL(CN(3,4))
            CN34I=AIMAG(CN(3,4))
            CN44R=REAL(CN(4,4))
            CN44I=AIMAG(CN(4,4))
            OCN22R=0.5*OMEGA*CN22R
            OCN22I=0.5*OMEGA*CN22I
            OCN32R=OMEGA*CN32R
            OCN32I=OMEGA*CN32I
            OCN33R=0.5*OMEGA*CN33R
            OCN33I=0.5*OMEGA*CN33I
            OCN24R=OMEGA*CN24R
            OCN24I=OMEGA*CN24I
            OCN34R=OMEGA*CN34R
            OCN34I=OMEGA*CN34I
            OCN44R=0.5*OMEGA*CN44R
            OCN44I=0.5*OMEGA*CN44I
            OW1=OMEGA*W1
            OW2=OMEGA*W2
            OW4=OMEGA*W4
C
C           Loop over ray history
            IX1=NINT((Y(5)-O1)/D1)+1
            IX2=NINT((Y(4)-O2)/D2)+1
            NHIST=IRAM(N1*(IX2-1)+IX1)
            DO 45 IHIST=1,NHIST
              R4=R4A+RAM(IMTT+IRAM(INUM+N1*(IX2-1)+IX1)+IHIST)
              IF(ABS(R4).LT.2.*SQRT(1./(OMEGA*AIMAG(CN(4,4))))) THEN
                GO TO 46
              END IF
   45       CONTINUE
            GO TO 79
   46       CONTINUE
C
C           Initial amplitudes of Gaussian packets
            IGRD=(((IXR-1)*NTR+ITR-1)*NP+IP-1)*NOMEGA+IOMEGA
            CGPAR=RAM(ICGPAR+IGRD)
            CGPAI=RAM(ICGPAI+IGRD)
CV          IF(((CGPAR.GT.-100.).AND.(CGPAR.LT.100.)).AND.
CV   *         ((CGPAI.GT.-100.).AND.(CGPAI.LT.100.))) THEN
CV            GO TO 79
CV          END IF
            IF(CGPAR.EQ.0..AND.CGPAI.EQ.0.) THEN
              GO TO 79
            END IF
C
C           Amplitudes of Gaussian packets
            CGPA=SQRT(CMPLX(0.,1.)*(Y(20)*Y(25)-Y(21)*Y(24))
     *                    /(ABS(Y(20)*Y(25)-Y(21)*Y(24))*CDET))
            CGPA=CMPLX(ABS(REAL(CGPA)),AIMAG(CGPA))
            CGPA=CGPA*CMPLX(CGPAR,CGPAI)*SQRT(VI/V)
     *         *CEXP(CMPLX(0.,1.)*PI/2.*(KMAH-0.5))
            CGPAR= REAL(CGPA)
            CGPAI=AIMAG(CGPA)
C
C           Limits of the loops over the grid points:
C           Imaginary part Y of the used submatrix of matrix N:
            Y22=AIMAG(CN(2,2))+2.*(KAPPA*Y(7)/DTRAY)**2
            Y23=AIMAG(CN(2,3))+2.*(KAPPA*Y(7)/DTRAY)*(KAPPA*Y(8)/DTRAY)
            Y33=AIMAG(CN(3,3))+2.*(KAPPA*Y(8)/DTRAY)**2
            Y24=AIMAG(CN(2,4))
            Y34=AIMAG(CN(3,4))
            Y44=AIMAG(CN(4,4))
C           Determinant of matrix Y
            YDET=Y22*Y33*Y44+2.*Y23*Y24*Y34
     *          -Y22*Y34*Y34-Y33*Y24*Y24-Y44*Y23*Y23
C           Diagonal elements of the inverse matrix
            X22=(Y33*Y44-Y34*Y34)/YDET
            X33=(Y22*Y44-Y24*Y24)/YDET
C           X44=(Y22*Y33-Y23*Y23)/YDET
C           Limits of the loops
            ID1=NINT(SQRT(X33*RELAMP*2./OMEGA))
            ID2=NINT(SQRT(X22*RELAMP*2./OMEGA))
            IX1=NINT((Y(5)-O1)/D1)+1
            IX2=NINT((Y(4)-O2)/D2)+1
            IX1U=IX1-ID1
            IX1D=IX1+ID1
            IX2L=IX2-ID2
            IX2R=IX2+ID2
            IF(IX1U.LT.IXTZU) THEN
              IX1U=IXTZU
            END IF
            IF(IX1D.GT.IXTZD) THEN
              IX1D=IXTZD
            END IF
            IF(IX2L.LT.IXTZL) THEN
              IX2L=IXTZL
            END IF
            IF(IX2R.GT.IXTZR) THEN
              IX2R=IXTZR
            END IF
C
            YDET=Y22*Y44-Y24*Y24
            X22=Y44/YDET
            YINVY3=(Y44*Y23-Y24*Y34)/YDET
            ID2=NINT(SQRT(X22*RELAMP*2./OMEGA))
            AUX1=KAPPA*Y(8)/DTRAY
            AUX2=KAPPA*Y(7)/DTRAY
            EXPR00=-OCN33I-AUX1*AUX1
            EXPR10=-OCN32I-AUX1*(AUX2+AUX2)
            EXPR2=-OCN22I-AUX2*AUX2
C
C           Loops over the target zone
            IX1=IX1U
   50       CONTINUE
              IX1=IX1+1
              IF(IX1.GT.IX1D) THEN
                GO TO 79
              END IF
              X1=O1+D1*(IX1-1)
              R1=X1-Y(5)
              EXPR0=EXPR00*R1*R1
              EXPR1=EXPR10*R1
              EXPI0=(OW1+OCN33R*R1)*R1
              EXPI1=OW2+OCN32R*R1
              EXPI2=OCN22R
              CAUXR0=OW4+OCN34R*R1
              CAUXI0=OCN34I*R1
C
C             Improving the limits of the inner loop:
C-C           R2=(Y44*Y23*R1-Y24*Y34*R1)/YDET
              R2=YINVY3*R1
              IX2=NINT((Y(4)-R2-O2)/D2)+1
              IX2L=IX2-ID2
              IX2R=IX2+ID2
              IF(IX2L.LT.IXTZL) THEN
                IX2L=IXTZL
              END IF
              IF(IX2R.GT.IXTZR) THEN
                IX2R=IXTZR
              END IF
C
              IX2=IX2L
   60         CONTINUE
                IX2=IX2+1
                IF(IX2.GT.IX2R) THEN
                  GO TO 50
                END IF
                IGRDB=N1*(IX2-1)+IX1
                IGRDC=N2*(N1-IX1)+IX2
                X2=O2+D2*(IX2-1)
                R2=X2-Y(4)
                GPEXPR=EXPR0+(EXPR1+EXPR2*R2)*R2
                GPEXPI=EXPI0+(EXPI1+EXPI2*R2)*R2
                CAUXR1=CAUXR0+OCN24R*R2
                CAUXI1=CAUXI0+OCN24I*R2
                CGP1R=0.
                CGP2R=0.
                NHIST=IRAM(IGRDB)
C
C               Loop over ray history
                K=IRAM(INUM+IGRDB)
                DO 70 J=K+1,K+NHIST
                  R4=R4A+RAM(IMTT+J)
                  CAUX2R=GPEXPR-(CAUXI1+OCN44I*R4)*R4
                  CAUX2I=GPEXPI+(CAUXR1+OCN44R*R4)*R4
C                 Weighting factor?
                  RFAR=RAM(IAPR+J)*1.E10
                  AUX=EXP(CAUX2R)
                  GPIMGR=COS(CAUX2I)
                  GPIMGI=SIN(CAUX2I)
                  CGP1R=CGP1R+RFAR*(CGPAR*GPIMGR-CGPAI*GPIMGI)*AUX
                  CGP2R=CGP2R+RFAR*RFAR
   69             CONTINUE
   70           CONTINUE
                IF(CGP2R.NE.0.) THEN
                  I=IGPSEC+IGRDC
                  RAM(I)=RAM(I)+CGP1R/CGP2R
                END IF
              GO TO 60
            GO TO 50
   79       CONTINUE
   80     CONTINUE
   90   CONTINUE
      GO TO 20
  100 CONTINUE
C
C     Null output array
      DO 105 I=1,N1N2MA
        RAM(IGPSEC+N1N2+I)=0.
  105 CONTINUE
      D= D1/4.
      ID=NINT(D)
      IF(ID.EQ.1) THEN
C       Without interpolation
        DO 107 I=1,N1N2MA
          RAM(IGPSEC+N1N2+I)=RAM(IGPSEC+I)
  107   CONTINUE
        GO TO 160
      END IF
C
C     Interpolating output grid values
      WRITE(*,'(A)')'+GPMIG: Interpolating... '
      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
          IGRDNE=N2MAX*(IX1NEW-1)+IX2NEW
          RAM(IGPSEC+N1N2+IGRDNE)=RAM(IGPSEC+IGRD)
  110   CONTINUE
  115 CONTINUE
      PID=PI/D
      DO 130 IX2=1,N2
        IX2NEW=(IX2-1)*ID+1
        DO 125 IX1=1,N1
          IX1NEW=(IX1-1)*ID+1
          IGRDNE=IGPSEC+N1N2+N2MAX*(IX1NEW-1)+IX2NEW
          DO 120 I=1-IX1NEW,N1MAX-IX1NEW
            IF(I.NE.0) THEN
              IGRD=IGRDNE+N2MAX*I
              AUX=PID*FLOAT(I)
              RAM(IGRD)=RAM(IGRD)+RAM(IGRDNE)*SIN(AUX)/AUX
            END IF
  120     CONTINUE
  125   CONTINUE
  130 CONTINUE
      DO 150 I=1,N1MAX
        DO 145 IX2=0,N2-1
          IX2NEW=IX2*ID+1
          IGRDNE=IGPSEC+N1N2+N2MAX*(I-1)+IX2NEW
          DO 140 J=1-IX2NEW,N2MAX-IX2NEW
            IF(J.NE.0) THEN
              IGRD=IGRDNE+J
              AUX=PID*FLOAT(J)
              RAM(IGRD)=RAM(IGRD)+RAM(IGRDNE)*SIN(AUX)/AUX
            END IF
  140     CONTINUE
  145   CONTINUE
  150 CONTINUE
C
  160 CONTINUE
C     Write migrated section to output file
      IF(GPMSEC.NE.' ') THEN
        CALL OARRAY(LU1,GPMSEC,'W')
        CALL WARRAY(LU1,' ',.FALSE.,0.,.FALSE.,0.,
     *                             N1N2MA,RAM(IGPSEC+N1N2+1))
        CLOSE(LU1)
      END IF
C
      WRITE(*,'(A)')'+GPMIG: Done. '
      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.for
C
C=======================================================================
C