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)
14 use machine ,
only : kind_phys
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)
24 real(kind=kind_phys),
intent(in) :: dtp, dxres
25 real(kind=kind_phys),
intent(in) :: taub(im)
27 real(kind=kind_phys),
intent(in) :: sinlat(im), xlatd(im)
28 real(kind=kind_phys),
intent(in),
dimension(im) :: sigma, &
31 real(kind=kind_phys),
intent(in),
dimension(im) :: xn, yn
33 real(kind=kind_phys),
intent(in),
dimension(im, levs) :: &
34 u1, v1, t1, bn2, rho, prsl, del
36 real(kind=kind_phys),
intent(in),
dimension(im, levs+1) :: prsi
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
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
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
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)
68 real(kind=kind_phys),
dimension(nworo, levs+1) :: wrms, akzw
70 real(kind=kind_phys) :: tauz(levs+1), rms_wind(levs+1)
71 real(kind=kind_phys) :: wave_act(nworo,levs+1)
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
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
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
91 taud(:, :) = 0.0 ; pkdis(:,:) = 0.0 ; taup(:,:) = 0.0
92 tau_sp(:,:) = 0.0 ; wrms(:,:) = 0.0
98 kxw = kxmin+(iw-1)*dkx
101 aspkx(iw) = kxw ** (kx_slope)
102 tau_kx(iw) = aspkx(iw)*dkx
105 tau_norm = sum(tau_kx)
106 tau_kx(:) = tau_kx(:)/tau_norm
109 write(6,771) maxval(tau_kx)*maxval(taub)*1.e3, minval(tau_kx), maxval(tau_kx)
111771
format(
' oro_spectral_solver ', 3(2x,f8.3))
121 if (taub(i) > 0.)
then
127 wave_act(1:nw, 1:levs+1) = 1.0
129 tauz(1:ksrc) = taub(i)
130 taub_kx(1:nw) = tau_kx(1:nw) * taub(i)
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, &
139 fcor2 = omega2*sinlat(j)*omega2*sinlat(j)*fc_flag
149 kzw = max(sqrt(bv2)/max(cxmin, uz), mkzmin)
156 c2f2(iw) = fcor2/akx2(iw)
157 wrms(iw,k)= taub_kx(iw)/rhoint*kzw/kxw
158 tau_sp(iw, k) = taub_kx(iw)
161 if (cxoro(iw) > cxmin)
then
162 wave_act(iw,k:levs+1) = 0.
164 cdf2(iw) = cxoro(iw)*cxoro(iw) -c2f2(iw)
165 if ( cdf2(iw) < cxmin2)
then
166 wave_act(iw,k:levs+1) = 0.
168 kzw2 = max(bv2/cdf2(iw) - akx2(iw), mkz2min)
171 wrms(iw,k)= taub_kx(iw)/rhoint * kzw/kxw
187 rhofac = rhoi(k-1)/rhoi(k)
191 if (wave_act(iw, k-1) <= 0.0) cycle
193 if ( cxoro(iw) > cxmin)
then
194 wave_act(iw,k:levs+1) = 0.0
196 cdf2(iw) = cxoro(iw)*cxoro(iw) -c2f2(iw)
197 if ( cdf2(iw) < cxmin2) wave_act(iw,k:levs+1) = 0.0
199 if ( wave_act(iw,k) <= 0.0) cycle
203 kzw2 = bv2/cdf2(iw) - akx2(iw)
205 if (kzw2 < mkz2min)
then
206 wave_act(iw,k:levs+1) = 0.0
217 betadis = cdf2(iw) / (cx*cx+c2f2(iw))
218 betam = 1.0 / (1.0+betadis)
222 etws = wrms(iw,k-1)*rhofac * kzw/akzw(iw,k-1)
224 kturb = ktur(k)+pkdis(j,k-1)
225 wfim = kturb*kzw2 +rayf
227 cdf1 = sqrt(cdf2(iw))
228 wfdm = wfim/(kxw*cdf1)*betam
229 wfdt = wfit/(kxw*cdf1)*betat
230 kzi = 2.*kzw*(wfdm+wfdt)*dzmet
234 cx2sat = linsat2*cdf2(iw)
236 if (etwk > cx2sat)
then
237 kds = kxw*cdf1*rhp2/kzw3
240 wfdm = wfim/(kxw*cdf1)
241 kzi = 2.*kzw*(wfdm + wfdm)*dzmet
242 etwk = cx2sat*exp(-kzi)
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)
257 tauz(k) = sum( tau_sp(:,k)*wave_act(:,k) )
258 rms_wind(k) = sum( wrms(:,k)*wave_act(:,k) )
260 pkdis(j,k) = sum(wkdis(:,k)*wave_act(:,k))+rms_wind(k)*rtau
262 if (pkdis(j,k) > kedmax) pkdis(j,k) = kedmax
267 tauz(k) = sum(tau_sp(:,k)*wave_act(:,k))
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)
276 taud(i,k) = ( tauz(k+1) - tauz(k))*grav/del(j,k)
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
297 real(kind=kind_phys),
dimension(nz ) :: u1, v1, t1, delp, rho, pmid
298 real(kind=kind_phys),
dimension(nz ) :: bn2
299 real(kind=kind_phys),
dimension(nz+1) :: pint
300 real(kind=kind_phys) :: xn, yn
304 real(kind=kind_phys),
dimension(nz+1) :: dzi, uzi, rhoi, ktur, kalp
308 real(kind=kind_phys) :: ui, vi, ti, uz, vz, shr2, rdz, kamp
309 real(kind=kind_phys) :: zgrow, zmet, rdpm, ritur, kmol, w1
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
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
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
346 rhoi(k) = rdi*pint(k)/t1(k+1)
347 dzi(k) = rgrav*delp(k)/rhoi(k)
352 rhoi(k) = rhoi(k-1)*.5