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: 6.70
C Date: 2012, November 7
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.htm
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 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 NFFT,INDRAM,MODF,NCHECK,ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,
     *RARRAY,WARRAY,FOURN,UARRAY
      REAL UARRAY
      INTEGER NFFT,INDRAM,MODF
C     NFFT,INDRAM,MODF,NCHECK... This file.
C     ERROR... File error.for.
C     RSEP1,RSEP3T,RSEP3R,RSEP3I... File sep.for.
C     RARRAY,WARRAY,UARRAY... File forms.for.
C     FOURN... File 'fourn.for' of Numerical Recipes.
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)
      DOUBLE PRECISION PID
      PARAMETER (PID=3.141592653589793D0)
      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
      REAL O1,O2,O3,D1,D2,D3,FFT,FFTFIL,
     * 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,M1,M2,M3
      INTEGER IRAM,I1,I2,I3,I4,I,J,K,L,NDIMFF,NFORFF,OFORFF
      REAL RRA,RRB,RR0,RRD,RIA,RIB,RI0,RID,RRK
      DOUBLE PRECISION  FFTD,O1D,O2D,O3D,D1D,D2D,D3D,M1D,M2D,M3D,
     *  D1OUTD,D2OUTD,D3OUTD,O1OUTD,O2OUTD,O3OUTD,
     *  AUXD,AUXRD,AUXID,RM0RD,RM0ID,RMTRD,RMTID,RMULTD
C
      UNDEF=UARRAY()
C
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)
      FFTD=DBLE(FFT)
      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
          FFTD=DBLE(I)*PID
        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.')
