NCEPLIBS-ip  5.1.0
splat.F
Go to the documentation of this file.
1 C> @file
2 C> @brief Computes cosines of colatitude and Gaussian weights
3 C> for sets of latitudes.
4 C>
5 C> ### Program History Log
6 C> Date | Programmer | Comments
7 C> -----|------------|---------
8 C> 96-02-20 | Iredell | Initial.
9 C> 97-10-20 | Iredell | Adjust precision.
10 C> 98-06-11 | Iredell | Generalize precision using FORTRAN 90 intrinsic.
11 C> 1998-12-03 | Iredell | Generalize precision further.
12 C> 1998-12-03 | Iredell | Uses AIX ESSL BLAS calls.
13 C> 2009-12-27 | D. Stark | Updated to switch between ESSL calls on an AIX platform, and Numerical Recipies calls elsewise.
14 C> 2010-12-30 | Slovacek | Update alignment so preprocessor does not cause compilation failure.
15 C> 2012-09-01 | E. Mirvis & M.Iredell | Merging & debugging linux errors of _d and _8 using generic LU factorization.
16 C> 2012-11-05 | E. Mirvis | Generic FFTPACK and LU lapack were removed.
17 C>
18 C> @author Iredell @date 96-02-20
19 
20 C> Computes cosines of colatitude and Gaussian weights
21 C> for one of the following specific global sets of latitudes.
22 C> - Gaussian latitudes (IDRT=4)
23 C> - Equally-spaced latitudes including poles (IDRT=0)
24 C> - Equally-spaced latitudes excluding poles (IDRT=256)
25 C>
26 C> The Gaussian latitudes are located at the zeroes of the
27 C> Legendre polynomial of the given order. These latitudes
28 C> are efficient for reversible transforms from spectral space.
29 C> (About twice as many equally-spaced latitudes are needed.)
30 C> The weights for the equally-spaced latitudes are based on
31 C> Ellsaesser (JAM,1966). (No weight is given the pole point.)
32 C> Note that when analyzing grid to spectral in latitude pairs,
33 C> if an equator point exists, its weight should be halved.
34 C> This version invokes the ibm essl matrix solver.
35 C>
36 C> @param[in] IDRT grid identifier
37 C> - 4 for Gaussian grid
38 C> - 0 for equally-spaced grid including poles
39 C> - 256 for equally-spaced grid excluding poles
40 C> @param[in] JMAX number of latitudes
41 C> @param[out] SLAT sines of latitude
42 C> @param[out] WLAT Gaussian weights
43 C>
44 C> @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)
70 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 C 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
102 CDIR$ 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
116 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117 C 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
149 CDIR$ 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
159 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160 C 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
192 CDIR$ 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
202 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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