C
C Program GPSTEP to calculate the discretization steps
C for (a) the decomposition of the wavefield into optimized
C 2-D Gaussian packets and (b) migration.
C
C Version: 7.40
C Date: 2017, May 19
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.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:
C     RYMM='string'... Name of the input file containing the
C             minimum and maximal values of the optimized and smoothed
C             initial parameters of Gaussian beams R0MIN,R0MAX,
C             Y0MIN,Y0MAX.
C             Default: RYMM='ryminmax.out'
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     DT=real... Grid interval along the time axis T of the decomposed
C             wavefield.
C             Default: DT=1.
C     OT=real, OX=real ... Coordinates of the origin of the grid.
C             Default: OT=0., OX=0.
C     ISRC=positive integer...
C             Default: ISRC=1
C     VSURF=real... Velocity determined according to maximum velocity
C             in the depth of the receivers.
C             Default: VSURF=2000.
C     FMIN=real... Minimum frequency of the Gaussian packets.
C             Default: FMIN=1.
C     FMAX=real... Maximum frequency of the Gaussian packets.
C             Default: FMAX=100.
C     FREF=real... Reference frequencies of the Gaussian packets.
C             Default: FREF=50.
C     PKSI=real... Parameter ksi with values in the interval (0,1).
C             Default: PKSI=0.5
C     PNY=real... Parameters ny with values in the interval (0,1).
C             Default: PNY=0.5
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 Name of the output file:
C     GPSTEP='string'... Name of the output 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='gpstep.out'
C     GBGRD='string'... Name of the output file containing new
C             grid parameters of the Hamiltonian hypersurface (XR,TR,P)
C             computed for smoothed and optimized Im(N0) and Re(N0).
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
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL RARRAY,LENGTH,RSEP1,RSEP3T,RSEP3I,RSEP3R,ERROR
      INTEGER  LENGTH
C
C     Filenames and parameters:
      CHARACTER*80 FILSEP,RYMM,GPSTEP,GBGRD
      CHARACTER*80 FGBY22,FGBR22,GPSTEPI
      INTEGER LU1
      REAL PI
      PARAMETER (LU1=1,PI=3.141592654)
C
C     Other variables:
      INTEGER N,NX,NT,NXR,NTR,NP,NOMEGA
      INTEGER N1,N2,N3,N1NEW,N2NEW,N3NEW,ISRC,IVER
      INTEGER IFGBY0,IFGBR0,IFGBY,IFGBR,IFGB
      REAL DX,DT,FMAX,FMIN,FREF,OMEMAX,OMEREF,VSURF
      REAL DXR,OXR,DTR,OTR,DP,OP,DOMEGA,OOMEGA,OOME1,D1,D2,D3,O1,O2,O3
      REAL R0MIN,R0MAX,Y0MIN,Y0MAX,PKSI,PNY,PKSIT,DSRC,OXRINI
      REAL ANIMIN,ANIMAX,ANRMIN,ANRMAX,AK0,AK0T,KAPPA
      COMPLEX CN,CNT,CN0,CN44,CN44T,CPT
      COMPLEX CK11,CK12,CK22,CKDET,CKI11,CKI12,CKI22
C
C.......................................................................
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPSTEP: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPSTEP-01
        CALL ERROR('GPSTEP-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)')'+GPSTEP: Reading, calculating... '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading filename of the input file
      CALL RSEP3T('RYMM',RYMM,'ryminmax.out')
      CALL RSEP3T('FGBY22',FGBY22,'y0m.out')
      CALL RSEP3T('FGBR22',FGBR22,'r0m.out')
      CALL RSEP3T('GPSTEPI',GPSTEPI,'gpstep.ini')
C
C     Reading input parameters
      CALL RSEP3I('NX',NX,116)
      CALL RSEP3R('DX',DX,25.)
      CALL RSEP3R('OX',OX,0.)
      CALL RSEP3I('NT',NT,750)
      CALL RSEP3R('DT',DT,0.004)
      CALL RSEP3R('OT',OT,0.)
C
      CALL RSEP3I('ISRC',ISRC,1)
      CALL RSEP3R('DSRC',DSRC,1.)
      CALL RSEP3R('OXRINI',OXRINI,1.)
