NCEPLIBS-sp  2.5.0
sptranfv.f
Go to the documentation of this file.
1 C> @file
2 C> @brief Perform a vector spherical transform
3 C>
4 C> ### Program History Log
5 C> Date | Programmer | Comments
6 C> -----|------------|---------
7 C> 96-02-29 | Iredell | Initial.
8 C> 1998-12-15 | Iredell | Generic fft used, openmp directives inserted
9 C> 2013-01-16 | Iredell & MIRVIS | Fixing afft negative sharing effect during omp loops
10 C>
11 C> @author Iredell @date 96-02-29
12 
13 C> This subprogram performs a spherical transform
14 C> between spectral coefficients of divergences and curls
15 C> and vector fields on a global cylindrical grid.
16 C>
17 C> The wave-space can be either triangular or rhomboidal.
18 C>
19 C> The grid-space can be either an equally-spaced grid
20 C> (with or without pole points) or a Gaussian grid.
21 C>
22 C> The wave and grid fields may have general indexing,
23 C> but each wave field is in sequential 'ibm order',
24 C> i.e. with zonal wavenumber as the slower index.
25 C>
26 C> Transforms are done in latitude pairs for efficiency;
27 C> thus grid arrays for each hemisphere must be passed.
28 C> if so requested, just a subset of the latitude pairs
29 C> may be transformed in each invocation of the subprogram.
30 C>
31 C> The transforms are all multiprocessed over latitude except
32 C> the transform from fourier to spectral is multiprocessed
33 C> over zonal wavenumber to ensure reproducibility.
34 C>
35 C> Transform several fields at a time to improve vectorization.
36 C> subprogram can be called from a multiprocessing environment.
37 C>
38 C> Minimum grid dimensions for unaliased transforms to spectral:
39 C> DIMENSION |LINEAR |QUADRATIC
40 C> ----------------------- |--------- |-------------
41 C> IMAX |2*MAXWV+2 |3*MAXWV/2*2+2
42 C> JMAX (IDRT=4,IROMB=0) |1*MAXWV+1 |3*MAXWV/2+1
43 C> JMAX (IDRT=4,IROMB=1) |2*MAXWV+1 |5*MAXWV/2+1
44 C> JMAX (IDRT=0,IROMB=0) |2*MAXWV+3 |3*MAXWV/2*2+3
45 C> JMAX (IDRT=0,IROMB=1) |4*MAXWV+3 |5*MAXWV/2*2+3
46 C> JMAX (IDRT=256,IROMB=0) |2*MAXWV+1 |3*MAXWV/2*2+1
47 C> JMAX (IDRT=256,IROMB=1) |4*MAXWV+1 |5*MAXWV/2*2+1
48 C>
49 C> @param IROMB spectral domain shape
50 C> (0 for triangular, 1 for rhomboidal)
51 C> @param MAXWV spectral truncation
52 C> @param IDRT grid identifier
53 C> - IDRT=4 for Gaussian grid
54 C> - IDRT=0 for equally-spaced grid including poles
55 C> - IDRT=256 for equally-spaced grid excluding poles
56 C> @param IMAX even number of longitudes.
57 C> @param JMAX number of latitudes.
58 C> @param KMAX number of fields to transform.
59 C> @param IP longitude index for the prime meridian
60 C> @param IS skip number between longitudes
61 C> @param JN skip number between n.h. latitudes from north
62 C> @param JS skip number between s.h. latitudes from south
63 C> @param KW skip number between wave fields
64 C> @param KG skip number between grid fields
65 C> @param JB latitude index (from pole) to begin transform
66 C> @param JE latitude index (from pole) to end transform
67 C> @param JC number of cpus over which to multiprocess
68 C> @param[out] WAVED wave divergence fields if IDIR>0
69 C> [WAVED=(D(GRIDU)/DLAM+D(CLAT*GRIDV)/DPHI)/(CLAT*RERTH)]
70 C> @param[out] WAVEZ wave vorticity fields if IDIR>0
71 C> [WAVEZ=(D(GRIDV)/DLAM-D(CLAT*GRIDU)/DPHI)/(CLAT*RERTH)]
72 C> @param[out] GRIDUN N.H. grid u-winds (starting at jb) if IDIR<0
73 C> @param[out] GRIDUS S.H. grid u-winds (starting at jb) if IDIR<0
74 C> @param[out] GRIDVN N.H. grid v-winds (starting at jb) if IDIR<0
75 C> @param[out] GRIDVS S.H. grid v-winds (starting at jb) if IDIR<0
76 C> @param IDIR transform flag
77 C> (IDIR>0 for wave to grid, IDIR<0 for grid to wave).
78 C>
79 C> @author Iredell @date 96-02-29
80  SUBROUTINE sptranfv(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
81  & IP,IS,JN,JS,KW,KG,JB,JE,JC,
82  & WAVED,WAVEZ,GRIDUN,GRIDUS,GRIDVN,GRIDVS,IDIR)
83 
84  REAL WAVED(*),WAVEZ(*),GRIDUN(*),GRIDUS(*),GRIDVN(*),GRIDVS(*)
85  REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
86  REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
87  REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
88  REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
89  REAL(8) AFFT(50000+4*IMAX), AFFT_TMP(50000+4*IMAX)
90  REAL CLAT(JB:JE),SLAT(JB:JE),WLAT(JB:JE)
91  REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2,JB:JE)
92  REAL PLNTOP(MAXWV+1,JB:JE)
93  INTEGER MP(2)
94  REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2)
95  REAL WTOP(2*(MAXWV+1),2)
96  REAL G(IMAX,2,2)
97  REAL WINC((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2)
98 
99 C SET PARAMETERS
100  mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
101  mp=1
102  CALL sptranf0(iromb,maxwv,idrt,imax,jmax,jb,je,
103  & eps,epstop,enn1,elonn1,eon,eontop,
104  & afft,clat,slat,wlat,pln,plntop)
105 
106 C TRANSFORM WAVE TO GRID
107  IF(idir.GT.0) THEN
108 C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS)
109  DO k=1,kmax
110  afft_tmp=afft
111  kws=(k-1)*kw
112  CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
113  & waved(kws+1),wavez(kws+1),
114  & w(1,1),w(1,2),wtop(1,1),wtop(1,2))
115  DO j=jb,je
116  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
117  & eps,epstop,enn1,elonn1,eon,eontop,
118  & afft_tmp,clat(j),slat(j),wlat(j),
119  & pln(1,j),plntop(1,j),mp,
120  & w(1,1),wtop(1,1),g(1,1,1),idir)
121  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
122  & eps,epstop,enn1,elonn1,eon,eontop,
123  & afft_tmp,clat(j),slat(j),wlat(j),
124  & pln(1,j),plntop(1,j),mp,
125  & w(1,2),wtop(1,2),g(1,1,2),idir)
126  IF(ip.EQ.1.AND.is.EQ.1) THEN
127  DO i=1,imax
128  ijkn=i+(j-jb)*jn+(k-1)*kg
129  ijks=i+(j-jb)*js+(k-1)*kg
130  gridun(ijkn)=g(i,1,1)
131  gridus(ijks)=g(i,2,1)
132  gridvn(ijkn)=g(i,1,2)
133  gridvs(ijks)=g(i,2,2)
134  ENDDO
135  ELSE
136  DO i=1,imax
137  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
138  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
139  gridun(ijkn)=g(i,1,1)
140  gridus(ijks)=g(i,2,1)
141  gridvn(ijkn)=g(i,1,2)
142  gridvs(ijks)=g(i,2,2)
143  ENDDO
144  ENDIF
145  ENDDO
146  ENDDO
147 
148 C TRANSFORM GRID TO WAVE
149  ELSE
150 C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS,WINC)
151  DO k=1,kmax
152  afft_tmp=afft
153  kws=(k-1)*kw
154  w=0
155  wtop=0
156  DO j=jb,je
157  IF(wlat(j).GT.0.) THEN
158  IF(ip.EQ.1.AND.is.EQ.1) THEN
159  DO i=1,imax
160  ijkn=i+(j-jb)*jn+(k-1)*kg
161  ijks=i+(j-jb)*js+(k-1)*kg
162  g(i,1,1)=gridun(ijkn)/clat(j)**2
163  g(i,2,1)=gridus(ijks)/clat(j)**2
164  g(i,1,2)=gridvn(ijkn)/clat(j)**2
165  g(i,2,2)=gridvs(ijks)/clat(j)**2
166  ENDDO
167  ELSE
168  DO i=1,imax
169  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
170  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
171  g(i,1,1)=gridun(ijkn)/clat(j)**2
172  g(i,2,1)=gridus(ijks)/clat(j)**2
173  g(i,1,2)=gridvn(ijkn)/clat(j)**2
174  g(i,2,2)=gridvs(ijks)/clat(j)**2
175  ENDDO
176  ENDIF
177  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
178  & eps,epstop,enn1,elonn1,eon,eontop,
179  & afft_tmp,clat(j),slat(j),wlat(j),
180  & pln(1,j),plntop(1,j),mp,
181  & w(1,1),wtop(1,1),g(1,1,1),idir)
182  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
183  & eps,epstop,enn1,elonn1,eon,eontop,
184  & afft_tmp,clat(j),slat(j),wlat(j),
185  & pln(1,j),plntop(1,j),mp,
186  & w(1,2),wtop(1,2),g(1,1,2),idir)
187  ENDIF
188  ENDDO
189  CALL spuv2dz(iromb,maxwv,enn1,elonn1,eon,eontop,
190  & w(1,1),w(1,2),wtop(1,1),wtop(1,2),
191  & winc(1,1),winc(1,2))
192  waved(kws+1:kws+2*mx)=waved(kws+1:kws+2*mx)+winc(1:2*mx,1)
193  wavez(kws+1:kws+2*mx)=wavez(kws+1:kws+2*mx)+winc(1:2*mx,2)
194  ENDDO
195  ENDIF
196  END
spdz2uv
subroutine spdz2uv(I, M, ENN1, ELONN1, EON, EONTOP, D, Z, U, V, UTOP, VTOP)
Computes the wind components from divergence and vorticity in spectral space.
Definition: spdz2uv.f:49
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:37
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:44
sptranfv
subroutine sptranfv(IROMB, MAXWV, IDRT, IMAX, JMAX, KMAX, IP, IS, JN, JS, KW, KG, JB, JE, JC, WAVED, WAVEZ, GRIDUN, GRIDUS, GRIDVN, GRIDVS, IDIR)
This subprogram performs a spherical transform between spectral coefficients of divergences and curls...
Definition: sptranfv.f:83
spuv2dz
subroutine spuv2dz(I, M, ENN1, ELONN1, EON, EONTOP, U, V, UTOP, VTOP, D, Z)
Computes the divergence and vorticity from wind components in spectral space.
Definition: spuv2dz.f:49