C 
C Program GRDFFT to compute the 1-D, 2-D or 3-D Fourier transformation
C of a complex function defined on 1-D, 2-D or 3-D grid of points.
C (The dimension of the FFT is less or equal to the grid dimension.)
C If the number of gridpoints in any direction of the input grid is not
C a power of 2, the input grid is enlarged to the nearest power of 2
C and the functional values in new gridpoints are completed according
C to input parameter FFTFIL.
C Subroutines from the Numerical Recipes are then used
C for the entire FFT.
C
C Version: 5.80
C Date: 2004, June 11
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Parameter describing the form of the FFT:
C     FFT=real ... The Fourier transform F(y) has the form of integral
C             of f(x)exp(i*FFT*x*y)dx, where f(x) is the input function
C             to be transformed. The integral is then multiplied by the
C             factor of (ABS(FFT)/(2.*PI))**(NDIM/2), where NDIM is the
C             dimension of the part of the input grid to be transformed
C             (and of the transformed grid).
C             The mostly used values of FFT are 6.28, -6.28, 1, -1.
C             If FFT is input as a multiple of PI plus minus 1%, the
C             value of FFT is rounded to the multiple of PI.
C             Default: FFT=6.28 (which means 2.*PI)
C Data specifying the parameters of the input grid:
C     O1=real, O2=real, O3=real ... Coordinates of the origin of the
C             grid.
C             Default: O1=0.  O2=0.  O3=0.
C     D1=real... Grid interval along the X1 axis.
C             Default: D1=0.
C     D2=real... Grid interval along the X2 axis.
C             Default: D2=0.
C     D3=real... Grid interval along the X3 axis.
C             Default: D3=0.
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     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     N4=positive integer... Number of space grids.  The Fourier
C             transform with respect to X1, X2 and X3 will be repeated
C             N4 times, for each space grid.
C             If N4=0, parameter M4 is used to specify the number of
C             space grids.
C             Note that each space grid must begin at a new line of
C             the input file, because each space grid is read by a
C             a single read statement in free format.
C             Default: N4=1
C     M4='string'... Name of the file containing a single integer number
C             specifying N4.
C             Default: M4=' ' means that N4=1.
C Data specifying the parameters of the grid for the FFT:
C     N1FFT=positive integer... Number of gridpoints along the X1 axis.
C     N2FFT=positive integer... Number of gridpoints along the X2 axis.
C     N3FFT=positive integer... Number of gridpoints along the X3 axis.
C             NiFFT must be an integer power of 2, and must be greater
C             than or equal to the corresponding Ni of the input data.
C             NiFFT=0 means, that the Fourier transform is not to be
C             done in the direction of the corresponding axis. In this
C             case the corresponding Ni number of gridpoints is
C             considered, this influences the default value of NiOUT.
C             Default: NiFFT equal to the lowest power of 2, which is
C             greater or equal to the corresponding Ni.
C     FFTFIL=real ... Value, which is used at the gridpoints of the grid
C             for FFT, at which the value is not given by input data.
C             In present version FFTFIL is used for real part, imaginary
C             part is zero, with exception of FFTFIL=-999999999. In such
C             case, input data are linearly interpolated from the values
C             in the first and in the last gridpoints of the input grid,
C             to the gridpoints of the grid for the FFT.
C             Default: FFTFIL=-999999999. (interpolation of the values)
C Data specifying the parameters of the output grid:
C     O1OUT=real, O2OUT=real, O3OUT=real ... Coordinates of the origin
C             of the output grid.
C             Default: OiOUT=-1./(2.*Di)*2.*PI/ABS(FFT)
C     D1OUT=real, D2OUT=real, D3OUT=real ... Grid intervals along the
C             first, second and third axes of the output grid. DiOUT
C             must not differ from the default value.
C             Default: DiOUT=1./(NiFFT*Di)*2.*PI/ABS(FFT)
C     N1OUT=positive integer, N2OUT=positive integer,
C     N3OUT=positive integer ... Numbers of gridpoints along the first,
C             second and third axes of the output grid.
C             Default: NiOUT=NiFFT
C Names of input and output formatted files with the functional values:
C     FFTINR='string', FFTINI='string' ... real and imaginary parts of
C             the input function.
C             Default: FFTINR=' ' FFTINI=' '
C     FFTOUTR='string', FFTOUTI='string' ... real and imaginary parts of
C             the output function.
C             Default: FFTOUTR=' ' FFTOUTI=' '
C     For general description of the files with gridded data refer to
C     file forms.htm.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL NFFT,INDRAM,MODF,NCHECK,ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,
     *RARRAY,WARRAY,FOURN
      INTEGER NFFT,INDRAM,MODF
