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