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