SABER
tools_sp.F90
Go to the documentation of this file.
1 ! (C) Copyright 2018-2019 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 !----------------------------------------------------------------------
6 ! Module: tools_sp
7 !> Spectral transforms
8 !----------------------------------------------------------------------
9 module tools_sp
10 
11 use tools_kinds, only: kind_real
12 use tools_const, only: pi
13 
14 implicit none
15 
16 private
17 real(kind=kind_real),parameter:: zero = 0.0_kind_real
18 real(kind=kind_real),parameter:: one = 1.0_kind_real
19 real(kind=kind_real),parameter:: two = 2.0_kind_real
20 real(kind=kind_real),parameter:: ten = 10.0_kind_real
21 real(kind=kind_real),parameter:: half = 0.5_kind_real
22 real(kind=kind_real),parameter:: quarter = 0.25_kind_real
23 
24 public :: splat
25 
26 contains
27 
28 !-----------------------------------------------------------------
29 ! Subroutine: splat
30 !> Compute latitude functions
31 !-----------------------------------------------------------------
32 subroutine splat(IDRT,JMAX,SLAT,WLAT)
33 
34 implicit none
35 
36 ! Passed variables
37 integer,intent(in ) :: idrt !< Grid identifier
38 integer,intent(in ) :: jmax !< Number of latitudes
39 real(kind=kind_real),intent(out) :: slat(jmax) !< Sines of latitude
40 real(kind=kind_real),intent(out) :: wlat(jmax) !< Gaussian weights-
41 
42 ! Local variables
43 real(kind=kind_real):: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2)
44 real(kind=kind_real):: slatd(jmax/2),sp,spmax,eps
45 integer,parameter:: jz=50
46 real(kind=kind_real):: bz(jz)
47 real(kind=kind_real):: dlt,d1,d,r
48 real(kind=kind_real):: awork((jmax+1)/2,((jmax+1)/2))
49 real(kind=kind_real):: bwork(((jmax+1)/2))
50 integer:: jh,jhe,jho,j0
51 integer:: ipvt((jmax+1)/2)
52 real(kind=kind_real),parameter:: c=(one-(two/pi)**2)*quarter
53 integer:: j,js,n
54 
55 data bz / 2.4048255577_kind_real, 5.5200781103_kind_real, &
56  & 8.6537279129_kind_real, 11.7915344391_kind_real, &
57  & 14.9309177086_kind_real, 18.0710639679_kind_real, &
58  & 21.2116366299_kind_real, 24.3524715308_kind_real, &
59  & 27.4934791320_kind_real, 30.6346064684_kind_real, &
60  & 33.7758202136_kind_real, 36.9170983537_kind_real, &
61  & 40.0584257646_kind_real, 43.1997917132_kind_real, &
62  & 46.3411883717_kind_real, 49.4826098974_kind_real, &
63  & 52.6240518411_kind_real, 55.7655107550_kind_real, &
64  & 58.9069839261_kind_real, 62.0484691902_kind_real, &
65  & 65.1899648002_kind_real, 68.3314693299_kind_real, &
66  & 71.4729816036_kind_real, 74.6145006437_kind_real, &
67  & 77.7560256304_kind_real, 80.8975558711_kind_real, &
68  & 84.0390907769_kind_real, 87.1806298436_kind_real, &
69  & 90.3221726372_kind_real, 93.4637187819_kind_real, &
70  & 96.6052679510_kind_real, 99.7468198587_kind_real, &
71  & 102.888374254_kind_real, 106.029930916_kind_real, &
72  & 109.171489649_kind_real, 112.313050280_kind_real, &
73  & 115.454612653_kind_real, 118.596176630_kind_real, &
74  & 121.737742088_kind_real, 124.879308913_kind_real, &
75  & 128.020877005_kind_real, 131.162446275_kind_real, &
76  & 134.304016638_kind_real, 137.445588020_kind_real, &
77  & 140.587160352_kind_real, 143.728733573_kind_real, &
78  & 146.870307625_kind_real, 150.011882457_kind_real, &
79  & 153.153458019_kind_real, 156.295034268_kind_real /
80 
81 ! Initialization
82 eps=ten*epsilon(sp)
83 d1=one
84 j0=0
85 
86 if(idrt==4) then
87  ! Gaussian latitudes
88  jh=jmax/2
89  jhe=(jmax+1)/2
90  r=one/sqrt((real(jmax,kind_real)+half)**2+c)
91  do j=1,min(jh,jz)
92  slatd(j)=cos(bz(j)*r)
93  end do
94  do j=jz+1,jh
95  slatd(j)=cos((bz(jz)+(j-jz)*pi)*r)
96  end do
97  spmax=one
98  do while(spmax>eps)
99  spmax=zero
100  do j=1,jh
101  pkm1(j)=one
102  pk(j)=slatd(j)
103  end do
104  do n=2,jmax
105  do j=1,jh
106  pkm2(j)=pkm1(j)
107  pkm1(j)=pk(j)
108  pk(j)=real((2*n-1)*slatd(j)*pkm1(j)-(n-1)*pkm2(j),kind_real)/real(n,kind_real)
109  end do
110  end do
111  do j=1,jh
112  sp=pk(j)*(one-slatd(j)**2)/(jmax*(pkm1(j)-slatd(j)*pk(j)))
113  slatd(j)=slatd(j)-sp
114  spmax=max(spmax,abs(sp))
115  end do
116  end do
117  do j=1,jh
118  slat(j)=slatd(j)
119  wlat(j)=(two*(one-slatd(j)**2))/(jmax*pkm1(j))**2
120  slat(jmax+1-j)=-slat(j)
121  wlat(jmax+1-j)=wlat(j)
122  end do
123  if(jhe>jh) then
124  slat(jhe)=zero
125  wlat(jhe)=two/jmax**2
126  do n=2,jmax,2
127  wlat(jhe)=wlat(jhe)*n**2/(n-1)**2
128  end do
129  end if
130 else if(idrt==0) then
131  ! Equally-spaced latitudes including poles
132  jh=jmax/2
133  jhe=(jmax+1)/2
134  jho=jhe-1
135  dlt=pi/(jmax-1)
136  slat(1)=one
137  do j=2,jh
138  slat(j)=cos(real((j-1),kind_real)*dlt)
139  end do
140  do js=1,jho
141  do j=1,jho
142  awork(js,j)=cos(real(2*(js-1)*j,kind_real)*dlt)
143  end do
144  end do
145  do js=1,jho
146  bwork(js)=-d1/real((4*(js-1)**2-1),kind_real)
147  end do
148  call ludcmp(awork,jho,jhe,ipvt)
149  call lubksb(awork,jho,jhe,ipvt,bwork)
150  wlat(1)=zero
151  do j=1,jho
152  wlat(j+1)=bwork(j)
153  end do
154  do j=1,jh
155  slat(jmax+1-j)=-slat(j)
156  wlat(jmax+1-j)=wlat(j)
157  end do
158  if(jhe>jh) then
159  slat(jhe)=zero
160  wlat(jhe)=two*wlat(jhe)
161  end if
162 else if(idrt==256) then
163  ! Equally-spaced latitudes excluding poles
164  jh=jmax/2
165  jhe=(jmax+1)/2
166  jho=jhe
167  dlt=pi/jmax
168  slat(1)=one
169  do j=1,jh
170  slat(j)=cos((real(j,kind_real)-half)*dlt)
171  end do
172  do js=1,jho
173  do j=1,jho
174  awork(js,j)=cos(real(2*(js-1),kind_real)*(real(j,kind_real)-half)*dlt)
175  end do
176  end do
177  do js=1,jho
178  bwork(js)=-d1/real(4*(js-1)**2-1,kind_real)
179  end do
180  call ludcmp(awork,jho,jhe,ipvt,d)
181  call lubksb(awork,jho,jhe,ipvt,bwork)
182  wlat(1)=zero
183  do j=1,jho
184  wlat(j)=bwork(j)
185  end do
186  do j=1,jh
187  slat(jmax+1-j)=-slat(j)
188  wlat(jmax+1-j)=wlat(j)
189  end do
190  if(jhe>jh) then
191  slat(jhe)=zero
192  wlat(jhe)=two*wlat(jhe)
193  end if
194 end if
195 
196 end subroutine splat
197 
198 !-----------------------------------------------------------------
199 ! Subroutine: lubksb
200 !> Solves a system of linear equations, follows call to LUDCMP
201 !-----------------------------------------------------------------
202 subroutine lubksb(A,N,NP,INDX,B)
203 
204 implicit none
205 
206 ! Passed varaibles
207 integer,intent(in):: NP !< ?
208 integer,intent(in):: N !< ?
209 real(kind=kind_real),intent(in):: a(np,np) !< ?
210 real(kind=kind_real),intent(inout):: b(n) !< ?
211 integer,intent(in):: INDX(N) !< ?
212 
213 ! Local variables
214 real(kind=kind_real):: sum
215 integer:: I,II,J,LL
216 
217 ii=0
218 do i=1,n
219  ll=indx(i)
220  sum=b(ll)
221  b(ll)=b(i)
222  if(ii/=0) then
223  do j=ii,i-1
224  sum=sum-a(i,j)*b(j)
225  end do
226  else if(abs(sum)>zero) then
227  ii=i
228  end if
229  b(i)=sum
230 end do
231 do i=n,1,-1
232  sum=b(i)
233  if(i<n) then
234  do j=i+1,n
235  sum=sum-a(i,j)*b(j)
236  end do
237  end if
238  b(i)=sum/a(i,i)
239 end do
240 
241 end subroutine lubksb
242 
243 !-----------------------------------------------------------------
244 ! Subroutine: ludcmp
245 !> Replaces an NxN matrix A with the LU decomposition
246 !-----------------------------------------------------------------
247 subroutine ludcmp(A,N,NP,INDX,D)
248 
249 implicit none
250 
251 ! Passed variables
252 integer,intent(in):: N !< ?
253 integer,intent(in):: NP !< ?
254 real(kind=kind_real),intent(inout):: a(np,np) !< ?
255 integer,intent(out):: INDX(N) !< ?
256 real(kind=kind_real),intent(out),optional:: d !< ?
257 
258 ! Local variables
259 real(kind=kind_real),parameter:: tiny=1.0e-20
260 real(kind=kind_real):: aamax,sum,dum
261 real(kind=kind_real):: vv(n)
262 integer:: IMAX
263 integer:: I,J,K
264 
265 d=one
266 do i=1,n
267  aamax=zero
268  do j=1,n
269  if(abs(a(i,j))>aamax) aamax=abs(a(i,j))
270  end do
271  if(.NOT.(abs(aamax)>zero)) print *, 'SINGULAR MATRIX.'
272  vv(i)=one/aamax
273 end do
274 do j=1,n
275  if(j>1) then
276  do i=1,j-1
277  sum=a(i,j)
278  if(i>1) then
279  do k=1,i-1
280  sum=sum-a(i,k)*a(k,j)
281  end do
282  a(i,j)=sum
283  end if
284  end do
285  end if
286  aamax=zero
287  do i=j,n
288  sum=a(i,j)
289  if(j>1) then
290  do k=1,j-1
291  sum=sum-a(i,k)*a(k,j)
292  end do
293  a(i,j)=sum
294  end if
295  dum=vv(i)*abs(sum)
296  if(dum>=aamax) then
297  imax=i
298  aamax=dum
299  end if
300  end do
301  if(j/=imax) then
302  do k=1,n
303  dum=a(imax,k)
304  a(imax,k)=a(j,k)
305  a(j,k)=dum
306  end do
307  d=-d
308  vv(imax)=vv(j)
309  end if
310  indx(j)=imax
311  if(j/=n) then
312  if(.NOT.(abs(a(j,j))>zero)) a(j,j)=tiny
313  dum=one/a(j,j)
314  do i=j+1,n
315  a(i,j)=a(i,j)*dum
316  end do
317  end if
318 end do
319 if(.NOT.(abs(a(n,n))>zero)) a(n,n)=tiny
320 
321 end subroutine ludcmp
322 
323 end module tools_sp
tools_const
Define usual constants and missing values.
Definition: tools_const.F90:8
tools_sp::ten
real(kind=kind_real), parameter ten
Definition: tools_sp.F90:20
tools_sp::quarter
real(kind=kind_real), parameter quarter
Definition: tools_sp.F90:22
tools_sp::lubksb
subroutine lubksb(A, N, NP, INDX, B)
Solves a system of linear equations, follows call to LUDCMP.
Definition: tools_sp.F90:203
tools_sp::zero
real(kind=kind_real), parameter zero
Definition: tools_sp.F90:17
tools_sp::ludcmp
subroutine ludcmp(A, N, NP, INDX, D)
Replaces an NxN matrix A with the LU decomposition.
Definition: tools_sp.F90:248
tools_sp::splat
subroutine, public splat(IDRT, JMAX, SLAT, WLAT)
Compute latitude functions.
Definition: tools_sp.F90:33
tools_sp::half
real(kind=kind_real), parameter half
Definition: tools_sp.F90:21
tools_sp::two
real(kind=kind_real), parameter two
Definition: tools_sp.F90:19
tools_kinds
Kinds definition.
Definition: tools_kinds.F90:8
tools_sp
Spectral transforms.
Definition: tools_sp.F90:9
tools_const::pi
real(kind_real), parameter, public pi
Definition: tools_const.F90:14
tools_kinds::kind_real
integer, parameter, public kind_real
Definition: tools_kinds.F90:18
tools_sp::one
real(kind=kind_real), parameter one
Definition: tools_sp.F90:18