C     NFFT,INDRAM,MODF,NCHECK ... This file.
C     ERROR ... File error.for.
C     RSEP1,RSEP3T,RSEP3R,RSEP3I ...
C             File sep.for.
C     RARRAY,WARRAY ... File forms.for.
C     FOURN ... File fourn.for.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C Common block /GRDFF/:
      INTEGER N1FFT,N1N2FF
      COMMON /GRDFF/ N1N2FF,N1FFT
      SAVE   /GRDFF/
C
      CHARACTER*80 FILSEP,FM4,FILINR,FILINI,FILOUR,FILOUI
      INTEGER LU1,LU2,LU3,LU4
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4)
      REAL PI
      PARAMETER (PI=3.141592653589793)
      INTEGER MODFFT
      INTEGER N1,N2,N3,N4,N1N2,NN
      INTEGER N2FFT,N3FFT,NDIFFT(3),NNFFT,NN2FFT,NT1FFT,NT2FFT,NT3FFT
      INTEGER N1TMP,N2TMP,N3TMP
      INTEGER N1OUT,N2OUT,N3OUT,N1N2OU,NNOUT
      REAL UNDEF
      PARAMETER (UNDEF=-999999999.)
      REAL O1,O2,O3,D1,D2,D3,FFT,FFTFIL
      REAL O1OUT,O2OUT,O3OUT,D1OUT,D2OUT,D3OUT,D1TMP,D2TMP,D3TMP
      INTEGER IR,II,I1MI,I1MA,I2MI,I2MA,I3MI,I3MA,IO1,IO2,IO3
      INTEGER K1MA,K2MA,K3MA,K1,K2,K3
      INTEGER IRAM,I1,I2,I3,I4,I,J,K,L,NDIMFF,NFORFF,OFORFF
      REAL RRA,RRB,RR0,RRD,RIA,RIB,RI0,RID,RRK,RMULT
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      FILSEP=' '
      WRITE(*,'(A)') '+GRDFFT: Enter input filename: '
      READ(*,*) FILSEP
      IF (FILSEP.EQ.' ') THEN
C       GRDFFT-19
        CALL ERROR('GRDFFT-19: 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)') '+GRDFFT: Working...            '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading the filenames of the files with the real and imaginary
C     parts of the input function:
      CALL RSEP3T('FFTINR',FILINR,' ')
      CALL RSEP3T('FFTINI',FILINI,' ')
      IF (FILINR.EQ.' ' .AND. FILINI.EQ.' ') THEN
C       GRDFFT-01
        CALL ERROR('GRDFFT-01: No input files specified.')
      ENDIF
C     Reading the filenames of the output files:
      CALL RSEP3T('FFTOUTR',FILOUR,' ')
      CALL RSEP3T('FFTOUTI',FILOUI,' ')
      IF (FILOUR.EQ.' ' .AND. FILOUI.EQ.' ') THEN
C       GRDFFT-02
        CALL ERROR('GRDFFT-02: No output files specified.')
      ENDIF
C     Reading the multiplicative constant FFT:
      CALL RSEP3R('FFT',FFT,2.*PI)
      I=NINT(FFT/PI)
      IF (I.NE.0) THEN
        IF (ABS((FFT/PI-FLOAT(I))/FLOAT(I)).LE.0.01) THEN
          FFT=FLOAT(I)*PI
        ENDIF
      ENDIF
C     Mode of the FFT:
      IF (FFT.GT.0.) THEN
        MODFFT=1
      ELSEIF (FFT.LT.0.) THEN
        MODFFT=-1
      ELSE
C       GRDFFT-03
        CALL ERROR('GRDFFT-03: Wrong value of FFT.')
      ENDIF
