C
C Program GPANAL to calculate the amplitudes of 2-D Gaussian packets
C
C Version: 6.10
C Date: 2007, June 14
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:
C     TIMESEC='string'... Name of the input file containing the time
C             section.
C             Default: TIMESEC='timesec.out'
C     FGBY22='string'... Name of the input file with the optimum initial
C             parameters Y0 of Gaussian beams - imaginary part.
C             Default: FGBY22='fgby22.out'
C     FGBR22='string'... Name of the input file with the optimum initial
C             parameters R0 of Gaussian beams - real part.
C             Default: FGBR22='fgbr22.out'
C     FGBYMIN='string'... Name of the input file with the minimum
C             value of Y0.
C             Default: FGBY22='fgbymin.out'
C     FGBYMAX='string'... Name of the input file with the maximum
C             value of Y0.
C             Default: FGBY22='fgbymax.out'
C     FGBRMIN='string'... Name of the input file with the minimum
C             value of R0.
C             Default: FGBY22='fgbrmin.out'
C     FGBRMAX='string'... Name of the input file with the maximum
C             value of R0.
C             Default: FGBY22='fgbrmax.out'
C Name of the output files:
C     GPOUTR='string'... Name of the output file containing
C             amplitudes of Gaussian packets - real part.
C             Default: GPOUTR='gpampr.out'
C     GPOUTI='string'... Name of the output file containing
C             amplitudes of Gaussian packets - imaginary part.
C             Default: GPOUTI='gpampi.out'
C     SYNTSEC='string'... Name of the output file containing synthetic
C             seismograms.
C             Default: SSC='synts.out'
C Names of input parameters:
C     PARK=real... Parameter ksi with values in the interval (0,1).
C             Default: PARK=0.5
C     PARN=real... Parameters ny with values in the interval (0,1).
C             Default: PARN=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     GWM=real... ????????????????????.
C             Default: GWM=3.
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=1.
C     FREF=real... Reference frequencies of the Gaussian packets.
C             Default: FREF=1.
C Data specifying the dimensions of the input grid of synthetic
C             seismograms to be decomposed:
C     OT=real, OX=real ... Coordinates of the origin
C             of the grid.
C             Default: OT=0., OX=0.
C     DT=real... Grid interval along the T axis.
C             Default: DT=1.
C     DX=real... Grid interval along the X axis.
C             Default: DX=1.
C     NT=positive integer... Number of gridpoints along the T axis.
C             Default: NT=1
C     NX=positive integer... Number of gridpoints along the X axis.
C             Default: NX=1
C Data specifying the dimensions of the input grid (Hamiltonian
C             hypersurface) of initial parameters of Gaussian beams:
C     OTR=real, OXR=real, OPR=real ... Coordinates of the origin
C             of the grid.
C             Default: OT=0., OX=0., OP=0.
C     DTR=real... Grid interval along the T axis.
C             Default: DT=1.
C     DXR=real... Grid interval along the X axis.
C             Default: DX=1.
C     DPR=real... Grid interval along the P axis.
C             Default: DP=1.
C     NTR=positive integer... Number of gridpoints along the T axis.
C             Default: NT=1
C     NXR=positive integer... Number of gridpoints along the X axis.
C             Default: NX=1
C     NPR=positive integer... Number of gridpoints along the P axis.
C             Default: NPR=1
C     NXTP=positive integer... Product NXR*NTR*NPR.
C             Default: NXTP=1
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             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             forms.for.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      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-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
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
      INCLUDE 'sep.inc'
C     sep.inc
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL OX,OT,OP,DX,DT,DP,FMAX,FMIN,OMEGAMAX,VSURF,FREF,OMEGAREF
      REAL XX,T,P,XR,TR
      INTEGER LU1,LU2,LU3,LU4,LU5
      INTEGER IRAM(MRAM)
      INTEGER NX,NT,NP,NXT,NXTP,N
      INTEGER IXR,ITR,IPR,IOMEGA,IX,IT
      INTEGER IGRD,IGRD0,IFGBY0,IFGBR0,IFGBY,IFGBR,IFGB
      INTEGER IGRDFR,IGRDFI,IGRDF0,IGRDF1,IGRDF2,IGRDF3