C     Grid parameters of the Hamiltonian hypersurface (NTR,NP,NXR)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,1.)
      CALL RSEP3R('O2',O2,1.)
      CALL RSEP3R('O3',O3,1.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
C
      CALL RSEP3R('VSURF',VSURF,2000.)
      CALL RSEP3R('FMAX',FMAX,1.)
      CALL RSEP3R('FMIN',FMIN,1.)
      CALL RSEP3R('FREF',FREF,1.)
      CALL RSEP3R('PKSI',PKSI,0.5)
      CALL RSEP3R('PNY',PNY,0.5)
      CALL RSEP3R('PKSIT',PKSIT,0.5)
      CALL RSEP3R('KAPPA',KAPPA,1.253314137)
      CALL RSEP3I('IVER',IVER,1)
C
C     Reading filenames of the output files
      CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.out')
      CALL RSEP3T('GBGRD',GBGRD,'gbgrd.out')
      CLOSE(LU1)
C
C     Read maximum and minimum values of the optimized and smoothed
C     initial parameters of Gaussian beams R0, Y0
      OPEN(LU1,FILE=RYMM,FORM='FORMATTED',STATUS='OLD')
      READ(LU1,*) R0MIN,R0MAX,Y0MIN,Y0MAX
      CLOSE(LU1)
      IF(Y0MIN.LT.0.) THEN
C       GPSTEP-02
        CALL ERROR('GPSTEP-02: Y0 contains negative values')
      ENDIF
C
C     N=R0+iY0
      ANRMIN=R0MIN
      ANRMAX=R0MAX
      ANIMIN=Y0MIN
      ANIMAX=Y0MAX
      WRITE(*,'(A,1(G16.6,1X))') 'ANRMIN =',ANRMIN
      WRITE(*,'(A,1(G16.6,1X))') 'ANRMAX =',ANRMAX
      WRITE(*,'(A,1(G16.6,1X))') 'ANIMIN =',ANIMIN
      WRITE(*,'(A,1(G16.6,1X))') 'ANIMAX =',ANIMAX
C
      OMEMAX=2.*PI*FMAX
      OOMEGA=2.*PI*FMIN
      OMEREF=2.*PI*FREF
      IF(OOMEGA.LT.(1./(DT*(NT-1)))) THEN
        OOMEGA= 1./(DT*(NT-1))
      END IF
C
C     Calculation of parameters N0 and K0 equations (33), (34)
C     (see Zacek, 2003)
      CN0=PNY*2.*CMPLX(0.,ANIMIN)
      AK0=PKSI*OOMEGA*ANIMIN*VSURF**2
C
C     Calculation of parameters N0, K0 and ~K0 according to ????
C     CN0=CMPLX(ANRMIN,ANIMIN)
C     AK0=0.5*OOMEGA*ANIMIN*VSURF**2
C     AK0T=AK0
C
C-----------------------------------------------------------------------
C     Version 1
      IF(IVER.EQ.1) THEN
C
C       Calculation of parameter ~K0, equation (35new),(Zacek,2003)
        AK0T=PKSIT*OOMEGA*ANIMIN*VSURF**2
C
C       Reading initial discretization steps
        OPEN(LU1,FILE=GPSTEPI,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,OOME1
C       Reading complex-valued constants
C       READ(LU1,*) AN0R,AN0I
C       READ(LU1,*) AK0,AK0T,Y0MIN
        CLOSE(LU1)
C
C       Writing discretization steps to the screen
        WRITE(*,*) ' '
        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)')          ' NXR   =',NXR
        WRITE(*,'(A,1(G16.6,1X))') ' DXR   =',DXR
        WRITE(*,'(A,1(G16.6,1X))') ' OXR   =',OXR
        WRITE(*,'(A,I6)')          ' NOMEGA=',NOMEGA
        WRITE(*,'(A,1(G16.6,1X))') ' DOMEGA=',DOMEGA
        WRITE(*,'(A,1(G16.6,1X))') ' OOME1=',OOME1
C
C       Writing complex-valued constants to the screen
        WRITE(*,'(A,1(G16.6,1X))') ' N0R   =',REAL(CN0)
        WRITE(*,'(A,1(G16.6,1X))') ' N0I   =',AIMAG(CN0)
        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))') ' PKSI=',PKSI
        WRITE(*,'(A,1(G16.6,1X))') ' PNY =',PNY
        WRITE(*,'(A,1(G16.6,1X))') ' PKSIT =',PKSIT
        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))') ' DT    =',DT
        WRITE(*,'(A,1(G16.6,1X))') ' OT    =',OT
        WRITE(*,*) ' '
C
        NXTP=NXR*NTR*NP
        IFGBY0=0
        IFGBR0=IFGBY0+NXTP
        DXRMIN=999999.
        DTRMIN=999999.
        AUXMIN=999999.
