C 
C Program GPSYNT to compose the synthetic time section
C from the amplitudes of 2-D Gaussian packets
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     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     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 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     GWM=real.. Gaussian beam decay at the halfwidth is EXP(-GWM*GWM/2)
C             Default: GWM=3.
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     RELERR=real... Relative root mean square error of the seismograms
C             caused by replacing the smallest amplitudes of Gaussian
C             packets by zeroes.
C             Default: RELERR=0.1
C Names of the output files:
C     GPAMPR='string'... Name of the output file containing
C             real part of Gaussian packet amplitudes.
C             Default: GPAMPR='gpampr.out'
C     GPAMPI='string'... Name of the output file containing
C             imaginary part of Gaussian packet amplitudes.
C             Default: GPAMPI='gpampi.out'
C     SSC='string'... Name of the output file with composed synthetic
C             seismograms.
C             Default: SSC='ssc.out'
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 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 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 OARRAY,RARRAY,WARRAY,ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R
C
C     Filenames and parameters:
      CHARACTER*80 FILSEP,SSC,FGBY22,FGBR22,GPAMPR,GPAMPI,GPSTEP
      CHARACTER*80 GRDNEW,GPA0R,GPA0I
      INTEGER LU1
      REAL PI
      PARAMETER (LU1=1,PI=3.141592654)
C
C     Other variables:
      INTEGER IRAM(MRAM)
      INTEGER NX,NT,NXT,NXTP,NXTPO
      INTEGER IXR,ITR,IP,IOMEGA,IX,IT,IT0
      INTEGER IGRD,IGRD0,IFGBY0,IFGBR0,IFGBY,IFGBR,IFGB
      INTEGER IGRDFR,IGRDFI,IGRDF0,IGRDF1,IGRDF2,IGRDF3
      INTEGER NXR,NTR,NP,NOMEGA
      INTEGER IXMAX,IXMIN,ITMAX,ITMIN,IXR0,IP0,MXR,MP
      REAL OX,OT,DX,DT,XX,T
      REAL DXR,DTR,DP,DOMEGA,OXR,OTR,OP,OOMEGA,XRX,TRT,AX,BT
      REAL XR,TR,P,OMEGA,FR,FI,W,WR,WI,GWM,DXTPO
      REAL DXT,PMAX,ANIMIN,Y0MIN
      REAL AK0,AK0T,AMPLR,AMPLI,AN0R,AN0I
      COMPLEX CN,CNT,CN0,CAMPL,CN44,CN44T,CPT
      REAL THETAR,THETAI
      COMPLEX THETA1,THETA4,CNT1,CPT4
      REAL THET2R,THET2I,THET3R,THET3I,THET4R,THET4I
      REAL CNT1R,CNT1I,CPT4R,CPT4I,AUX,COSWW,SINWW
      REAL RELAMP,RELERR,ERRSUM,TMIN,TMAX,WR0,WI0,WR3,WI3,WR4,WI4
      CHARACTER*80 TEXT
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPSYNT: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPANAL-01
        CALL ERROR('GPSYNT-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)')'+GPSYNT: 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('FGBY22',FGBY22,'y0m-res.out')
      CALL RSEP3T('FGBR22',FGBR22,'r0m-res.out')
      CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.out')
      CALL RSEP3T('GPAMPR',GPAMPR,'gpampr.out')
      CALL RSEP3T('GPAMPI',GPAMPI,'gpampi.out')
      CALL RSEP3T('GRDNEW',GRDNEW,'grdnew.out')
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.1)
      CALL RSEP3R('OT',OT,0.)
      CALL RSEP3R('GWM',GWM,3.)
      CALL RSEP3R('RELAMP',RELAMP,0.1)
      CALL RSEP3R('RELERR',RELERR,0.1)
      RELAMP=ALOG(RELAMP)
C
C     Reading filenames of the output files
      CALL RSEP3T('SSC',SSC,'ssc.out')
      CALL RSEP3T('GPA0R',GPA0R,'gpa0r.out')
      CALL RSEP3T('GPA0I',GPA0I,'gpa0i.out')
      CLOSE(LU1)
C
C     Reading discretization steps
      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,Y0MIN
      CLOSE(LU1)
