NCEPLIBS-sp  2.3.3
sptranfv.f
Go to the documentation of this file.
1 C> @file
2 C>
3 C> Perform a vector spherical transform
4 C> @author IREDELL @date 96-02-29
5 
6 C> THIS SUBPROGRAM PERFORMS A SPHERICAL TRANSFORM
7 C> BETWEEN SPECTRAL COEFFICIENTS OF DIVERGENCES AND CURLS
8 C> AND VECTOR FIELDS ON A GLOBAL CYLINDRICAL GRID.
9 C> THE WAVE-SPACE CAN BE EITHER TRIANGULAR OR RHOMBOIDAL.
10 C> THE GRID-SPACE CAN BE EITHER AN EQUALLY-SPACED GRID
11 C> (WITH OR WITHOUT POLE POINTS) OR A GAUSSIAN GRID.
12 C> THE WAVE AND GRID FIELDS MAY HAVE GENERAL INDEXING,
13 C> BUT EACH WAVE FIELD IS IN SEQUENTIAL 'IBM ORDER',
14 C> I.E. WITH ZONAL WAVENUMBER AS THE SLOWER INDEX.
15 C> TRANSFORMS ARE DONE IN LATITUDE PAIRS FOR EFFICIENCY;
16 C> THUS GRID ARRAYS FOR EACH HEMISPHERE MUST BE PASSED.
17 C> IF SO REQUESTED, JUST A SUBSET OF THE LATITUDE PAIRS
18 C> MAY BE TRANSFORMED IN EACH INVOCATION OF THE SUBPROGRAM.
19 C> THE TRANSFORMS ARE ALL MULTIPROCESSED OVER LATITUDE EXCEPT
20 C> THE TRANSFORM FROM FOURIER TO SPECTRAL IS MULTIPROCESSED
21 C> OVER ZONAL WAVENUMBER TO ENSURE REPRODUCIBILITY.
22 C> TRANSFORM SEVERAL FIELDS AT A TIME TO IMPROVE VECTORIZATION.
23 C> SUBPROGRAM CAN BE CALLED FROM A MULTIPROCESSING ENVIRONMENT.
24 C>
25 C> PROGRAM HISTORY LOG:
26 C> - 96-02-29 IREDELL
27 C> - 1998-12-15 IREDELL GENERIC FFT USED, OPENMP DIRECTIVES INSERTED
28 C> - 2013-01-16 IREDELL & MIRVIS FIXING AFFT NEGATIVE SHARING EFFECT DURING
29 C> OMP LOOPS BY CREATING TMP AFFT COPY (AFFT_TMP)
30 C> TO BE PRIVATE DURING OMP LOOP THREADING
31 C>
32 C> @param IROMB - INTEGER SPECTRAL DOMAIN SHAPE
33 C> (0 FOR TRIANGULAR, 1 FOR RHOMBOIDAL)
34 C> @param MAXWV - INTEGER SPECTRAL TRUNCATION
35 C> @param IDRT - INTEGER GRID IDENTIFIER
36 C> (IDRT=4 FOR GAUSSIAN GRID,
37 C> IDRT=0 FOR EQUALLY-SPACED GRID INCLUDING POLES,
38 C> IDRT=256 FOR EQUALLY-SPACED GRID EXCLUDING POLES)
39 C> @param IMAX - INTEGER EVEN NUMBER OF LONGITUDES.
40 C> @param JMAX - INTEGER NUMBER OF LATITUDES.
41 C> @param KMAX - INTEGER NUMBER OF FIELDS TO TRANSFORM.
42 C> @param IP - INTEGER LONGITUDE INDEX FOR THE PRIME MERIDIAN
43 C> @param IS - INTEGER SKIP NUMBER BETWEEN LONGITUDES
44 C> @param JN - INTEGER SKIP NUMBER BETWEEN N.H. LATITUDES FROM NORTH
45 C> @param JS - INTEGER SKIP NUMBER BETWEEN S.H. LATITUDES FROM SOUTH
46 C> @param KW - INTEGER SKIP NUMBER BETWEEN WAVE FIELDS
47 C> @param KG - INTEGER SKIP NUMBER BETWEEN GRID FIELDS
48 C> @param JB - INTEGER LATITUDE INDEX (FROM POLE) TO BEGIN TRANSFORM
49 C> @param JE - INTEGER LATITUDE INDEX (FROM POLE) TO END TRANSFORM
50 C> @param JC - INTEGER NUMBER OF CPUS OVER WHICH TO MULTIPROCESS
51 C> @param[out] WAVED - REAL (*) WAVE DIVERGENCE FIELDS IF IDIR>0
52 C> [WAVED=(D(GRIDU)/DLAM+D(CLAT*GRIDV)/DPHI)/(CLAT*RERTH)]
53 C> @param[out] WAVEZ - REAL (*) WAVE VORTICITY FIELDS IF IDIR>0
54 C> [WAVEZ=(D(GRIDV)/DLAM-D(CLAT*GRIDU)/DPHI)/(CLAT*RERTH)]
55 C> @param[out] GRIDUN - REAL (*) N.H. GRID U-WINDS (STARTING AT JB) IF IDIR<0
56 C> @param[out] GRIDUS - REAL (*) S.H. GRID U-WINDS (STARTING AT JB) IF IDIR<0
57 C> @param[out] GRIDVN - REAL (*) N.H. GRID V-WINDS (STARTING AT JB) IF IDIR<0
58 C> @param[out] GRIDVS - REAL (*) S.H. GRID V-WINDS (STARTING AT JB) IF IDIR<0
59 C> @param IDIR - INTEGER TRANSFORM FLAG
60 C> (IDIR>0 FOR WAVE TO GRID, IDIR<0 FOR GRID TO WAVE)
61 C>
62 C> SUBPROGRAMS CALLED:
63 C> - SPTRANF0 SPTRANF SPECTRAL INITIALIZATION
64 C> - SPTRANF1 SPTRANF SPECTRAL TRANSFORM
65 C> - SPDZ2UV COMPUTE WINDS FROM DIVERGENCE AND VORTICITY
66 C> - SPUV2DZ COMPUTE DIVERGENCE AND VORTICITY FROM WINDS
67 C>
68 C> REMARKS: MINIMUM GRID DIMENSIONS FOR UNALIASED TRANSFORMS TO SPECTRAL:
69 C> DIMENSION |LINEAR |QUADRATIC
70 C> ----------------------- |--------- |-------------
71 C> IMAX |2*MAXWV+2 |3*MAXWV/2*2+2
72 C> JMAX (IDRT=4,IROMB=0) |1*MAXWV+1 |3*MAXWV/2+1
73 C> JMAX (IDRT=4,IROMB=1) |2*MAXWV+1 |5*MAXWV/2+1
74 C> JMAX (IDRT=0,IROMB=0) |2*MAXWV+3 |3*MAXWV/2*2+3
75 C> JMAX (IDRT=0,IROMB=1) |4*MAXWV+3 |5*MAXWV/2*2+3
76 C> JMAX (IDRT=256,IROMB=0) |2*MAXWV+1 |3*MAXWV/2*2+1
77 C> JMAX (IDRT=256,IROMB=1) |4*MAXWV+1 |5*MAXWV/2*2+1
78  SUBROUTINE sptranfv(IROMB,MAXWV,IDRT,IMAX,JMAX,KMAX,
79  & IP,IS,JN,JS,KW,KG,JB,JE,JC,
80  & WAVED,WAVEZ,GRIDUN,GRIDUS,GRIDVN,GRIDVS,IDIR)
81 
82  REAL WAVED(*),WAVEZ(*),GRIDUN(*),GRIDUS(*),GRIDVN(*),GRIDVS(*)
83  REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
84  REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
85  REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
86  REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
87  REAL(8) AFFT(50000+4*IMAX), AFFT_TMP(50000+4*IMAX)
88  REAL CLAT(JB:JE),SLAT(JB:JE),WLAT(JB:JE)
89  REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2,JB:JE)
90  REAL PLNTOP(MAXWV+1,JB:JE)
91  INTEGER MP(2)
92  REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2)
93  REAL WTOP(2*(MAXWV+1),2)
94  REAL G(IMAX,2,2)
95  REAL WINC((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2,2)
96 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 C SET PARAMETERS
98  mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
99  mp=1
100  CALL sptranf0(iromb,maxwv,idrt,imax,jmax,jb,je,
101  & eps,epstop,enn1,elonn1,eon,eontop,
102  & afft,clat,slat,wlat,pln,plntop)
103 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 C TRANSFORM WAVE TO GRID
105  IF(idir.GT.0) THEN
106 C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS)
107  DO k=1,kmax
108  afft_tmp=afft
109  kws=(k-1)*kw
110  CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
111  & waved(kws+1),wavez(kws+1),
112  & w(1,1),w(1,2),wtop(1,1),wtop(1,2))
113  DO j=jb,je
114  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
115  & eps,epstop,enn1,elonn1,eon,eontop,
116  & afft_tmp,clat(j),slat(j),wlat(j),
117  & pln(1,j),plntop(1,j),mp,
118  & w(1,1),wtop(1,1),g(1,1,1),idir)
119  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
120  & eps,epstop,enn1,elonn1,eon,eontop,
121  & afft_tmp,clat(j),slat(j),wlat(j),
122  & pln(1,j),plntop(1,j),mp,
123  & w(1,2),wtop(1,2),g(1,1,2),idir)
124  IF(ip.EQ.1.AND.is.EQ.1) THEN
125  DO i=1,imax
126  ijkn=i+(j-jb)*jn+(k-1)*kg
127  ijks=i+(j-jb)*js+(k-1)*kg
128  gridun(ijkn)=g(i,1,1)
129  gridus(ijks)=g(i,2,1)
130  gridvn(ijkn)=g(i,1,2)
131  gridvs(ijks)=g(i,2,2)
132  ENDDO
133  ELSE
134  DO i=1,imax
135  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
136  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
137  gridun(ijkn)=g(i,1,1)
138  gridus(ijks)=g(i,2,1)
139  gridvn(ijkn)=g(i,1,2)
140  gridvs(ijks)=g(i,2,2)
141  ENDDO
142  ENDIF
143  ENDDO
144  ENDDO
145 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146 C TRANSFORM GRID TO WAVE
147  ELSE
148 C$OMP PARALLEL DO PRIVATE(AFFT_TMP,KWS,W,WTOP,G,IJKN,IJKS,WINC)
149  DO k=1,kmax
150  afft_tmp=afft
151  kws=(k-1)*kw
152  w=0
153  wtop=0
154  DO j=jb,je
155  IF(wlat(j).GT.0.) THEN
156  IF(ip.EQ.1.AND.is.EQ.1) THEN
157  DO i=1,imax
158  ijkn=i+(j-jb)*jn+(k-1)*kg
159  ijks=i+(j-jb)*js+(k-1)*kg
160  g(i,1,1)=gridun(ijkn)/clat(j)**2
161  g(i,2,1)=gridus(ijks)/clat(j)**2
162  g(i,1,2)=gridvn(ijkn)/clat(j)**2
163  g(i,2,2)=gridvs(ijks)/clat(j)**2
164  ENDDO
165  ELSE
166  DO i=1,imax
167  ijkn=mod(i+ip-2,imax)*is+(j-jb)*jn+(k-1)*kg+1
168  ijks=mod(i+ip-2,imax)*is+(j-jb)*js+(k-1)*kg+1
169  g(i,1,1)=gridun(ijkn)/clat(j)**2
170  g(i,2,1)=gridus(ijks)/clat(j)**2
171  g(i,1,2)=gridvn(ijkn)/clat(j)**2
172  g(i,2,2)=gridvs(ijks)/clat(j)**2
173  ENDDO
174  ENDIF
175  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
176  & eps,epstop,enn1,elonn1,eon,eontop,
177  & afft_tmp,clat(j),slat(j),wlat(j),
178  & pln(1,j),plntop(1,j),mp,
179  & w(1,1),wtop(1,1),g(1,1,1),idir)
180  CALL sptranf1(iromb,maxwv,idrt,imax,jmax,j,j,
181  & eps,epstop,enn1,elonn1,eon,eontop,
182  & afft_tmp,clat(j),slat(j),wlat(j),
183  & pln(1,j),plntop(1,j),mp,
184  & w(1,2),wtop(1,2),g(1,1,2),idir)
185  ENDIF
186  ENDDO
187  CALL spuv2dz(iromb,maxwv,enn1,elonn1,eon,eontop,
188  & w(1,1),w(1,2),wtop(1,1),wtop(1,2),
189  & winc(1,1),winc(1,2))
190  waved(kws+1:kws+2*mx)=waved(kws+1:kws+2*mx)+winc(1:2*mx,1)
191  wavez(kws+1:kws+2*mx)=wavez(kws+1:kws+2*mx)+winc(1:2*mx,2)
192  ENDDO
193  ENDIF
194 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195  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
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:41
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
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:81
spdz2uv
subroutine spdz2uv(I, M, ENN1, ELONN1, EON, EONTOP, D, Z, U, V, UTOP, VTOP)
Compute winds from divergence and vorticity.
Definition: spdz2uv.f:45