C
      INTEGER NXR,NTR,NPR,NOMEGA
      REAL DXR,DTR,DPR,DOMEGA,OXR,OTR,OPR,OOMEGA,XRX,TRT,AX,BT
C
      REAL PR,OMEGA,DOMEGA1,KAPPA,PI,FR,FI,W,WW,WR,WI,GWM,DXTPO
C
      REAL DXT,PRMIN,PRMAX,ANIMIN,ANIMAX,ANRMIN,ANRMAX,PARK,PARN
C
      INTEGER IXMAX,IXMIN,ITMAX,ITMIN,IXR0,IPR0,MXR,MPR
C
      REAL AK0,AK0D,AMPLR,AMPLI,DXR1,DPR1
C
      COMPLEX AN,ANDD,AN0,THETA,AMPL,N44,N44D,PRD
      COMPLEX THETA1,THETA2,THETA3,THETA4
C
      PARAMETER (LU1=21,LU2=22,LU3=23,LU4=24,LU5=25,LU6=26,LU7=27)
      PARAMETER (LU8=28,LU9=29,LU10=30,LU11=31)
      PARAMETER (PI=3.141592654)
      CHARACTER*80 FILSEP,TIMESEC,FGBY22,FGBR22,GPOUTR,GPOUTI,SYNTSEC
      CHARACTER*80 FGBRMIN,FGBRMAX,FGBYMIN,FGBYMAX
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPANAL: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPANAL-01
        CALL ERROR('GPANAL-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)')'+GPANAL: 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('TIMESEC',TIMESEC,'times.out')
      CALL RSEP3T('FGBY22',FGBY22,'fgby22.out')
      CALL RSEP3T('FGBR22',FGBR22,'fgbr22.out')
      CALL RSEP3T('FGBYMIN',FGBYMIN,'fgbymin.out')
      CALL RSEP3T('FGBYMAX',FGBYMAX,'fgbymax.out')
      CALL RSEP3T('FGBRMIN',FGBRMIN,'fgbrmin.out')
      CALL RSEP3T('FGBRMAX',FGBRMAX,'fgbrmax.out')
C
C     Reading filenames of the output files
      CALL RSEP3T('GPOUTR',GPOUTR,'gpampr.out')
      CALL RSEP3T('GPOUTI',GPOUTI,'gpampi.out')
      CALL RSEP3T('SYNTSEC',SYNTSEC,'synts.out')
C
      CALL RSEP3I('NX',NX,1)
      CALL RSEP3I('NT',NT,1)
      CALL RSEP3I('NXTP',NXTP,1)
C
      CALL RSEP3R('DX',DX,1.)
      CALL RSEP3R('DT',DT,1.)
      CALL RSEP3R('OX',OX,0.)
      CALL RSEP3R('OT',OT,0.)
C
      CALL RSEP3R('PARK',PARK,0.5)
      CALL RSEP3R('PARN',PARN,0.5)
      CALL RSEP3R('KAPPA',KAPPA,1.253314137)
      CALL RSEP3R('GWM',GWM,3.)
      CALL RSEP3R('VSURF',VSURF,2000.)
      CALL RSEP3R('FMAX',FMAX,1.)
      CALL RSEP3R('FMIN',FMIN,1.)
      CALL RSEP3R('FREF',FREF,1.)
C
      OPEN(LU2,FILE=TIMESEC,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FGBY22,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU4,FILE=FGBR22,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU5,FILE=GPOUTR,FORM='FORMATTED')
      OPEN(LU6,FILE=GPOUTI,FORM='FORMATTED')
      OPEN(LU7,FILE=SYNTSEC,FORM='FORMATTED')
      OPEN(LU8,FILE=FGBYMIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU9,FILE=FGBYMAX,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU10,FILE=FGBRMIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU11,FILE=FGBRMAX,FORM='FORMATTED',STATUS='OLD')
C
      NXT=NX*NT
C      NXTP=NXT*NP
C
      CALL RARRAY(LU2,' ','FORMATTED',.FALSE.,0.,NXT,RAM)