C
      CN0=CMPLX(AN0R,AN0I)
      ANIMIN=Y0MIN
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))') ' OOMEGA=',OOMEGA
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,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
      NXT=NX*NT
      NXTPO=NXR*NTR*NP*NOMEGA
      DXTPO=DXR*DTR*DP*DOMEGA
      DXT=DX*DT
      IFGBY0=MAX0(NXT,2*NXTPO)
      IFGBR0=IFGBY0+NXTP
      IGRDF0=IFGBR0+NXTP
C
      IF((IGRDF0+NXTPO).GT.MRAM) THEN
C       GPANAL-02
        CALL ERROR('GPSYNT-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
C     Read initial parameters of Gaussian beams - imaginary part (Y0)
      CALL OARRAY(LU1,FGBY22,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXTP,RAM(IFGBY0+1))
      CLOSE(LU1)
C
C     Read initial parameters of Gaussian beams - real part (R0)
      CALL OARRAY(LU1,FGBR22,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXTP,RAM(IFGBR0+1))
      CLOSE(LU1)
C
C     Read amplitudes of 2-D Gaussian packets
      CALL OARRAY(LU1,GPAMPR,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXTPO,RAM(IGRDF0+1))
      CLOSE(LU1)
      CALL OARRAY(LU1,GPAMPI,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,
     *                                     NXTPO,RAM(IGRDF0+NXTPO+1))
      CLOSE(LU1)
      CALL OARRAY(LU1,GRDNEW,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,1,RAM(1))
      CLOSE(LU1)
C
      ERCOEF=2.*PI*PI*PI*DXTPO/FLOAT(NXT)/DXT/RAM(1)**2
C
C     Replacing the smallest amplitudes of Gaussian packets by zeroes
C
C     Loop over coordinate of intersection of central ray
C     of Gaussian packet with the profile
      WRITE(*,'(A)') '+GPSYNT: Zeroing amplitudes'
      I=0
      DO 100 IXR=1,NXR
        IGRDF3=(IXR-1)*NTR
C
C       Loop over corresponding arrival time of Gaussian packet
        DO 95 ITR=1,NTR
C         WRITE(*,'(A,I6)') '  ITR   =',ITR
          IGRDF2=(IGRDF3+ITR-1)*NP
C
C         Loop over component of slowness vector of Gaussian packet
C         along profile
          DO 90 IP=1,NP
            P=OP+DP*(IP-1)
            IGRDF1=IGRDF0+(IGRDF2+IP-1)*NOMEGA
            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           Loop over positive circular frequency of Gaussian packet
            DO 85 IOMEGA=1,NOMEGA
              I=I+1
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO
C
              AMPLR=RAM(IGRDFR)
              AMPLI=RAM(IGRDFI)
              PMAX=PI/(OMEGA*DX)
              IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN
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)
                AUX=(AMPLR**2+AMPLI**2)*ERCOEF/OMEGA**2
                RAM(I)=AUX/SQRT(AIMAG(CN)*AIMAG(CN44))
C               WRITE(*,*) AMPLR,AMPLI,ERCOEF,RAM(I)
C               WRITE(*,*) AUX,AIMAG(CN),AIMAG(CN44),RAM(I)
              ELSE
                IF(AMPLR.NE.0..OR.AMPLI.NE.0.) THEN
C               GPSYNT-03
                  CALL ERROR('GPSYNT-03: Non zero amplitude')
                END IF
                RAM(I)=0.
              END IF
   85       CONTINUE
   90     CONTINUE
   95   CONTINUE
  100 CONTINUE
C     WRITE(*,*) (RAM(I),I=1,NXTPO)
      CALL INDEXX(NXTPO,RAM(1),IRAM(NXTPO+1))
      J=0
      ERRSUM=0.
      DO 200 I=1,NXTPO
        II=IRAM(NXTPO+I)
        IF(RAM(II).EQ.0.) J=I
        ERRSUM=ERRSUM+RAM(II)
        IF(ERRSUM.GT.RELERR**2) THEN
          GO TO 201
        END IF
        RAM(IGRDF0+II)=0.
        RAM(IGRDF0+NXTPO+II)=0.
  200 CONTINUE
  201 CONTINUE
      WRITE(*,*) 'Number of zero amplitudes=',J
      WRITE(*,*) 'Number of zero amplitudes after zeroing=',I-1
      WRITE(*,*) 'Number of amplitudes=',NXTPO