C
C       Read initial parameters of Gaussian beams - imaginary part (Y0)
        OPEN(LU1,FILE=FGBY22,FORM='FORMATTED',STATUS='OLD')
        CALL RARRAY(LU1,' ','FORMATTED',.FALSE.,0.,NXTP,
     *              RAM(1))
        CLOSE(LU1)
C
C       Read initial parameters of Gaussian beams - real part (R0)
        OPEN(LU1,FILE=FGBR22,FORM='FORMATTED',STATUS='OLD')
        CALL RARRAY(LU1,' ','FORMATTED',.FALSE.,0.,NXTP,
     *              RAM(NXTP+1))
        CLOSE(LU1)
        DO 50 IXR=1,NXR
C         WRITE(*,'(2(A,I3))') '+GPSTEP: NXR=',IXR,'/',NXR
C
C         Loop over corresponding arrival time of Gaussian packet
          DO 45 ITR=1,NTR
C
C           Loop over component of slowness vector of Gaussian packet
C           along profile
            DO 40 IP=1,NP
              P=OP+DP*(IP-1)
              DO 35 IOMEGA=1,NOMEGA
                OMEGA=OOME1+DOMEGA*(IOMEGA-1)
                IF(OMEGA.GT.OMEREF) THEN
                  OMEGA=OMEREF
                END IF
                PMAX=PI/(OMEGA*DX)
                IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN
                  IFGB=((IXR-1)*NP+IP-1)*NTR+ITR
                  IFGBY=IFGBY0+IFGB
                  IFGBR=IFGBR0+IFGB
C
C                 Complex valued parameter N determines initial shape
C                 of Gaussian packets
                  CN=CMPLX(RAM(IFGBR),RAM(IFGBY))
C
C                 Computation of ~N, equation (16),(Zacek,2003)
C                                             (28),(Zacek,2006)
                  CNT=CN*CN0/(2.*CN-CN0)
C
C                 Computation of ~p, equation (17),(Zacek,2003)
C                                             (29),(Zacek,2006)
                  CPT=-P*CNT/CN
C
C                 Computation of ~N44, equation (18),(Zacek,2003)
C                                               (30),(Zacek,2006)
                  CN44T=CNT*CMPLX(0.,AK0T)
     *                /(OMEGA*CNT-CMPLX(0.,AK0T)*CPT**2)
C
C                 Computation of N44, equation (19),(Zacek,2003)
C                                              (31),(Zacek,2006)
                  CN44=CN*CMPLX(0.,AK0)
     *               /(OMEGA*CN-CMPLX(0.,AK0)*P**2)
C
C                 Computation of equations (6,10),(Zacek,2003)
C                                          (14,17),(Zacek,2006)
                  CK11=CN+CN44*P*P+CNT+CN44T*CPT*CPT
                  CK12=-CN44*P-CN44T*CPT
                  CK22=CN44+CN44T
                  CKDET=CK11*CK22-CK12*CK12
                  IF(CKDET.EQ.0.) THEN
                    CALL ERROR('GPSTEP-XX: CKDET equal zero')
                  END IF
                  CKI11=CK22/CKDET
                  CKI12=-CK12/CKDET
                  CKI22=CK11/CKDET
C
C                 Computation of equations (25,26),(Zacek,2003)
C                                          (45,46),(Zacek,2006)
                  DELTAX=KAPPA*SQRT(-AIMAG(CKI11)*2./OMEGA)
                  DELTAT=KAPPA*SQRT(-AIMAG(CKI22)*2./OMEGA)
C
C                 Computation of equation (47),(Zacek,2006)
                  AUX=0.5*SQRT(AIMAG(CKI11)*AIMAG(CKI22))
     &                -ABS(AIMAG(CKI12))
                  DXRMIN=AMIN1(DXRMIN,DELTAX)
                  DTRMIN=AMIN1(DTRMIN,DELTAT)
                  AUXMIN=AMIN1(AUXMIN,AUX)
C                 WRITE(2,*)DELTAX,DELTAT,AUX
                END IF
                IF(OMEGA.EQ.OMEREF) THEN
                  GO TO 36
                END IF
   35         CONTINUE
   36         CONTINUE
   40       CONTINUE
   45     CONTINUE
   50   CONTINUE
C
        N=NINT((NX-1)*DX/DXRMIN)
        NXR=N+2
        DXR=(NX-1)*DX/(N+1)
C
        IF(OOMEGA.LT.OMEMAX) THEN
C         Computation of maximum discretization step DOMEGA according to
C         equation (30), (Zacek, 2003)
          DOMEGA=KAPPA*SQRT(2.*AK0*AK0T/(AK0+AK0T))
          N=NINT((OMEMAX-OOMEGA)/DOMEGA)
          NOMEGA=N+2
          DOMEGA=(OMEMAX-OOMEGA)/(N+1)
        ELSE
          OOMEGA=OMEMAX
          DOMEGA=1.
          NOMEGA=1
        END IF