C     Reading the values describing the input grid:
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('N4',N4,1)
      N1N2=N1*N2
      NN=N1N2*N3
      IF (N4.LE.0) THEN
        CALL RSEP3T('M4',FM4,' ')
        IF (FM4.EQ.' ') THEN
          N4=1
        ELSE
          OPEN(LU1,FILE=FM4,STATUS='OLD')
          READ(LU1,*) N4
          CLOSE(LU1)
        ENDIF
      ENDIF
      CALL RSEP3R('D1',D1,0.)
      CALL RSEP3R('D2',D2,0.)
      CALL RSEP3R('D3',D3,0.)
      IF (((D1.EQ.0.).AND.(N1.NE.1)).OR.
     *    ((D2.EQ.0.).AND.(N2.NE.1)).OR.
     *    ((D3.EQ.0.).AND.(N3.NE.1))) THEN
C       GRDFFT-04
        CALL ERROR('GRDFFT-04: Wrong input grid.')
      ENDIF
      IF ((D1.EQ.0.).AND.(N1.EQ.1)) D1=1.
      IF ((D2.EQ.0.).AND.(N2.EQ.1)) D2=1.
      IF ((D3.EQ.0.).AND.(N3.EQ.1)) D3=1.
      N1TMP=NFFT(N1)
      N2TMP=NFFT(N2)
      N3TMP=NFFT(N3)
      CALL RSEP3I('N1FFT',NT1FFT,N1TMP)
      CALL RSEP3I('N2FFT',NT2FFT,N2TMP)
      CALL RSEP3I('N3FFT',NT3FFT,N3TMP)
      IF (NT1FFT.EQ.0) THEN
        N1FFT=N1
      ELSE
        N1FFT=NT1FFT
        CALL NCHECK(N1FFT,N1TMP)
      ENDIF
      IF (NT2FFT.EQ.0) THEN
        N2FFT=N2
      ELSE
        N2FFT=NT2FFT
        CALL NCHECK(N2FFT,N2TMP)
      ENDIF
      IF (NT3FFT.EQ.0) THEN
        N3FFT=N3
      ELSE
        N3FFT=NT3FFT
        CALL NCHECK(N3FFT,N3TMP)
      ENDIF
      N1N2FF=N1FFT*N2FFT
      NNFFT=N1N2FF*N3FFT
      NN2FFT=2*NNFFT
C     Reading the constant FFTFIL for values missing in input grid:
      CALL RSEP3R('FFTFIL',FFTFIL,UNDEF)
C     Dimension of the temporary grid:
      NFORFF=NNFFT
      IF (NT1FFT.EQ.0) NFORFF=NFORFF/N1FFT
      IF (NT2FFT.EQ.0) NFORFF=NFORFF/N2FFT
      IF (NT3FFT.EQ.0) NFORFF=NFORFF/N3FFT
      OFORFF=MRAM-2*NFORFF
C     Reading the values describing the output grid:
      CALL RSEP3I('N1OUT',N1OUT,N1FFT)
      CALL RSEP3I('N2OUT',N2OUT,N2FFT)
      CALL RSEP3I('N3OUT',N3OUT,N3FFT)
      N1N2OU=N1OUT*N2OUT
      NNOUT=N1N2OU*N3OUT
      CALL RSEP3R('O1OUT',O1OUT,-1./(2.*D1)*2.*PI/ABS(FFT))
      CALL RSEP3R('O2OUT',O2OUT,-1./(2.*D2)*2.*PI/ABS(FFT))
      CALL RSEP3R('O3OUT',O3OUT,-1./(2.*D3)*2.*PI/ABS(FFT))
      D1TMP=1./(N1FFT*D1)*2.*PI/ABS(FFT)
      D2TMP=1./(N2FFT*D2)*2.*PI/ABS(FFT)
      D3TMP=1./(N3FFT*D3)*2.*PI/ABS(FFT)
      CALL RSEP3R('D1OUT',D1OUT,D1TMP)
      CALL RSEP3R('D2OUT',D2OUT,D2TMP)
      CALL RSEP3R('D3OUT',D3OUT,D3TMP)
      D1TMP=D1OUT/D1TMP
      D2TMP=D2OUT/D2TMP
      D3TMP=D3OUT/D3TMP
      IF ((ABS(D1TMP-1.).GT.0.001).AND.(N1OUT.NE.1)) THEN
        IF ((NT1FFT.NE.0).OR.(D1OUT.NE.D1)) THEN
