diff --git a/.github/workflows/developer.yml b/.github/workflows/developer.yml index bbfc776..4d3119d 100644 --- a/.github/workflows/developer.yml +++ b/.github/workflows/developer.yml @@ -39,6 +39,12 @@ jobs: cmake -DOPENMP=ON -DCMAKE_PREFIX_PATH="~/" -DENABLE_DOCS=ON -DCMAKE_Fortran_FLAGS="-g -fprofile-arcs -ftest-coverage -O0 -fsanitize=address" -DCMAKE_BUILD_TYPE=Debug .. make -j2 VERBOSE=1 + - uses: actions/upload-artifact@v2 + with: + name: docs + path: | + sp/build/docs/html + - name: test-sp run: | cd sp/build diff --git a/src/fftpack.F b/src/fftpack.F index 4f67e83..bb415a6 100644 --- a/src/fftpack.F +++ b/src/fftpack.F @@ -1,5 +1,15 @@ C> @file -C> @brief fftpack. +C> @brief A concatination of the (FFTPACK)[https://netlib.org/fftpack/] library code. +C> +C> FFTPACK is a package of Fortran subprograms for the fast Fourier +C> transform of periodic and other symmetric sequences. It includes +C> complex, real, sine, cosine, and quarter-wave transforms. +C> +C>Reference: +C>- P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations +C>(G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83. +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO C> dcrft C> @@ -18,6 +28,7 @@ C> @param n2 C> @param z C> @param nz +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale, & table,n1,wrk,n2,z,nz) @@ -61,6 +72,9 @@ SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale, C> @param n2 C> @param z C> @param nz +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO + SUBROUTINE scrft(init,x,ldx,y,ldy,n,m,isign,scale, & table,n1,wrk,n2,z,nz) @@ -97,6 +111,8 @@ SUBROUTINE scrft(init,x,ldx,y,ldy,n,m,isign,scale, C> @param table C> @param work C> @param isys +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys) implicit none @@ -137,6 +153,8 @@ SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys) C> @param n2 C> @param z C> @param nz +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale, & table,n1,wrk,n2,z,nz) @@ -184,6 +202,8 @@ SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale, C> @param n2 C> @param z C> @param nz +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE srcft(init,x,ldx,y,ldy,n,m,isign,scale, & table,n1,wrk,n2,z,nz) @@ -224,6 +244,8 @@ SUBROUTINE srcft(init,x,ldx,y,ldy,n,m,isign,scale, C> @param table C> @param work C> @param isys +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys) implicit none @@ -255,6 +277,8 @@ SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys) C> @param N C> @param R C> @param WSAVE +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RFFTF (N,R,WSAVE) DIMENSION R(1) ,WSAVE(1) IF (N .EQ. 1) RETURN @@ -267,6 +291,8 @@ SUBROUTINE RFFTF (N,R,WSAVE) C> @param N C> @param R C> @param WSAVE +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RFFTB (N,R,WSAVE) DIMENSION R(1) ,WSAVE(1) IF (N .EQ. 1) RETURN @@ -278,6 +304,8 @@ SUBROUTINE RFFTB (N,R,WSAVE) C> C> @param N C> @param WSAVE +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RFFTI (N,WSAVE) DIMENSION WSAVE(1) IF (N .EQ. 1) RETURN @@ -292,6 +320,8 @@ SUBROUTINE RFFTI (N,WSAVE) C> @param CH C> @param WA C> @param IFAC +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC) DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(*) NF = IFAC(2) @@ -359,6 +389,8 @@ SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC) C> @param CH C> @param WA C> @param IFAC +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC) DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(*) NF = IFAC(2) @@ -424,6 +456,8 @@ SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC) C> @param N C> @param WA C> @param IFAC +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RFFTI1 (N,WA,IFAC) DIMENSION WA(1) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ @@ -491,6 +525,8 @@ SUBROUTINE RFFTI1 (N,WA,IFAC) C> @param CC C> @param CH C> @param WA1 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(1) @@ -528,6 +564,8 @@ SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1) C> @param CH C> @param WA1 C> @param WA2 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(1) ,WA2(1) @@ -576,6 +614,8 @@ SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2) C> @param WA1 C> @param WA2 C> @param WA3 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(1) ,WA2(1) ,WA3(1) @@ -645,6 +685,8 @@ SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3) C> @param WA2 C> @param WA3 C> @param WA4 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) @@ -721,6 +763,8 @@ SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) C> @param CH C> @param CH2 C> @param WA +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), @@ -897,6 +941,8 @@ SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) C> @param CC C> @param CH C> @param WA1 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1) DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) , 1 WA1(1) @@ -933,6 +979,8 @@ SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1) C> @param CH C> @param WA1 C> @param WA2 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2) DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) , 1 WA1(1) ,WA2(1) @@ -978,6 +1026,8 @@ SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2) C> @param WA1 C> @param WA2 C> @param WA3 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3) DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) , 1 WA1(1) ,WA2(1) ,WA3(1) @@ -1043,6 +1093,8 @@ SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3) C> @param WA2 C> @param WA3 C> @param WA4 +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) , 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1) @@ -1115,6 +1167,8 @@ SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4) C> @param CH C> @param CH2 C> @param WA +C> +C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), diff --git a/src/spanaly.f b/src/spanaly.f index cb8a025..417524f 100644 --- a/src/spanaly.f +++ b/src/spanaly.f @@ -1,5 +1,5 @@ C> @file -C> Analyze spectral from fourier +C> @brief Analyze spectral from Fourier. C> C> ### Program History Log C> Date | Programmer | Comments @@ -10,25 +10,26 @@ C> C> @author IREDELL @date 91-10-31 -C> Analyzes spectral coefficients from fourier coefficients -C> for a latitude pair (northern and southern hemispheres). +C> Analyzes spectral coefficients from Fourier coefficients +C> for a latitude pair (Northern and Southern hemispheres). +C> C> Vector components are multiplied by cosine of latitude. C> -C> @param I SPECTRAL DOMAIN SHAPE (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param M SPECTRAL TRUNCATION -C> @param IM EVEN NUMBER OF FOURIER COEFFICIENTS -C> @param IX DIMENSION OF FOURIER COEFFICIENTS (IX>=IM+2) -C> @param NC DIMENSION OF SPECTRAL COEFFICIENTS (NC>=(M+1)*((I+1)*M+2)) -C> @param NCTOP DIMENSION OF SPECTRAL COEFFICIENTS OVER TOP (NCTOP>=2*(M+1)) -C> @param KM NUMBER OF FIELDS -C> @param WGT GAUSSIAN WEIGHT -C> @param CLAT COSINE OF LATITUDE -C> @param PLN ((M+1)*((I+1)*M+2)/2) LEGENDRE POLYNOMIALS -C> @param PLNTOP (M+1) LEGENDRE POLYNOMIAL OVER TOP -C> @param MP (KM) IDENTIFIERS (0 FOR SCALAR, 1 FOR VECTOR) -C> @param F (IX,2,KM) FOURIER COEFFICIENTS COMBINED -C> @param SPC (NC,KM) SPECTRAL COEFFICIENTS -C> @param SPCTOP (NCTOP,KM) SPECTRAL COEFFICIENTS OVER TOP +C> @param I spectral domain shape (0 for triangular, 1 for rhomboidal) +C> @param M spectral truncation +C> @param IM even number of Fourier coefficients +C> @param IX dimension of Fourier coefficients (IX>=IM+2) +C> @param NC dimension of spectral coefficients (NC>=(M+1)*((I+1)*M+2)) +C> @param NCTOP dimension of spectral coefficients over top (NCTOP>=2*(M+1)) +C> @param KM number of fields +C> @param WGT Gaussian weight +C> @param CLAT cosine of latitude +C> @param PLN Legendre polynomials +C> @param PLNTOP Legendre polynomial over top +C> @param MP identifiers (0 for scalar, 1 for vector) +C> @param F Fourier coefficients combined +C> @param SPC spectral coefficients +C> @param SPCTOP spectral coefficients over top C> C> @author IREDELL @date 91-10-31 SUBROUTINE SPANALY(I,M,IM,IX,NC,NCTOP,KM,WGT,CLAT,PLN,PLNTOP,MP, diff --git a/src/spfft.f b/src/spfft.f index de5af98..672bb02 100644 --- a/src/spfft.f +++ b/src/spfft.f @@ -1,44 +1,42 @@ C> @file -C> -C> Perform multiple fast fourier transforms. -C> @author IREDELL ORG: W/NMC23 @date 96-02-20 +C> @brief Perform multiple fast fourier transforms. +C> @author Iredell @date 96-02-20 C> This subprogram performs multiple fast fourier transforms C> between complex amplitudes in fourier space and real values C> in cyclic physical space. +C> C> Subprogram spfft must be invoked first with idir=0 C> to initialize trigonemetric data. Use subprogram spfft1 C> to perform an fft without previous initialization. C> This version invokes the ibm essl fft. C> -C> the restrictions on imax are that it must be a multiple -C> of 1 to 25 factors of two, up to 2 factors of three, -C> and up to 1 factor of five, seven and eleven. +C> The restrictions on imax are that it must be a multiple +C> of 1 to 25 factors of two, up to 2 factors of three, +C> and up to 1 factor of five, seven and eleven. C> -C> If IDIR=0, then W and G need not contain any valid data. -C> the other parameters must be supplied and cannot change -C> in succeeding calls until the next time it is called with IDIR=0. +C> If IDIR=0, then W and G need not contain any valid data. +C> the other parameters must be supplied and cannot change +C> in succeeding calls until the next time it is called with IDIR=0. C> -C> This subprogram is not thread-safe when IDIR=0. On the other hand, -C> when IDIR is not zero, it can be called from a threaded region. +C> This subprogram is not thread-safe when IDIR=0. On the other hand, +C> when IDIR is not zero, it can be called from a threaded region. C> -C> @param IMAX - INTEGER NUMBER OF VALUES IN THE CYCLIC PHYSICAL SPACE -C> (SEE LIMITATIONS ON IMAX IN REMARKS BELOW.) -C> @param INCW - INTEGER FIRST DIMENSION OF THE COMPLEX AMPLITUDE ARRAY -C> (INCW >= IMAX/2+1) -C> @param INCG - INTEGER FIRST DIMENSION OF THE REAL VALUE ARRAY -C> (INCG >= IMAX) -C> @param KMAX - INTEGER NUMBER OF TRANSFORMS TO PERFORM -C> @param[out] W - COMPLEX(INCW,KMAX) COMPLEX AMPLITUDES IF IDIR>0 -C> @param[out] G - REAL(INCG,KMAX) REAL VALUES IF IDIR<0 -C> @param IDIR - INTEGER DIRECTION FLAG -C> IDIR=0 TO INITIALIZE INTERNAL TRIGONOMETRIC DATA -C> IDIR>0 TO TRANSFORM FROM FOURIER TO PHYSICAL SPACE -C> IDIR<0 TO TRANSFORM FROM PHYSICAL TO FOURIER SPACE +C> @param IMAX number of values in the cyclic physical space +C> (see limitations on imax in remarks below.) +C> @param INCW first dimension of the complex amplitude array +C> (INCW >= IMAX/2+1) +C> @param INCG first dimension of the real value array +C> (INCG >= IMAX) +C> @param KMAX number of transforms to perform +C> @param[out] W complex amplitudes if IDIR>0 +C> @param[out] G real values if IDIR<0 +C> @param IDIR direction flag +C> - IDIR=0 to initialize internal trigonometric data +C> - IDIR>0 TO transform from Fourier to physical space +C> - IDIR<0 TO transform from physical to fourier space C> -C> CALLED: -C> - SCRFT() IBM ESSL COMPLEX TO REAL FOURIER TRANSFORM -C> - SRCFT() IBM ESSL REAL TO COMPLEX FOURIER TRANSFORM +C> @author Iredell @date 96-02-20 SUBROUTINE SPFFT(IMAX,INCW,INCG,KMAX,W,G,IDIR) IMPLICIT NONE @@ -49,9 +47,9 @@ SUBROUTINE SPFFT(IMAX,INCW,INCG,KMAX,W,G,IDIR) REAL,SAVE,ALLOCATABLE:: AUX1CR(:),AUX1RC(:) INTEGER:: NAUX2 REAL:: AUX2(20000+INT(0.57*IMAX)) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + NAUX2=20000+INT(0.57*IMAX) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C INITIALIZATION. C ALLOCATE AND FILL AUXILIARY ARRAYS WITH TRIGONOMETRIC DATA SELECT CASE(IDIR) @@ -63,12 +61,12 @@ SUBROUTINE SPFFT(IMAX,INCW,INCG,KMAX,W,G,IDIR) & AUX1CR,NAUX1,AUX2,NAUX2,0.,0) CALL SRCFT(1,G,INCG,W,INCW,IMAX,KMAX,+1,1./IMAX, & AUX1RC,NAUX1,AUX2,NAUX2,0.,0) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C FOURIER TO PHYSICAL TRANSFORM. CASE(1:) CALL SCRFT(0,W,INCW,G,INCG,IMAX,KMAX,-1,1., & AUX1CR,NAUX1,AUX2,NAUX2,0.,0) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C PHYSICAL TO FOURIER TRANSFORM. CASE(:-1) CALL SRCFT(0,G,INCG,W,INCW,IMAX,KMAX,+1,1./IMAX, diff --git a/src/spgradq.f b/src/spgradq.f index e4f1f00..c667822 100644 --- a/src/spgradq.f +++ b/src/spgradq.f @@ -1,7 +1,6 @@ C> @file -C> -C> Compute gradient in spectral space. -C> @author IREDELL @date 92-10-31 +C> @brief Compute gradient in spectral space. +C> @author Iredell @date 92-10-31 C> Computes the horizontal vector gradient of a scalar field c> in spectral space. @@ -19,16 +18,18 @@ C> Advantage is taken of the fact that eps(l,l)=0 c> in order to vectorize over the entire spectral domain. C> -C> @param I SPECTRAL DOMAIN SHAPE (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param M SPECTRAL TRUNCATION -C> @param ENN1 ((M+1)*((I+1)*M+2)/2) N*(N+1)/A**2 -C> @param ELONN1 ((M+1)*((I+1)*M+2)/2) L/(N*(N+1))*A -C> @param EON ((M+1)*((I+1)*M+2)/2) EPSILON/N*A -C> @param EONTOP (M+1) EPSILON/N*A OVER TOP -C> @param Q ((M+1)*((I+1)*M+2)) SCALAR FIELD -C> @param QDX ((M+1)*((I+1)*M+2)) ZONAL GRADIENT (TIMES COSLAT) -C> @param QDY ((M+1)*((I+1)*M+2)) MERID GRADIENT (TIMES COSLAT) -C> @param QDYTOP (2*(M+1)) MERID GRADIENT (TIMES COSLAT) OVER TOP +C> @param I spectral domain shape (0 for triangular, 1 for rhomboidal) +C> @param M spectral truncation +C> @param ENN1 +C> @param ELONN1 +C> @param EON EPSILON/N*A +C> @param EONTOP EPSILON/N*A over top +C> @param Q scalar field +C> @param QDX zonal gradient (times coslat) +C> @param QDY merid gradient (times coslat) +C> @param QDYTOP merid gradient (times coslat) over top +C> +C> @author IREDELL @date 92-10-31 SUBROUTINE SPGRADQ(I,M,ENN1,ELONN1,EON,EONTOP,Q,QDX,QDY,QDYTOP) REAL ENN1((M+1)*((I+1)*M+2)/2),ELONN1((M+1)*((I+1)*M+2)/2) @@ -36,7 +37,7 @@ SUBROUTINE SPGRADQ(I,M,ENN1,ELONN1,EON,EONTOP,Q,QDX,QDY,QDYTOP) REAL Q((M+1)*((I+1)*M+2)) REAL QDX((M+1)*((I+1)*M+2)),QDY((M+1)*((I+1)*M+2)) REAL QDYTOP(2*(M+1)) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C TAKE ZONAL AND MERIDIONAL GRADIENTS K=1 QDX(2*K-1)=0. @@ -54,13 +55,12 @@ SUBROUTINE SPGRADQ(I,M,ENN1,ELONN1,EON,EONTOP,Q,QDX,QDY,QDYTOP) QDX(2*K)=ELONN1(K)*ENN1(K)*Q(2*K-1) QDY(2*K-1)=-EON(K)*ENN1(K-1)*Q(2*K-3) QDY(2*K)=-EON(K)*ENN1(K-1)*Q(2*K-2) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C TAKE MERIDIONAL GRADIENT OVER TOP DO L=0,M K=L*(2*M+(I-1)*(L-1))/2+I*L+M+1 QDYTOP(2*L+1)=-EONTOP(L+1)*ENN1(K)*Q(2*K-1) QDYTOP(2*L+2)=-EONTOP(L+1)*ENN1(K)*Q(2*K) ENDDO -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - RETURN END diff --git a/src/sptgptv.f b/src/sptgptv.f index b0726d2..37cf3b0 100644 --- a/src/sptgptv.f +++ b/src/sptgptv.f @@ -1,51 +1,53 @@ C> @file +C> @brief Transform spectral vector to station points. +C> +C> ### Program History Log +C> Date | Programmer | Comments +C> -----|------------|--------- +C> 96-02-29 | IREDELL | Initial +C> 1998-12-15 | IREDELL | Openmp directives inserted +C> 1999-08-18 | IREDELL | Openmp directive typo fixed +C> 2003-06-30 | IREDELL | use spfftpt() C> -C> Transform spectral vector to station points C> @author IREDELL @date 96-02-29 -C> THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM -C> FROM SPECTRAL COEFFICIENTS OF DIVERGENCES AND CURLS -C> TO SPECIFIED SETS OF STATION POINT VECTORS ON THE GLOBE. -C> THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL. -C> THE WAVE AND POINT FIELDS MAY HAVE GENERAL INDEXING, -C> BUT EACH WAVE FIELD IS IN SEQUENTIAL 'IBM ORDER', -C> I.E. WITH ZONAL WAVENUMBER AS THE SLOWER INDEX. -C> THE TRANSFORMS ARE ALL MULTIPROCESSED OVER STATIONS. -C> TRANSFORM SEVERAL FIELDS AT A TIME TO IMPROVE VECTORIZATION. -C> SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. +C> This subprogram performs a spherical transform +C> from spectral coefficients of divergences and curls +C> to specified sets of station point vectors on the globe. +C> +C> The wave-space can be either triangular or rhomboidal. +C> +C> The wave and point fields may have general indexing, +C> but each wave field is in sequential 'IBM order', +C> i.e. with zonal wavenumber as the slower index. +C> +C> The transforms are all multiprocessed over stations. +C> +C> Transform several fields at a time to improve vectorization. C> -C> PROGRAM HISTORY LOG: -C> - 96-02-29 IREDELL -C> - 1998-12-15 IREDELL OPENMP DIRECTIVES INSERTED -C> - 1999-08-18 IREDELL OPENMP DIRECTIVE TYPO FIXED -C> - 2003-06-30 IREDELL USE SPFFTPT +C> Subprogram can be called from a multiprocessing environment. C> -C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE -C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param MAXWV - INTEGER SPECTRAL TRUNCATION -C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM. -C> @param NMAX - INTEGER NUMBER OF STATION POINTS TO RETURN -C> @param KWSKIP - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS -C> (DEFAULTS TO (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0) -C> @param KGSKIP - INTEGER SKIP NUMBER BETWEEN STATION POINT SETS -C> (DEFAULTS TO NMAX IF KGSKIP=0) -C> @param NRSKIP - INTEGER SKIP NUMBER BETWEEN STATION LATS AND LONS -C> (DEFAULTS TO 1 IF NRSKIP=0) -C> @param NGSKIP - INTEGER SKIP NUMBER BETWEEN STATION POINTS -C> (DEFAULTS TO 1 IF NGSKIP=0) -C> @param RLAT - REAL (*) STATION LATITUDES IN DEGREES -C> @param RLON - REAL (*) STATION LONGITUDES IN DEGREES -C> @param WAVED - REAL (*) WAVE DIVERGENCE FIELDS -C> @param WAVEZ - REAL (*) WAVE VORTICITY FIELDS -C> @param UP - REAL (*) STATION POINT U-WIND SETS -C> @param VP - REAL (*) STATION POINT V-WIND SETS +C> @param IROMB spectral domain shape +c> (0 for triangular, 1 for rhomboidal) +C> @param MAXWV spectral truncation +C> @param KMAX number of fields to transform. +C> @param NMAX number of station points to return +C> @param KWSKIP skip number between wave fields +c> (defaults to (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0) +C> @param KGSKIP skip number between station point sets +c> (defaults to NMAX IF KGSKIP=0) +C> @param NRSKIP skip number between station lats and lons +c> (defaults to 1 if NRSKIP=0) +C> @param NGSKIP skip number between station points +c> (defaults to 1 if NGSKIP=0) +C> @param RLAT station latitudes in degrees +C> @param RLON station longitudes in degrees +C> @param WAVED wave divergence fields +C> @param WAVEZ wave vorticity fields +C> @param UP station point u-wind sets +C> @param VP station point v-wind sets C> -C> SUBPROGRAMS CALLED: -C> - SPWGET GET WAVE-SPACE CONSTANTS -C> - SPLEGEND COMPUTE LEGENDRE POLYNOMIALS -C> - SPSYNTH SYNTHESIZE FOURIER FROM SPECTRAL -C> - SPDZ2UV COMPUTE WINDS FROM DIVERGENCE AND VORTICITY -C> - SPFFTPT POINTWISE FOURIER TRANSFORM +C> @author IREDELL @date 96-02-29 SUBROUTINE SPTGPTV(IROMB,MAXWV,KMAX,NMAX, & KWSKIP,KGSKIP,NRSKIP,NGSKIP, & RLAT,RLON,WAVED,WAVEZ,UP,VP) @@ -62,7 +64,7 @@ SUBROUTINE SPTGPTV(IROMB,MAXWV,KMAX,NMAX, REAL F(2*MAXWV+3,2,2*KMAX) REAL G(2*KMAX) PARAMETER(PI=3.14159265358979) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C CALCULATE PRELIMINARY CONSTANTS CALL SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP) MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 @@ -78,7 +80,7 @@ SUBROUTINE SPTGPTV(IROMB,MAXWV,KMAX,NMAX, IF(NR.EQ.0) NR=1 IF(NG.EQ.0) NG=1 MP=1 -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C CALCULATE SPECTRAL WINDS C$OMP PARALLEL DO PRIVATE(KWS) DO K=1,KMAX @@ -87,7 +89,7 @@ SUBROUTINE SPTGPTV(IROMB,MAXWV,KMAX,NMAX, & WAVED(KWS+1),WAVEZ(KWS+1), & W(1,K),W(1,KMAX+K),WTOP(1,K),WTOP(1,KMAX+K)) ENDDO -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C CALCULATE STATION FIELDS C$OMP PARALLEL DO PRIVATE(KU,KV,RADLAT,SLAT1,CLAT1) C$OMP& PRIVATE(PLN,PLNTOP,F,G,NK) @@ -116,5 +118,5 @@ SUBROUTINE SPTGPTV(IROMB,MAXWV,KMAX,NMAX, VP(NK)=G(KV) ENDDO ENDDO -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + END diff --git a/src/sptran.f b/src/sptran.f index 4171486..1d52a45 100644 --- a/src/sptran.f +++ b/src/sptran.f @@ -1,88 +1,93 @@ C> @file +C> @brief Perform a scalar spherical transform. +C> +C> ### Program History Log +C> Date | Programmer | Comments +C> -----|------------|--------- +C> 96-02-29 | IREDELL | Initial +C> 1998-12-15 | IREDELL | Generic fft used, openmp directives inserted C> -C> Perform a scalar spherical transform C> @author IREDELL @date 96-02-29 C> This subprogram performs a spherical transform between spectral C> coefficients of scalar quantities and fields on a global -C> cylindrical grid. The wave-space can be either triangular or -C> rhomboidal. The grid-space can be either an equally-spaced grid -C> (with or without pole points) or a gaussian grid. +C> cylindrical grid. +C> +C> The wave-space can be either triangular or +C> rhomboidal. +C> +C> The grid-space can be either an equally-spaced grid +C> (with or without pole points) or a Gaussian grid. +C> C> The wave and grid fields may have general indexing, -C> but each wave field is in sequential 'ibm order', +C> but each wave field is in sequential 'IBM order', C> i.e. with zonal wavenumber as the slower index. +C> C> Transforms are done in latitude pairs for efficiency; C> thus grid arrays for each hemisphere must be passed. C> If so requested, just a subset of the latitude pairs C> may be transformed in each invocation of the subprogram. +C> C> The transforms are all multiprocessed over latitude except -C> the transform from fourier to spectral is multiprocessed +C> the transform from Fourier to spectral is multiprocessed C> over zonal wavenumber to ensure reproducibility. +C> C> Transform several fields at a time to improve vectorization. C> Subprogram can be called from a multiprocessing environment. C> -C> PROGRAM HISTORY LOG: -C> - 96-02-29 IREDELL -C> - 1998-12-15 IREDELL GENERIC FFT USED, OPENMP DIRECTIVES INSERTED -C> -C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE -C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param MAXWV - INTEGER SPECTRAL TRUNCATION -C> @param IDRT - INTEGER GRID IDENTIFIER -C> (IDRT=4 FOR GAUSSIAN GRID, -C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES, -C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES) -C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES. -C> @param JMAX - INTEGER NUMBER OF LATITUDES. -C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM. -C> @param IPRIME - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN. -C> (DEFAULTS TO 1 IF IPRIME=0) -C> @param ISKIP - INTEGER SKIP NUMBER BETWEEN LONGITUDES -C> (DEFAULTS TO 1 IF ISKIP=0) -C> @param JNSKIP - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH -C> (DEFAULTS TO IMAX IF JNSKIP=0) -C> @param JSSKIP - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH -C> (DEFAULTS TO -IMAX IF JSSKIP=0) -C> @param KWSKIP - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS -C> (DEFAULTS TO (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0) -C> @param KGSKIP - INTEGER SKIP NUMBER BETWEEN GRID FIELDS -C> (DEFAULTS TO IMAX*JMAX IF KGSKIP=0) -C> @param JBEG - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM -C> (DEFAULTS TO 1 IF JBEG=0) -C> (IF JBEG=0 AND IDIR<0, WAVE IS ZEROED BEFORE TRANSFORM) -C> @param JEND - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM -C> (DEFAULTS TO (JMAX+1)/2 IF JEND=0) -C> @param JCPU - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS -C> @param[out] WAVE - REAL (*) WAVE FIELDS IF IDIR>0 -C> @param[out] GRIDN - REAL (*) N.H. GRID FIELDS (STARTING AT JBEG) IF IDIR<0 -C> @param[out] GRIDS - REAL (*) S.H. GRID FIELDS (STARTING AT JBEG) IF IDIR<0 -C> @param IDIR - INTEGER TRANSFORM FLAG -C> (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE) -C> -C> SUBPROGRAMS CALLED: -C> - sptranf() Perform a scalar spherical transform -C> C> Minimum grid dimensions for unaliased transforms to spectral: -C> DIMENSION |LINEAR |QUADRATIC -C> ----------------------- |--------- |------------- -C> IMAX | 2*MAXWV+2 | 3*MAXWV/2*2+2 -C> JMAX (IDRT=4,IROMB=0) | 1*MAXWV+1 | 3*MAXWV/2+1 -C> JMAX (IDRT=4,IROMB=1) | 2*MAXWV+1 | 5*MAXWV/2+1 -C> JMAX (IDRT=0,IROMB=0) | 2*MAXWV+3 | 3*MAXWV/2*2+3 -C> JMAX (IDRT=0,IROMB=1) | 4*MAXWV+3 | 5*MAXWV/2*2+3 -C> JMAX (IDRT=256,IROMB=0) | 2*MAXWV+1 | 3*MAXWV/2*2+1 -C> JMAX (IDRT=256,IROMB=1) | 4*MAXWV+1 | 5*MAXWV/2*2+1 +C> DIMENSION |LINEAR |QUADRATIC +C> ----------------------- |--------- |------------- +C> IMAX | 2*MAXWV+2 | 3*MAXWV/2*2+2 +C> JMAX (IDRT=4,IROMB=0) | 1*MAXWV+1 | 3*MAXWV/2+1 +C> JMAX (IDRT=4,IROMB=1) | 2*MAXWV+1 | 5*MAXWV/2+1 +C> JMAX (IDRT=0,IROMB=0) | 2*MAXWV+3 | 3*MAXWV/2*2+3 +C> JMAX (IDRT=0,IROMB=1) | 4*MAXWV+3 | 5*MAXWV/2*2+3 +C> JMAX (IDRT=256,IROMB=0) | 2*MAXWV+1 | 3*MAXWV/2*2+1 +C> JMAX (IDRT=256,IROMB=1) | 4*MAXWV+1 | 5*MAXWV/2*2+1 C> +C> @param IROMB spectral domain shape +c> (0 for triangular, 1 for rhomboidal) +C> @param MAXWV spectral truncation +C> @param IDRT grid identifier +C> - IDRT=4 for Gaussian grid, +C> - IDRT=0 for equally-spaced grid including poles, +C> - IDRT=256 for equally-spaced grid excluding poles +C> @param IMAX even number of longitudes. +C> @param JMAX number of latitudes. +C> @param KMAX number of fields to transform. +C> @param IPRIME longitude index for the prime meridian. +C> (defaults to 1 if IPRIME=0) +C> @param ISKIP skip number between longitudes +C> (defaults to 1 if ISKIP=0) +C> @param JNSKIP skip number between n.h. latitudes from north +C> (defaults to imax if JNSKIP=0) +C> @param JSSKIP skip number between s.h. latitudes from south +c> (defaults to -imax if JSSKIP=0) +C> @param KWSKIP skip number between wave fields +c> (defaults to (MAXWV+1)*((IROMB+1)*MAXWV+2) IF KWSKIP=0) +C> @param KGSKIP skip number between grid fields +c> (defaults to IMAX*JMAX IF KGSKIP=0) +C> @param JBEG latitude index (from pole) to begin transform +c> (defaults to 1 if JBEG=0) +C> (if JBEG=0 and IDIR<0, wave is zeroed before transform) +C> @param JEND latitude index (from pole) to end transform +c> (defaults to (JMAX+1)/2 IF JEND=0) +C> @param JCPU number of cpus over which to multiprocess +C> @param[out] WAVE wave fields if IDIR>0 +C> @param[out] gridn n.h. grid fields (starting at jbeg) if IDIR<0 +C> @param[out] grids s.h. grid fields (starting at jbeg) if IDIR<0 +C> @param IDIR transform flag +C> (idir>0 for wave to grid, idir<0 for grid to wave) +C> +C> @author IREDELL @date 96-02-29 SUBROUTINE SPTRAN(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, & IPRIME,ISKIP,JNSKIP,JSSKIP,KWSKIP,KGSKIP, & JBEG,JEND,JCPU, & WAVE,GRIDN,GRIDS,IDIR) REAL WAVE(*),GRIDN(*),GRIDS(*) -!==EM== integer JC, JCPU -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!==EM== -! print *, "jjjjjjjcccccccccc from SPTRAN- before NCPUS", JC + MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 IP=IPRIME IS=ISKIP @@ -102,21 +107,16 @@ SUBROUTINE SPTRAN(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, IF(JB.EQ.0) JB=1 IF(JE.EQ.0) JE=(JMAX+1)/2 IF(JC.EQ.0) JC=NCPUS() -! JC=NCPUS() -!==EM== -! print *, "jjjjjjjjjjjjjjjjjjjjcccccccccc from SPTRAN", JC -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + IF(IDIR.LT.0.AND.JBEG.EQ.0) THEN DO K=1,KMAX KWS=(K-1)*KW WAVE(KWS+1:KWS+2*MX)=0 ENDDO ENDIF -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -!==EM== -! print *, "jjjjjjjcccccccccc from SPTRANF- before NCPUS", JC + CALL SPTRANF(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, & IP,IS,JN,JS,KW,KG,JB,JE,JC, & WAVE,GRIDN,GRIDS,IDIR) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + END diff --git a/src/sptranf.f b/src/sptranf.f index 25eb081..99219f2 100644 --- a/src/sptranf.f +++ b/src/sptranf.f @@ -1,73 +1,76 @@ C> @file +C> @brief Perform a scalar spherical transform C> -C> Perform a scalar spherical transform -C> @author IREDELL @date 96-02-29 +C> ### Program History Log +C> Date | Programmer | Comments +C> -----|------------|--------- +C> 96-02-29 | Iredell | Initial. +C> 1998-12-15 | Iredell | Generic fft used, openmp directives inserted +C> 2013-01-16 | Iredell, Mirvis | Fixing afft negative sharing effect +C> +C> @author Iredell @date 96-02-29 C> This subprogram performs a spherical transform between spectral C> coefficients of scalar quantities and fields on a global -C> cylindrical grid. The wave-space can be either triangular or +C> cylindrical grid. +C> +C> The wave-space can be either triangular or C> rhomboidal. The grid-space can be either an equally-spaced grid C> (with or without pole points) or a gaussian grid. +C> C> The wave and grid fields may have general indexing, C> but each wave field is in sequential 'ibm order', C> i.e. with zonal wavenumber as the slower index. +C> C> Transforms are done in latitude pairs for efficiency; C> thus grid arrays for each hemisphere must be passed. +C> C> If so requested, just a subset of the latitude pairs C> may be transformed in each invocation of the subprogram. C> The transforms are all multiprocessed over latitude except C> the transform from fourier to spectral is multiprocessed C> over zonal wavenumber to ensure reproducibility. +C> C> Transform several fields at a time to improve vectorization. C> Subprogram can be called from a multiprocessing environment. C> -C> PROGRAM HISTORY LOG: -C> - 96-02-29 IREDELL -C> - 1998-12-15 IREDELL GENERIC FFT USED, OPENMP DIRECTIVES INSERTED -C> - 2013-01-16 IREDELL MIRVIS FIXING AFFT NEGATIVE SHARING EFFECT DURING -C> OMP LOOPS BY CREATING TMP AFFT COPY (AFFT_TMP) -C> TO BE PRIVATE DURING OMP LOOP THREADING -C> -C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE -C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param MAXWV - INTEGER SPECTRAL TRUNCATION -C> @param IDRT - INTEGER GRID IDENTIFIER -C> (IDRT=4 FOR GAUSSIAN GRID, -C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES, -C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES) -C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES. -C> @param JMAX - INTEGER NUMBER OF LATITUDES. -C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM. -C> @param IP - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN -C> @param IS - INTEGER SKIP NUMBER BETWEEN LONGITUDES -C> @param JN - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH -C> @param JS - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH -C> @param KW - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS -C> @param KG - INTEGER SKIP NUMBER BETWEEN GRID FIELDS -C> @param JB - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM -C> @param JE - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM -C> @param JC - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS -C> @param[out] WAVE - REAL (*) WAVE FIELDS IF IDIR>0 -C> @param[out] GRIDN - REAL (*) N.H. GRID FIELDS (STARTING AT JB) IF IDIR<0 -C> @param[out] GRIDS - REAL (*) S.H. GRID FIELDS (STARTING AT JB) IF IDIR<0 -C> @param IDIR - INTEGER TRANSFORM FLAG -C> (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE) -C> -C> SUBPROGRAMS CALLED: -C> - sptranf0() sptranf() spectral initialization -C> - sptranf1() sptranf() spectral transform -C> C> Minimum grid dimensions for unaliased transforms to spectral: -C> DIMENSION |LINEAR |QUADRATIC -C> ----------------------- |--------- |------------- -C> IMAX | 2*MAXWV+2 | 3*MAXWV/2*2+2 -C> JMAX (IDRT=4,IROMB=0) | 1*MAXWV+1 | 3*MAXWV/2+1 -C> JMAX (IDRT=4,IROMB=1) | 2*MAXWV+1 | 5*MAXWV/2+1 -C> JMAX (IDRT=0,IROMB=0) | 2*MAXWV+3 | 3*MAXWV/2*2+3 -C> JMAX (IDRT=0,IROMB=1) | 4*MAXWV+3 | 5*MAXWV/2*2+3 -C> JMAX (IDRT=256,IROMB=0) | 2*MAXWV+1 | 3*MAXWV/2*2+1 -C> JMAX (IDRT=256,IROMB=1) | 4*MAXWV+1 | 5*MAXWV/2*2+1 +C> DIMENSION |LINEAR |QUADRATIC +C> ----------------------- |--------- |------------- +C> IMAX | 2*MAXWV+2 | 3*MAXWV/2*2+2 +C> JMAX (IDRT=4,IROMB=0) | 1*MAXWV+1 | 3*MAXWV/2+1 +C> JMAX (IDRT=4,IROMB=1) | 2*MAXWV+1 | 5*MAXWV/2+1 +C> JMAX (IDRT=0,IROMB=0) | 2*MAXWV+3 | 3*MAXWV/2*2+3 +C> JMAX (IDRT=0,IROMB=1) | 4*MAXWV+3 | 5*MAXWV/2*2+3 +C> JMAX (IDRT=256,IROMB=0) | 2*MAXWV+1 | 3*MAXWV/2*2+1 +C> JMAX (IDRT=256,IROMB=1) | 4*MAXWV+1 | 5*MAXWV/2*2+1 +C> +C> @param IROMB spectral domain shape +c> (0 for triangular, 1 for rhomboidal) +C> @param MAXWV spectral truncation +C> @param IDRT grid identifier +C> - IDRT=4 for Gaussian grid, +C> - IDRT=0 for equally-spaced grid including poles +C> - IDRT=256 for equally-spaced grid excluding poles +C> @param IMAX even number of longitudes. +C> @param JMAX number of latitudes. +C> @param KMAX number of fields to transform. +C> @param IP longitude index for the prime meridian +C> @param IS skip number between longitudes +C> @param JN skip number between n.h. latitudes from north +C> @param JS skip number between s.h. latitudes from south +C> @param KW skip number between wave fields +C> @param KG skip number between grid fields +C> @param JB latitude index (from pole) to begin transform +C> @param JE latitude index (from pole) to end transform +C> @param JC number of cpus over which to multiprocess +C> @param[out] WAVE wave fields if IDIR>0 +C> @param[out] GRIDN n.h. grid fields (starting at JB) if IDIR<0 +C> @param[out] GRIDS s.h. grid fields (starting at JB) if IDIR<0 +C> @param IDIR transform flag +C> (IDIR>0 for wave to grid, IDIR<0 for grid to wave) C> +C> @author Iredell @date 96-02-29 SUBROUTINE SPTRANF(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, & IP,IS,JN,JS,KW,KG,JB,JE,JC, & WAVE,GRIDN,GRIDS,IDIR) @@ -84,7 +87,7 @@ SUBROUTINE SPTRANF(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, REAL WTOP(2*(MAXWV+1)) REAL G(IMAX,2) ! write(0,*) 'sptranf top' -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C SET PARAMETERS MP=0 CALL SPTRANF0(IROMB,MAXWV,IDRT,IMAX,JMAX,JB,JE, @@ -121,7 +124,7 @@ SUBROUTINE SPTRANF(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, ENDIF ENDDO ENDDO -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C TRANSFORM GRID TO WAVE ELSE C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,WTOP,G,IJKN,IJKS) @@ -155,5 +158,4 @@ SUBROUTINE SPTRANF(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, ENDDO ENDDO ENDIF -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END diff --git a/src/sptranf0.f b/src/sptranf0.f index 51cbac2..4feb90d 100644 --- a/src/sptranf0.f +++ b/src/sptranf0.f @@ -1,43 +1,36 @@ C> @file -C> -C> Sptranf spectral initialization +C> @brief Sptranf spectral initialization. C> @author IREDELL @date 96-02-29 C> This subprogram performs an initialization for C> subprogram sptranf(). Use this subprogram outside -C> the sptranf family context at your own risk. -C> -C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE -C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param MAXWV - INTEGER SPECTRAL TRUNCATION -C> @param IDRT - INTEGER GRID IDENTIFIER -C> (IDRT=4 FOR GAUSSIAN GRID, -C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES, -C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES) -C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES -C> @param JMAX - INTEGER NUMBER OF LATITUDES -C> @param JB - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM -C> @param JE - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM -C> @param EPS - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param EPSTOP - REAL (MAXWV+1) -C> @param ENN1 - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param ELONN1 - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param EON - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param EONTOP - REAL (MAXWV+1) -C> @param AFFT - REAL(8) (50000+4*IMAX) AUXILIARY ARRAY IF IDIR=0 -C> @param CLAT - REAL (JB:JE) COSINES OF LATITUDE -C> @param SLAT - REAL (JB:JE) SINES OF LATITUDE -C> @param WLAT - REAL (JB:JE) GAUSSIAN WEIGHTS -C> @param PLN - REAL ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2,JB:JE) -C> LEGENDRE POLYNOMIALS -C> @param PLNTOP - REAL (MAXWV+1,JB:JE) LEGENDRE POLYNOMIAL OVER TOP +C> the sptranf() family context at your own risk. C> -C> SUBPROGRAMS CALLED: -C> - spwget() GET WAVE-SPACE CONSTANTS -C> - spffte() PERFORM FAST FOURIER TRANSFORM -C> - splat() COMPUTE LATITUDE FUNCTIONS -C> - splegend() COMPUTE LEGENDRE POLYNOMIALS +C> @param IROMB spectral domain shape +c> (0 for triangular, 1 for rhomboidal) +C> @param MAXWV spectral truncation +C> @param IDRT grid identifier +C> - IDRT=4 for Gaussian grid, +C> - IDRT=0 for equally-spaced grid including poles, +C> - IDRT=256 for equally-spaced grid excluding poles +C> @param IMAX even number of longitudes +C> @param JMAX number of latitudes +C> @param JB latitude index (from pole) to begin transform +C> @param JE latitude index (from pole) to end transform +C> @param EPS +C> @param EPSTOP +C> @param ENN1 +C> @param ELONN1 +C> @param EON +C> @param EONTOP +C> @param AFFT auxiliary array if IDIR=0 +C> @param CLAT cosines of latitude +C> @param SLAT sines of latitude +C> @param WLAT Gaussian weights +C> @param PLN Legendre polynomials +C> @param PLNTOP Legendre polynomial over top C> +C> @author IREDELL @date 96-02-29 SUBROUTINE SPTRANF0(IROMB,MAXWV,IDRT,IMAX,JMAX,JB,JE, & EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP, & AFFT,CLAT,SLAT,WLAT,PLN,PLNTOP) @@ -51,7 +44,7 @@ SUBROUTINE SPTRANF0(IROMB,MAXWV,IDRT,IMAX,JMAX,JB,JE, REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2,JB:JE) REAL PLNTOP(MAXWV+1,JB:JE) REAL SLATX(JMAX),WLATX(JMAX) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + CALL SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP) CALL SPFFTE(IMAX,(IMAX+2)/2,IMAX,2,0.,0.,0,AFFT) CALL SPLAT(IDRT,JMAX,SLATX,WLATX) @@ -67,5 +60,5 @@ SUBROUTINE SPTRANF0(IROMB,MAXWV,IDRT,IMAX,JMAX,JB,JE, CALL SPLEGEND(IROMB,MAXWV,SLAT(J),CLAT(J),EPS,EPSTOP, & PLN(1,J),PLNTOP(1,J)) ENDDO -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + END diff --git a/src/sptranfv.f b/src/sptranfv.f index 7025fdc..addc132 100644 --- a/src/sptranfv.f +++ b/src/sptranfv.f @@ -1,80 +1,82 @@ C> @file +C> @brief Perform a vector spherical transform C> -C> Perform a vector spherical transform -C> @author IREDELL @date 96-02-29 +C> ### Program History Log +C> Date | Programmer | Comments +C> -----|------------|--------- +C> 96-02-29 | Iredell | Initial. +C> 1998-12-15 | Iredell | Generic fft used, openmp directives inserted +C> 2013-01-16 | Iredell & MIRVIS | Fixing afft negative sharing effect during omp loops +C> +C> @author Iredell @date 96-02-29 -C> THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM -C> BETWEEN SPECTRAL COEFFICIENTS OF DIVERGENCES AND CURLS -C> AND VECTOR FIELDS ON A GLOBAL CYLINDRICAL GRID. -C> THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL. -C> THE GRID-SPACE CAN BE EITHER AN EQUALLY-SPACED GRID -C> (WITH OR WITHOUT POLE POINTS) OR A GAUSSIAN GRID. -C> THE WAVE AND GRID FIELDS MAY HAVE GENERAL INDEXING, -C> BUT EACH WAVE FIELD IS IN SEQUENTIAL 'IBM ORDER', -C> I.E. WITH ZONAL WAVENUMBER AS THE SLOWER INDEX. -C> TRANSFORMS ARE DONE IN LATITUDE PAIRS FOR EFFICIENCY; -C> THUS GRID ARRAYS FOR EACH HEMISPHERE MUST BE PASSED. -C> IF SO REQUESTED, JUST A SUBSET OF THE LATITUDE PAIRS -C> MAY BE TRANSFORMED IN EACH INVOCATION OF THE SUBPROGRAM. -C> THE TRANSFORMS ARE ALL MULTIPROCESSED OVER LATITUDE EXCEPT -C> THE TRANSFORM FROM FOURIER TO SPECTRAL IS MULTIPROCESSED -C> OVER ZONAL WAVENUMBER TO ENSURE REPRODUCIBILITY. -C> TRANSFORM SEVERAL FIELDS AT A TIME TO IMPROVE VECTORIZATION. -C> SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT. +C> This subprogram performs a spherical transform +c> between spectral coefficients of divergences and curls +c> and vector fields on a global cylindrical grid. +C> +C> The wave-space can be either triangular or rhomboidal. +C> +C> The grid-space can be either an equally-spaced grid +c> (with or without pole points) or a gaussian grid. +C> +C> The wave and grid fields may have general indexing, +c> but each wave field is in sequential 'ibm order', +c> i.e. with zonal wavenumber as the slower index. +C> +C> Transforms are done in latitude pairs for efficiency; +c> thus grid arrays for each hemisphere must be passed. +c> if so requested, just a subset of the latitude pairs +c> may be transformed in each invocation of the subprogram. +C> +C> The transforms are all multiprocessed over latitude except +c> the transform from fourier to spectral is multiprocessed +c> over zonal wavenumber to ensure reproducibility. C> -C> PROGRAM HISTORY LOG: -C> - 96-02-29 IREDELL -C> - 1998-12-15 IREDELL GENERIC FFT USED, OPENMP DIRECTIVES INSERTED -C> - 2013-01-16 IREDELL & MIRVIS FIXING AFFT NEGATIVE SHARING EFFECT DURING -C> OMP LOOPS BY CREATING TMP AFFT COPY (AFFT_TMP) -C> TO BE PRIVATE DURING OMP LOOP THREADING +C> Transform several fields at a time to improve vectorization. +c> subprogram can be called from a multiprocessing environment. C> -C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE -C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) -C> @param MAXWV - INTEGER SPECTRAL TRUNCATION -C> @param IDRT - INTEGER GRID IDENTIFIER -C> (IDRT=4 FOR GAUSSIAN GRID, -C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES, -C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES) -C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES. -C> @param JMAX - INTEGER NUMBER OF LATITUDES. -C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM. -C> @param IP - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN -C> @param IS - INTEGER SKIP NUMBER BETWEEN LONGITUDES -C> @param JN - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH -C> @param JS - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH -C> @param KW - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS -C> @param KG - INTEGER SKIP NUMBER BETWEEN GRID FIELDS -C> @param JB - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM -C> @param JE - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM -C> @param JC - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS -C> @param[out] WAVED - REAL (*) WAVE DIVERGENCE FIELDS IF IDIR>0 +C> Minimum grid dimensions for unaliased transforms to spectral: +C> DIMENSION |LINEAR |QUADRATIC +C> ----------------------- |--------- |------------- +C> IMAX |2*MAXWV+2 |3*MAXWV/2*2+2 +C> JMAX (IDRT=4,IROMB=0) |1*MAXWV+1 |3*MAXWV/2+1 +C> JMAX (IDRT=4,IROMB=1) |2*MAXWV+1 |5*MAXWV/2+1 +C> JMAX (IDRT=0,IROMB=0) |2*MAXWV+3 |3*MAXWV/2*2+3 +C> JMAX (IDRT=0,IROMB=1) |4*MAXWV+3 |5*MAXWV/2*2+3 +C> JMAX (IDRT=256,IROMB=0) |2*MAXWV+1 |3*MAXWV/2*2+1 +C> JMAX (IDRT=256,IROMB=1) |4*MAXWV+1 |5*MAXWV/2*2+1 +C> +C> @param IROMB spectral domain shape +c> (0 for triangular, 1 for rhomboidal) +C> @param MAXWV spectral truncation +C> @param IDRT grid identifier +C> - IDRT=4 for gaussian grid +C> - IDRT=0 for equally-spaced grid including poles +C> - IDRT=256 for equally-spaced grid excluding poles +C> @param IMAX even number of longitudes. +C> @param JMAX number of latitudes. +C> @param KMAX number of fields to transform. +C> @param IP longitude index for the prime meridian +C> @param IS skip number between longitudes +C> @param JN skip number between n.h. latitudes from north +C> @param JS skip number between s.h. latitudes from south +C> @param KW skip number between wave fields +C> @param KG skip number between grid fields +C> @param JB latitude index (from pole) to begin transform +C> @param JE latitude index (from pole) to end transform +C> @param JC number of cpus over which to multiprocess +C> @param[out] WAVED wave divergence fields if IDIR>0 C> [WAVED=(D(GRIDU)/DLAM+D(CLAT*GRIDV)/DPHI)/(CLAT*RERTH)] -C> @param[out] WAVEZ - REAL (*) WAVE VORTICITY FIELDS IF IDIR>0 +C> @param[out] WAVEZ wave vorticity fields if IDIR>0 C> [WAVEZ=(D(GRIDV)/DLAM-D(CLAT*GRIDU)/DPHI)/(CLAT*RERTH)] -C> @param[out] GRIDUN - REAL (*) N.H. GRID U-WINDS (STARTING AT JB) IF IDIR<0 -C> @param[out] GRIDUS - REAL (*) S.H. GRID U-WINDS (STARTING AT JB) IF IDIR<0 -C> @param[out] GRIDVN - REAL (*) N.H. GRID V-WINDS (STARTING AT JB) IF IDIR<0 -C> @param[out] GRIDVS - REAL (*) S.H. GRID V-WINDS (STARTING AT JB) IF IDIR<0 -C> @param IDIR - INTEGER TRANSFORM FLAG -C> (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE) -C> -C> SUBPROGRAMS CALLED: -C> - SPTRANF0 SPTRANF SPECTRAL INITIALIZATION -C> - SPTRANF1 SPTRANF SPECTRAL TRANSFORM -C> - SPDZ2UV COMPUTE WINDS FROM DIVERGENCE AND VORTICITY -C> - SPUV2DZ COMPUTE DIVERGENCE AND VORTICITY FROM WINDS +C> @param[out] GRIDUN N.H. grid u-winds (starting at jb) if IDIR<0 +C> @param[out] GRIDUS S.H. grid u-winds (starting at jb) if IDIR<0 +C> @param[out] GRIDVN N.H. grid v-winds (starting at jb) if IDIR<0 +C> @param[out] GRIDVS S.H. grid v-winds (starting at jb) if IDIR<0 +C> @param IDIR transform flag +C> (IDIR>0 for wave to grid, IDIR<0 for grid to wave). C> -C> REMARKS: MINIMUM GRID DIMENSIONS FOR UNALIASED TRANSFORMS TO SPECTRAL: -C> DIMENSION |LINEAR |QUADRATIC -C> ----------------------- |--------- |------------- -C> IMAX |2*MAXWV+2 |3*MAXWV/2*2+2 -C> JMAX (IDRT=4,IROMB=0) |1*MAXWV+1 |3*MAXWV/2+1 -C> JMAX (IDRT=4,IROMB=1) |2*MAXWV+1 |5*MAXWV/2+1 -C> JMAX (IDRT=0,IROMB=0) |2*MAXWV+3 |3*MAXWV/2*2+3 -C> JMAX (IDRT=0,IROMB=1) |4*MAXWV+3 |5*MAXWV/2*2+3 -C> JMAX (IDRT=256,IROMB=0) |2*MAXWV+1 |3*MAXWV/2*2+1 -C> JMAX (IDRT=256,IROMB=1) |4*MAXWV+1 |5*MAXWV/2*2+1 +C> @author Iredell @date 96-02-29 SUBROUTINE SPTRANFV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, & IP,IS,JN,JS,KW,KG,JB,JE,JC, & WAVED,WAVEZ,GRIDUN,GRIDUS,GRIDVN,GRIDVS,IDIR) @@ -93,14 +95,14 @@ SUBROUTINE SPTRANFV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, REAL WTOP(2*(MAXWV+1),2) REAL G(IMAX,2,2) REAL WINC((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C SET PARAMETERS MX=(MAXWV+1)*((IROMB+1)*MAXWV+2)/2 MP=1 CALL SPTRANF0(IROMB,MAXWV,IDRT,IMAX,JMAX,JB,JE, & EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP, & AFFT,CLAT,SLAT,WLAT,PLN,PLNTOP) -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C TRANSFORM WAVE TO GRID IF(IDIR.GT.0) THEN C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS) @@ -142,7 +144,7 @@ SUBROUTINE SPTRANFV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, ENDIF ENDDO ENDDO -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + C TRANSFORM GRID TO WAVE ELSE C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS,WINC) @@ -191,5 +193,4 @@ SUBROUTINE SPTRANFV(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX, WAVEZ(KWS+1:KWS+2*MX)=WAVEZ(KWS+1:KWS+2*MX)+WINC(1:2*MX,2) ENDDO ENDIF -C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END diff --git a/src/spwget.f b/src/spwget.f index 7990c72..c7641ac 100644 --- a/src/spwget.f +++ b/src/spwget.f @@ -1,18 +1,17 @@ C> @file -C> -C> Get wave-space constants +C> @brief Get wave-space constants. C> @author IREDELL @date 96-02-29 C> This subprogram gets wave-space constants. C> C> @param IROMB SPECTRAL DOMAIN SHAPE (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL) C> @param MAXWV SPECTRAL TRUNCATION -C> @param EPS ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param EPSTOP (MAXWV+1) -C> @param ENN1 ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param ELONN1 ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param EON ((MAXWV+1)*((IROMB+1)*MAXWV+2)/2) -C> @param EONTOP (MAXWV+1) +C> @param EPS +C> @param EPSTOP +C> @param ENN1 +C> @param ELONN1 +C> @param EON +C> @param EONTOP C> C> @author IREDELL @date 96-02-29 SUBROUTINE SPWGET(IROMB,MAXWV,EPS,EPSTOP,ENN1,ELONN1,EON,EONTOP)