forked from NOAA-PSL/stochastic_physics
-
Notifications
You must be signed in to change notification settings - Fork 0
/
pln2eo_stochy.f
290 lines (275 loc) · 8.68 KB
/
pln2eo_stochy.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
!>@brief The module 'pln2eo_a_stochy_mod' contains the subroutine pln2eo_a_stochy
module pln2eo_a_stochy_mod
implicit none
contains
!>@brief The subroutine 'pln2eo_a_stochy' calculates the assoicate legendre polynomials
!>@details This code is taken from the legacy spectral GFS
subroutine pln2eo_a_stochy(plnev_a,plnod_a,epse,epso,colrad_a,
& ls_node,num_lat)
!
! use x-number method to archieve accuracy due to recursive to avoid
! underflow and overflow if necessary by henry juang 2012 july
!
use stochy_resol_def
use spectral_layout_mod
use kinddef
implicit none
!
! define x number constant for real8 start
integer, parameter :: in_f = 960 , in_h = in_f/2
real(kind=kind_dbl_prec), parameter :: bb_f = 2.d0 ** ( in_f )
real(kind=kind_dbl_prec), parameter :: bs_f = 2.d0 ** (-in_f )
real(kind=kind_dbl_prec), parameter :: bb_h = 2.d0 ** ( in_h )
real(kind=kind_dbl_prec), parameter :: bs_h = 2.d0 ** (-in_h )
! define x number constant end
cc
real(kind=kind_dbl_prec) plnev_a(len_trie_ls,latg2)
real(kind=kind_dbl_prec) plnod_a(len_trio_ls,latg2)
cc
real(kind=kind_dbl_prec) epse(len_trie_ls)
real(kind=kind_dbl_prec) epso(len_trio_ls)
cc
real(kind=kind_dbl_prec) colrad_a(latg2)
cc
integer ls_node(ls_dim,3)
cc
integer num_lat
cc
!cmr ls_node(1,1) ... ls_node(ls_max_node,1) : values of L
!cmr ls_node(1,2) ... ls_node(ls_max_node,2) : values of jbasev
!cmr ls_node(1,3) ... ls_node(ls_max_node,3) : values of jbasod
cc
integer l,lat,locl,max_l,n
cc
integer indev
integer indod
cc
! need index for alp to be x-number
integer id, ialp1, ialp2, ialp3, iprod
integer ialp10(0:jcap)
real(kind=kind_dbl_prec) aa, bb, w
real(kind=kind_dbl_prec) a,alp1,alp2,alp3,b
real(kind=kind_dbl_prec) cos2,fl,prod,sinlat,coslat
cc
real(kind=kind_dbl_prec) alp10(0:jcap)
cc
real(kind=kind_dbl_prec) cons0,cons0p5,cons1,cons2,cons3 !constant
cc
cc
integer indlsev,jbasev
integer indlsod,jbasod
cc
include 'function2'
cc
cc
cons0=0.0d0 !constant
cons0p5=0.5d0 !constant
cons1=1.0d0 !constant
cons2=2.0d0 !constant
cons3=3.0d0 !constant
cc
cc
max_l=-1
do locl=1,ls_max_node
max_l = max ( max_l, ls_node(locl,1) )
enddo
cc
cc
do lat=1,num_lat
cc
sinlat = cos(colrad_a(lat))
cos2=cons1-sinlat*sinlat !constant
coslat = sqrt(cos2)
! use x number for alp10
alp10(0) = sqrt(0.5)
ialp10(0) = 0
do l=1,max_l
fl = l
prod=coslat*sqrt(cons1+cons1/(cons2*fl))
iprod=0
w = abs(prod)
if( w.ge.bb_h ) then
prod = prod * bs_f
iprod = iprod + 1
elseif( w.lt.bs_h ) then
prod = prod * bb_f
iprod = iprod - 1
endif
alp10(l)=alp10(l-1)*prod
ialp10(l)=ialp10(l-1)+iprod
w = abs(alp10(l))
if( w.ge.bb_h ) then
alp10(l) = alp10(l) * bs_f
ialp10(l) = ialp10(l) + 1
elseif( w.lt.bs_h ) then
alp10(l) = alp10(l) * bb_f
ialp10(l) = ialp10(l) - 1
endif
enddo
cc
do locl=1,ls_max_node
l=ls_node(locl,1)
jbasev=ls_node(locl,2)
jbasod=ls_node(locl,3)
n=l
fl=l
! get m=normalized x number for alp1 start
alp1=alp10(l)
ialp1=ialp10(l)
indev=indlsev(n ,l)
indod=indlsod(n+1,l)
! x2f plnev_a(indev ,lat)=alp1
! x2f start
if( ialp1.eq.0 ) then
plnev_a(indev ,lat)=alp1
elseif( ialp1.eq.-1 ) then
plnev_a(indev ,lat)=alp1 * bs_f
elseif( ialp1.lt.-1 ) then
plnev_a(indev ,lat)=0.0
!! plnev_a(indev ,lat)=alp1 * bs_f * bs_f
else
plnev_a(indev ,lat)=alp1 * bb_f
endif
! x2f end
! xltime alp2=sqrt(cons2*fl+cons3)*sinlat*alp1 !constant
! xltime start
prod=sqrt(cons2*fl+cons3)*sinlat
iprod=0
w = abs(prod)
if( w.ge.bb_h ) then
prod = prod * bs_f
iprod = iprod + 1
elseif( w.lt.bs_h ) then
prod = prod * bb_f
iprod = iprod - 1
endif
alp2=alp1*prod
ialp2 = ialp1 + iprod
! xltime end
! norm alp2 start
w = abs(alp2)
if( w.ge.bb_h ) then
alp2 = alp2*bs_f
ialp2 = ialp2 + 1
elseif( w.lt.bs_h ) then
alp2 = alp2*bb_f
ialp2 = ialp2 - 1
endif
! norm alp2 end
! x2f plnod_a(indod ,lat)=alp2
! x2f start
if( ialp2.eq.0 ) then
plnod_a(indod ,lat)=alp2
elseif( ialp2.eq.-1 ) then
plnod_a(indod ,lat)=alp2 * bs_f
elseif( ialp2.lt.-1 ) then
plnod_a(indod ,lat)=0.0
!! plnod_a(indod ,lat)=alp2 * bs_f * bs_f
else
plnod_a(indod ,lat)=alp2 * bb_f
endif
! x2f end
cc
do n=l+2,jcap+1
if(mod(n+l,2).eq.0) then
indev=indev+1
! xlsum2 start
aa = sinlat / epse(indev)
bb = epso(indod) / epse(indev)
id = ialp2 - ialp1
if( id.eq.0 ) then
alp3 = aa*alp2 - bb*alp1
ialp3 = ialp1
elseif( id.eq.1 ) then
alp3 = aa*alp2 - bb*alp1*bs_f
ialp3 = ialp2
elseif( id.eq.-1 ) then
alp3 = aa*alp2*bs_f - bb*alp1
ialp3 = ialp1
elseif( id.gt.1 ) then
alp3 = aa*alp2
ialp3 = ialp2
else
alp3 = - bb*alp1
ialp3 = ialp1
endif
! xlsum2 end
! xnorm alp3 start
w = abs(alp3)
if( w.ge.bb_h ) then
alp3 = alp3*bs_f
ialp3 = ialp3 + 1
elseif( w.lt.bs_h ) then
alp3 = alp3*bb_f
ialp3 = ialp3 - 1
endif
! xnorm alp3 end
! x2f alp3 start
if( ialp3.eq.0 ) then
plnev_a(indev,lat)=alp3
elseif( ialp3.eq.-1 ) then
plnev_a(indev,lat)=alp3 * bs_f
elseif( ialp3.lt.-1 ) then
plnev_a(indev,lat)=0.0
else
plnev_a(indev,lat)=alp3 * bb_f
endif
! x2f alp3 end
else
indod=indod+1
! xlsum2 start
aa = sinlat / epso(indod)
bb = epse(indev) / epso(indod)
id = ialp2 - ialp1
if( id.eq.0 ) then
alp3 = aa*alp2 - bb*alp1
ialp3 = ialp1
elseif( id.eq.1 ) then
alp3 = aa*alp2 - bb*alp1*bs_f
ialp3 = ialp2
elseif( id.eq.-1 ) then
alp3 = aa*alp2*bs_f - bb*alp1
ialp3 = ialp1
elseif( id.gt.1 ) then
alp3 = aa*alp2
ialp3 = ialp2
else
alp3 = - bb*alp1
ialp3 = ialp1
endif
! xlsum2 end
! xnorm alp3 start
w = abs(alp3)
if( w.ge.bb_h ) then
alp3 = alp3*bs_f
ialp3 = ialp3 + 1
elseif( w.lt.bs_h ) then
alp3 = alp3*bb_f
ialp3 = ialp3 - 1
endif
! xnorm alp3 end
! x2f alp3 start
if( ialp3.eq.0 ) then
plnod_a(indod,lat)=alp3
elseif( ialp3.eq.-1 ) then
plnod_a(indod,lat)=alp3 * bs_f
elseif( ialp3.lt.-1 ) then
plnod_a(indod,lat)=0.0
else
plnod_a(indod,lat)=alp3 * bb_f
endif
! x2f alp3 end
endif
alp1=alp2
alp2=alp3
ialp1 = ialp2
ialp2 = ialp3
enddo
cc
enddo
cc
enddo
cc
return
end
end module pln2eo_a_stochy_mod