NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
splat.F
Go to the documentation of this file.
1C> @file
2C> @brief Computes cosines of colatitude and Gaussian weights
3C> for sets of latitudes.
4C>
5C> ### Program History Log
6C> Date | Programmer | Comments
7C> -----|------------|---------
8C> 96-02-20 | Iredell | Initial.
9C> 97-10-20 | Iredell | Adjust precision.
10C> 98-06-11 | Iredell | Generalize precision using FORTRAN 90 intrinsic.
11C> 1998-12-03 | Iredell | Generalize precision further.
12C> 1998-12-03 | Iredell | Uses AIX ESSL BLAS calls.
13C> 2009-12-27 | D. Stark | Updated to switch between ESSL calls on an AIX platform, and Numerical Recipies calls elsewise.
14C> 2010-12-30 | Slovacek | Update alignment so preprocessor does not cause compilation failure.
15C> 2012-09-01 | E. Mirvis & M.Iredell | Merging & debugging linux errors of _d and _8 using generic LU factorization.
16C> 2012-11-05 | E. Mirvis | Generic FFTPACK and LU lapack were removed.
17C>
18C> @author Iredell @date 96-02-20
19
20C> Computes cosines of colatitude and Gaussian weights
21C> for one of the following specific global sets of latitudes.
22C> - Gaussian latitudes (IDRT=4)
23C> - Equally-spaced latitudes including poles (IDRT=0)
24C> - Equally-spaced latitudes excluding poles (IDRT=256)
25C>
26C> The Gaussian latitudes are located at the zeroes of the
27C> Legendre polynomial of the given order. These latitudes
28C> are efficient for reversible transforms from spectral space.
29C> (About twice as many equally-spaced latitudes are needed.)
30C> The weights for the equally-spaced latitudes are based on
31C> Ellsaesser (JAM,1966). (No weight is given the pole point.)
32C> Note that when analyzing grid to spectral in latitude pairs,
33C> if an equator point exists, its weight should be halved.
34C> This version invokes the ibm essl matrix solver.
35C>
36C> @param[in] IDRT grid identifier
37C> - 4 for Gaussian grid
38C> - 0 for equally-spaced grid including poles
39C> - 256 for equally-spaced grid excluding poles
40C> @param[in] JMAX number of latitudes
41C> @param[out] SLAT sines of latitude
42C> @param[out] WLAT Gaussian weights
43C>
44C> @author Iredell @date 96-02-20
45 SUBROUTINE splat(IDRT,JMAX,SLAT,WLAT)
46 REAL SLAT(JMAX),WLAT(JMAX)
47 INTEGER,PARAMETER:: KD=selected_real_kind(15,45)
48 REAL(KIND=kd):: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2)
49 REAL(KIND=kd):: slatd(jmax/2),sp,spmax,eps=10.*epsilon(sp)
50 parameter(jz=50)
51 REAL BZ(JZ)
52 DATA bz / 2.4048255577, 5.5200781103,
53 $ 8.6537279129, 11.7915344391, 14.9309177086, 18.0710639679,
54 $ 21.2116366299, 24.3524715308, 27.4934791320, 30.6346064684,
55 $ 33.7758202136, 36.9170983537, 40.0584257646, 43.1997917132,
56 $ 46.3411883717, 49.4826098974, 52.6240518411, 55.7655107550,
57 $ 58.9069839261, 62.0484691902, 65.1899648002, 68.3314693299,
58 $ 71.4729816036, 74.6145006437, 77.7560256304, 80.8975558711,
59 $ 84.0390907769, 87.1806298436, 90.3221726372, 93.4637187819,
60 $ 96.6052679510, 99.7468198587, 102.888374254, 106.029930916,
61 $ 109.171489649, 112.313050280, 115.454612653, 118.596176630,
62 $ 121.737742088, 124.879308913, 128.020877005, 131.162446275,
63 $ 134.304016638, 137.445588020, 140.587160352, 143.728733573,
64 $ 146.870307625, 150.011882457, 153.153458019, 156.295034268 /
65 REAL:: DLT,D1=1.
66 REAL AWORK((JMAX+1)/2,((JMAX+1)/2)),BWORK(((JMAX+1)/2))
67 INTEGER:: JHE,JHO,INFO
68 INTEGER IPVT((JMAX+1)/2)
69 parameter(pi=3.14159265358979,c=(1.-(2./pi)**2)*0.25)
70C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71C GAUSSIAN LATITUDES
72 IF(idrt.EQ.4) THEN
73 jh=jmax/2
74 jhe=(jmax+1)/2
75 r=1./sqrt((jmax+0.5)**2+c)
76 DO j=1,min(jh,jz)
77 slatd(j)=cos(bz(j)*r)
78 ENDDO
79 DO j=jz+1,jh
80 slatd(j)=cos((bz(jz)+(j-jz)*pi)*r)
81 ENDDO
82 spmax=1.
83 DO WHILE(spmax.GT.eps)
84 spmax=0.
85 DO j=1,jh
86 pkm1(j)=1.
87 pk(j)=slatd(j)
88 ENDDO
89 DO n=2,jmax
90 DO j=1,jh
91 pkm2(j)=pkm1(j)
92 pkm1(j)=pk(j)
93 pk(j)=((2*n-1)*slatd(j)*pkm1(j)-(n-1)*pkm2(j))/n
94 ENDDO
95 ENDDO
96 DO j=1,jh
97 sp=pk(j)*(1.-slatd(j)**2)/(jmax*(pkm1(j)-slatd(j)*pk(j)))
98 slatd(j)=slatd(j)-sp
99 spmax=max(spmax,abs(sp))
100 ENDDO
101 ENDDO
102CDIR$ IVDEP
103 DO j=1,jh
104 slat(j)=real(slatd(j))
105 wlat(j)=real((2.*(1.-slatd(j)**2))/(jmax*pkm1(j))**2)
106 slat(jmax+1-j)=-slat(j)
107 wlat(jmax+1-j)=wlat(j)
108 ENDDO
109 IF(jhe.GT.jh) THEN
110 slat(jhe)=0.
111 wlat(jhe)=2./jmax**2
112 DO n=2,jmax,2
113 wlat(jhe)=wlat(jhe)*n**2/(n-1)**2
114 ENDDO
115 ENDIF
116C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117C EQUALLY-SPACED LATITUDES INCLUDING POLES
118 ELSEIF(idrt.EQ.0) THEN
119 jh=jmax/2
120 jhe=(jmax+1)/2
121 jho=jhe-1
122 dlt=pi/(jmax-1)
123 slat(1)=1.
124 DO j=2,jh
125 slat(j)=cos((j-1)*dlt)
126 ENDDO
127 DO js=1,jho
128 DO j=1,jho
129 awork(js,j)=cos(2*(js-1)*j*dlt)
130 ENDDO
131 ENDDO
132 DO js=1,jho
133 bwork(js)=-d1/(4*(js-1)**2-1)
134 ENDDO
135
136 ! Call LAPACK routines
137#if (LSIZE==4)
138 CALL sgetrf(jho, jho, awork, jhe, ipvt, info)
139 CALL sgetrs('N', jho, 1, awork, jhe, ipvt, bwork, jho, info)
140#else
141 CALL dgetrf(jho, jho, awork, jhe, ipvt, info)
142 CALL dgetrs('N', jho, 1, awork, jhe, ipvt, bwork, jho, info)
143#endif
144
145 wlat(1)=0.
146 DO j=1,jho
147 wlat(j+1)=bwork(j)
148 ENDDO
149CDIR$ IVDEP
150 DO j=1,jh
151 print *, j, jmax, jmax+1-j
152 slat(jmax+1-j)=-slat(j)
153 wlat(jmax+1-j)=wlat(j)
154 ENDDO
155 IF(jhe.GT.jh) THEN
156 slat(jhe)=0.
157 wlat(jhe)=2.*wlat(jhe)
158 ENDIF
159C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160C EQUALLY-SPACED LATITUDES EXCLUDING POLES
161 ELSEIF(idrt.EQ.256) THEN
162 jh=jmax/2
163 jhe=(jmax+1)/2
164 jho=jhe
165 dlt=pi/jmax
166 slat(1)=1.
167 DO j=1,jh
168 slat(j)=cos((j-0.5)*dlt)
169 ENDDO
170 DO js=1,jho
171 DO j=1,jho
172 awork(js,j)=cos(2*(js-1)*(j-0.5)*dlt)
173 ENDDO
174 ENDDO
175 DO js=1,jho
176 bwork(js)=-d1/(4*(js-1)**2-1)
177 ENDDO
178
179 ! Call LAPACK routines
180#if (LSIZE==4)
181 CALL sgetrf(jho, jho, awork, jhe, ipvt, info)
182 CALL sgetrs('N', jho, 1, awork, jhe, ipvt, bwork, jho, info)
183#else
184 CALL dgetrf(jho, jho, awork, jhe, ipvt, info)
185 CALL dgetrs('N', jho, 1, awork, jhe, ipvt, bwork, jho, info)
186#endif
187
188 wlat(1)=0.
189 DO j=1,jho
190 wlat(j)=bwork(j)
191 ENDDO
192CDIR$ IVDEP
193 DO j=1,jh
194 slat(jmax+1-j)=-slat(j)
195 wlat(jmax+1-j)=wlat(j)
196 ENDDO
197 IF(jhe.GT.jh) THEN
198 slat(jhe)=0.
199 wlat(jhe)=2.*wlat(jhe)
200 ENDIF
201 ENDIF
202C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203 RETURN
204 END
splat
subroutine splat(idrt, jmax, slat, wlat)
Computes cosines of colatitude and Gaussian weights for one of the following specific global sets of ...
Definition splat.F:46