C
C       Computation of maximum discretization step DP according to
C       equation (29), (Zacek, 2003)
        DP=KAPPA*SQRT(AIMAG(CN0)/OMEREF)
C
        N=NINT(2./(VSURF*DP))
        NP=N+2
        DP=2./(VSURF*(N+1))
C
        N=NINT((NT-1)*DT/DTRMIN)
        NTR=N+2
        DTR=(NT-1)*DT/(N+1)
C
        OXR=OXRINI+(ISRC-1)*DSRC
C
C       OXR=1587.5-FLOAT(NXR-1)*DXR/2.+ISRC*25.
C       OXR=OX+FLOAT(NX-1)*DX/2.+FLOAT(ISRC-1)*DSRC
        OP=-1./VSURF
        OTR=0.
C
        AUX=DXRMIN*DTRMIN*DP*DOMEGA
        WRITE(*,*)DXRMIN,DTRMIN,DP,DOMEGA,AUXMIN,AUX
        WRITE(*,*)NXR,NTR,NP
C
C       Computing new grid parameters of the Hamiltonian hypersurface
        N1NEW=NTR
        N2NEW=NP
        N3NEW=NXR
        D1=FLOAT(N1-1)/FLOAT(N1NEW-1)
        D2=FLOAT(N2-1)/FLOAT(N2NEW-1)
        D3=FLOAT(N3-1)/FLOAT(N3NEW-1)
        OPEN(LU1,FILE=GBGRD,FORM='FORMATTED')
        CALL WSEP3I(1,'N1',N1NEW)
        CALL WSEP3I(1,'N2',N2NEW)
        CALL WSEP3I(1,'N3',N3NEW)
        CALL WSEP3R(1,'D1',D1)
        CALL WSEP3R(1,'D2',D2)
        CALL WSEP3R(1,'D3',D3)
        CALL WSEP3R(1,'O1',O1)
        CALL WSEP3R(1,'O2',O2)
        CALL WSEP3R(1,'O3',O3)
        CLOSE(LU1)
C
C       Writing discretization steps to the output file GPSTEP
        WRITE(*,'(A)')'+GPSTEP: Writing output file.'
        OPEN(LU1,FILE=GPSTEP,FORM='FORMATTED')
        WRITE(LU1,'(3A)') '''',GPSTEP(1:LENGTH(GPSTEP)),''''
        WRITE(LU1,*) NXR,DXR,OXR,' / (NXR,DXR,OXR)'
        WRITE(LU1,*) NTR,DTR,OTR,' / (NTR,DTR,OTR)'
        WRITE(LU1,*) NP, DP, OP ,' / (NP,DP,OP)'
        WRITE(LU1,*) NOMEGA,DOMEGA,OOMEGA,' / (NOMEGA,DOMEGA,OOMEGA)'
C
C       Writing complex-valued constants to the output file GPSTEP
        WRITE(LU1,*) REAL(CN0),AIMAG(CN0),' / (Re(N0),Im(N0))'
        WRITE(LU1,*) AK0,AK0T,ANIMIN,' / (K0,~K0,Im(N)min)'
        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))') ' Re(N0)  =',REAL(CN0)
        WRITE(*,'(A,1(G16.6,1X))') ' Im(N0)  =',AIMAG(CN0)
        WRITE(*,'(A,1(G16.6,1X))') ' K0      =',AK0
        WRITE(*,'(A,1(G16.6,1X))') ' ~K0     =',AK0T
        WRITE(*,'(A,1(G16.6,1X))') ' Im(N)min=',ANIMIN
C
C-----------------------------------------------------------------------
C     Version 2 (old)
      ELSE
C
C       Calculation of parameter ~K0, equation (35),(see Zacek,2003)
        AK0T=AK0*(1.-PNY)/PNY
C
        IF(OOMEGA.LT.OMEMAX) THEN
C         Computation of maximum discretization step DOMEGA according to
C         equation (30), (Zacek, 2003)
          DOMEGA=KAPPA*SQRT(2.*AK0*AK0T/(AK0+AK0T))
          N=NINT((OMEMAX-OOMEGA)/DOMEGA)
          NOMEGA=N+2
          DOMEGA=(OMEMAX-OOMEGA)/(N+1)
        ELSE
          OOMEGA=OMEMAX
          DOMEGA=1.
          NOMEGA=1
        END IF