C     Writing output grid values
      IF(GPA0R.NE.' ') THEN
        CALL OARRAY(LU1,GPA0R,'W')
        CALL WARRAY(LU1,' ',.FALSE.,0.,.FALSE.,0.,
     *                                    NXTPO,RAM(IGRDF0+1))
        CLOSE(LU1)
      END IF
      IF(GPA0I.NE.' ') THEN
        CALL OARRAY(LU1,GPA0I,'W')
        CALL WARRAY(LU1,' ',.FALSE.,0.,.FALSE.,0.,
     *                              NXTPO,RAM(IGRDF0+NXTPO+1))
        CLOSE(LU1)
      END IF
C
C-----------------------------------------------------------------------
C     Composition of the synthetic time section
      WRITE(*,'(A)')'+GPSYNT: Composition '
C
      BT=GWM/SQRT(AK0)
      OXR=OX
      WRITE(*,'(A,F7.2)') ' OXR=',OXR
C
      DO 600 IGRD=1,NXT
        RAM(IGRD)=0.
  600 CONTINUE
C
C     Loop over coordinate of intersection of central ray
C     of Gaussian packet with the profile
      DO 1000 IXR=1,NXR
        WRITE(*,'(2(A,I3))') '+GPSYNT: Composition NXR=',IXR,'/',NXR
        XR=OXR+DXR*(IXR-1)
        IGRDF3=(IXR-1)*NTR
C
C       Loop over corresponding arrival time of Gaussian packet
        DO 950 ITR=1,NTR
C         WRITE(*,'(A,I6)') '  ITR   =',ITR
          TR=OTR+DTR*(ITR-1)
          IGRDF2=(IGRDF3+ITR-1)*NP
C
C         Loop over component of slowness vector of Gaussian packet
C         along profile
          DO 900 IP=1,NP
            P=OP+DP*(IP-1)
            IGRDF1=IGRDF0+(IGRDF2+IP-1)*NOMEGA
            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           Loop over positive circular frequency of Gaussian packet
            DO 850 IOMEGA=1,NOMEGA
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO
C
C             IF(((RAM(IGRDFR).GT.900.).OR.(RAM(IGRDFI).GT.900.)).OR.
C    *        ((RAM(IGRDFR).LT.-900.).OR.(RAM(IGRDFI).LT.-900.)))THEN
              IF(RAM(IGRDFR).NE.0..OR.RAM(IGRDFI).NE.0.) THEN
C
                AMPLR=RAM(IGRDFR)
                AMPLI=RAM(IGRDFI)
                PMAX=PI/(OMEGA*DX)
                IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN
                  THETA1=0.5*CN/(CN-CMPLX(0.,AK0)/OMEGA*P**2)
                  THETA4=THETA1*CMPLX(0.,AK0)
                  THET4R=REAL (THETA4)
                  THET4I=AIMAG(THETA4)
                  CNT1=CN*THETA1
                  CNT1R=OMEGA*REAL (CNT1)
                  CNT1I=OMEGA*AIMAG(CNT1)
C
C                 Gaussian beam horizontal halfwidth
                  AX=GWM/SQRT(OMEGA*AIMAG(CN))
                  IXMIN=INT((XR-AX)/DX)+1
                  IF(IXMIN.LT.1) THEN
                    IXMIN=1
                  END IF
                  IXMAX=INT((XR+AX)/DX)+2
                  IF(IXMAX.GT.NX) THEN
                    IXMAX=NX
                  END IF
                  ITMIN=INT((TR-BT)/DT)+1
                  IF(ITMIN.LT.1) THEN
                    ITMIN=1
                  END IF
                  ITMAX=INT((TR+BT)/DT)+2
                  IF(ITMAX.GT.NT) THEN
                    ITMAX=NT
                  END IF
