C 
C Program GPANAL to calculate the amplitudes of 2-D Gaussian packets
C
C Version: 6.00
C Date: 2006, May 30
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: zacek@karel.troja.mff.cuni.cz
C
C.......................................................................
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings, read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'SEP'/
C     'SEP'...Name of the file with input parameters.
C             Description of file SEP
C             If blank, default values of the corresponding data are
C             considered.
C Default: 'SEP'=' '.
C
C                                                     
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file
C     'sep.for'.
C Names of input files:
C     TIMESEC='string'... Name of the input file containing the time
C             section.
C             Default: TIMESEC='timesec.out' (usually sufficient)
C     FGBY22='string'... Name of the input file with the initial
C             parameters of Gaussian beams - imaginary part.
C             Default: FGBY22='fgby22.out' (usually sufficient)
C     FGBR22='string'... Name of the input file with the initial
C             parameters of Gaussian beams - real part.
C             Default: FGBR22='fgbr22.out' (usually sufficient)
C Name of input parameters:
C     N44=real... ????????????????????.
C             Default: N44=?
C     FMIN=real... ????????????????????.
C             Default: FMIN=0.
C     FMAX=real... ????????????????????.
C             Default: FMAX=100.
C Name of the output file:
C     GPOUT='string'... Name of the output file containing
C             amplitudes of Gaussian packets.
C             Default: GPOUT='gpampl.out'
C Data specifying the dimensions of the input grid:
C     OT=real, OX=real, OP=real ... Coordinates of the origin
C             of the grid.
C             Default: OT=0., OX=0., OP=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     DP=real... Grid interval along the P axis.
C             Default: DP=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     NP=positive integer... Number of gridpoints along the P axis.
C             Default: NP=1
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 Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc

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

      INCLUDE 'sep.inc'
C     sep.inc

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

C     Auxiliary storage locations:
      REAL 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

      INTEGER NXR,NTR,NPR,NOMEGA
      REAL DXR,DTR,DPR,DOMEGA,OXR,OTR,OPR,OOMEGA,XRX,TRT,AX,BT

      REAL PR,OMEGA,DOMEGA1,KAPPA,PI,FR,FI,W,WW,WR,WI,GWM,DXTPO

      REAL DXT,PRMIN,PRMAX,ANIMIN,ANIMAX,ANRMIN,ANRMAX,PARK,PARN

      INTEGER IXMAX,IXMIN,ITMAX,ITMIN,IXR0,IPR0,MXR,MPR

      REAL AK0,AK0D,AMPLR,AMPLI,DXR1,DPR1

      COMPLEX AN,ANDD,AN0,THETA,AMPL,N44,N44D,PRD
      COMPLEX THETA1,THETA2,THETA3,THETA4

      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     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     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)

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     Reading filenames of the output files
      CALL RSEP3T('GPOUTR',GPOUTR,'gpampr.out')
      CALL RSEP3T('GPOUTI',GPOUTI,'gpampi.out')
      CALL RSEP3T('SYNTSEC',SYNTSEC,'synts.out')

      CALL RSEP3I('NX',NX,1)
      CALL RSEP3I('NT',NT,1)
      CALL RSEP3I('NXTP',NXTP,1)

      CALL RSEP3R('DX',DX,1.)
      CALL RSEP3R('DT',DT,1.)
      CALL RSEP3R('OX',OX,0.)
      CALL RSEP3R('OT',OT,0.)

      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.)

      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')

      NXT=NX*NT
C      NXTP=NXT*NP

      CALL RARRAY(LU2,' ','FORMATTED',.FALSE.,0.,NXT,RAM)

      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

      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      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)

      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

      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

      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

      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

      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

      NXTPO=NXR*NTR*NPR*NOMEGA
      DXTPO=DXR*DTR*DPR*DOMEGA
      DXT=DX*DT

      WRITE(*,'(A,I6)') 'Expected time is',
     *                   NINT(NXTPO/237006.0*NXT/87000.0*1.5*7200.0)

ccc   IFGBY0=NXT+1
      IFGBY0=NXT
      IFGBR0=IFGBY0+NXTP
      IGRDF0=NXT+2*NXTP

      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     Decomposition of the time section into 2-D Gaussian packets