C       The dimension MRAM in the common block RAMC defined in the file
C       ram.inc should be
C       at least  2*(N1FFT*N2FFT*N3FFT +
C          + MAX(N1*N2*N3,N1OUT*N2OUT*N3OUT,N1FFT*N2FFT*N3FFT).
      ENDIF
C
C     Preparing the quantities needed in double precision:
      O1D=DBLE(O1)
      O2D=DBLE(O2)
      O3D=DBLE(O3)
      D1D=DBLE(D1)
      D2D=DBLE(D2)
      D3D=DBLE(D3)
      D1OUTD=DBLE(D1OUT)
      D2OUTD=DBLE(D2OUT)
      D3OUTD=DBLE(D3OUT)
      O1OUTD=DBLE(O1OUT)
      O2OUTD=DBLE(O2OUT)
      O3OUTD=DBLE(O3OUT)
C     Preparing number NDIMFF describing the dimension
C     of the part of the input grid to be transformed.
      NDIMFF=3
      IF ((N1FFT.EQ.1).OR.(NT1FFT.EQ.0))  NDIMFF=NDIMFF-1
      IF ((N2FFT.EQ.1).OR.(NT2FFT.EQ.0))  NDIMFF=NDIMFF-1
      IF ((N3FFT.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:
      RMULTD=1.D0
      IF ((N1FFT.NE.1).AND.(NT1FFT.NE.0)) RMULTD=RMULTD*D1D
      IF ((N2FFT.NE.1).AND.(NT2FFT.NE.0)) RMULTD=RMULTD*D2D
      IF ((N3FFT.NE.1).AND.(NT3FFT.NE.0)) RMULTD=RMULTD*D3D
      RMULTD=RMULTD*DSQRT(DABS(FFTD)/(2.D0*PID))**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+1
        II=MRAM-NN+1
        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
        WRITE(*,'(A)') '+GRDFFT: Working...            '
C
        I=N1FFT-N1
        M1=I/2
        M1D=DBLE(M1)
        I1MI=1+M1
        I1MA=M1+N1
        I=N2FFT-N2
        M2=I/2
        M2D=DBLE(M2)
        I2MI=1+M2
        I2MA=M2+N2
        I=N3FFT-N3
        M3=I/2
        M3D=DBLE(M3)
        I3MI=1+M3
        I3MA=M3+N3
        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
C       Recording the known values:
        IRAM=1
        DO 24, I3=1,N3FFT
          DO 23, I2=1,N2FFT
            DO 22, I1=1,N1FFT
              IF (IRAM+1.GE.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=IR + I1-I1MI + (I2-I2MI)*N1 + (I3-I3MI)*N1N2
                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       Adding the before-FFT multiplicative factor:
        DO 58, I3=1,N3FFT
          DO 57, I2=1,N2FFT
            DO 56, I1=1,N1FFT
              RMTRD=1.D0
              RMTID=0.D0
              IF ((N1FFT.NE.1).AND.(NT1FFT.NE.0)) THEN
                AUXD=FFTD*
     *          (DBLE(I1-1)*D1D*(O1OUTD-DBLE(NINT(O1OUT/D1OUT))*D1OUTD))
                AUXRD=DCOS(AUXD)
                AUXID=DSIN(AUXD)
                RM0RD=RMTRD
                RM0ID=RMTID
                RMTRD=RM0RD*AUXRD - RM0ID*AUXID
                RMTID=RM0RD*AUXID + RM0ID*AUXRD
              ENDIF
              IF ((N2FFT.NE.1).AND.(NT2FFT.NE.0)) THEN
                AUXD=FFTD*
     *          (DBLE(I2-1)*D2D*(O2OUTD-DBLE(NINT(O2OUT/D2OUT))*D2OUTD))
                AUXRD=DCOS(AUXD)
                AUXID=DSIN(AUXD)
                RM0RD=RMTRD
                RM0ID=RMTID
                RMTRD=RM0RD*AUXRD - RM0ID*AUXID
                RMTID=RM0RD*AUXID + RM0ID*AUXRD
              ENDIF
              IF ((N3FFT.NE.1).AND.(NT3FFT.NE.0)) THEN
                AUXD=FFTD*
     *          (DBLE(I3-1)*D3D*(O3OUTD-DBLE(NINT(O3OUT/D3OUT))*D3OUTD))
                AUXRD=DCOS(AUXD)
                AUXID=DSIN(AUXD)
                RM0RD=RMTRD
                RM0ID=RMTID
                RMTRD=RM0RD*AUXRD - RM0ID*AUXID
                RMTID=RM0RD*AUXID + RM0ID*AUXRD
              ENDIF
              RM0RD=DBLE(RAM(INDRAM(I1,I2,I3)  ))
              RM0ID=DBLE(RAM(INDRAM(I1,I2,I3)+1))
              RAM(INDRAM(I1,I2,I3)  )=SNGL(RM0RD*RMTRD - RM0ID*RMTID)
              RAM(INDRAM(I1,I2,I3)+1)=SNGL(RM0ID*RMTRD + RM0RD*RMTID)
  56        CONTINUE
  57      CONTINUE
  58    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       FFT 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 FFT loop over subgrids.
C
C       Reordering the computed values to the output field,
C       adding the multiplicative factor:
        IO1=MODF(NINT(O1OUT/D1OUT)+1,N1FFT)
        IF (IO1.LT.0) IO1=IO1+N1FFT
        IO2=MODF(NINT(O2OUT/D2OUT)+1,N2FFT)
        IF (IO2.LT.0) IO2=IO2+N2FFT
        IO3=MODF(NINT(O3OUT/D3OUT)+1,N3FFT)
        IF (IO3.LT.0) IO3=IO3+N3FFT
        IF ((IO1.LE.0).OR.(IO2.LE.0).OR.(IO3.LE.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
        IR=MRAM-2*NNOUT+1
        II=MRAM-NNOUT+1
        DO 74, I3=1,N3OUT
          DO 73, I2=1,N2OUT
            DO 72, I1=1,N1OUT
              I=INDRAM(MODF(IO1+I1-1,N1FFT),MODF(IO2+I2-1,N2FFT),
     *                 MODF(IO3+I3-1,N3FFT))
              RMTRD=RMULTD
              RMTID=DBLE(0.)
              IF ((N1FFT.NE.1).AND.(NT1FFT.NE.0)) THEN
                AUXD=FFTD*(O1OUTD+DBLE(I1-1)*D1OUTD)*(O1D-M1D*D1D)
                AUXRD=DCOS(AUXD)
                AUXID=DSIN(AUXD)
                RM0RD=RMTRD
                RM0ID=RMTID
                RMTRD=RM0RD*AUXRD - RM0ID*AUXID
                RMTID=RM0RD*AUXID + RM0ID*AUXRD
              ENDIF
              IF ((N2FFT.NE.1).AND.(NT2FFT.NE.0)) THEN
                AUXD=FFTD*(O2OUTD+DBLE(I2-1)*D2OUTD)*(O2D-M2D*D2D)
                AUXRD=DCOS(AUXD)
                AUXID=DSIN(AUXD)
                RM0RD=RMTRD
                RM0ID=RMTID
                RMTRD=RM0RD*AUXRD - RM0ID*AUXID
                RMTID=RM0RD*AUXID + RM0ID*AUXRD
              ENDIF
              IF ((N3FFT.NE.1).AND.(NT3FFT.NE.0)) THEN
                AUXD=FFTD*(O3OUTD+DBLE(I3-1)*D3OUTD)*(O3D-M3D*D3D)
                AUXRD=DCOS(AUXD)
                AUXID=DSIN(AUXD)
                RM0RD=RMTRD
                RM0ID=RMTID
                RMTRD=RM0RD*AUXRD - RM0ID*AUXID
                RMTID=RM0RD*AUXID + RM0ID*AUXRD
              ENDIF
              RM0RD=DBLE(RAM(I))
              RM0ID=DBLE(RAM(I+1))
              RAM(IR)=SNGL(RM0RD*RMTRD - RM0ID*RMTID)
              RAM(II)=SNGL(RM0ID*RMTRD + RM0RD*RMTID)
              IR=IR+1
              II=II+1
  72        CONTINUE
  73      CONTINUE
  74    CONTINUE
C
C
C       Writing the results of the FFT:
        IR=MRAM-2*NNOUT+1
        II=MRAM-NNOUT+1
        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
      MODF=MOD(I,J)
      IF (MODF.EQ.0) MODF=J
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION NFFT(N)
      INTEGER N
      IF     (N.LE.         1) THEN
              NFFT=         1
      ELSEIF (N.LE.         2) THEN
              NFFT=         2
      ELSEIF (N.LE.         4) THEN
              NFFT=         4
      ELSEIF (N.LE.         8) THEN
              NFFT=         8
      ELSEIF (N.LE.        16) THEN
              NFFT=        16
      ELSEIF (N.LE.        32) THEN
              NFFT=        32
      ELSEIF (N.LE.        64) THEN
              NFFT=        64
      ELSEIF (N.LE.       128) THEN
              NFFT=       128
      ELSEIF (N.LE.       256) THEN
              NFFT=       256
      ELSEIF (N.LE.       512) THEN
              NFFT=       512
      ELSEIF (N.LE.      1024) THEN
              NFFT=      1024
      ELSEIF (N.LE.      2048) THEN
              NFFT=      2048
      ELSEIF (N.LE.      4096) THEN
              NFFT=      4096
      ELSEIF (N.LE.      8192) THEN
              NFFT=      8192
      ELSEIF (N.LE.     16384) THEN
              NFFT=     16384
      ELSEIF (N.LE.     32768) THEN
              NFFT=     32768
      ELSEIF (N.LE.     65536) THEN
              NFFT=     65536
      ELSEIF (N.LE.    131072) THEN
              NFFT=    131072
      ELSEIF (N.LE.    262144) THEN
              NFFT=    262144
      ELSEIF (N.LE.    524288) THEN
              NFFT=    524288
      ELSEIF (N.LE.   1048576) THEN
              NFFT=   1048576
      ELSEIF (N.LE.   2097152) THEN
              NFFT=   2097152
      ELSEIF (N.LE.   4194304) THEN
              NFFT=   4194304
      ELSEIF (N.LE.   8388608) THEN
              NFFT=   8388608
      ELSEIF (N.LE.  16777216) THEN
              NFFT=  16777216
      ELSEIF (N.LE.  33554432) THEN
              NFFT=  33554432
      ELSEIF (N.LE.  67108864) THEN
              NFFT=  67108864
      ELSEIF (N.LE. 134217728) THEN
              NFFT= 134217728
      ELSEIF (N.LE. 268435456) THEN
              NFFT= 268435456
      ELSEIF (N.LE. 536870912) THEN
              NFFT= 536870912
      ELSEIF (N.LE.1073741824) THEN
              NFFT=1073741824
      ELSE
C       GRDFFT-09
        CALL ERROR ('GRDFFT-09: Too large N')
C       One of the N1, N2, N3 or N4 is greater than 2**31.
      ENDIF
cc      Old:
cc      REAL AUX
cc      AUX=LOG10(FLOAT(N))/LOG10(2.)
cc      NFFT=2**INT(AUX+0.999999)
      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 grid for FFT')
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 grid for FFT')
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 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'fourn.for'
C     fourn.for of Numerical Recipes
C
C=======================================================================
C