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