C
      CALL RARRAY(LU3,' ','FORMATTED',.FALSE.,0.,NXTP,
     *            RAM(NXT+1))
      CALL RARRAY(LU4,' ','FORMATTED',.FALSE.,0.,NXTP,
     *            RAM(NXT+NXTP+1))
      READ(LU8,*)  ANIMIN
      READ(LU9,*)  ANIMAX
      READ(LU10,*) ANRMIN
      READ(LU11,*) ANRMAX
C
      OMEGAMAX=2.0*PI*FMAX
      OOMEGA=2.0*PI*FMIN
      OMEGAREF=2.0*PI*FREF
      IF(OOMEGA.LT.(1.0/(DT*(NT-1)))) THEN
        OOMEGA= 1.0/(DT*(NT-1))
      END IF
C
C      AK0=PARK*OOMEGA*ANIMIN*VSURF**2
C      AK0D=AK0*(1.0-PARN)/PARN
C      AN0=PARN*2.0*CMPLX(ANRMIN,ANIMIN)
      AK0=0.5*OOMEGA*ANIMIN*VSURF**2
      AK0D=AK0
      AN0=CMPLX(ANRMIN,ANIMIN)
C
      WRITE(*,'(A,1(G16.6,1X))') 'K0    =',AK0
      WRITE(*,'(A,1(G16.6,1X))') 'K0D   =',AK0D
      WRITE(*,'(A,1(G16.6,1X))') 'N0R   =',REAL(AN0)
      WRITE(*,'(A,1(G16.6,1X))') 'N0I   =',AIMAG(AN0)
      WRITE(*,'(A,1(G16.6,1X))') 'NMIN  =',ANIMIN
      WRITE(*,'(A,1(G16.6,1X))') 'NMAX  =',ANIMAX
C
      IF(OOMEGA.LT.OMEGAMAX) THEN
        DOMEGA=KAPPA*SQRT(2.0*AK0*AK0D/(AK0+AK0D))
C        DOMEGA1=2.0*PI/(NT*DT)
C        IF(DOMEGA1.LT.DOMEGA) THEN
C          WRITE(*,'(A)')'+GPANAL: Aliasing '
C        END IF
        N=NINT((OMEGAMAX-OOMEGA)/DOMEGA)
        NOMEGA=N+2
        DOMEGA=(OMEGAMAX-OOMEGA)/(N+1)
      ELSE
        OOMEGA=OMEGAMAX
        DOMEGA=1.0
        NOMEGA=1
      END IF
C      WRITE(*,'(A,1(G16.6,1X))') 'FMIN  =',OOMEGA/(2.0*PI)
      WRITE(*,'(A,1(G16.6,1X))') 'DF    =',DOMEGA/(2.0*PI)
      WRITE(*,'(A,I6)')          'NF    =',NOMEGA
C
      OXR=OX
C      DXR=KAPPA*SQRT(-2.0*(1.0-PARN)/OMEGAREF*AIMAG(
C     *   (CMPLX(1.0,0.0)-OOMEGA/OMEGAREF*
C     *    CMPLX(0.0,ANIMIN)/CMPLX(ANRMAX,ANIMAX)*PARK)/
C     *    CMPLX(ANRMAX,ANIMAX)))
      DXR=KAPPA*SQRT(-AIMAG(0.5/(OMEGAREF*AN0)))
      N=NINT((NX-1)*DX/DXR)
      NXR=N+2
      DXR=(NX-1)*DX/(N+1)
      WRITE(*,'(A,1(G16.6,1X))') 'DXR   =',DXR
      WRITE(*,'(A,I6)')          'NXR   =',NXR
C
      OTR=OT
C      DTR=KAPPA*SQRT(-2.0*PARN*AIMAG(
C     *   (CMPLX(1.0,0.0)-CMPLX(0.0,ANIMIN)/CMPLX(ANRMIN,ANIMIN)*PARK)
C     *   /CMPLX(0.0,AK0)))
      DTR=KAPPA*SQRT(0.5/AK0)
      N=NINT((NT-1)*DT/DTR)
      NTR=N+2
      DTR=(NT-1)*DT/(N+1)
      WRITE(*,'(A,1(G16.6,1X))') 'DTR   =',DTR
      WRITE(*,'(A,I6)')          'NTR   =',NTR