C
C       Computation of maximum discretization step DXR according to
C       equation (38), (Zacek, 2003)
        DXR=KAPPA*SQRT(-2.*(1.-PNY)/OMEREF*AIMAG(
     *     (CMPLX(1.,0.)-OOMEGA/OMEREF*
     *      CMPLX(0.,ANIMIN)/CMPLX(0.,ANIMAX)*PKSI)/
     *      CMPLX(0.,ANIMAX)))
C
C       Computation of DXR according to ?????
C       DXR=KAPPA*SQRT(-AIMAG(0.5/(OMEREF*CN0)))
C
        N=NINT((NX-1)*DX/DXR)
        NXR=N+2
        DXR=(NX-1)*DX/(N+1)
C
C       Computation of maximum discretization step DP according to
C       equation (29), (Zacek, 2003)
        DP=KAPPA*SQRT(AIMAG(CN0)/OMEREF)
C
        N=NINT(2./(VSURF*DP))
        NP=N+2
        DP=2./(VSURF*(N+1))
C
C       Computation of maximum discretization step DTR according to
C       equation (40), (Zacek, 2003)
        DTR=KAPPA*SQRT(-2.*PNY*AIMAG(
     *     (CMPLX(1.,0.)-CMPLX(0.,ANIMIN)/CMPLX(0.,ANIMIN)*PKSI)
     *     /CMPLX(0.,AK0)))
C
C       Computation of DTR according to ?????
C       DTR=KAPPA*SQRT(0.5/AK0)
C
        N=NINT((NT-1)*DT/DTR)
        NTR=N+2
        DTR=(NT-1)*DT/(N+1)
C
        OXR=OXRINI+(ISRC-1)*DSRC
C
C       OXR=1587.5-FLOAT(NXR-1)*DXR/2.+ISRC*25.
C       OXR=OX+FLOAT(NX-1)*DX/2.+FLOAT(ISRC-1)*DSRC
        OP=-1./VSURF
        OTR=0.
C
C       Computing new grid parameters of the Hamiltonian hypersurface
        N1NEW=NTR
        N2NEW=NP
        N3NEW=NXR
        D1=FLOAT(N1-1)/FLOAT(N1NEW-1)
        D2=FLOAT(N2-1)/FLOAT(N2NEW-1)
        D3=FLOAT(N3-1)/FLOAT(N3NEW-1)
        OPEN(LU1,FILE=GBGRD,FORM='FORMATTED')
        CALL WSEP3I(1,'N1',N1NEW)
        CALL WSEP3I(1,'N2',N2NEW)
        CALL WSEP3I(1,'N3',N3NEW)
        CALL WSEP3R(1,'D1',D1)
        CALL WSEP3R(1,'D2',D2)
        CALL WSEP3R(1,'D3',D3)
        CALL WSEP3R(1,'O1',O1)
        CALL WSEP3R(1,'O2',O2)
        CALL WSEP3R(1,'O3',O3)
        CLOSE(LU1)
C
C       Writing discretization steps to the output file GPSTEP
        WRITE(*,'(A)')'+GPSTEP: Writing output file.'
        OPEN(LU1,FILE=GPSTEP,FORM='FORMATTED')
        WRITE(LU1,'(3A)') '''',GPSTEP(1:LENGTH(GPSTEP)),''''
        WRITE(LU1,*) NXR,DXR,OXR,' / (NXR,DXR,OXR)'
        WRITE(LU1,*) NTR,DTR,OTR,' / (NTR,DTR,OTR)'
        WRITE(LU1,*) NP, DP, OP ,' / (NP,DP,OP)'
        WRITE(LU1,*) NOMEGA,DOMEGA,OOMEGA,' / (NOMEGA,DOMEGA,OOMEGA)'
C
C       Writing complex-valued constants to the output file GPSTEP
        WRITE(LU1,*) REAL(CN0),AIMAG(CN0),' / (Re(N0),Im(N0))'
        WRITE(LU1,*) AK0,AK0T,ANIMIN,' / (K0,~K0,Im(N)min)'
        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))') ' Re(N0)  =',REAL(CN0)
        WRITE(*,'(A,1(G16.6,1X))') ' Im(N0)  =',AIMAG(CN0)
        WRITE(*,'(A,1(G16.6,1X))') ' K0      =',AK0
        WRITE(*,'(A,1(G16.6,1X))') ' ~K0     =',AK0T
        WRITE(*,'(A,1(G16.6,1X))') ' Im(N)min=',ANIMIN
      END IF
C
      WRITE(*,'(A)') '+GPSTEP: 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 'forms.for'
C     forms.for
C
C=======================================================================
C