C         GRDFFT-05
          CALL ERROR('GRDFFT-05: Wrong D1OUT.')
        ENDIF
      ENDIF
      IF ((ABS(D2TMP-1.).GT.0.001).AND.(N2OUT.NE.1)) THEN
        IF ((NT2FFT.NE.0).OR.(D2OUT.NE.D2)) THEN
C         GRDFFT-16
          CALL ERROR('GRDFFT-16: Wrong D2OUT.')
        ENDIF
      ENDIF
      IF ((ABS(D3TMP-1.).GT.0.001).AND.(N3OUT.NE.1)) THEN
        IF ((NT3FFT.NE.0).OR.(D3OUT.NE.D3)) THEN
C         GRDFFT-17
          CALL ERROR('GRDFFT-17: Wrong D3OUT.')
        ENDIF
      ENDIF
C
      IF ((NN2FFT+2*MAX0(NN,NNOUT,NFORFF)).GT.MRAM) THEN
C       GRDFFT-06
        CALL ERROR('GRDFFT-06: Small array RAM.')
      ENDIF
C
C     Preparing number NDIMFF describing the dimension
C     of the part of the input grid to be transformed.
      NDIMFF=3
      IF ((N1.EQ.1).OR.(NT1FFT.EQ.0))  NDIMFF=NDIMFF-1
      IF ((N2.EQ.1).OR.(NT2FFT.EQ.0))  NDIMFF=NDIMFF-1
      IF ((N3.EQ.1).OR.(NT3FFT.EQ.0))  NDIMFF=NDIMFF-1
      IF (NDIMFF.EQ.0) THEN
C       GRDFFT-07
        CALL ERROR('GRDFFT-07: Input grid is 1*1*1.')
      ENDIF
C     Computing the multiplicative factor:
      RMULT=1.
      IF ((N1.NE.1).AND.(NT1FFT.NE.0)) RMULT=RMULT*D1
      IF ((N2.NE.1).AND.(NT2FFT.NE.0)) RMULT=RMULT*D2
      IF ((N3.NE.1).AND.(NT3FFT.NE.0)) RMULT=RMULT*D3
      RMULT=RMULT*SQRT(ABS(FFT)/(2.*PI))**NDIMFF
C
C     Opening files with input and output grids:
      IF (N4.GT.1) THEN
        IF (FILINR.NE.' ') THEN
          OPEN(LU1,FILE=FILINR,FORM='FORMATTED',STATUS='OLD')
        ENDIF
        IF (FILINI.NE.' ') THEN
          OPEN(LU2,FILE=FILINI,FORM='FORMATTED',STATUS='OLD')
        ENDIF
        IF (FILOUR.NE.' ') THEN
          OPEN(LU3,FILE=FILOUR,FORM='FORMATTED')
        ENDIF
        IF (FILOUI.NE.' ') THEN
          OPEN(LU4,FILE=FILOUI,FORM='FORMATTED')
        ENDIF
      ENDIF
C
C     Loop over N4 space grids:
      DO 90, I4=1,N4
C       Reading the input function:
        IR=MRAM-2*NN
        II=MRAM-NN
        IF (FILINR.NE.' ') THEN
          IF (N4.GT.1) THEN
            CALL RARRAY(LU1,' '   ,'FORMATTED',.TRUE.,0.,NN,RAM(IR))
          ELSE
            CALL RARRAY(LU1,FILINR,'FORMATTED',.TRUE.,0.,NN,RAM(IR))
          ENDIF
        ELSE
          DO 10, I1=IR,II-1
            RAM(I1)=0.
  10      CONTINUE
        ENDIF
        IF (FILINI.NE.' ') THEN
          IF (N4.GT.1) THEN
            CALL RARRAY(LU2,' '   ,'FORMATTED',.TRUE.,0.,NN,RAM(II))
          ELSE
            CALL RARRAY(LU2,FILINI,'FORMATTED',.TRUE.,0.,NN,RAM(II))
          ENDIF
        ELSE
          DO 11, I1=II,MRAM
            RAM(I1)=0.
  11      CONTINUE
        ENDIF
