CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cires_ugwpv1_sporo.F90
1
3
6 contains
7
9 subroutine oro_spectral_solver(im, levs,npt,ipt, kref,kdt,me,master, &
10 dtp,dxres, taub, u1, v1, t1, xn, yn, bn2, rho, prsi, prsL, &
11 del, sigma, hprime, gamma, theta, &
12 sinlat, xlatd, taup, taud, pkdis)
13!
14 use machine , only : kind_phys
15 use ugwp_common, only : grav, omega2, rd
16!
17 implicit none
18
19 integer, intent(in) :: im, levs
20 integer, intent(in) :: npt
21 integer, intent(in) :: kdt, me, master
22 integer, intent(in) :: kref(im), ipt(im)
23
24 real(kind=kind_phys), intent(in) :: dtp, dxres
25 real(kind=kind_phys), intent(in) :: taub(im)
26
27 real(kind=kind_phys), intent(in) :: sinlat(im), xlatd(im)
28 real(kind=kind_phys), intent(in), dimension(im) :: sigma, &
29 hprime, gamma, theta
30
31 real(kind=kind_phys), intent(in), dimension(im) :: xn, yn
32
33 real(kind=kind_phys), intent(in), dimension(im, levs) :: &
34 u1, v1, t1, bn2, rho, prsl, del
35
36 real(kind=kind_phys), intent(in), dimension(im, levs+1) :: prsi
37!
38! out : taup, taud, pkdis
39!
40 real(kind=kind_phys), intent(inout), dimension(im, levs+1) :: taup
41 real(kind=kind_phys), intent(inout), dimension(im, levs) :: taud
42 real(kind=kind_phys), intent(inout), dimension(im, levs) :: pkdis
43!
44! multiwave oro-spectra
45! locals
46!
47
48 integer, parameter :: nworo = 30
49 real(kind=kind_phys), parameter :: fc_flag = 0.0
50 real(kind=kind_phys), parameter :: mkzmin = 6.28e-3/50.0
51 real(kind=kind_phys), parameter :: mkz2min = mkzmin* mkzmin
52 real(kind=kind_phys), parameter :: kedmin = 1.e-3
53 real(kind=kind_phys), parameter :: kedmax = 350.,axmax=250.e-5
54 real(kind=kind_phys), parameter :: rtau = 0.01 ! nonlin-OGW scale 1/10sec
55 real(kind=kind_phys), parameter :: linsat2 =0.5
56 real(kind=kind_phys), parameter :: kxmin = 6.28e-3/100.
57 real(kind=kind_phys), parameter :: kxmax = 6.28e-3/5.0
58 real(kind=kind_phys), parameter :: dkx = (kxmax -kxmin)/(nworo-1)
59 real(kind=kind_phys), parameter :: kx_slope= -5./3.
60 real(kind=kind_phys), parameter :: hps =7000., rhp2 = .5/hps
61 real(kind=kind_phys), parameter :: cxmin=0.5, cxmin2=cxmin*cxmin
62
63 real(kind=kind_phys) :: akx(nworo), cxoro(nworo), akx2(nworo)
64 real(kind=kind_phys) :: aspkx(nworo), c2f2(nworo), cdf2(nworo)
65 real(kind=kind_phys) :: tau_sp(nworo,levs+1), wkdis(nworo, levs+1)
66 real(kind=kind_phys) :: tau_kx(nworo),taub_kx(nworo)
67
68 real(kind=kind_phys), dimension(nworo, levs+1) :: wrms, akzw
69
70 real(kind=kind_phys) :: tauz(levs+1), rms_wind(levs+1)
71 real(kind=kind_phys) :: wave_act(nworo,levs+1)
72
73 real(kind=kind_phys) :: kxw, kzw, kzw2, kzw3, kzi, dzmet, rhoint
74 real(kind=kind_phys) :: rayf, kturb
75 real(kind=kind_phys) :: uz, bv, bv2,kxsp, fcor2, cf2
76
77 real(kind=kind_phys) :: fdis
78 real(kind=kind_phys) :: wfdm, wfdt, wfim, wfit
79 real(kind=kind_phys) :: betadis, betam, betat, kds, cx, rhofac
80 real(kind=kind_phys) :: etwk, etws, tauk, cx2sat
81 real(kind=kind_phys) :: cdf1, tau_norm
82!
83! mean flow
84!
85 real(kind=kind_phys), dimension(levs+1) :: uzi,rhoi,ktur, kalp, dzi
86 real(kind=kind_phys) :: belps, aelps, nhills, selps
87 integer :: i, j, k, isp, iw
88 integer :: nw, nzi, ksrc
89
90
91 taud(:, :) = 0.0 ; pkdis(:,:) = 0.0 ; taup(:,:) = 0.0
92 tau_sp(:,:) = 0.0 ; wrms(:,:) = 0.0
93 nw = nworo
94 nzi = levs+1
95
96 do iw = 1, nw
97! !kxw = 0.25/(dxres)*iw
98 kxw = kxmin+(iw-1)*dkx
99 akx(iw) = kxw
100 akx2(iw) = kxw*kxw
101 aspkx(iw) = kxw ** (kx_slope)
102 tau_kx(iw) = aspkx(iw)*dkx
103 enddo
104
105 tau_norm = sum(tau_kx)
106 tau_kx(:) = tau_kx(:)/tau_norm
107
108 if (kdt == 1) then
109 write(6,771) maxval(tau_kx)*maxval(taub)*1.e3, minval(tau_kx), maxval(tau_kx)
110 endif
111771 format( ' oro_spectral_solver ', 3(2x,f8.3))
112!
113! main loop over oro-points
114!
115 do i =1, npt
116 j = ipt(i)
117
118!
119! estimate "nhills" => stochastic choices for OGWs
120!
121 if (taub(i) > 0.) then
122!
123! max_kxridge =min( .5*sigma(j)/hprime(j), kmax)
124! ridge-dependent dkx = (max_kxridge -kxmin)/(nw-1)
125! option to make grid-box variable kx-spectra kxw = kxmin+(iw-1)*dkx
126!
127 wave_act(1:nw, 1:levs+1) = 1.0
128 ksrc = kref(i)
129 tauz(1:ksrc) = taub(i)
130 taub_kx(1:nw) = tau_kx(1:nw) * taub(i)
131 wkdis(:,:) = kedmin
132
133 call oro_meanflow(levs, nzi, u1(j,:), v1(j,:), t1(j,:), &
134 & prsi(j,:), prsl(j,:), &
135 & del(j,:), rho(i,:), &
136 & bn2(i,:), uzi, rhoi,ktur, kalp,dzi, &
137 & xn(i), yn(i))
138
139 fcor2 = omega2*sinlat(j)*omega2*sinlat(j)*fc_flag
140
141 k = ksrc
142
143 bv2 = bn2(i,k)
144 uz = uzi(k) !u1(j,ksrc)*xn(i)+v1(j,ksrc)*yn(i)!
145 kturb = ktur(k)
146 rayf = kalp(k)
147 rhoint = rhoi(k)
148 dzmet = dzi(k)
149 kzw = max(sqrt(bv2)/max(cxmin, uz), mkzmin)
150!
151! specify oro-kx spectra and related variables k=ksrc
152!
153 do iw = 1, nw
154 kxw = akx(iw)
155 cxoro(iw) = 0.0 - uz
156 c2f2(iw) = fcor2/akx2(iw)
157 wrms(iw,k)= taub_kx(iw)/rhoint*kzw/kxw
158 tau_sp(iw, k) = taub_kx(iw)
159!
160!
161 if (cxoro(iw) > cxmin) then
162 wave_act(iw,k:levs+1) = 0. ! crit-level
163 else
164 cdf2(iw) = cxoro(iw)*cxoro(iw) -c2f2(iw)
165 if ( cdf2(iw) < cxmin2) then
166 wave_act(iw,k:levs+1) = 0. ! coriolis cut-off
167 else
168 kzw2 = max(bv2/cdf2(iw) - akx2(iw), mkz2min)
169 kzw = sqrt(kzw2)
170 akzw(iw,k)= kzw
171 wrms(iw,k)= taub_kx(iw)/rhoint * kzw/kxw
172 endif
173 endif
174 enddo ! nw-spectral loop
175!
176! defined abobe, k = ksrc: akx(nworo), cxoro(nworo), tau_sp(ksrc, nworo)
177! propagate upward multiwave-spectra are filtered by dissipation & instability
178!
179! tau_sp(:,ksrc+1:levs+1) = tau_sp(:, ksrc)
180 do k= ksrc+1, levs
181 uz = uzi(k)
182 bv2 =bn2(i,k)
183 bv = sqrt(bv2)
184 rayf = kalp(k)
185 rhoint= rhoi(k)
186 dzmet = dzi(k)
187 rhofac = rhoi(k-1)/rhoi(k)
188
189 do iw = 1, nworo
190!
191 if (wave_act(iw, k-1) <= 0.0) cycle
192 cxoro(iw)= 0.0 - uz
193 if ( cxoro(iw) > cxmin) then
194 wave_act(iw,k:levs+1) = 0.0 ! crit-level
195 else
196 cdf2(iw) = cxoro(iw)*cxoro(iw) -c2f2(iw)
197 if ( cdf2(iw) < cxmin2) wave_act(iw,k:levs+1) = 0.0
198 endif
199 if ( wave_act(iw,k) <= 0.0) cycle
200!
201! upward propagation
202!
203 kzw2 = bv2/cdf2(iw) - akx2(iw)
204
205 if (kzw2 < mkz2min) then
206 wave_act(iw,k:levs+1) = 0.0
207 else
208!
209! upward propagation w/o reflection effects
210!
211 kxw = akx(iw)
212 kzw = sqrt(kzw2)
213 akzw(iw,k) = kzw
214 kzw3 = kzw2*kzw
215
216 cx = cxoro(iw)
217 betadis = cdf2(iw) / (cx*cx+c2f2(iw))
218 betam = 1.0 / (1.0+betadis)
219 betat = 1.0 - betam
220 kds = wkdis(iw,k-1)
221
222 etws = wrms(iw,k-1)*rhofac * kzw/akzw(iw,k-1)
223
224 kturb = ktur(k)+pkdis(j,k-1)
225 wfim = kturb*kzw2 +rayf
226 wfit = wfim ! do updates with Pr-numbers Kv/Kt
227 cdf1 = sqrt(cdf2(iw))
228 wfdm = wfim/(kxw*cdf1)*betam
229 wfdt = wfit/(kxw*cdf1)*betat
230 kzi = 2.*kzw*(wfdm+wfdt)*dzmet
231 fdis = exp(-kzi)
232
233 etwk = etws*fdis
234 cx2sat = linsat2*cdf2(iw)
235
236 if (etwk > cx2sat) then
237 kds = kxw*cdf1*rhp2/kzw3
238 etwk = cx2sat
239 wfim = kds*kzw2
240 wfdm = wfim/(kxw*cdf1)
241 kzi = 2.*kzw*(wfdm + wfdm)*dzmet
242 etwk = cx2sat*exp(-kzi)
243 endif
244! if( lat(j) eq 40.5 ) then stop
245 wkdis(iw,k) = kds
246 wrms(iw,k) = etwk
247 tauk = etwk*kxw/kzw
248 tau_sp(iw,k) = tauk *rhoint
249 if ( tau_sp(iw,k) > tau_sp(iw,k-1)) &
250 tau_sp(iw,k) = tau_sp(iw,k-1)
251
252 ENDIF ! upward
253 ENDDO ! spectral
254
255!......... do spectral sum of rms, wkdis, tau
256
257 tauz(k) = sum( tau_sp(:,k)*wave_act(:,k) )
258 rms_wind(k) = sum( wrms(:,k)*wave_act(:,k) )
259
260 pkdis(j,k) = sum(wkdis(:,k)*wave_act(:,k))+rms_wind(k)*rtau
261
262 if (pkdis(j,k) > kedmax) pkdis(j,k) = kedmax
263
264 ENDDO ! k=ksrc+1, levs
265
266 k = ksrc
267 tauz(k) = sum(tau_sp(:,k)*wave_act(:,k))
268 tauz(k) = tauz(k+1) ! zero momentum dep-n at k=ksrc
269
270 pkdis(j,k) = sum(wkdis(:,k)*wave_act(:,k))
271 rms_wind(k) = sum(wrms(:,k)*wave_act(:,k))
272 tauz(levs+1) = tauz(levs)
273 taup(i, 1:levs+1) = tauz(1:levs+1)
274
275 do k=ksrc, levs
276 taud(i,k) = ( tauz(k+1) - tauz(k))*grav/del(j,k)
277!
278! limiters can be applied to avoid "large" wave accelerations
279!
280! if (taud(i,k) .gt. 0)taud(i,k)=taud(i,k)*.01
281! if (abs(taud(i,k)).ge.axmax)taud(i,k)=sign(taud(i,k),axmax)
282 enddo
283 endif ! taub > 0
284 enddo ! oro-points (i, j, ipt)
285!
286 end subroutine oro_spectral_solver
287!-------------------------------------------------------------
288!
290 subroutine oro_meanflow(nz, nzi, u1, v1, t1, pint, pmid, &
291 & delp, rho, bn2, uzi, rhoi, ktur, kalp, dzi, xn, yn)
292 use machine , only : kind_phys
293 use ugwp_common , only : velmin, dw2min, rdi, grav, rgrav, hpscale, rhp, rh4
294 implicit none
295
296 integer :: nz, nzi
297 real(kind=kind_phys), dimension(nz ) :: u1, v1, t1, delp, rho, pmid
298 real(kind=kind_phys), dimension(nz ) :: bn2 ! define at the interfaces
299 real(kind=kind_phys), dimension(nz+1) :: pint
300 real(kind=kind_phys) :: xn, yn
301
302! output
303
304 real(kind=kind_phys), dimension(nz+1) :: dzi, uzi, rhoi, ktur, kalp
305
306! locals
307 integer :: i, j, k
308 real(kind=kind_phys) :: ui, vi, ti, uz, vz, shr2, rdz, kamp
309 real(kind=kind_phys) :: zgrow, zmet, rdpm, ritur, kmol, w1
310
311! paremeters
312! real(kind=kind_phys), parameter :: hps = 7000., rpspa = 1.e-5
313! real(kind=kind_phys), parameter :: rhps=1.0/hps
314! real(kind=kind_phys), parameter :: h4= 0.25/hps
315
316 real(kind=kind_phys), parameter :: rimin = 0.125, kedmin = 0.01
317 real(kind=kind_phys), parameter :: lturb = 30. , uturb = 150.0
318 real(kind=kind_phys), parameter :: lsc2 = lturb*lturb,usc2 = uturb*uturb
319
320 kalp(1:nzi) = 2.e-7 ! radiative damping scale
321
322 do k=2, nz
323 rdpm = grav/(pmid(k-1)-pmid(k))
324 ui = .5*(u1(k-1)+u1(k))
325 vi = .5*(v1(k-1)+v1(k))
326 uzi(k) = ui*xn + vi*yn
327 ti = .5*(t1(k-1)+t1(k))
328 rhoi(k) = rdi*pint(k)/ti
329 rdz = rdpm *rhoi(k)
330 dzi(k) = 1./rdz
331 uz = u1(k)-u1(k-1)
332 vz = v1(k)-v1(k-1)
333 shr2 = rdz*rdz*(max(uz*uz+vz*vz, dw2min))
334 zmet = -hpscale*alog(pint(k)*1.e-5)
335 zgrow = exp(zmet*rh4)
336 kmol = 2.e-5*exp(zmet*rhp) + kedmin
337 ritur = max(bn2(k)/shr2, rimin)
338 kamp = sqrt(shr2)*lsc2 *zgrow
339 w1 = 1./(1. + 5*ritur)
340 ktur(k) = kamp * w1 * w1 + kmol
341 enddo
342
343 k = 1
344 uzi(k) = uzi(k+1)
345 ktur(k) = ktur(k+1)
346 rhoi(k) = rdi*pint(k)/t1(k+1)
347 dzi(k) = rgrav*delp(k)/rhoi(k)
348
349 k = nzi
350 uzi(k) = uzi(k-1)
351 ktur(k) = ktur(k-1)
352 rhoi(k) = rhoi(k-1)*.5
353 dzi(k) = dzi(k-1)
354
355 end subroutine oro_meanflow
356
357 end module cires_ugwpv1_sporo