C
      OPR=-1.0/VSURF
      DPR=KAPPA*SQRT(AIMAG(AN0)/OMEGAREF)
      N=NINT(2.0/(VSURF*DPR))
      NPR=N+2
      DPR=2.0/(VSURF*(N+1))
      WRITE(*,'(A,1(G16.6,1X))') 'DP    =',DPR
      WRITE(*,'(A,I6)')          'NP    =',NPR
C
      NXTPO=NXR*NTR*NPR*NOMEGA
      DXTPO=DXR*DTR*DPR*DOMEGA
      DXT=DX*DT
C
      WRITE(*,'(A,I6)') 'Expected time is',
     *                   NINT(NXTPO/237006.0*NXT/87000.0*1.5*7200.0)
C
ccc   IFGBY0=NXT+1
      IFGBY0=NXT
      IFGBR0=IFGBY0+NXTP
      IGRDF0=NXT+2*NXTP
C
      IF((IGRDF0+NXTPO).GT.MRAM) THEN
C       GPANAL-02
        CALL ERROR('GPANAL-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
C
C     Decomposition of the time section into 2-D Gaussian packets
C     -----------------------------------------------------------
C
      WRITE(*,'(A)')'+GPANAL: Decomposition '
C
      BT=GWM/SQRT(AK0D)
C
      DO 500 IXR=1,NXR
        XR=OXR+DXR*(IXR-1)
C
        WRITE(*,'(A,I6)')          'IXR   =',IXR
C
        IGRDF3=(IXR-1)*NTR
C
        DO 450 ITR=1,NTR
          TR=OTR+DTR*(ITR-1)
C
          IGRDF2=(IGRDF3+ITR-1)*NPR
C
          DO 400 IPR=1,NPR
            PR=OPR+DPR*(IPR-1)
C
            IGRDF1=IGRDF0+(IGRDF2+IPR-1)*NOMEGA
C
            IFGB=((IXR-1)*NPR+IPR-1)*NTR+ITR
            IFGBY=IFGBY0+IFGB
            IFGBR=IFGBR0+IFGB
C
            AN=CMPLX(RAM(IFGBR),RAM(IFGBY))
            ANDD=AN*AN0/(2.0*AN-AN0)
C
            PRD=-PR*ANDD/AN
C
            DO 350 IOMEGA=1,NOMEGA
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
C
              DXR1=KAPPA*SQRT(-2.0*(1.0-PARN)/OMEGA*AIMAG(
     *            (CMPLX(1.0,0.0)-OOMEGA/OMEGA
     *            *CMPLX(0.0,ANIMIN)/AN*PARK)/AN))
              MXR=INT(DXR1/DXR)
              IF(MXR.LT.1) THEN
                MXR=1
              END IF
              IXR0=INT(FLOAT(NXR-MXR*INT(FLOAT(NXR-1)
     *            /FLOAT(MXR))-1)/2.0)+1
C
              DPR1=KAPPA*SQRT(AIMAG(AN0)/OMEGA)
              MPR=INT(DPR1/DPR)
              IF(MPR.LT.1) THEN
                MPR=1
              END IF
              IPR0=INT(FLOAT(NPR-MPR*INT(FLOAT(NPR-1)
     *            /FLOAT(MPR))-1)/2.0)+1
C
              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO
C
              IF((MOD((IXR-IXR0),MXR).EQ.0).AND.
     *           (MOD((IPR-IPR0),MPR).EQ.0)) THEN
C
                FR=0.0
                FI=0.0
                PRMAX=PI/(OMEGA*DX)
C
                IF((PR.GT.(-PRMAX)).AND.(PR.LT.PRMAX)) THEN
C
                  N44=AN*CMPLX(0.0,AK0)
     *               /(OMEGA*AN-CMPLX(0.0,AK0)*PR**2)
                  N44D=ANDD*CMPLX(0.0,AK0D)
     *                /(OMEGA*ANDD-CMPLX(0.0,AK0D)*PRD**2)
C
                  AMPL=0.5*CSQRT(CMPLX(0.0,-1.0)*(AN+ANDD))
     *                *CSQRT(CMPLX(0.0,-1.0)*N44*N44D
     *                *(OMEGA/CMPLX(0.0,AK0)+OMEGA/CMPLX(0.0,AK0D)))
C
                  AMPLR=REAL(AMPL)
                  AMPLI=AIMAG(AMPL)
C
                  THETA1=0.5*ANDD/(ANDD-CMPLX(0.0,AK0D)/OMEGA*PRD**2)
                  THETA4=CMPLX(0.0,AK0D)/OMEGA
C
                  AX=GWM/SQRT(OMEGA*AIMAG(ANDD))
C
                  IXMIN=INT((XR-AX)/DX)+1
                  IF(IXMIN.LT.1) THEN
                    IXMIN=1
                  END IF
C
                  IXMAX=INT((XR+AX)/DX)+2
                  IF(IXMAX.GT.NX) THEN
                    IXMAX=NX
                  END IF
C
                  ITMIN=INT((TR-BT)/DT)+1
                  IF(ITMIN.LT.1) THEN
                    ITMIN=1
                  END IF
C
                  ITMAX=INT((TR+BT)/DT)+2
                  IF(ITMAX.GT.NT) THEN
                    ITMAX=NT
                  END IF
C
                  DO 300 IX=IXMIN,IXMAX
                    XX=OX+DX*(IX-1)
                    XRX=XX-XR
                    IGRD0=(IX-1)*NT
C
                    THETA2=ANDD*XRX**2
                    THETA3=2.0*PRD*THETA4*XRX
C
                    DO 250 IT=ITMIN,ITMAX
                      T=OT+DT*(IT-1)
                      TRT=T-TR
C
                      IGRD=IGRD0+IT
C
C                      THETA=0.5*ANDD/(ANDD-CMPLX(0.0,AK0D)
C     *                     /OMEGA*PRD**2)*(ANDD*XRX**2-2.0*PRD
C     *                     *CMPLX(0.0,AK0D)/OMEGA*XRX*TRT
C     *                     +CMPLX(0.0,AK0D)/OMEGA*TRT**2)
C
                      THETA=THETA1*(THETA2-THETA3*TRT+THETA4*TRT**2)
C
                      W=RAM(IGRD)*EXP(-OMEGA*AIMAG(THETA))
                      WW=OMEGA*(REAL(THETA)+TRT-PR*XRX)
                      WR=W*(AMPLR*COS(WW)-AMPLI*SIN(WW))
                      WI=W*(AMPLR*SIN(WW)+AMPLI*COS(WW))
C
                      FR=FR+WR
                      FI=FI+WI
C
  250               CONTINUE
C
  300             CONTINUE
C
                END IF
C
C               Complex valued amplitudes of GPs
                RAM(IGRDFR)=MXR*MPR*DXT*FR*OMEGA**2/(2.0*PI**3)
                RAM(IGRDFI)=MXR*MPR*DXT*FI*OMEGA**2/(2.0*PI**3)
C
              ELSE
C
                RAM(IGRDFR)=0.0
                RAM(IGRDFI)=0.0
C                WRITE(*,'(A,4(G16.6,1X))') '0 - ',
C     *                                      XR,TR,PR,OMEGA/(2.0*PI)
C
              END IF
C
  350       CONTINUE
C
  400     CONTINUE
C
  450   CONTINUE
C
  500 CONTINUE
C     Writing output grid values
      IF(GPOUTR.NE.' ') THEN
        CALL WARRAY(LU5,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                    NXTPO,RAM(IGRDF0+1))
      END IF
      IF(GPOUTI.NE.' ') THEN
        CALL WARRAY(LU6,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                              NXTPO,RAM(IGRDF0+NXTPO+1))
      END IF
C
C     Composition of the synthetic time section
C     -----------------------------------------
C
      WRITE(*,'(A)')'+GPANAL: Composition '
C
      BT=GWM/SQRT(AK0)
C
      DO 600 IGRD=1,NXT
        RAM(IGRD)=0.0
  600 CONTINUE
C
      DO 1000 IXR=1,NXR
        XR=OXR+DXR*(IXR-1)
C
        WRITE(*,'(A,I6)')          'IXR   =',IXR
C
        IGRDF3=(IXR-1)*NTR
C
        DO 950 ITR=1,NTR
          TR=OTR+DTR*(ITR-1)
C
          IGRDF2=(IGRDF3+ITR-1)*NPR
C
          DO 900 IPR=1,NPR
            PR=OPR+DPR*(IPR-1)
C
            IGRDF1=IGRDF0+(IGRDF2+IPR-1)*NOMEGA
C
            IFGB=((IXR-1)*NPR+IPR-1)*NTR+ITR
            IFGBY=IFGBY0+IFGB
            IFGBR=IFGBR0+IFGB
C
            AN=CMPLX(RAM(IFGBR),RAM(IFGBY))
C
            DO 850 IOMEGA=1,NOMEGA
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
C
              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO
C
              IF(((RAM(IGRDFR).GT.0.0).OR.(RAM(IGRDFI).GT.0.0)).OR.
     *        ((RAM(IGRDFR).LT.0.0).OR.(RAM(IGRDFI).LT.0.0)))THEN
C
                PRMAX=PI/(OMEGA*DX)
C
                IF((PR.GT.(-PRMAX)).AND.(PR.LT.PRMAX)) THEN
C
                  THETA1=0.5*AN/(AN-CMPLX(0.0,AK0)/OMEGA*PR**2)
                  THETA4=CMPLX(0.0,AK0)/OMEGA
C
                  AX=GWM/SQRT(OMEGA*AIMAG(AN))
C
                  IXMIN=INT((XR-AX)/DX)+1
                  IF(IXMIN.LT.1) THEN
                    IXMIN=1
                  END IF
C
                  IXMAX=INT((XR+AX)/DX)+2
                  IF(IXMAX.GT.NX) THEN
                    IXMAX=NX
                  END IF
C
                  ITMIN=INT((TR-BT)/DT)+1
                  IF(ITMIN.LT.1) THEN
                    ITMIN=1
                  END IF
C
                  ITMAX=INT((TR+BT)/DT)+2
                  IF(ITMAX.GT.NT) THEN
                    ITMAX=NT
                  END IF
C
                  DO 800 IX=IXMIN,IXMAX
                    XX=OX+DX*(IX-1)
                    XRX=XX-XR
                    IGRD0=(IX-1)*NT
C
                    THETA2=AN*XRX**2
                    THETA3=2.0*PR*THETA4*XRX
C
                    DO 750 IT=ITMIN,ITMAX
                      T=OT+DT*(IT-1)
                      TRT=T-TR
C
                      IGRD=IGRD0+IT
C
C                      THETA=0.5*AN/(AN-CMPLX(0.0,AK0)/OMEGA*PR**2)
C     *                     *(AN*XRX**2-2.0*PR*CMPLX(0.0,AK0)/OMEGA
C     *                     *XRX*TRT+CMPLX(0.0,AK0)/OMEGA*TRT**2)
C
                      THETA=THETA1*(THETA2-THETA3*TRT+THETA4*TRT**2)
C
                      WW=OMEGA*(REAL(THETA)+PR*XRX-TRT)
                      WR=RAM(IGRDFR)*COS(WW)
     *                  -RAM(IGRDFI)*SIN(WW)
                      WR=WR*EXP(-OMEGA*AIMAG(THETA))
C
C                     Synthetic time section
                      RAM(IGRD)=RAM(IGRD)+WR*DXTPO
C
  750               CONTINUE
C
  800             CONTINUE
C
                END IF
C
              END IF
C
  850       CONTINUE
C
  900     CONTINUE
C
  950   CONTINUE
C
 1000 CONTINUE
C     Writing output grid values
      IF(SYNTSEC.NE.' ') THEN
        CALL WARRAY(LU7,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                                NXT,RAM)
      END IF
C
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
      CLOSE(LU7)
      CLOSE(LU8)
      CLOSE(LU9)
      CLOSE(LU10)
      CLOSE(LU11)
      WRITE(*,'(A)') '+GPANAL: 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