!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module pft_module #ifndef USE_LIMA use ecmfft implicit none private public pft2d, pft2d_phys, pft_init, pft_end ! integer ifax(13) ! FFT internal information real, allocatable :: trigs(:) real, allocatable :: sc(:), se(:) real, allocatable :: damp(:,:,:) integer js2g0 ! if "SGI_FFT" is defined SGI's FFT library will be used (instead of the ! default FFT f77 code from ECMWF) ! (Need to link to SGI library with "-lcomplib.sgimath" ) ! if "ALT_PFT" is defined the prognostic winds on the D grid will be directly ! filtered (instead of the tendencies), and there will be 2 less calls ! to pft2d. However, the faster algebraic 3-pt filter will not be used. ! Programmer: S.-J. Lin, NOAA/GFDL, Princeton, NJ CONTAINS !----------------------------------------------------------------------- !BOP ! !IROUTINE: pft_init --- Calculate algebraic and FFT polar filters ! ! !INTERFACE: subroutine pft_init(im, jm, beglat, endlat, ighost, cosp, cose, print_msg) ! !USES: ! !INPUT PARAMETERS: integer, intent(in):: im ! Global X dimension integer, intent(in):: jm ! Global Y dimension integer, intent(in):: beglat, endlat integer, intent(in):: ighost ! Normally 0 real, intent(in):: cosp(jm) ! cosine array real, intent(in):: cose(jm) ! cosine array logical, intent(in):: print_msg ! !DESCRIPTION: ! ! Compute coefficients for the 3-point algebraic and the FFT ! polar filters. ! !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: integer ipft1, ipft2 integer jn1g1 ! j north limit ghosted 1 (starts jm) integer i, j real ycrit ! critical value real lat_s real dl, coszc, cutoff, phi, dmp, rat real pi cutoff = 1.e-20 pi = 4.d0 * datan(1.d0) ipft1 = max(2, beglat-ighost) js2g0 = max(2, beglat) ipft2 = min(jm-1, endlat+ighost) jn1g1 = min(jm, endlat+1) allocate( sc(ipft1:ipft2) ) allocate( se(js2g0:jn1g1) ) allocate( damp(im,ipft1:jn1g1,2) ) allocate( trigs(3*im/2+1) ) call fftfax(im, trigs) !--------------------------------------------- ! Determine ycrit such that effective DX >= DY !--------------------------------------------- rat = float(im)/float(2*(jm-1)) ycrit = acos( min(0.81, rat) ) * (180./pi) if(print_msg) then write(6,*) 'Initializing polar filter .........' #ifdef ALT_PFT write(6,*) 'Starting latitude for FFT pft=', ycrit #else write(6,*) 'Starting latitude for algebraic pft=', ycrit #endif endif coszc = cos(ycrit*pi/180.) ! INIT fft polar coefficients: dl = pi/dble(im) do j=ipft1, ipft2 do i=1,im damp(i,j,1) = 1. enddo enddo do j=js2g0,jn1g1 do i=1,im damp(i,j,2) = 1. enddo enddo !************ ! Cell center !************ do j=ipft1, ipft2 sc(j) = (coszc/cosp(j))**2 #ifdef ALT_PFT if(sc(j) > 1.) then #else if(sc(j) > 1. .and. sc(j) <= 2.0) then sc(j) = 1. + (sc(j)-1.)/(sc(j)+1.) elseif(sc(j) > 2. .and. sc(j) <= 4.) then sc(j) = 1. + sc(j)/(8.-sc(j)) sc(j) = min(2., sc(j)) lat_s = acos(cosp(j+1))*180/pi elseif(sc(j) > 4. ) then #endif ! FFT filter do i=1,im/2 phi = dl * i dmp = min((cosp(j)/coszc)/sin(phi), 1.)**2 if(dmp < cutoff) dmp = 0. damp(2*i-1,j,1) = dmp damp(2*i ,j,1) = dmp enddo endif enddo #if ( !defined ALT_PFT ) ! if(print_msg) write(6,*) 'Switching to FFT at latitude =', lat_s #endif !************ ! Cell edges !************ do j=js2g0,jn1g1 se(j) = (coszc/cose(j))**2 #ifdef ALT_PFT if( se(j) > 1. ) then #else if(se(j) > 1. .and. se(j) <= 2.0 ) then se(j) = 1. + (se(j)-1.)/(se(j)+1.) elseif(se(j) > 2. .and. se(j) <= 4.) then se(j) = 1. + se(j)/(8.-se(j)) se(j) = min(2., se(j)) elseif(se(j) > 4. ) then #endif ! FFT do i=1,im/2 phi = dl * i dmp = min((cose(j)/coszc)/sin(phi), 1.)**2 if(dmp < cutoff) dmp = 0. damp(2*i-1,j,2) = dmp damp(2*i ,j,2) = dmp enddo endif enddo !EOC end subroutine pft_init !----------------------------------------------------------------------- subroutine pft_end deallocate( trigs ) deallocate( sc ) deallocate( se ) deallocate( damp ) end subroutine pft_end !----------------------------------------------------------------------- !BOP ! !IROUTINE: pft2d --- Two-dimensional fast Fourier transform ! ! !INTERFACE: subroutine pft2d(p, im, jp, q1, q2, igrid) ! !USES: ! !INPUT PARAMETERS: integer, intent(in):: im ! Total X dimension integer, intent(in):: jp ! Total Y dimension integer, intent(in):: igrid ! 1: cell center; 2: edge ! !INPUT/OUTPUT real, intent(inout):: p(im,jp) ! Array to be polar filtered real, intent(inout):: q1( im+2, *) ! Work array real, intent(inout):: q2(*) ! Work array ! !DESCRIPTION: ! ! Perform a two-dimensional fast Fourier transformation. ! ! !REVISION HISTORY: !EOP !----------------------------------------------------------------------- !BOC ! !LOCAL VARIABLES: real s(jp) ! 3-point algebraic filter real ptmp(0:im+1) integer jf(jp) real rsc, bt integer i, j, n, nj if ( igrid == 1 ) then do j=1,jp s(j) = sc(j+js2g0-1) enddo else do j=1,jp s(j) = se(j+js2g0-1) enddo endif nj = 0 ! Number of latitudes to be filtered by zonal FFT do j=1,jp #ifdef ALT_PFT if(s(j) > 1.01) then #else if(s(j) > 1.01 .and. s(j) <= 2.) then rsc = 1./s(j) bt = 0.5*(s(j)-1.) do i=1,im ptmp(i) = p(i,j) enddo ptmp( 0) = ptmp(im) ptmp(im+1) = ptmp(1) do i=1,im p(i,j) = rsc * ( ptmp(i) + bt*(ptmp(i-1)+ptmp(i+1)) ) enddo elseif(s(j) > 2.) then #endif ! Packing for FFT nj = nj + 1 jf(nj) = j do i=1,im q1(i,nj) = p(i,j) enddo q1(im+1,nj) = 0. q1(im+2,nj) = 0. endif enddo if( nj == 0) return call rfftmlt(q1, q2, 1, im+2, im, nj, -1) do n=1,nj do i=5,im+2 q1(i,n) = q1(i,n) * damp(i-2,jf(n)+js2g0-1, igrid) enddo enddo call rfftmlt(q1, q2, 1, im+2, im, nj, 1) do n=1,nj do i=1,im p(i,jf(n)) = q1(i,n) enddo enddo !EOC end subroutine pft2d subroutine pft2d_phys(p, im, jp) ! This is a very weak zonal filter for filtering the ! tendencies from the physics ! !INPUT PARAMETERS: integer, intent(in):: im ! Total X dimension integer, intent(in):: jp ! Total Y dimension ! !INPUT/OUTPUT PARAMETERS: real, intent(inout):: p(im,jp) ! Array to be polar filtered ! !LOCAL VARIABLES: real ptmp(0:im+1) real rsc, bt integer i, j do j=1,jp if(sc(js2g0+j-1) > 1.01 ) then rsc = 1./ min(2., sc(js2g0+j-1)) bt = 0.5*(min(2., sc(js2g0+j-1))-1.) do i=1,im ptmp(i) = p(i,j) enddo ptmp( 0) = ptmp(im) ptmp(im+1) = ptmp(1) do i=1,im p(i,j) = rsc * ( ptmp(i) + bt*(ptmp(i-1)+ptmp(i+1)) ) enddo endif enddo end subroutine pft2d_phys !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !IROUTINE: rfftmlt --- Apply fast Fourier transform ! ! !INTERFACE: subroutine rfftmlt(a,work,inc,jump,n,lot,isign) ! !DESCRIPTION: ! ! Apply the fast Fourier transform. If CPP token SGI_FFT is ! set, SGI libraries will be used. Otherwise the Fortran code ! is inlined. ! ! !REVISION HISTORY: ! ! 99.11.24 Sawyer Added wrappers for SGI ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- !BOC integer jump, inc, n, lot, isign #ifdef SGI_FFT ! Need to link to SGI library with "-lcomplib.sgimath" ! real a(jump,lot) real work(1) ! Not used; here for plug reason ! Local integer*4 iisign,in,iinc,ijump,ilot integer i, j real scale !-----convert to i4 iisign = isign iinc = inc ijump = jump in = n ilot = lot if( iisign < 0 ) then !-----forward call dzfftm1du (iisign,in,ilot,a,iinc,ijump,trigs) elseif( iisign > 0 ) then !-----backward call zdfftm1du (iisign,in,ilot,a,iinc,ijump,trigs) scale = 1.0/float(n) do j=1,lot do i=1,jump a(i,j) = scale*a(i,j) enddo enddo endif #else ! ! Default f77 version ! real a(jump*lot) real work((n+1)*lot) integer nfax, nx, ink, ibase, i, j, k, l, ia, ib, nb, m integer la, nh, jbase,igo ! SUBROUTINE "FFT991" - MULTIPLE REAL/HALF-COMPLEX PERIODIC ! FAST FOURIER TRANSFORM ! SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO ! THAT IN MRFFT2 ! PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM ! IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 ! (1970), 315-337) ! A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA ! WORK IS AN AREA OF SIZE (N+1)*LOT ! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES ! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 ! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' ! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) ! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR ! N IS THE LENGTH OF THE DATA VECTORS ! LOT IS THE NUMBER OF DATA VECTORS ! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT ! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL ! ORDERING OF COEFFICIENTS: ! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) ! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED ! ! ORDERING OF DATA: ! X(0),X(1),X(2),...,X(N-1) ! ! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN ! PARALLEL ! *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER ! ! DEFINITION OF TRANSFORMS: ! ------------------------- ! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) ! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) ! ! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) ! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) nfax=ifax(1) nx=n+1 nh=n/2 ink=inc+inc if (isign==+1) go to 30 ! IF NECESSARY, TRANSFER DATA TO WORK AREA igo=50 if (mod(nfax,2)==1) goto 40 ibase=1 jbase=1 do 20 l=1,lot i=ibase j=jbase !DIR$ IVDEP do 10 m=1,n work(j)=a(i) i=i+inc j=j+1 10 continue ibase=ibase+jump jbase=jbase+nx 20 continue igo=60 go to 40 ! PREPROCESSING (ISIGN=+1) 30 continue call fft99a(a,work,trigs,inc,jump,n,lot) igo=60 ! COMPLEX TRANSFORM 40 continue ia=1 la=1 do 80 k=1,nfax if (igo==60) go to 60 50 continue call vpassm(a(ia),a(ia+inc),work(1),work(2),trigs, & ink,2,jump,nx,lot,nh,ifax(k+1),la) igo=60 go to 70 60 continue call vpassm(work(1),work(2),a(ia),a(ia+inc),trigs, & 2,ink,nx,jump,lot,nh,ifax(k+1),la) igo=50 70 continue la=la*ifax(k+1) 80 continue if (isign==-1) go to 130 ! IF NECESSARY, TRANSFER DATA FROM WORK AREA if (mod(nfax,2)==1) go to 110 ibase=1 jbase=1 do 100 l=1,lot i=ibase j=jbase !DIR$ IVDEP do 90 m=1,n a(j)=work(i) i=i+1 j=j+inc 90 continue ibase=ibase+nx jbase=jbase+jump 100 continue ! FILL IN ZEROS AT END 110 continue ib=n*inc+1 !DIR$ IVDEP do 120 l=1,lot a(ib)=0.0 a(ib+inc)=0.0 ib=ib+jump 120 continue go to 140 ! POSTPROCESSING (ISIGN=-1): 130 continue call fft99b(work,a,trigs,inc,jump,n,lot) 140 continue #endif !EOC end subroutine rfftmlt subroutine fftfax (n, trigs) integer n #ifdef SGI_FFT real trigs(1) ! local integer*4 nn nn=n call dzfftm1dui (nn,trigs) #else real trigs(3*n/2+1) integer:: mode=3 integer i call fax (ifax, n, mode) i = ifax(1) if (ifax(i+1) > 5 .or. n <= 4) ifax(1) = -99 if (ifax(1) <= 0 ) then write(6,*) ' set99 -- invalid n' stop 'set99' endif call fftrig (trigs, n, mode) #endif end subroutine fftfax ! ! Add Lima code here. ! #else use ecmfft !BOP ! ! !PUBLIC MEMBER FUNCTIONS: public pft2d, pft2d_phys, pft_init, rfftmlt, fftfax ! ! !DESCRIPTION: ! ! This module provides fast-Fourier transforms ! ! \begin{tabular}{|l|l|} \hline \hline ! pft2d & \\ \hline ! pft\_cf & \\ \hline ! rfftmlt & \\ \hline ! fftfax & \\ \hline ! \hline ! \end{tabular} ! ! !REVISION HISTORY: ! 01.01.30 Lin Integrated into this module ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- CONTAINS !----------------------------------------------------------------------- !BOP ! !IROUTINE: pft2d --- Two-dimensional fast Fourier transform ! ! !INTERFACE: subroutine pft2d(p, s, damp, im, jp, ifax, trigs, q1, q2) implicit none ! !USES: ! !INPUT PARAMETERS: integer im ! Total X dimension integer jp ! Total Y dimension integer ifax(13) ! FFT internal information real s(jp) ! 3-point algebraic filter real damp(im,jp) ! FFT damping coefficients real trigs(3*im/2+1) ! !INPUT/OUTPUT PARAMETERS: real q1( im+2, *) ! Work array real q2(*) ! Work array real p(im,jp) ! Array to be polar filtered ! !DESCRIPTION: ! ! Perform a two-dimensional fast Fourier transformation. ! ! !REVISION HISTORY: ! 01.01.30 Lin Put into this module ! 01.03.26 Sawyer Added ProTeX documentation ! 02.04.05 Sawyer Integrated newest FVGCM version ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: real rsc, bt integer i, j, n, nj !Local Auto arrays: real ptmp(0:im+1) !!! real q1( im+2, jp) !!! real q2( (im+1)*jp ) integer jf(jp) nj = 0 do 200 j=1,jp #ifdef ALT_PFT if(s(j) > 1.01) then #else if(s(j) > 1.01 .and. s(j) <= 2.) then rsc = 1./s(j) bt = 0.5*(s(j)-1.) do i=1,im ptmp(i) = p(i,j) enddo ptmp( 0) = p(im,j) ptmp(im+1) = p(1 ,j) do i=1,im p(i,j) = rsc * ( ptmp(i) + bt*(ptmp(i-1)+ptmp(i+1)) ) enddo elseif(s(j) > 4.) then #endif ! Packing for FFT nj = nj + 1 jf(nj) = j do i=1,im q1(i,nj) = p(i,j) enddo q1(im+1,nj) = 0. q1(im+2,nj) = 0. endif 200 continue if( nj == 0) return call rfftmlt(q1, q2, trigs, ifax, 1, im+2, im, nj, -1) do n=1,nj do i=5,im+2 q1(i,n) = q1(i,n) * damp(i-2,jf(n)) enddo enddo call rfftmlt(q1, q2, trigs, ifax, 1, im+2, im, nj, 1) do n=1,nj do i=1,im p(i,jf(n)) = q1(i,n) enddo enddo !EOC end subroutine pft2d subroutine pft2d_phys(p, s, im, jp) implicit none ! This is a very weak zonal filter for filtering the ! tendencies from the physics ! !INPUT PARAMETERS: integer, intent(in):: im ! Total X dimension integer, intent(in):: jp ! Total Y dimension real, intent(in):: s(jp) ! 3-point algebraic filter ! !INPUT/OUTPUT PARAMETERS: real, intent(inout):: p(im,jp) ! Array to be polar filtered ! !LOCAL VARIABLES: real ptmp(0:im+1) real rsc, bt integer i, j do j=1,jp if(s(j) > 1.01 ) then rsc = 1./ min(2., s(j)) bt = 0.5*(min(2.,s(j))-1.) do i=1,im ptmp(i) = p(i,j) enddo ptmp( 0) = p(im,j) ptmp(im+1) = p(1 ,j) do i=1,im p(i,j) = rsc * ( ptmp(i) + bt*(ptmp(i-1)+ptmp(i+1)) ) enddo endif enddo end subroutine pft2d_phys !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !IROUTINE: pft_init --- Calculate algebraic and FFT polar filters ! ! !INTERFACE: subroutine pft_init(im, jm, js2g0, jn1g1, sc, se, dc, de, & cosp, cose, ycrit, ipft1, ipft2) implicit none ! !USES: ! !INPUT PARAMETERS: integer im ! Total X dimension integer jm ! Total Y dimension integer js2g0 ! j south limit ghosted 0 (SP: from 2) integer jn1g1 ! j north limit ghosted 1 (starts jm) integer ipft1, ipft2 real cosp(jm) ! cosine array real cose(jm) ! cosine array real ycrit ! critical value ! !OUTPUT PARAMETERS: real sc(ipft1:ipft2) ! Algebric filter at center real se(js2g0:jn1g1) ! Algebric filter at edge real dc(im,ipft1:ipft2) ! FFT filter at center real de(im,js2g0:jn1g1) ! FFT filter at edge ! !DESCRIPTION: ! ! Compute coefficients for the 3-point algebraic and the FFT ! polar filters. ! ! !REVISION HISTORY: ! ! 99.01.01 Lin Creation ! 99.08.20 Sawyer/Lin Changes for SPMD mode ! 01.01.30 Lin Put into this module ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- !BOC ! ! !LOCAL VARIABLES: integer i, j real dl, coszc, cutoff, phi, damp real pi pi = 4.d0 * datan(1.d0) coszc = cos(ycrit*pi/180.) ! INIT fft polar coefficients: dl = pi/dble(im) cutoff = 1.e-20 do j=ipft1, ipft2 do i=1,im dc(i,j) = 1. enddo enddo do j=js2g0,jn1g1 do i=1,im de(i,j) = 1. enddo enddo ! write(6,*) 'Initializing 3-point polar filter coefficients:' !************ ! Cell center !************ do j=ipft1, ipft2 sc(j) = (coszc/cosp(j))**2 #ifdef ALT_PFT if(sc(j) > 1.) then #else if(sc(j) > 1. .and. sc(j) <= 2.0) then sc(j) = 1. + (sc(j)-1.)/(sc(j)+1.) elseif(sc(j) > 2. .and. sc(j) <= 4.) then sc(j) = 1. + sc(j)/(8.-sc(j)) sc(j) = min(2., sc(j)) elseif(sc(j) > 4. ) then #endif ! FFT filter do i=1,im/2 phi = dl * i damp = min((cosp(j)/coszc)/sin(phi), 1.)**2 if(damp < cutoff) damp = 0. dc(2*i-1,j) = damp dc(2*i ,j) = damp enddo endif enddo !************ ! Cell edges !************ do j=js2g0,jn1g1 se(j) = (coszc/cose(j))**2 #ifdef ALT_PFT if( se(j) > 1. ) then #else if(se(j) > 1. .and. se(j) <= 2.0 ) then se(j) = 1. + (se(j)-1.)/(se(j)+1.) elseif(se(j) > 2. .and. se(j) <= 4.) then se(j) = 1. + se(j)/(8.-se(j)) se(j) = min(2., se(j)) elseif(se(j) > 4. ) then #endif ! FFT do i=1,im/2 phi = dl * i damp = min((cose(j)/coszc)/sin(phi), 1.)**2 if(damp < cutoff) damp = 0. de(2*i-1,j) = damp de(2*i ,j) = damp enddo endif enddo !EOC end subroutine pft_init !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !IROUTINE: rfftmlt --- Apply fast Fourier transform ! ! !INTERFACE: subroutine rfftmlt(a,work,trigs,ifax,inc,jump,n,lot,isign) ! !DESCRIPTION: ! ! Apply the fast Fourier transform. If CPP token SGI_FFT is ! set, SGI libraries will be used. Otherwise the Fortran code ! is inlined. ! ! !REVISION HISTORY: ! ! 99.11.24 Sawyer Added wrappers for SGI ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- !BOC #ifdef SGI_FFT ! ! WS 99.11.24 : Added SGI wrappers ! integer inc, jump, n, lot, isign real a(jump,lot) real trigs(1) real work(1) ! Not used; here for plug reason integer ifax(*) ! Local integer*4 iisign,in,iinc,ijump,ilot integer i, j real scale !-----convert to i4 iisign = isign iinc = inc ijump = jump in = n ilot = lot if( iisign < 0 ) then !-----forward call dzfftm1du (iisign,in,ilot,a,iinc,ijump,trigs) endif if( iisign > 0 ) then !-----backward call zdfftm1du (iisign,in,ilot,a,iinc,ijump,trigs) scale = 1.0/float(n) do j=1,lot do i=1,jump a(i,j) = scale*a(i,j) enddo enddo endif #else ! ! Default f77 version ! real a(jump*lot) real work((n+1)*lot) real trigs(3*n/2+1) integer ifax(13) ! SUBROUTINE "FFT991" - MULTIPLE REAL/HALF-COMPLEX PERIODIC ! FAST FOURIER TRANSFORM ! SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO ! THAT IN MRFFT2 ! PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM ! IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12 ! (1970), 315-337) ! A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA ! WORK IS AN AREA OF SIZE (N+1)*LOT ! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES ! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2 ! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' ! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) ! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR ! N IS THE LENGTH OF THE DATA VECTORS ! LOT IS THE NUMBER OF DATA VECTORS ! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT ! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL ! ORDERING OF COEFFICIENTS: ! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) ! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED ! ! ORDERING OF DATA: ! X(0),X(1),X(2),...,X(N-1) ! ! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN ! PARALLEL ! *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER ! ! DEFINITION OF TRANSFORMS: ! ------------------------- ! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) ! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) ! ! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) ! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) nfax=ifax(1) nx=n+1 nh=n/2 ink=inc+inc if (isign==+1) go to 30 ! IF NECESSARY, TRANSFER DATA TO WORK AREA igo=50 if (mod(nfax,2)==1) goto 40 ibase=1 jbase=1 do 20 l=1,lot i=ibase j=jbase !DIR$ IVDEP do 10 m=1,n work(j)=a(i) i=i+inc j=j+1 10 continue ibase=ibase+jump jbase=jbase+nx 20 continue igo=60 go to 40 ! PREPROCESSING (ISIGN=+1) 30 continue call fft99a(a,work,trigs,inc,jump,n,lot) igo=60 ! COMPLEX TRANSFORM 40 continue ia=1 la=1 do 80 k=1,nfax if (igo==60) go to 60 50 continue call vpassm(a(ia),a(ia+inc),work(1),work(2),trigs, & ink,2,jump,nx,lot,nh,ifax(k+1),la) igo=60 go to 70 60 continue call vpassm(work(1),work(2),a(ia),a(ia+inc),trigs, & 2,ink,nx,jump,lot,nh,ifax(k+1),la) igo=50 70 continue la=la*ifax(k+1) 80 continue if (isign==-1) go to 130 ! IF NECESSARY, TRANSFER DATA FROM WORK AREA if (mod(nfax,2)==1) go to 110 ibase=1 jbase=1 do 100 l=1,lot i=ibase j=jbase !DIR$ IVDEP do 90 m=1,n a(j)=work(i) i=i+1 j=j+inc 90 continue ibase=ibase+nx jbase=jbase+jump 100 continue ! FILL IN ZEROS AT END 110 continue ib=n*inc+1 !DIR$ IVDEP do 120 l=1,lot a(ib)=0.0 a(ib+inc)=0.0 ib=ib+jump 120 continue go to 140 ! POSTPROCESSING (ISIGN=-1): 130 continue call fft99b(work,a,trigs,inc,jump,n,lot) 140 continue #endif !EOC end subroutine rfftmlt !----------------------------------------------------------------------- !----------------------------------------------------------------------- !BOP ! !IROUTINE: fftfax --- Initialize FFT ! ! !INTERFACE: subroutine fftfax (n, ifax, trigs) ! !USES: ! !DESCRIPTION: ! ! Initialize the fast Fourier transform. If CPP token SGI_FFT is ! set, SGI libraries will be used. Otherwise the Fortran code ! is inlined. ! ! !REVISION HISTORY: ! ! 99.11.24 Sawyer Added wrappers for SGI ! 01.03.26 Sawyer Added ProTeX documentation ! !EOP !----------------------------------------------------------------------- !BOC #ifdef SGI_FFT real trigs(1) integer ifax(*) integer n ! local integer*4 nn nn=n call dzfftm1dui (nn,trigs) #else integer ifax(13) real trigs(3*n/2+1) ! MODE 3 IS USED FOR REAL/HALF-COMPLEX TRANSFORMS. IT IS POSSIBLE ! TO DO COMPLEX/COMPLEX TRANSFORMS WITH OTHER VALUES OF MODE, BUT ! DOCUMENTATION OF THE DETAILS WERE NOT AVAILABLE WHEN THIS ROUTINE ! WAS WRITTEN. data mode /3/ call fax (ifax, n, mode) i = ifax(1) if (ifax(i+1) > 5 .or. n <= 4) ifax(1) = -99 if (ifax(1) <= 0 ) then write(6,*) ' set99 -- invalid n' stop 'set99' endif call fftrig (trigs, n, mode) #endif !EOC end subroutine fftfax !----------------------------------------------------------------------- #endif end module pft_module ! End of Lima pft_module