NCEPLIBS-sp  2.3.3
sptranf.f
Go to the documentation of this file.
1 C> @file
2 C>
3 C> Perform a scalar spherical transform
4 C> @author IREDELL @date 96-02-29
5 
6 C> This subprogram performs a spherical transform between spectral
7 C> coefficients of scalar quantities and fields on a global
8 C> cylindrical grid. The wave-space can be either triangular or
9 C> rhomboidal. The grid-space can be either an equally-spaced grid
10 C> (with or without pole points) or a gaussian grid.
11 C> The wave and grid fields may have general indexing,
12 C> but each wave field is in sequential 'ibm order',
13 C> i.e. with zonal wavenumber as the slower index.
14 C> Transforms are done in latitude pairs for efficiency;
15 C> thus grid arrays for each hemisphere must be passed.
16 C> If so requested, just a subset of the latitude pairs
17 C> may be transformed in each invocation of the subprogram.
18 C> The transforms are all multiprocessed over latitude except
19 C> the transform from fourier to spectral is multiprocessed
20 C> over zonal wavenumber to ensure reproducibility.
21 C> Transform several fields at a time to improve vectorization.
22 C> Subprogram can be called from a multiprocessing environment.
23 C>
24 C> PROGRAM HISTORY LOG:
25 C> - 96-02-29 IREDELL
26 C> - 1998-12-15 IREDELL GENERIC FFT USED, OPENMP DIRECTIVES INSERTED
27 C> - 2013-01-16 IREDELL MIRVIS FIXING AFFT NEGATIVE SHARING EFFECT DURING
28 C> OMP LOOPS BY CREATING TMP AFFT COPY (AFFT_TMP)
29 C> TO BE PRIVATE DURING OMP LOOP THREADING
30 C>
31 C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE
32 C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
33 C> @param MAXWV - INTEGER SPECTRAL TRUNCATION
34 C> @param IDRT - INTEGER GRID IDENTIFIER
35 C> (IDRT=4 FOR GAUSSIAN GRID,
36 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
37 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
38 C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES.
39 C> @param JMAX - INTEGER NUMBER OF LATITUDES.
40 C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM.
41 C> @param IP - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN
42 C> @param IS - INTEGER SKIP NUMBER BETWEEN LONGITUDES
43 C> @param JN - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH
44 C> @param JS - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH
45 C> @param KW - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS
46 C> @param KG - INTEGER SKIP NUMBER BETWEEN GRID FIELDS
47 C> @param JB - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM
48 C> @param JE - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM
49 C> @param JC - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS
50 C> @param[out] WAVE - REAL (*) WAVE FIELDS IF IDIR>0
51 C> @param[out] GRIDN - REAL (*) N.H. GRID FIELDS (STARTING AT JB) IF IDIR<0
52 C> @param[out] GRIDS - REAL (*) S.H. GRID FIELDS (STARTING AT JB) IF IDIR<0
53 C> @param IDIR - INTEGER TRANSFORM FLAG
54 C> (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE)
55 C>
56 C> SUBPROGRAMS CALLED:
57 C> - sptranf0() sptranf() spectral initialization
58 C> - sptranf1() sptranf() spectral transform
59 C>
60 C> Minimum grid dimensions for unaliased transforms to spectral:
61 C> DIMENSION |LINEAR |QUADRATIC
62 C> ----------------------- |--------- |-------------
63 C> IMAX | 2*MAXWV+2 | 3*MAXWV/2*2+2
64 C> JMAX (IDRT=4,IROMB=0) | 1*MAXWV+1 | 3*MAXWV/2+1
65 C> JMAX (IDRT=4,IROMB=1) | 2*MAXWV+1 | 5*MAXWV/2+1
66 C> JMAX (IDRT=0,IROMB=0) | 2*MAXWV+3 | 3*MAXWV/2*2+3
67 C> JMAX (IDRT=0,IROMB=1) | 4*MAXWV+3 | 5*MAXWV/2*2+3
68 C> JMAX (IDRT=256,IROMB=0) | 2*MAXWV+1 | 3*MAXWV/2*2+1
69 C> JMAX (IDRT=256,IROMB=1) | 4*MAXWV+1 | 5*MAXWV/2*2+1
70 C>
71  SUBROUTINE sptranf(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
72  & IP,IS,JN,JS,KW,KG,JB,JE,JC,
73  & WAVE,GRIDN,GRIDS,IDIR)
74 
75  REAL WAVE(*),GRIDN(*),GRIDS(*)
76  REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
77  REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
78  REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
79  REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
80  REAL(8) AFFT(50000+4*IMAX), AFFT_TMP(50000+4*IMAX)
81  REAL CLAT(JB:JE),SLAT(JB:JE),WLAT(JB:JE)
82  REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2,JB:JE)
83  REAL PLNTOP(MAXWV+1,JB:JE)
84  REAL WTOP(2*(MAXWV+1))
85  REAL G(IMAX,2)
86 ! write(0,*) 'sptranf top'
87 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88 C SET PARAMETERS
89  mp=0
90  CALL sptranf0(iromb,maxwv,idrt,imax,jmax,jb,je,
91  & eps,epstop,enn1,elonn1,eon,eontop,
92  & afft,clat,slat,wlat,pln,plntop)
93 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94 C TRANSFORM WAVE TO GRID
95  IF(idir.GT.0) THEN
96 C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,WTOP,G,IJKN,IJKS)
97  DO k=1,kmax
98  afft_tmp=afft
99  kws=(k-1)*kw
100  wtop=0
101  DO j=jb,je
102  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
103  & eps,epstop,enn1,elonn1,eon,eontop,
104  & afft_tmp,clat(j),slat(j),wlat(j),
105  & pln(1,j),plntop(1,j),mp,
106  & wave(kws+1),wtop,g,idir)
107  IF(ip.EQ.1.AND.is.EQ.1) THEN
108  DO i=1,imax
109  ijkn=i+(j-jb)*jn+(k-1)*kg
110  ijks=i+(j-jb)*js+(k-1)*kg
111  gridn(ijkn)=g(i,1)
112  grids(ijks)=g(i,2)
113  ENDDO
114  ELSE
115  DO i=1,imax
116  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
117  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
118  gridn(ijkn)=g(i,1)
119  grids(ijks)=g(i,2)
120  ENDDO
121  ENDIF
122  ENDDO
123  ENDDO
124 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125 C TRANSFORM GRID TO WAVE
126  ELSE
127 C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,WTOP,G,IJKN,IJKS)
128  DO k=1,kmax
129  afft_tmp=afft
130  kws=(k-1)*kw
131  wtop=0
132  DO j=jb,je
133  IF(wlat(j).GT.0.) THEN
134  IF(ip.EQ.1.AND.is.EQ.1) THEN
135  DO i=1,imax
136  ijkn=i+(j-jb)*jn+(k-1)*kg
137  ijks=i+(j-jb)*js+(k-1)*kg
138  g(i,1)=gridn(ijkn)
139  g(i,2)=grids(ijks)
140  ENDDO
141  ELSE
142  DO i=1,imax
143  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
144  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
145  g(i,1)=gridn(ijkn)
146  g(i,2)=grids(ijks)
147  ENDDO
148  ENDIF
149  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
150  & eps,epstop,enn1,elonn1,eon,eontop,
151  & afft_tmp,clat(j),slat(j),wlat(j),
152  & pln(1,j),plntop(1,j),mp,
153  & wave(kws+1),wtop,g,idir)
154  ENDIF
155  ENDDO
156  ENDDO
157  ENDIF
158 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159  END
sptranf0
subroutine sptranf0(IROMB, MAXWV, IDRT, IMAX, JMAX, JB, JE, EPS, EPSTOP, ENN1, ELONN1, EON, EONTOP, AFFT, CLAT, SLAT, WLAT, PLN, PLNTOP)
This subprogram performs an initialization for subprogram sptranf().
Definition: sptranf0.f:44
sptranf
subroutine sptranf(IROMB, MAXWV, IDRT, IMAX, JMAX, KMAX, IP, IS, JN, JS, KW, KG, JB, JE, JC, WAVE, GRIDN, GRIDS, IDIR)
This subprogram performs a spherical transform between spectral coefficients of scalar quantities and...
Definition: sptranf.f:74
sptranf1
subroutine sptranf1(IROMB, MAXWV, IDRT, IMAX, JMAX, JB, JE, EPS, EPSTOP, ENN1, ELONN1, EON, EONTOP, AFFT, CLAT, SLAT, WLAT, PLN, PLNTOP, MP, W, WTOP, G, IDIR)
THIS SUBPROGRAM PERFORMS AN SINGLE LATITUDE TRANSFORM FOR SUBPROGRAM SPTRANF.
Definition: sptranf1.f:48