C
        I=N1FFT-N1
        I1MI=1+I/2
        I1MA=N1FFT-I/2-MOD(I,2)
        I=N2FFT-N2
        I2MI=1+I/2
        I2MA=N2FFT-I/2-MOD(I,2)
        I=N3FFT-N3
        I3MI=1+I/2
        I3MA=N3FFT-I/2-MOD(I,2)
        IF ((I1MI.LT.1).OR.(I2MI.LT.1).OR.(I3MI.LT.1).OR.
     *      (I1MI.GT.N1FFT).OR.(I2MI.GT.N2FFT).OR.(I3MI.GT.N3FFT).OR.
     *      (I1MA.LT.1).OR.(I2MA.LT.1).OR.(I3MA.LT.1).OR.
     *      (I1MA.GT.N1FFT).OR.(I2MA.GT.N2FFT).OR.(I3MA.GT.N3FFT)) THEN
C         GRDFFT-08
          CALL ERROR('GRDFFT-08: Wrong value of IiMA or IiMI.')
C         This error should not appear.
        ENDIF
        IO1=MODF(-NINT(O1/D1),N1)
        IF (IO1.LT.0) IO1=IO1+N1
        IO2=MODF(-NINT(O2/D2),N2)
        IF (IO2.LT.0) IO2=IO2+N2
        IO3=MODF(-NINT(O3/D3),N3)
        IF (IO3.LT.0) IO3=IO3+N3
        IF ((IO1.LT.0).OR.(IO2.LT.0).OR.(IO3.LT.0).OR.
     *      (IO1.GT.N1).OR.(IO2.GT.N2).OR.(IO3.GT.N3).OR.
     *      (IO1.LT.0).OR.(IO2.LT.0).OR.(IO3.LT.0).OR.
     *      (IO1.GT.N1).OR.(IO2.GT.N2).OR.(IO3.GT.N3)) THEN
C         GRDFFT-09
          CALL ERROR('GRDFFT-09: Wrong value of IOi.')
C         This error should not appear.
        ENDIF
C       Recording the known values:
        IRAM=1
        DO 24, I3=1,N3FFT
          DO 23, I2=1,N2FFT
            DO 22, I1=1,N1FFT
              IF (IRAM.GT.IR) THEN
C               GRDFFT-10
                CALL ERROR('GRDFFT-10: Wrong index of array RAM.')
C               This error should not appear.
              ENDIF
              IF ((I1.GE.I1MI).AND.(I1.LE.I1MA).AND.
     *            (I2.GE.I2MI).AND.(I2.LE.I2MA).AND.
     *            (I3.GE.I3MI).AND.(I3.LE.I3MA)) THEN
                I=MODF(IO1+I1,N1) + (MODF(IO2+I2,N2)-1)*N1 +
     *            (MODF(IO3+I3,N3)-1)*N1N2 + IR-1
                IF ((I.LT.IR).OR.(I.GE.II)) THEN
C                 
C                 GRDFFT-11
                  CALL ERROR('GRDFFT-11: Wrong index in the data array')
C                 This error should not appear.
                ENDIF
                RAM(IRAM)=RAM(I)
                RAM(IRAM+1)=RAM(I+NN)
              ENDIF
              IRAM=IRAM+2
  22        CONTINUE
  23      CONTINUE
  24    CONTINUE
C       Interpolating the unknown values in the first axis direction:
        DO 34, I3=I3MI,I3MA
          DO 33, I2=I2MI,I2MA
            RRA=RAM(INDRAM(I1MI,I2,I3))
            RRB=RAM(INDRAM(I1MA,I2,I3))
            RR0=(RRA+RRB)/2.
            RRD=(RRA-RRB)/2.
            RIA=RAM(INDRAM(I1MI,I2,I3)+1)
            RIB=RAM(INDRAM(I1MA,I2,I3)+1)
            RI0=(RIA+RIB)/2.
            RID=(RIA-RIB)/2.
            DO 31, I1=1,I1MI-1
              IF (FFTFIL.EQ.UNDEF) THEN
                RRK=FLOAT(I1-1)/FLOAT(I1MI-1)
                RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD
                RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID
              ELSE
                RAM(INDRAM(I1,I2,I3))=FFTFIL
                RAM(INDRAM(I1,I2,I3)+1)=0.
              ENDIF
  31        CONTINUE
            DO 32, I1=I1MA+1,N1FFT
              IF (FFTFIL.EQ.UNDEF) THEN
                RRK=FLOAT(I1-N1FFT)/FLOAT(N1FFT-I1MA)
                RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD
                RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID
              ELSE
                RAM(INDRAM(I1,I2,I3))=FFTFIL
                RAM(INDRAM(I1,I2,I3)+1)=0.
              ENDIF
  32        CONTINUE
  33      CONTINUE
  34    CONTINUE
