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: '
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')
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
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,
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                                                                 ```