C
C                 Loop over traces of wavefield (common shot gather)
C                 to be decomposed
                  DO 800 IX=IXMIN,IXMAX
                    XX=OX+DX*(IX-1)
                    XRX=XX-XR
                    IGRD0=(IX-1)*NT
                    AUX=XRX**2
                    THET2R=AUX*CNT1R+OMEGA*P*XRX
                    THET2I=AUX*CNT1I
                    AUX=2.*P*XRX
                    THET3R=AUX*THET4R+OMEGA
                    THET3I=AUX*THET4I
                    AUX=0.5*THET3I/THET4I
                    TMAX=AUX*AUX-(RELAMP+THET2I)/THET4I
                    IF(TMAX.GT.0.) THEN
                      TMAX=SQRT(TMAX)
                      TMIN=(AUX-TMAX)
                      TMAX=(AUX+TMAX)
C                     ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),ITMIN)
C                     ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),ITMAX)
                      ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),1)
                      ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),NT)
                    ELSE
                      ITMIN=1
                      ITMAX=-1
                    END IF
C
C                   Loop over samples of each trace of wavefield
                    IF(ITMIN.LT.ITMAX) THEN
                      IT=(ITMIN+ITMAX)/2
                      T=OT+DT*(IT-1)
                      TRT=T-TR
                      THETAR=THET2R-(THET3R-THET4R*TRT)*TRT
                      THETAI=THET2I-(THET3I-THET4I*TRT)*TRT
                      W=DXTPO*EXP(-THETAI)
                      COSWW=W*COS(THETAR)
                      SINWW=W*SIN(THETAR)
                      WR=AMPLR*COSWW-AMPLI*SINWW
                      WI=AMPLR*SINWW+AMPLI*COSWW
C
C                     Synthetic time section
                      IGRD=IGRD0+IT
                      RAM(IGRD)=RAM(IGRD)+WR
                      IT0=IT
                      WR0=WR
                      WI0=WI
                      W=EXP( (THET3I-THET4I*(2.*TRT-DT))*DT)
                      AUX=  -(THET3R-THET4R*(2.*TRT-DT))*DT
                      WR3=W*COS(AUX)
                      WI3=W*SIN(AUX)
                      W=EXP(-2.*THET4I*DT*DT)
                      AUX=2.*THET4R*DT*DT
                      WR4=W*COS(AUX)
                      WI4=W*SIN(AUX)
C
                      DO 751 IGRD=IGRD0+IT0+1,IGRD0+ITMAX
                        AUX=WR3*WR4-WI3*WI4
                        WI3=WR3*WI4+WI3*WR4
                        WR3=AUX
                        AUX=WR*WR3-WI*WI3
                        WI =WR*WI3+WI*WR3
                        WR =AUX
C
C                       Synthetic time section
                        RAM(IGRD)=RAM(IGRD)+WR
  751                 CONTINUE
                      WR=WR0
                      WI=WI0
                      W=EXP(-(THET3I-THET4I*(2.*TRT+DT))*DT)
                      AUX=   (THET3R-THET4R*(2.*TRT+DT))*DT
                      WR3=W*COS(AUX)
                      WI3=W*SIN(AUX)
                      DO 752 IGRD=IGRD0+IT0-1,IGRD0+ITMIN,-1
                        AUX=WR3*WR4-WI3*WI4
                        WI3=WR3*WI4+WI3*WR4
                        WR3=AUX
                        AUX=WR*WR3-WI*WI3
                        WI =WR*WI3+WI*WR3
                        WR =AUX
C
C                       Synthetic time section
                        RAM(IGRD)=RAM(IGRD)+WR
  752                 CONTINUE
                    END IF
  800             CONTINUE
                END IF
              END IF
  850       CONTINUE
  900     CONTINUE
  950   CONTINUE
 1000 CONTINUE
C     Writing output grid values
      IF(SSC.NE.' ') THEN
        CALL OARRAY(LU1,SSC,'W')
        CALL WARRAY(LU1,' ',.FALSE.,0.,.FALSE.,0.,NXT,RAM)
        CLOSE(LU1)
      END IF
C
      WRITE(*,'(A)') '+GPSYNT: 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
      INCLUDE 'indexx.for'
C     indexx.for
C
C=======================================================================
C