C       Interpolating the unknown values in the second axis direction:
        DO 44, I3=I3MI,I3MA
          DO 43, I1=1,N1FFT
            RRA=RAM(INDRAM(I1,I2MI,I3))
            RRB=RAM(INDRAM(I1,I2MA,I3))
            RR0=(RRA+RRB)/2.
            RRD=(RRA-RRB)/2.
            RIA=RAM(INDRAM(I1,I2MI,I3)+1)
            RIB=RAM(INDRAM(I1,I2MA,I3)+1)
            RI0=(RIA+RIB)/2.
            RID=(RIA-RIB)/2.
            DO 41, I2=1,I2MI-1
              IF (FFTFIL.EQ.UNDEF) THEN
                RRK=FLOAT(I2-1)/FLOAT(I2MI-1)
                RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD
                RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID
              ELSE
                RAM(INDRAM(I1,I2,I3))=FFTFIL
                RAM(INDRAM(I1,I2,I3)+1)=0.
              ENDIF
  41        CONTINUE
            DO 42, I2=I2MA+1,N2FFT
              IF (FFTFIL.EQ.UNDEF) THEN
                RRK=FLOAT(I2-N2FFT)/FLOAT(N2FFT-I2MA)
                RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD
                RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID
              ELSE
                RAM(INDRAM(I1,I2,I3))=FFTFIL
                RAM(INDRAM(I1,I2,I3)+1)=0.
              ENDIF
  42        CONTINUE
  43      CONTINUE
  44    CONTINUE
C       Interpolating the unknown values in the third axis direction:
        DO 54, I1=1,N1FFT
          DO 53, I2=1,N2FFT
            RRA=RAM(INDRAM(I1,I2,I3MI))
            RRB=RAM(INDRAM(I1,I2,I3MA))
            RR0=(RRA+RRB)/2.
            RRD=(RRA-RRB)/2.
            RIA=RAM(INDRAM(I1,I2,I3MI)+1)
            RIB=RAM(INDRAM(I1,I2,I3MA)+1)
            RI0=(RIA+RIB)/2.
            RID=(RIA-RIB)/2.
            DO 51, I3=1,I3MI-1
              IF (FFTFIL.EQ.UNDEF) THEN
                RRK=FLOAT(I3-1)/FLOAT(I3MI-1)
                RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD
                RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID
              ELSE
                RAM(INDRAM(I1,I2,I3))=FFTFIL
                RAM(INDRAM(I1,I2,I3)+1)=0.
              ENDIF
  51        CONTINUE
            DO 52, I3=I3MA+1,N3FFT
              IF (FFTFIL.EQ.UNDEF) THEN
                RRK=FLOAT(I3-N3FFT)/FLOAT(N3FFT-I3MA)
                RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD
                RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID
              ELSE
                RAM(INDRAM(I1,I2,I3))=FFTFIL
                RAM(INDRAM(I1,I2,I3)+1)=0.
              ENDIF
  52        CONTINUE
  53      CONTINUE
  54    CONTINUE
