!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
module wind_fft
!$$$ module documentation block
! . . . .
! module: wind_fft
! prgmmr: pondeca org: np23 date: 2006-08-01
!
! abstract: contains subroutines for 2-dimensional FFTs on the plane.
! adapted from jeff whitaker's August 1, 1991 two-layer
! qg-model (qg2l)
!
! program history log:
! 2006-08-01 pondeca
!
! subroutines included:
! sub divvort_to_psichi
! sub fft2d
! sub rfft
! sub cfft
! sub FFT2
! sub PASS2
! sub PREFFT
! sub FACTOR
!
! variable definitions
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind,r_single,r_kind
use constants, only: zero,half,one,two,four
implicit none
! set default to private
private
! set subroutines to public
public :: divvort_to_psichi
public :: fft2d
public :: rfft
public :: cfft
public :: FFT2
public :: PASS2
public :: PREFFT
public :: FACTOR
integer(i_kind) nx,ny,mmax,nwavesx,nmax,nwavesy
integer(i_kind) nfax0,nfay0
real(r_kind),allocatable::rk(:),rl(:)
real(r_kind),allocatable::indxy(:),trigx(:,:),trigy(:,:)
real(r_kind),allocatable::delimn(:,:)
complex(r_kind),allocatable::scrmn(:,:,:,:),scr(:,:),qmn(:,:,:)
integer(i_kind) ifax0(20),ifay0(20)
!-------------------------------------------------------------------------
contains
!-------------------------------------------------------------------------
! NOAA/NCEP, National Centers for Environmental Prediction GSI !
!-------------------------------------------------------------------------
subroutine divvort_to_psichi(nx0,ny0,mmax0,nmax0,rld0,qg)
!$$$ subprogram documentation block
! . . . .
! subprogram: divvort_to_psichi
! prgmmr: pondeca org: np22 date: 2006-08-01
!
!
! abstract: use two-dimensional FFTs to convert gridded
! fields of divergence and vorticity into
! stream-function and velocity-potential
!
! program history log:
! 2005-02-08 pondeca
!
! input argument list:
! nx0,ny0 - number of grid points in x/y
! mmax0,nmax0
! rld0
!
! output argument list:
! qg
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
integer(i_kind), intent(in ) :: nx0,ny0
integer(i_kind), intent(in ) :: mmax0,nmax0
real(r_single) , intent(in ) :: rld0(nx0,ny0)
real(r_single) , intent( out) :: qg(nx0,ny0,2)
real(r_kind) xmax,ymax,pi,wavey,delmn
integer(i_kind) i,j,k,ny2,m,mm,indx
nx=nx0
ny=ny0
mmax=mmax0
nmax=nmax0
nwavesx = mmax+1
nwavesy = 2*nmax+1
! write(6,*) ' in divvort_to_psichi: nx,ny,nwavesx,nwavesy=',nx,ny,nwavesx,nwavesy
allocate(scrmn(nwavesx,nwavesy,2,4))
allocate(scr(nwavesx,ny))
allocate(qmn(nwavesx,nwavesy,2))
allocate(delimn(nwavesx,nwavesy))
allocate(rk(nwavesx))
allocate(rl(nwavesy))
allocate(indxy(nwavesy))
allocate(trigx(2,nx))
allocate(trigy(2,ny))
pi = four*atan(one)
xmax = float(nx-1)
ymax = float(ny-1)
! write(6,*) ' in divvort_to_psichi: xmax,ymax,=',xmax,ymax
!
!==> compute trig tables for fft routines.
!
call prefft(nx,nfax0,ifax0,trigx)
call prefft(ny,nfay0,ifay0,trigy)
!
!==> set up wavenumbers used in fourier differentiation.
do 100 i=1,nwavesx
rk(i) = two*pi*float(i-1)/xmax
100 continue
ny2 = (ny/2)+1
do 200 j=1,ny
mm = j/ny2
m = mm*ny+1
wavey = two*pi*float(j-m)/ymax
if (j<=nmax+1 .or. j>=ny-nmax+1) then
indx = j-m+nmax+1
indxy(indx) = j
rl(indx) = wavey
end if
200 continue
!
do 300 j=1,nwavesy
do 300 i=1,nwavesx
delmn = -(rk(i)*rk(i) + rl(j)*rl(j))
delimn(i,j) = zero
if (delmn/=zero) delimn(i,j)=one/delmn
300 continue
! write(6,*) 'divvort_to_psichi: delimn,min,max=',minval(delimn),maxval(delimn)
do 600 k=1,2
scr=(zero,zero)
call fft2d(qg(1,1,k),qmn(1,1,k),scr,+1)
do 500 j=1,nwavesy
do 500 i=1,nwavesx
qmn(i,j,k)=delimn(i,j)*qmn(i,j,k)
500 continue
scr=(zero,zero)
call fft2d(qg(1,1,k),qmn(1,1,k),scr,-1)
qg(:,:,k)=qg(:,:,k)*rld0(:,:)**2
600 continue
! write(6,*) ' in divvort_to_psichi: all done'
deallocate(scrmn)
deallocate(scr)
deallocate(qmn)
deallocate(delimn)
deallocate(rk)
deallocate(rl)
deallocate(indxy)
deallocate(trigx)
deallocate(trigy)
return
end subroutine divvort_to_psichi
!
subroutine fft2d(data,coeff,scr,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: fft2d
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! idir
! data
! scr,coeff
!
! output argument list:
! data
! scr,coeff
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind), intent(in ) :: idir
real(r_single) , intent(inout) :: data(nx,ny)
complex(r_kind), intent(inout) :: scr(nwavesx,ny),coeff(nwavesx,nwavesy)
! write(6,*) 'IN fft2d,nx,ny,nwavesx,nwavesy=',nx,ny,nwavesx,nwavesy
!-------------------------------------------------
!==> forward transform.
!-------------------------------------------------
if (idir==1) then
call rfft(data,scr,1)
call cfft(scr,coeff,1)
!-------------------------------------------------
!==> inverse transform.
!-------------------------------------------------
else if (idir==-1) then
call cfft(scr,coeff,-1)
call rfft(data,scr,-1)
!
else
write(6,*) ' idir must be +1 or -1 in fft2d!'
write(6,*) ' execution terminated.'
stop
!
end if
!
return
end subroutine fft2d
!
subroutine rfft(data,coeff,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: rfft
! prgmmr:
!
! abstract: performs fast multiple real fft's pairwise using a
! multiple complex fft routine. the number of fft's
! to be performed must be even.
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! idir
!
! output argument list:
! data
! coeff
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind),intent(in ) :: idir
real(r_single) ,intent( out) :: data(nx,ny)
complex(r_kind),intent( out) :: coeff(nwavesx,ny)
integer(i_kind) i,j,n,npts,ndata,ndatah
real(r_kind) a(ny/2,nx,2),c(ny/2,nx,2)
npts = nx
ndata = ny
ndatah = ndata/2
!
!----------------------
!==> forward transform.
!----------------------
if (idir==+1) then
!
!==> copy the data into the work array.
! transforms are computed pairwise using a complex fft.
!
do 10 j=1,npts
do 10 i=1,ndatah
a(i,j,1) = data(j,i)
a(i,j,2) = data(j,i+ndatah)
10 continue
!
call fft2(a,c,npts,npts-1,nfax0,ifax0,-1,trigx,ndatah)
!
do 20 i=1,ndatah
coeff(1,i) = c(i,1,1)
coeff(1,i+ndatah) = c(i,1,2)
do 20 n=2,nwavesx
coeff(n,i) = half*cmplx(c(i,n,1),c(i,n,2)) + &
half*cmplx(c(i,npts-n+2,1),-c(i,npts-n+2,2))
coeff(n,i+ndatah) = half*cmplx(c(i,n,2),-c(i,n,1)) + &
half*cmplx(c(i,npts-n+2,2),c(i,npts-n+2,1))
20 continue
!----------------------
!==> inverse transform.
!----------------------
else if (idir==-1) then
!
do 25 j=1,npts
do 25 i=1,ndatah
c(i,j,1) = zero
c(i,j,2) = zero
25 continue
do 30 i=1,ndatah
c(i,1,1) = real(coeff(1,i))
c(i,1,2) = real(coeff(1,i+ndatah))
do 30 n=2,nwavesx
c(i,npts-n+2,1) = real(coeff(n,i))+aimag(coeff(n,i+ndatah))
c(i,n,1) = real(coeff(n,i))-aimag(coeff(n,i+ndatah))
c(i,n,2) = aimag(coeff(n,i))+real(coeff(n,i+ndatah))
c(i,npts-n+2,2) = real(coeff(n,i+ndatah))-aimag(coeff(n,i))
30 continue
!
call fft2(c,a,npts,npts-1,nfax0,ifax0,+1,trigx,ndatah)
!
do 40 j=1,npts
do 40 i=1,ndatah
data(j,i) = a(i,j,1)
data(j,i+ndatah) = a(i,j,2)
40 continue
!
else
write(6,*) ' idir must be +1 or -1 in rfft!'
write(6,*) ' execution terminated.'
stop
end if
!
return
end subroutine rfft
!
subroutine cfft(data,coeff,idir)
!$$$ subprogram documentation block
! . . . .
! subprogram: cfft
! prgmmr:
!
! abstract: performs fast multiple complex fft's.
! array data(ndata,npts) contains ndata distinct complex
! data sets of length npts to be transformed.
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! idir
!
! output argument list:
! data
! coeff
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind),intent(in ) :: idir
complex(r_kind),intent( out) :: data(nwavesx,ny),coeff(nwavesx,nwavesy)
integer(i_kind) i,j,jj,npts,ndata
real(r_kind) a(nwavesx,ny,2),c(nwavesx,ny,2)
npts = ny
ndata = nwavesx
!----------------------
!==> forward transform.
!----------------------
if (idir==+1) then
!
!==> copy the data into the work array.
!
do 10 j=1,npts
do 10 i=1,ndata
a(i,j,1) = real(data(i,j))
a(i,j,2) = aimag(data(i,j))
10 continue
!
call fft2(a,c,npts,npts-1,nfax0,ifax0,-1,trigy,ndata)
!
do 20 j=1,nwavesy
do 20 i=1,ndata
jj = indxy(j)
coeff(i,j) = cmplx(c(i,jj,1),c(i,jj,2))
20 continue
!----------------------
!==> inverse transform.
!----------------------
else if (idir==-1) then
!
do 25 j=1,npts
do 25 i=1,ndata
c(i,j,1) = zero
c(i,j,2) = zero
25 continue
do 30 j=1,nwavesy
do 30 i=1,ndata
jj = indxy(j)
c(i,jj,1) = real(coeff(i,j))
c(i,jj,2) = aimag(coeff(i,j))
30 continue
!
call fft2(c,a,npts,npts-1,nfax0,ifax0,+1,trigy,ndata)
!
do 40 j=1,npts
do 40 i=1,ndata
data(i,j) = cmplx(a(i,j,1),a(i,j,2))
40 continue
!
else
write(6,*) ' idir must be +1 or -1 in cfft!'
write(6,*) ' execution terminated.'
stop
end if
!
return
end subroutine cfft
!
SUBROUTINE FFT2(A,C,N,NDIM,NFAX,IFAX,ISIGN,TRIG,LEN)
!$$$ subprogram documentation block
! . . . .
! subprogram: FFT2
! prgmmr:
!
! abstract: PERFORMS A COMPLEX FFT ON MULTIPLE DATA IN A
! AND RETURNS THE RESULT IN C
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! A - INPUT ARRAY (DESTROYED DURING CALCULATION).
! DIMENSIONED AS A(LEN,0;NDIM,2), WHERE THE 1ST INDEX
! LABELS DISTINCT DATA TO BE TRANSFORMED, & 2ND INDEX LABELS
! THE N POINTS TO BE TRANSFORMED, & 3RD INDEX LABELS
! REAL (1) OR IMAGINARY (2) PARTS
! N - NO. OF POINTS IN TRANSFORM DIRECTION, MUST HAVE ONLY
! PRIME FACTORS OF 2 & 3
! NDIM - 2ND DIMENSION OF A & C
! NFAX - NO. OF PRIME FACTORS OF N
! IFAX
! ISIGN - SET ISIGN = -1 TO COMPUTE FOURIER COEFFICIENTS,
! = +1 TO COMPUTE GRID VALUES.
! TRIG - ARRAY CONTAINING TRIGONOMETRIC FACTORS:
! DIMENSIONED TRIG(2,0:N-1)
! TRIG(1,J) = COS (2.*PI*J/N)
! TRIG(2,J) = SIN (2.*PI*J/N)
! LEN - NO. OF DISTINCT TRANSFORMS TO BE PERFORMED.
!
! output argument list:
! C - OUTPUT ARRAY (MUST BE DISTINCT FROM INPUT ARRAY),
! DIMENSIONED THE SAME AS ARRAY A
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
INTEGER(i_kind),intent(in ) :: IFAX(*),LEN,NDIM,NFAX,ISIGN,N
REAL(r_kind) ,intent(inout) :: A(LEN,0:NDIM,2),TRIG(2,0:N-1)
REAL(r_kind) ,intent(inout) :: C(LEN,0:NDIM,2)
INTEGER(i_kind) I,IJ,LA,IFAC
REAL(r_kind) XNI
LOGICAL ODD
REAL(r_kind),PARAMETER::PI=3.1415926535898_r_kind
!
LA=1
ODD=.TRUE.
DO 10 I=1,NFAX
IFAC=IFAX(I)
IF (ODD) THEN
CALL PASS2(A,C,N,NDIM,ISIGN,IFAC,LA,TRIG,LEN)
ELSE
CALL PASS2(C,A,N,NDIM,ISIGN,IFAC,LA,TRIG,LEN)
END IF
ODD=.NOT.ODD
LA=LA*IFAC
10 CONTINUE
!
IF (ODD) THEN
DO 30 I=0,N-1
DO 20 IJ=1,LEN
C(IJ,I,1) = A(IJ,I,1)
C(IJ,I,2) = A(IJ,I,2)
20 CONTINUE
30 CONTINUE
END IF
IF (ISIGN==-1) THEN
XNI=one/N
DO 50 I=0,N-1
DO 40 IJ=1,LEN
C(IJ,I,1) = XNI * C(IJ,I,1)
C(IJ,I,2) = XNI * C(IJ,I,2)
40 CONTINUE
50 CONTINUE
END IF
RETURN
END SUBROUTINE FFT2
! ------------------------------------------------------------
SUBROUTINE PASS2(A,C,N,NDIM,ISIGN,IFAC,LA,TRIG,LEN)
!$$$ subprogram documentation block
! . . . .
! subprogram: PASS2
! prgmmr:
!
! abstract: PERFORMS ONE PASS OF FFT2
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! A
! N
! NDIM
! ISIGN
! IFAC
! LA
! TRIG
! LEN
!
! output argument list:
! C
!
! note:
!
! THIS ROUTINE IS NEVER CALLED DIRECTLY BY THE USER.
! THE ARGUMENTS ARE SIMILAR TO THOSE OF FFT2.
!
! THE INNER LOOPS IN THIS SUBROUTINE (THOSE OVER THE INDEX IJ)
! EXTEND OVER PART OF THE 2ND DIMENSIONS OF THE ARRAYS A & C.
! THIS PRODUCES LONGER VECTOR LENGTHS.
!
! THE FOLLOWING CONSTANTS PRESUME A 64-BIT MACHINE. THEY ARE
! SINES OF 45 AND 60 DEGREES.
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
INTEGER(i_kind),intent(in ) :: N,NDIM,ISIGN,IFAC,LA,LEN
REAL(r_kind) ,intent(in ) :: TRIG(2,0:N-1)
REAL(r_kind) ,intent(inout) :: A(LEN,0:NDIM,2),C(LEN,0:NDIM,2)
INTEGER(i_kind) I,I0,I01,I1,I11,J,J0,J01,J1,J11,K,M,LLA,IND(0:20),JND(0:20),JUMP,IJ,IJ1
INTEGER(i_kind) I2,I21,J2,J21
REAL(r_kind) SN60,CC,SS,AM1,AM2,AP1,AP2,TA1,TA2,C1,C2,S1,S2,T1,T2
REAL(r_kind),PARAMETER::ASN60 = half*1.732050807569_r_kind
SN60=ISIGN * ASN60
M = N/IFAC
!
!
! SET UP INDEXING
!
DO 10 K=0,IFAC-1
IND(K) = K*M
JND(K) = K*LA
10 CONTINUE
LLA =LA * LEN
!
! PERFORM THE ARITHMETIC
!
I = 0
J = 0
JUMP = (IFAC-1) * LA
DO 130 K = 0,M-LA,LA
IF (IFAC==2) THEN
I0 = IND(0) + I
I1 = IND(1) + I
J0 = JND(0) + J
J1 = JND(1) + J
CC = TRIG(1,K)
SS = ISIGN * TRIG(2,K)
IF (K==0) THEN
DO 20 IJ = 1,LLA
if(ij>len) then
ij1=mod(ij-1,len)+1
i01=(ij-1)/len+i0; i11=(ij-1)/len+i1
j01=(ij-1)/len+j0; j11=(ij-1)/len+j1
if(i01<=ndim.and.i11<=ndim.and.&
j01<=ndim.and.j11<=ndim) then
C(IJ1,J01,1) = A(IJ1,I01,1) + A(IJ1,I11,1)
C(IJ1,J01,2) = A(IJ1,I01,2) + A(IJ1,I11,2)
C(IJ1,J11,1) = A(IJ1,I01,1) - A(IJ1,I11,1)
C(IJ1,J11,2) = A(IJ1,I01,2) - A(IJ1,I11,2)
end if
else
C(IJ,J0,1) = A(IJ,I0,1) + A(IJ,I1,1)
C(IJ,J0,2) = A(IJ,I0,2) + A(IJ,I1,2)
C(IJ,J1,1) = A(IJ,I0,1) - A(IJ,I1,1)
C(IJ,J1,2) = A(IJ,I0,2) - A(IJ,I1,2)
end if
20 CONTINUE
ELSE
DO 50 IJ = 1,LLA
if(ij>len) then
ij1=mod(ij-1,len)+1
i01=(ij-1)/len+i0; i11=(ij-1)/len+i1
j01=(ij-1)/len+j0; j11=(ij-1)/len+j1
if(i01<=ndim.and.i11<=ndim.and.&
j01<=ndim.and.j11<=ndim) then
C(IJ1,J01,1) = A(IJ1,I01,1) + A(IJ1,I11,1)
C(IJ1,J01,2) = A(IJ1,I01,2) + A(IJ1,I11,2)
AM1 = A(IJ1,I01,1) - A(IJ1,I11,1)
AM2 = A(IJ1,I01,2) - A(IJ1,I11,2)
C(IJ1,J11,1) = CC * AM1 - SS * AM2
C(IJ1,J11,2) = SS * AM1 + CC * AM2
end if
else
C(IJ,J0,1) = A(IJ,I0,1) + A(IJ,I1,1)
C(IJ,J0,2) = A(IJ,I0,2) + A(IJ,I1,2)
AM1 = A(IJ,I0,1) - A(IJ,I1,1)
AM2 = A(IJ,I0,2) - A(IJ,I1,2)
C(IJ,J1,1) = CC * AM1 - SS * AM2
C(IJ,J1,2) = SS * AM1 + CC * AM2
end if
50 CONTINUE
END IF
ELSEIF (IFAC==3) THEN
I0 = IND(0) + I
I1 = IND(1) + I
I2 = IND(2) + I
J0 = JND(0) + J
J1 = JND(1) + J
J2 = JND(2) + J
IF (K==0) THEN
DO 60 IJ = 1,LLA
if(ij>len) then
ij1=mod(ij-1,len)+1
i01=(ij-1)/len+i0; i11=(ij-1)/len+i1; i21=(ij-1)/len+i2
j01=(ij-1)/len+j0; j11=(ij-1)/len+j1; j21=(ij-1)/len+j2
if(i01<=ndim.and.i11<=ndim.and.i21<=ndim.and. &
j01<=ndim.and.j11<=ndim.and.i21<=ndim) then
AP1 = A(IJ1,I11,1) + A(IJ1,I21,1)
AP2 = A(IJ1,I11,2) + A(IJ1,I21,2)
C(IJ1,J01,1) = A(IJ1,I01,1) + AP1
C(IJ1,J01,2) = A(IJ1,I01,2) + AP2
TA1 = A(IJ1,I01,1) - half * AP1
TA2 = A(IJ1,I01,2) - half * AP2
AM1 = SN60 * (A(IJ1,I11,1) - A(IJ1,I21,1))
AM2 = SN60 * (A(IJ1,I11,2) - A(IJ1,I21,2))
C(IJ1,J11,1) = TA1 - AM2
C(IJ1,J11,2) = TA2 + AM1
C(IJ1,J21,1) = TA1 + AM2
C(IJ1,J21,2) = TA2 - AM1
endif
else
AP1 = A(IJ,I1,1) + A(IJ,I2,1)
AP2 = A(IJ,I1,2) + A(IJ,I2,2)
C(IJ,J0,1) = A(IJ,I0,1) + AP1
C(IJ,J0,2) = A(IJ,I0,2) + AP2
TA1 = A(IJ,I0,1) - half * AP1
TA2 = A(IJ,I0,2) - half * AP2
AM1 = SN60 * (A(IJ,I1,1) - A(IJ,I2,1))
AM2 = SN60 * (A(IJ,I1,2) - A(IJ,I2,2))
C(IJ,J1,1) = TA1 - AM2
C(IJ,J1,2) = TA2 + AM1
C(IJ,J2,1) = TA1 + AM2
C(IJ,J2,2) = TA2 - AM1
end if
60 CONTINUE
ELSE
C1 = TRIG(1,K)
C2 = TRIG(1,2*K)
S1 = ISIGN * TRIG(2,K)
S2 = ISIGN * TRIG(2,2*K)
DO 70 IJ = 1,LLA
if(ij>len) then
ij1=mod(ij-1,len)+1
i01=(ij-1)/len+i0; i11=(ij-1)/len+i1; i21=(ij-1)/len+i2
j01=(ij-1)/len+j0; j11=(ij-1)/len+j1; j21=(ij-1)/len+j2
if(i01<=ndim.and.i11<=ndim.and.i21<=ndim.and. &
j01<=ndim.and.j11<=ndim.and.i21<=ndim) then
AP1 = A(IJ1,I11,1) +A(IJ1,I21,1)
AP2 = A(IJ1,I11,2) +A(IJ1,I21,2)
C(IJ1,J01,1) = A(IJ1,I01,1) +AP1
C(IJ1,J01,2) = A(IJ1,I01,2) +AP2
TA1 = A(IJ1,I01,1) - half*AP1
TA2 = A(IJ1,I01,2) - half*AP2
AM1 = SN60 * (A(IJ1,I11,1) - A(IJ1,I21,1))
AM2 = SN60 * (A(IJ1,I11,2) - A(IJ1,I21,2))
T1 = TA1 - AM2
T2 = TA2 + AM1
C(IJ1,J11,1) = C1 * T1 - S1 * T2
C(IJ1,J11,2) = S1 * T1 + C1 *T2
T1 = TA1 + AM2
T2 = TA2 - AM1
C(IJ1,J21,1) = C2 * T1 - S2 * T2
C(IJ1,J21,2) = S2 * T1 + C2 * T2
endif
else
AP1 = A(IJ,I1,1) +A(IJ,I2,1)
AP2 = A(IJ,I1,2) +A(IJ,I2,2)
C(IJ,J0,1) = A(IJ,I0,1) +AP1
C(IJ,J0,2) = A(IJ,I0,2) +AP2
TA1 = A(IJ,I0,1) - half*AP1
TA2 = A(IJ,I0,2) - half*AP2
AM1 = SN60 * (A(IJ,I1,1) - A(IJ,I2,1))
AM2 = SN60 * (A(IJ,I1,2) - A(IJ,I2,2))
T1 = TA1 - AM2
T2 = TA2 + AM1
C(IJ,J1,1) = C1 * T1 - S1 * T2
C(IJ,J1,2) = S1 * T1 + C1 *T2
T1 = TA1 + AM2
T2 = TA2 - AM1
C(IJ,J2,1) = C2 * T1 - S2 * T2
C(IJ,J2,2) = S2 * T1 + C2 * T2
end if
70 CONTINUE
END IF
END IF
I = I+LA
J = J+LA
J = J+JUMP
130 CONTINUE
RETURN
END SUBROUTINE PASS2
!-------------------------------------------------------------
SUBROUTINE PREFFT (N,NFAX,IFAX,TRIG)
!$$$ subprogram documentation block
! . . . .
! subprogram: PREFFT
! prgmmr:
!
! abstract: C0MPUTES PRELIMINARY QUANTITIES FOR FFT ROUTINES
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! N - NO. OF POINTS IN DATA, MUST HAVE PRIME FACTORS 2 & 3
! TRIG - ARRAY CONTAINING TRIGONOMETRIC FACTORS, DIMNSN (2,0;N-1)
! TRIG(1,J) = COS (2*PI*J/N)
! TRIG(2,J) = SIN (2.*PI*J/N)
!
! output argument list:
! NFAX - NO. OF PRIME FACTORS OF N (OUTPUT)
! IFAX - ARRAY CONTAINING PRIME FACTORS OF N (OUTPUT)
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
INTEGER(i_kind),intent(in ) :: N
REAL(r_kind) ,intent(inout) :: TRIG(2,0:N-1)
INTEGER(i_kind),intent( out) :: IFAX(*),NFAX
INTEGER(i_kind) K
REAL(r_kind) ARG
REAL(r_kind),PARAMETER::PI=3.1415926535898_r_kind
CALL FACTOR (N,NFAX,IFAX)
DO 10 K=0,N-1
ARG = two*PI*K/N
TRIG(1,K) = COS(ARG)
TRIG(2,K) = SIN(ARG)
10 CONTINUE
RETURN
END SUBROUTINE PREFFT
! --------------------------------------------------------------
SUBROUTINE FACTOR (N,NFAX,IFAX)
!$$$ subprogram documentation block
! . . . .
! subprogram: FACTOR
! prgmmr:
!
! abstract: COMPUTES THE FACTORS OF N
!
! program history log:
! 2009-10-02 lueken - added subprogram doc block
!
! input argument list:
! N
!
! output argument list:
! NFAX
! IFAX
!
! note:
!
! THIS ROUTINE IS NOT CALLED DIRECTLY BY THE USER
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
INTEGER(i_kind),intent(in ) :: N
INTEGER(i_kind),intent( out) :: IFAX(*),NFAX
INTEGER(i_kind) NN,II
NFAX = 0
NN= N
! EXTRACT FACTORS OF 3
DO 10 II = 1,20
IF (NN==3*(NN/3)) THEN
NFAX = NFAX+1
IFAX(NFAX) = 3
NN = NN/3
ELSE
EXIT
END IF
10 CONTINUE
! EXTRACT FACTORS OF 2
DO 30 II = NFAX+1,20
IF (NN==2*(NN/2)) THEN
NFAX = NFAX +1
IFAX(NFAX) =2
NN = NN/2
ELSE
EXIT
END IF
30 CONTINUE
IF (NN/=1) THEN
write(6,*) 'PORRA 4'
STOP
END IF
RETURN
END SUBROUTINE FACTOR
!=======================================================================
!=======================================================================
end module wind_fft