C     -----------------------------------------------------------

      WRITE(*,'(A)')'+GPANAL: Decomposition '

      BT=GWM/SQRT(AK0D)

      DO 500 IXR=1,NXR
        XR=OXR+DXR*(IXR-1)

        WRITE(*,'(A,I6)')          'IXR   =',IXR

        IGRDF3=(IXR-1)*NTR

        DO 450 ITR=1,NTR
          TR=OTR+DTR*(ITR-1)

          IGRDF2=(IGRDF3+ITR-1)*NPR

          DO 400 IPR=1,NPR
            PR=OPR+DPR*(IPR-1)

            IGRDF1=IGRDF0+(IGRDF2+IPR-1)*NOMEGA

            IFGB=((IXR-1)*NPR+IPR-1)*NTR+ITR
            IFGBY=IFGBY0+IFGB
            IFGBR=IFGBR0+IFGB

            AN=CMPLX(RAM(IFGBR),RAM(IFGBY))
            ANDD=AN*AN0/(2.0*AN-AN0)

            PRD=-PR*ANDD/AN

            DO 350 IOMEGA=1,NOMEGA
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)

              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

              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

              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO

              IF((MOD((IXR-IXR0),MXR).EQ.0).AND.
     *           (MOD((IPR-IPR0),MPR).EQ.0)) THEN

                FR=0.0
                FI=0.0
                PRMAX=PI/(OMEGA*DX)

                IF((PR.GT.(-PRMAX)).AND.(PR.LT.PRMAX)) THEN

                  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)

                  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)))

                  AMPLR=REAL(AMPL)
                  AMPLI=AIMAG(AMPL)

                  THETA1=0.5*ANDD/(ANDD-CMPLX(0.0,AK0D)/OMEGA*PRD**2)
                  THETA4=CMPLX(0.0,AK0D)/OMEGA

                  AX=GWM/SQRT(OMEGA*AIMAG(ANDD))

                  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

                  DO 300 IX=IXMIN,IXMAX
                    XX=OX+DX*(IX-1)
                    XRX=XX-XR
                    IGRD0=(IX-1)*NT

                    THETA2=ANDD*XRX**2
                    THETA3=2.0*PRD*THETA4*XRX

                    DO 250 IT=ITMIN,ITMAX
                      T=OT+DT*(IT-1)
                      TRT=T-TR

                      IGRD=IGRD0+IT

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)

                      THETA=THETA1*(THETA2-THETA3*TRT+THETA4*TRT**2)

                      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))

                      FR=FR+WR
                      FI=FI+WI

  250               CONTINUE

  300             CONTINUE

                END IF

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)

              ELSE

                RAM(IGRDFR)=0.0
                RAM(IGRDFI)=0.0
C                WRITE(*,'(A,4(G16.6,1X))') '0 - ',
C     *                                      XR,TR,PR,OMEGA/(2.0*PI)

              END IF

  350       CONTINUE

  400     CONTINUE

  450   CONTINUE

  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     -----------------------------------------

      WRITE(*,'(A)')'+GPANAL: Composition '

      BT=GWM/SQRT(AK0)

      DO 600 IGRD=1,NXT
        RAM(IGRD)=0.0
  600 CONTINUE

      DO 1000 IXR=1,NXR
        XR=OXR+DXR*(IXR-1)

        WRITE(*,'(A,I6)')          'IXR   =',IXR

        IGRDF3=(IXR-1)*NTR

        DO 950 ITR=1,NTR
          TR=OTR+DTR*(ITR-1)

          IGRDF2=(IGRDF3+ITR-1)*NPR

          DO 900 IPR=1,NPR
            PR=OPR+DPR*(IPR-1)

            IGRDF1=IGRDF0+(IGRDF2+IPR-1)*NOMEGA

            IFGB=((IXR-1)*NPR+IPR-1)*NTR+ITR
            IFGBY=IFGBY0+IFGB
            IFGBR=IFGBR0+IFGB

            AN=CMPLX(RAM(IFGBR),RAM(IFGBY))

            DO 850 IOMEGA=1,NOMEGA
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)

              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO

              IF(((RAM(IGRDFR).GT.900.0).OR.(RAM(IGRDFI).GT.900.0)).OR.
     *        ((RAM(IGRDFR).LT.-900.0).OR.(RAM(IGRDFI).LT.-900.0)))THEN

                PRMAX=PI/(OMEGA*DX)

                IF((PR.GT.(-PRMAX)).AND.(PR.LT.PRMAX)) THEN

                  THETA1=0.5*AN/(AN-CMPLX(0.0,AK0)/OMEGA*PR**2)
                  THETA4=CMPLX(0.0,AK0)/OMEGA

                  AX=GWM/SQRT(OMEGA*AIMAG(AN))

                  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

                  DO 800 IX=IXMIN,IXMAX
                    XX=OX+DX*(IX-1)
                    XRX=XX-XR
                    IGRD0=(IX-1)*NT

                    THETA2=AN*XRX**2
                    THETA3=2.0*PR*THETA4*XRX

                    DO 750 IT=ITMIN,ITMAX
                      T=OT+DT*(IT-1)
                      TRT=T-TR

                      IGRD=IGRD0+IT

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)

                      THETA=THETA1*(THETA2-THETA3*TRT+THETA4*TRT**2)

                      WW=OMEGA*(REAL(THETA)+PR*XRX-TRT)
                      WR=RAM(IGRDFR)*COS(WW)
     *                  -RAM(IGRDFI)*SIN(WW)
                      WR=WR*EXP(-OMEGA*AIMAG(THETA))

C                     Synthetic time section
                      RAM(IGRD)=RAM(IGRD)+WR*DXTPO

  750               CONTINUE

  800             CONTINUE

                END IF

              END IF

  850       CONTINUE

  900     CONTINUE

  950   CONTINUE

 1000 CONTINUE
C     Writing output grid values
      IF(SYNTSEC.NE.' ') THEN
        CALL WARRAY(LU7,' ','FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                                                NXT,RAM)
      END IF

      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
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.forC
C
C=======================================================================
C