C
C
C       Computing the FFT:
C       Quantities describing the parts of the grid where FFT will
C       be done:
        NDIFFT(1)=N1FFT
        NDIFFT(2)=N2FFT
        NDIFFT(3)=N3FFT
        IF (NT1FFT.EQ.0) NDIFFT(1)=1
        IF (NT2FFT.EQ.0) NDIFFT(2)=1
        IF (NT3FFT.EQ.0) NDIFFT(3)=1
        IF (NDIFFT(2).EQ.1) THEN
          NDIFFT(2)=NDIFFT(3)
        ENDIF
        IF (NDIFFT(1).EQ.1) THEN
          NDIFFT(1)=NDIFFT(2)
          NDIFFT(2)=NDIFFT(3)
        ENDIF
        NFORFF=NNFFT
        IF (NT1FFT.EQ.0) NFORFF=NFORFF/N1FFT
        IF (NT2FFT.EQ.0) NFORFF=NFORFF/N2FFT
        IF (NT3FFT.EQ.0) NFORFF=NFORFF/N3FFT
        OFORFF=MRAM-2*NFORFF
        K1MA=1
        K2MA=1
        K3MA=1
        IF (NT1FFT.EQ.0) K1MA=N1FFT
        IF (NT2FFT.EQ.0) K2MA=N2FFT
        IF (NT3FFT.EQ.0) K3MA=N3FFT
        I1MI=1
        I1MA=N1FFT
        I2MI=1
        I2MA=N2FFT
        I3MI=1
        I3MA=N3FFT
C
C       Loop over subgrids:
        DO 69, K3=1,K3MA
        DO 68, K2=1,K2MA
        DO 67, K1=1,K1MA
          IF (NT1FFT.EQ.0) THEN
            I1MI=K1
            I1MA=K1
          ENDIF
          IF (NT2FFT.EQ.0) THEN
            I2MI=K2
            I2MA=K2
          ENDIF
          IF (NT3FFT.EQ.0) THEN
            I3MI=K3
            I3MA=K3
          ENDIF
C         Moving the subgrid to the temporary location:
          K=OFORFF-2
          DO 63, I3=I3MI,I3MA
          DO 62, I2=I2MI,I2MA
          DO 61, I1=I1MI,I1MA
            I=INDRAM(I1,I2,I3)
            J=I+1
            K=K+2
            L=K+1
            IF ((I.LE.0).OR.(I.GT.NN2FFT).OR.
     *          (J.LE.0).OR.(J.GT.NN2FFT).OR.
     *          (K.GT.MRAM).OR.(L.GT.MRAM)) THEN
C             GRDFFT-15
              CALL ERROR('GRDFFT-15: Wrong index of array RAM.')
C             This error should not appear.
            ENDIF
            RAM(K)=RAM(I)
            RAM(L)=RAM(J)
  61      CONTINUE
  62      CONTINUE
  63      CONTINUE
C         FFT:
          CALL FOURN(RAM(OFORFF),NDIFFT,NDIMFF,MODFFT)
C         Moving the subgrid back from the temporary location:
          K=OFORFF-2
          DO 66, I3=I3MI,I3MA
          DO 65, I2=I2MI,I2MA
          DO 64, I1=I1MI,I1MA
            I=INDRAM(I1,I2,I3)
            J=I+1
            K=K+2
            L=K+1
            IF ((I.LE.0).OR.(I.GT.NN2FFT).OR.
     *          (J.LE.0).OR.(J.GT.NN2FFT).OR.
     *          (K.GT.MRAM).OR.(L.GT.MRAM)) THEN
C             GRDFFT-18
              CALL ERROR('GRDFFT-18: Wrong index of array RAM.')
C             This error should not appear.
            ENDIF
            RAM(I)=RAM(K)
            RAM(J)=RAM(L)
  64      CONTINUE
  65      CONTINUE
  66      CONTINUE
  67    CONTINUE
  68    CONTINUE
  69    CONTINUE
C       End of the loop over subgrids.
C
C
C       Adding the multiplicative factor:
        DO 70, I1=1,NN2FFT
          RAM(I1)=RAM(I1)*RMULT
  70    CONTINUE
C
C
C       Writing the results of the FFT:
        IO1=MODF(NINT(O1OUT/D1OUT),N1FFT)
        IF (IO1.LT.0) IO1=IO1+N1FFT
        IO2=MODF(NINT(O2OUT/D2OUT),N2FFT)
        IF (IO2.LT.0) IO2=IO2+N2FFT
        IO3=MODF(NINT(O3OUT/D3OUT),N3FFT)
        IF (IO3.LT.0) IO3=IO3+N3FFT
        IF ((IO1.LT.0).OR.(IO2.LT.0).OR.(IO3.LT.0).OR.
     *      (IO1.GT.N1FFT).OR.(IO2.GT.N2FFT).OR.(IO3.GT.N3FFT).OR.
     *      (IO1.LT.0).OR.(IO2.LT.0).OR.(IO3.LT.0).OR.
     *      (IO1.GT.N1FFT).OR.(IO2.GT.N2FFT).OR.(IO3.GT.N3FFT)) THEN
