NCEPLIBS-sp 2.4.0
sptgpmv.f
Go to the documentation of this file.
1C> @file
2C> @brief Transform spectral vector to Mercator
3C> ### Program history log:
4C> Date | Programmer | Comments
5C> -----|------------|----------
6C> 96-02-29 | IREDELL | Initial.
7C> 1998-12-15 | IREDELL | OpenMP directives inserted.
8C> @author IREDELL @date 96-02-29
9
10C> This subprogram performs a spherical transform
11C> from spectral coefficients of divergences and curls
12C> to vector fields on a Mercator grid.
13C>
14C> The wave-space can be either triangular or rhomboidal.
15C>
16C> The wave and grid fields may have general indexing,
17C> but each wave field is in sequential 'ibm order',
18C> i.e., with zonal wavenumber as the slower index.
19C>
20C> The Mercator grid is identified by the location
21C> of its first point and by its respective increments.
22C>
23C> The transforms are all multiprocessed over sector points.
24C> Transform several fields at a time to improve vectorization.
25C>
26C> Subprogram can be called from a multiprocessing environment.
27C>
28C> @param IROMB Spectral domain shape
29C> (0 for triangular, 1 for rhomboidal)
30C> @param MAXWV Spectral truncation
31C> @param KMAX Number of fields to transform
32C> @param MI Number of points in the faster zonal direction
33C> @param MJ Number of points in the slower merid direction
34C> @param KWSKIP Skip number between wave fields
35C> (defaults to (MAXWV+1)*((IROMB+1)*MAXWV+2) if KWSKIP=0)
36C> @param KGSKIP Skip number between grid fields
37C> (defaults to MI*MJ if KGSKIP=0)
38C> @param NISKIP Skip number between grid i-points
39C> (defaults to 1 if NISKIP=0)
40C> @param NJSKIP Skip number between grid j-points
41C> (defaults to MI if NJSKIP=0)
42C> @param RLAT1 Latitude of the first grid point in degrees
43C> @param RLON1 Longitude of the first grid point in degrees
44C> @param DLAT Latitude increment in degrees such that
45C> D(PHI)/D(J)=DLAT*COS(PHI) where J is meridional index.
46C> DLAT is negative for grids indexed southward.
47C> (in terms of grid increment dy valid at latitude RLATI,
48C> The latitude increment DLAT is determined as
49C> DLAT=DPR*DY/(RERTH*COS(RLATI/DPR))
50C> where DPR=180/PI and RERTH is Earth's radius)
51C> @param DLON longitude increment in degrees such that
52C> D(LAMBDA)/D(I)=DLON where I is zonal index.
53C> DLON is negative for grids indexed westward.
54C> @param WAVED Wave divergence fields
55C> @param WAVEZ Wave vorticity fields
56C> @param UM Mercator u-winds
57C> @param VM Mercator v-winds
58C>
59C> @author IREDELL @date 96-02-29
60 SUBROUTINE sptgpmv(IROMB,MAXWV,KMAX,MI,MJ,
61 & KWSKIP,KGSKIP,NISKIP,NJSKIP,
62 & RLAT1,RLON1,DLAT,DLON,WAVED,WAVEZ,UM,VM)
63
64 REAL WAVED(*),WAVEZ(*),UM(*),VM(*)
65 REAL EPS((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EPSTOP(MAXWV+1)
66 REAL ENN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
67 REAL ELONN1((MAXWV+1)*((IROMB+1)*MAXWV+2)/2)
68 REAL EON((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),EONTOP(MAXWV+1)
69 INTEGER MP(2*KMAX)
70 REAL W((MAXWV+1)*((IROMB+1)*MAXWV+2)/2*2+1,2*KMAX)
71 REAL WTOP(2*(MAXWV+1),2*KMAX)
72 REAL PLN((MAXWV+1)*((IROMB+1)*MAXWV+2)/2),PLNTOP(MAXWV+1)
73 REAL F(2*MAXWV+3,2,2*KMAX)
74 REAL CLAT(MJ),SLAT(MJ),CLON(MAXWV,MI),SLON(MAXWV,MI)
75 parameter(rerth=6.3712e6)
76 parameter(pi=3.14159265358979,dpr=180./pi)
77C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78C CALCULATE PRELIMINARY CONSTANTS
79 CALL spwget(iromb,maxwv,eps,epstop,enn1,elonn1,eon,eontop)
80 mx=(maxwv+1)*((iromb+1)*maxwv+2)/2
81 mxtop=maxwv+1
82 mdim=2*mx+1
83 idim=2*maxwv+3
84 kw=kwskip
85 kg=kgskip
86 ni=niskip
87 nj=njskip
88 IF(kw.EQ.0) kw=2*mx
89 IF(kg.EQ.0) kg=mi*mj
90 IF(ni.EQ.0) ni=1
91 IF(nj.EQ.0) nj=mi
92 DO i=1,mi
93 rlon=mod(rlon1+dlon*(i-1)+3600,360.)
94 DO l=1,maxwv
95 clon(l,i)=cos(l*rlon/dpr)
96 slon(l,i)=sin(l*rlon/dpr)
97 ENDDO
98 ENDDO
99 ye=1-log(tan((rlat1+90)/2/dpr))*dpr/dlat
100 DO j=1,mj
101 rlat=atan(exp(dlat/dpr*(j-ye)))*2*dpr-90
102 clat(j)=cos(rlat/dpr)
103 slat(j)=sin(rlat/dpr)
104 ENDDO
105 mp=1
106C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107C CALCULATE SPECTRAL WINDS
108C$OMP PARALLEL DO PRIVATE(KWS)
109 DO k=1,kmax
110 kws=(k-1)*kw
111 CALL spdz2uv(iromb,maxwv,enn1,elonn1,eon,eontop,
112 & waved(kws+1),wavez(kws+1),
113 & w(1,k),w(1,kmax+k),wtop(1,k),wtop(1,kmax+k))
114 ENDDO
115C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116C TRANSFORM TO GRID
117C$OMP PARALLEL DO PRIVATE(PLN,PLNTOP,F,KU,KV,IJK)
118 DO j=1,mj
119 CALL splegend(iromb,maxwv,slat(j),clat(j),eps,epstop,
120 & pln,plntop)
121 CALL spsynth(iromb,maxwv,2*maxwv,idim,mdim,2*mxtop,2*kmax,
122 & clat(j),pln,plntop,mp,w,wtop,f)
123 DO k=1,kmax
124 ku=k
125 kv=k+kmax
126 DO i=1,mi
127 ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
128 um(ijk)=f(1,1,ku)
129 vm(ijk)=f(1,1,kv)
130 ENDDO
131 DO l=1,maxwv
132 DO i=1,mi
133 ijk=(i-1)*ni+(j-1)*nj+(k-1)*kg+1
134 um(ijk)=um(ijk)+2.*(f(2*l+1,1,ku)*clon(l,i)
135 & -f(2*l+2,1,ku)*slon(l,i))
136 vm(ijk)=vm(ijk)+2.*(f(2*l+1,1,kv)*clon(l,i)
137 & -f(2*l+2,1,kv)*slon(l,i))
138 ENDDO
139 ENDDO
140 ENDDO
141 ENDDO
142C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143 END
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
subroutine splegend(I, M, SLAT, CLAT, EPS, EPSTOP, PLN, PLNTOP)
Evaluates the orthonormal associated Legendre polynomials in the spectral domain at a given latitude.
Definition: splegend.f:45
subroutine spsynth(I, M, IM, IX, NC, NCTOP, KM, CLAT, PLN, PLNTOP, MP, SPC, SPCTOP, F)
Synthesizes Fourier coefficients from spectral coefficients for a latitude pair (Northern and Souther...
Definition: spsynth.f:39
subroutine sptgpmv(IROMB, MAXWV, KMAX, MI, MJ, KWSKIP, KGSKIP, NISKIP, NJSKIP, RLAT1, RLON1, DLAT, DLON, WAVED, WAVEZ, UM, VM)
This subprogram performs a spherical transform from spectral coefficients of divergences and curls to...
Definition: sptgpmv.f:63
subroutine spwget(IROMB, MAXWV, EPS, EPSTOP, ENN1, ELONN1, EON, EONTOP)
This subprogram gets wave-space constants.
Definition: spwget.f:18