C         GRDFFT-12
          CALL ERROR('GRDFFT-12: Wrong value of IOi.')
C         This error should not appear.
        ENDIF
C       Reordering the computed values:
        IR=MRAM-2*NNOUT
        II=MRAM-NNOUT
        DO 74, I3=1,N3OUT
          DO 73, I2=1,N2OUT
            DO 72, I1=1,N1OUT
              I=INDRAM(MODF(IO1+I1,N1FFT),MODF(IO2+I2,N2FFT),
     *                 MODF(IO3+I3,N3FFT))
              RAM(IR)=RAM(I)
              RAM(II)=RAM(I+1)
              IR=IR+1
              II=II+1
  72        CONTINUE
  73      CONTINUE
  74    CONTINUE
        IR=MRAM-2*NNOUT
        II=MRAM-NNOUT
        IF (N4.GT.1) THEN
          IF (FILOUR.NE.' ') THEN
            CALL WARRAY(LU3,' '   ,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                  NNOUT,RAM(IR))
          ENDIF
          IF (FILOUI.NE.' ') THEN
            CALL WARRAY(LU4,' '   ,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                  NNOUT,RAM(II))
          ENDIF
        ELSE
          IF (FILOUR.NE.' ') THEN
            CALL WARRAY(LU3,FILOUR,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                  NNOUT,RAM(IR))
          ENDIF
          IF (FILOUI.NE.' ') THEN
            CALL WARRAY(LU4,FILOUI,'FORMATTED',.FALSE.,0.,.FALSE.,0.,
     *                  NNOUT,RAM(II))
          ENDIF
        ENDIF
  90  CONTINUE
C
C     Closing files with input and output grids:
      IF (N4.GT.1) THEN
        IF (FILINR.NE.' ') THEN
          CLOSE(LU1)
        ENDIF
        IF (FILINI.NE.' ') THEN
          CLOSE(LU2)
        ENDIF
        IF (FILOUR.NE.' ') THEN
          CLOSE(LU3)
        ENDIF
        IF (FILOUI.NE.' ') THEN
          CLOSE(LU4)
        ENDIF
      ENDIF
      WRITE(*,'(A)') '+GRDFFT: Finished.             '
      STOP
      END
C
C=======================================================================
C
      INTEGER FUNCTION INDRAM(I,J,K)
C Common block /GRDFF/:
      INTEGER N1FFT,N1N2FF
      COMMON /GRDFF/ N1N2FF,N1FFT
      SAVE   /GRDFF/
      INTEGER I,J,K
      INDRAM=2*((K-1)*N1N2FF + (J-1)*N1FFT + (I-1)) + 1
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION MODF(I,J)
      INTEGER I,J
      IF (J.EQ.1) THEN
        MODF=1
      ELSE
        MODF=MOD(I,J)
        IF ((MODF.EQ.0).AND.(I.NE.0)) MODF=J
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION NFFT(N)
      INTEGER N
      REAL AUX
      AUX=LOG10(FLOAT(N))/LOG10(2.)
      IF (AUX.GT.AINT(AUX)) AUX=AUX+1.
      NFFT=2**AINT(AUX)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE NCHECK(N1,N2)
      INTEGER N1,N2
      EXTERNAL NFFT
      INTEGER NFFT
      IF (N1.LT.N2) THEN
C       GRDFFT-13
        CALL ERROR
     *       ('GRDFFT-13: Wrong specification of the output grid.')
C       Number of gridpoints for FFT must be greater than or equal to
C       the number of gridpoints in data.
      ENDIF
      IF (N1.NE.NFFT(N1)) THEN
C       GRDFFT-14
        CALL ERROR
     *       ('GRDFFT-14: Wrong specification of the output grid.')
C       Number of gridpoints for FFT must be an integer power of 2.
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'fourn.for'
C     fourn.for of Numerical Recipes
C
C=======================================================================
C