CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cires_ugwpv1_oro.F90
1
3
6contains
7
8 subroutine orogw_v1 (im, km, imx, me, master, dtp, kdt, do_tofd, &
9 xlatd, sinlat, coslat, sparea, &
10 cdmbgwd, hprime, oc, oa4, clx4, theta, sigmad, &
11 gammad, elvmaxd, sgh30, kpbl, &
12 u1 ,v1, t1, q1, prsi,del,prsl,prslk, zmeti, zmet, &
13 pdvdt, pdudt, pdtdt, pkdis, dusfc, dvsfc,rdxzb , &
14 zobl, zlwb, zogw, tau_ogw, dudt_ogw, dvdt_ogw, &
15 dudt_obl, dvdt_obl,dudt_ofd, dvdt_ofd, &
16 du_ogwcol, dv_ogwcol, du_oblcol, dv_oblcol, &
17 du_ofdcol, dv_ofdcol, errmsg,errflg )
18
19!---------------------------------------------------------------------------
20! ugwp_v1: orogw_v1 following recent updates of Lott & Miller 1997
21! eventually will be replaced with more "advanced" LLWB
22! and multi-wave solver that produce competitive FV3GFS-skills
23!
24! computation of kref for ogw + coorde diagnostics
25! all constants/parameters inside cires_ugwp_initialize.f90
26!
27! 10/2020 main updates
28! (a) introduce extra diagnostics of x-y obl-ofd-ogw as in the GSL-drag
29! for intercomparisons
30!
31! (b) quit with cdmbgwd(1:2)
32! cdmbgwd(1) = 1 for all resolutions, number of hills control SA-effects
33! cdmbgwd(2) = 1 ...............number of hills control SA-effects
34!
35! (c) cleff = pi2/(nk*dx) lheff = nk*dx (nk = 6,4,2, 1)
36! alternative lheff = min( dogw=hprime/sigma*gamma, dx)
37! we still not use the "broad spectral solver"
38!
39! (d) hefff = (nsig * hprime -znlk)/nsig, orchestrating MB and OGW
40!
41! (e) for linsat-solver the total "eddy" damping Ked = Ked * Nhills,
42! scale-aware amplification of the momentum deposition for low-res runs
43!----------------------------------------
44
45 use machine , only : kind_phys
46 use ugwp_common, only : dw2min, velmin, grav, omega1, rd, cpd, rv, pi, arad, fv
47 use ugwp_common, only : rcpdt, grav2, rgrav, rcpd, rcpd2
48 use ugwp_common, only : rad_to_deg, deg_to_rad, pi2, pih, rdi, gor, grcp, gocp, gr2, bnv2min
49
50 use ugwp_oro_init, only : rimin, ric, efmin, efmax, &
51 hpmax, hpmin, sigfaci => sigfac, &
52 dpmin, minwnd, hminmt, hncrit, &
53 rlolev, gmax, veleps, factop, &
54 frc, ce, ceofrc, frmax, cg, &
55 fdir, mdir, nwdir, &
56 cdmb, cleff, fcrit_v1, &
57 n_tofd, ze_tofd, ztop_tofd
58
59 use cires_ugwpv1_module, only : kxw, max_kdis, max_axyz
60
61!----------------------------------------
62 implicit none
63!----------------------------------------
64! internal parameters
65!----------------------------------------
66 real(kind=kind_phys), parameter :: sigfac = 3 ! N*hprime height of Subgrid Hill over which SSO-flo
67 real(kind=kind_phys), parameter :: sigfacs = 0.25 ! M*hprime height is the low boundary of the hill
68
69 real(kind=kind_phys), parameter :: dbmax = 1./3600./12. ! max-Krmtb in hours for u=10 m/s => 20 m/s/day
70 character(len=8) :: strsolver='pss-1986' ! current operational Ri-solver or 'spect_2020'
71
72
73 real(kind=kind_phys) :: gammin = 0.00999999 ! a/b = gammma_min =1% <====>
74 real(kind=kind_phys), parameter :: nhilmax = 15. ! max number of SSO-hills in grid-box
75 real(kind=kind_phys), parameter :: sso_min = 3000. ! min-lenghth of the hill, GTOP30 ~dx~1 km
76
77 real(kind=kind_phys), parameter :: nfr = 2.+1. ! power in the emprical Function(Fr/Frc)
78 real(kind=kind_phys), parameter :: afr = 1. ! (Fr/Frc)^2/(afr +[Fr/Frc]^nfr), Fr = h*mkz
79 real(kind=kind_phys), parameter :: frnorm =afr+1.0 ! to get cont-ous taulin(Fr=Frc) = tau_nonlin(Fr=Frc) !
80 real(kind=kind_phys), parameter :: max_frf =2.0 ! max-value of non-lin flux over the linear at Fr=Frc
81
82 logical, parameter :: do_adjoro = .false. !
83!----------------------------------------
84
85 integer, intent(in) :: im, km, imx, kdt
86 integer, intent(in) :: me, master
87 logical, intent(in) :: do_tofd
88
89 integer, intent(in) :: kpbl(im) ! index for the pbl top layer!
90 real(kind=kind_phys), intent(in) :: dtp ! time step
91 real(kind=kind_phys), intent(in) :: cdmbgwd(2)
92
93 real(kind=kind_phys), intent(in) :: hprime(im), oc(im), oa4(im,4), &
94 clx4(im,4), theta(im), &
95 sigmad(im), gammad(im), elvmaxd(im)
96!
97 real(kind=kind_phys), intent(in) :: sgh30(im)
98
99 real(kind=kind_phys), intent(in), dimension(im,km) :: &
100 u1, v1, t1, q1,del, prsl, prslk, zmet
101
102 real(kind=kind_phys), intent(in),dimension(im,km+1):: prsi, zmeti
103
104 real(kind=kind_phys), intent(in) :: xlatd(im),sinlat(im), coslat(im)
105 real(kind=kind_phys), intent(in) :: sparea(im)
106
107!
108!output -phys-tend
109!
110 real(kind=kind_phys),dimension(im,km),intent(out) :: &
111 pdvdt, pdudt, pkdis, pdtdt
112! output - diag-coorde
113 real(kind=kind_phys),dimension(im,km),intent(out) :: &
114 dudt_ogw,dvdt_ogw, dudt_obl,dvdt_obl, dudt_ofd,dvdt_ofd
115
116 real(kind=kind_phys),dimension(im),intent(out) :: dusfc, dvsfc, &
117 du_ogwcol,dv_ogwcol, du_oblcol,dv_oblcol, du_ofdcol,dv_ofdcol
118!
119 real(kind=kind_phys),dimension(im),intent(out) :: rdxzb
120 real(kind=kind_phys),dimension(im),intent(out) :: zobl, zogw, zlwb, tau_ogw
121
122 character(len=*), intent(out) :: errmsg
123 integer, intent(out) :: errflg
124!
125!
126! locals vars for SSO
127!
128
129 real(kind=kind_phys), dimension(im) :: oa, clx
130 real(kind=kind_phys), dimension(im) :: sigma, gamma, elvmax ! corrected sigmaD, gammaD, elvmaxD
131
132 real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
133 real(kind=kind_phys) :: belpmin, dsmin, dsmax
134
135 real(kind=kind_phys) :: arhills(im), mkd05_hills(im) ! number of hills in the grid
136 real(kind=kind_phys) :: taub_kd05(im)
137!
138! locals mean flow ...etc
139!
140 real(kind=kind_phys), dimension(im,km) :: ri_n, bnv2, ro
141 real(kind=kind_phys), dimension(im,km) :: vtk, vtj, velco
142!==================
143!mtb
144!==================
145 real(kind=kind_phys) :: ztoph,zlowh,ph_blk, dz_blk
146 real(kind=kind_phys), dimension(im) :: wk, pe, ek, up
147
148 real(kind=kind_phys), dimension(im,km) :: db, ang, uds
149
150 real(kind=kind_phys) :: zlen, dbtmp, r, phiang, dbim, zr
151 real(kind=kind_phys) :: eng0, eng1, cosang2, sinang2
152 real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
153
154!==================
155! tofd
156! some constants now in "use ugwp_oro_init" + "use ugwp_common"
157!
158!==================
159
160 real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
161 real(kind=kind_phys), dimension(km) :: utofd1, vtofd1
162 real(kind=kind_phys), dimension(km) :: epstofd1, krf_tofd1
163 real(kind=kind_phys), dimension(km) :: up1, vp1, zpm
164!==================
165! ogw
166!==================
167 real(kind=kind_phys) :: xlingfs
168 logical :: icrilv(im)
169!
170 real(kind=kind_phys), dimension(im) :: xn, yn, ubar, vbar, ulow, &
171 roll, bnv2bar, scor, dtfac, xlinv, delks, delks1
172!
173 real(kind=kind_phys) :: taup(im,km+1), taud(im,km)
174 real(kind=kind_phys) :: taub(im), taulin(im), tausat(im), ahdxres(im)
175 real(kind=kind_phys) :: heff, hsat, hdis
176
177 integer, dimension(im) :: kref, idxzb, ipt, khtop, iwk, izlow
178!
179! local real scalars
180!
181 real(kind=kind_phys) :: bnv, fr, ri_gw, brvf, fr2
182 real(kind=kind_phys) :: tem, tem1, tem2, temc, temv
183 real(kind=kind_phys) :: ti, rdz, dw2, shr2, bvf2
184 real(kind=kind_phys) :: rdelks, efact, coefm, gfobnv
185 real(kind=kind_phys) :: scork, rscor, hd, fro, sira
186 real(kind=kind_phys) :: dtaux, dtauy, zmetp, zmetk
187
188 real(kind=kind_phys) :: windik, wdir
189 real(kind=kind_phys) :: sigmin, dxres,sigres,hdxres, cdmb4, mtbridge
190
191 real(kind=kind_phys) :: kxridge, inv_b2eff, zw1, zw2
192 real(kind=kind_phys) :: belps, aelps, nhills, selps
193
194! real(kind=kind_phys) :: rgrav, rcpd, rcpd2, rad_to_deg, deg_to_rad
195! real(kind=kind_phys) :: pi2, pih, rdi, gor, grcp, gocp, gr2, bnv2min
196
197 real(kind=kind_phys) :: cleff_max ! resolution-aware max-wn
198 real(kind=kind_phys) :: nonh_fact ! non-hydroststic factor 1.-(kx/kz_hh)**2
199 real(kind=kind_phys) :: fcrit2
200 real(kind=kind_phys) :: fr_func, frnd
201!
202!
203! local integers
204!
205 integer :: kmm1, kmm2, lcap, lcapp1
206 integer :: npt, kbps
207 integer :: kmps, idir, nwd, klcap, kp1, kmpbl, kmll
208 integer :: k_mtb, k_zlow, ktrial, klevm1
209 integer :: i, j, k
210!
211! Initialize CCPP error handling variables
212 errmsg = ''
213 errflg = 0
214!===========================
215! First step Check do we have sub-grid hills
216!
217!
218! out-arrays are zreoed in unified_ugwp.F90
219!
220 do i=1,im
221 rdxzb(i) = 0.0
222 dusfc(i) = 0.0
223 dvsfc(i) = 0.0
224 ipt(i) = 0
225 enddo
226
227! ---- for lm and gwd calculation points
228!cires_ugwp_initialize.F90: real, parameter :: hpmax=2400.0, hpmin=25.0
229!cires_ugwp_initialize.F90: real, parameter :: hminmt=50. ! min mtn height (*j*)
230!---- for lm and gwd calculation points
231! ccpp-gwdps.f PARAMETER (hpmax=2400.0, hpmin=1.0) parameter (elvmax > hminmt=50.)
232
233 npt = 0
234 do i = 1,im
235 if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin ) then
236 npt = npt + 1
237 ipt(npt) = i
238 endif
239 enddo
240
241 if (npt == 0) then
242! print *, 'orogw_v1 npt = 0 elvmax ', maxval(elvmaxd), hminmt
243! print *, 'orogw_v1 npt = 0 hprime ', maxval(hprime), hpmin
244 return ! no ogw/mbl calculation done
245 endif
246
247
248!=================================
249! Start if npt >= 1
250! initialize gamma and sigma for
251! performing the QC of SSO inputs
252!=================================
253 gamma(:) = gammad(:)
254 sigma(:) = sigmad(:)
255!
256!=======================================================================
257! mtb-blocking sigma_min and dxres => cires_initialize (best way ....)
258!
259 sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
260 dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
261
262! ! GTOP30-arc dx~1Km res-n so sso_hill ~ (2-4)*dx
263 cleff_max = pi2/max(dsmin/5.,sso_min) ! maxval for kx = 6.28/(dx_min/5. ~2.5 km) for C768
264 cleff_max = pi2/dsmin
265
266 hdxres = 0.5*dsmax
267
268 gammin = min(sso_min/hdxres, 1.)
269 gammin = max(0.1, gammin)
270 ! sigma-degined as tan(angle) = h/2: L/2= h/L
271 sigmin = hpmin/hdxres ! min-slope Hmin= 2*hpmin, dxres=Lmax
272
273! if ( kdt == 1 .and. me == master) then
274! print *, ' orogw_v1 scale2 ', cdmbgwd(2)
275! print *, ' orogw_v1 imx ', imx
276! print *, ' orogw_v1 gam_min ', gammin
277! print *, ' orogw_v1 sso_min ', sso_min
278! print *, ' orogw_v1 gam_min ', gammin
279! print *, ' orogw_v1 npt number of GRID-cells with hills ', npt
280! endif
281
282!============================================================
283! Purpose to adjust oro-specification on the fly
284! needs to be done 1-time during init-n for each block
285! hprime sigma gamma and grid-length must be "related"
286! width_mount_a = hprime/sigma < dxres cannot access dxres
287! width_mount_b = width_mount_a * gamma
288!
289! Sellipse= pi * a*b = (width_mount_a)^2 *gamma <= Sarea
290! Limiters on "elongated" hills gamma= a/b < gam_min
291! Limiters on "longest" hills (b, a) <= sqrt(area)
292!
293! 0.01=gammin < gamma=a_hill/b_hill < 1
294! hpmin/(dx/2)=sigmin < sigma= hprime/a_ell < 1.
295! Nhills = (dx*dy=Sarea)/(pi* a_hill *b_hill)
296!=============================================================
297
298 arhills(:) =0.
299 mkd05_hills(:) =0.
300
301 do j = 1,npt
302 i = ipt(j)
303 dxres = sqrt(sparea(i))
304 ahdxres(j) = dxres
305 if (gamma(i) > 1.0) gamma(i) = 1.0
306
307 gamma(i) = max(gammin, gamma(i))
308!
309! min-adjustment: 1) abs(gamma(i)) ; 2) sigres = max(sigmin, sigma(i))
310!
311 sigres = max(sigmin, sigma(i))
312 sigma(i) =sigres
313 aelps = min( hprime(i)/sigres, dxres)
314 belps = min(aelps/abs(gamma(i)), dxres)
315 gamma(i) = aelps/belps
316
317 if (do_adjoro ) then
318!
319! more adjustments "lengths", gamma and sigma, assuminng H_hill=2*hprime
320!
321 if (hprime(i) > hdxres*sigres) sigres= hprime(i)/dxres
322 aelps = min( hprime(i)/sigres, hdxres)
323 sigma(i) = sigres
324 if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i), dxres)
325!
326! small-scale "turbulent" oro-hills < sso_min, sso_min_dx = 3km
327! will be treated as "circular" elevations
328!
329 if( aelps < sso_min ) then
330!
331! a, b > sso_min upscale ellipse a/b > 0.1 a>sso_min & h/b=>new_sigm
332!
333 aelps = sso_min
334 if (belps < sso_min ) then
335 gamma(i) = 1.0
336 belps = aelps*gamma(i)
337 else
338 gamma(i) = min(aelps/belps, 1.0)
339 endif
340
341 sigma(i) = hprime(i)/aelps
342 gamma(i) = min(aelps/belps, 1.0)
343
344 endif !aelps < sso_min
345 endif ! if (do_adjoro )
346
347 selps = belps*belps*gamma(i)*pi ! area of the elliptical hill
348 nhills = min(nhilmax, sparea(i)/selps)
349 arhills(j) = max(nhills, 1.0)
350
351! if (kdt==1 ) write(6,333) nint(nhills)+1,xlatd(i), hprime(i),aelps*1.e-3, belps*1.e-3, sigma(i),gamma(i)
352! 333 format( ' nhil: ', i6, 4(2x, f9.3), 2(2x, e9.3))
353
354 enddo
355
356!=======================================================================
357! mtb-blocking : LM-1997; Zadra et al. 2004 ;metoffice dec 2010 H Wells
358!=======================================================================
359
360 do i=1,npt
361 khtop(i) = 2
362 idxzb(i) = 0
363 izlow(i) = 1
364 enddo
365
366 do k=1,km
367 do i=1,im
368 db(i,k) = 0.0
369 ang(i,k) = 0.0
370 uds(i,k) = 0.0
371 enddo
372 enddo
373
374 kmm1 = km - 1 ; kmm2 = km - 2 ; kmll = kmm1
375 lcap = km ; lcapp1 = lcap + 1
376 cdmb4 = 0.25*cdmb
377
378 do i = 1, npt
379 j = ipt(i)
380 elvmax(j) = min( sigfac * hprime(j), hncrit)
381!
382!gfsv15/16: ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit=8000.)
383! SSO-effects from the surface to "ELVMAX" =4*hprime + ELVMAX
384 enddo
385
386
387!===================================================================
388! below khtop-level H= 3*hp, and izlow = 0.5*Hp or the "first" layer
389! are used tp estimate "Mean" Flow that interact with SG-HILL
390! if sig*HP < Hpbl => GWs-> above PBL
391! WRF: ( 1 to max(2*Hp or H_pbl)
392! GFS-15/16: OGWs (1 to max(Kpbl+1, or K_dPs=(Ps-Pk=50hPa) ~ 950 mb)
393! excitation above Kref
394! BLOCKING: ZDOMAIN (1 - Kaver => ELVMAX(J) + sigfac * hp)
395!===================================================================
396
397
398 do k = 1, kmm1
399 do i = 1, npt
400 j = ipt(i)
401 ztoph = sigfac * hprime(j)
402 zlowh = sigfacs* hprime(j)
403 zmetp = zmet(j,k+1)
404 zmetk = zmet(j,k)
405!
406! GFSv15/16: izlow=1
407! elvmax(j)=elvmaxd(J) + sig*hp: if (( elvmax(j) <= zmetp) .and. (elvmax(j).ge.zmetk) ) khtop(i) = max(khtop(i), k+1 )
408!
409
410 if (( ztoph <= zmetp) .and. (ztoph >= zmetk) ) khtop(i) = max(khtop(i), k+1 )
411 if (zlowh <= zmetp .and. zlowh >= zmetk) izlow(i) = max(izlow(i),k)
412 enddo
413 enddo
414!
415 do k = 1,km
416 do i =1,npt
417 j = ipt(i)
418 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
419 vtk(i,k) = vtj(i,k) / prslk(j,k)
420 ro(i,k) = rdi * prsl(j,k) / vtj(i,k) ! density mid
421 taup(i,k) = 0.0
422 enddo
423 enddo
424!
425! perform ri_n or ri_mf computation for both OGW and OBL
426!23456
427 do k = 1,kmm1
428 do i =1,npt
429 j = ipt(i)
430 rdz = 1. / (zmet(j,k+1) - zmet(j,k))
431 tem1 = u1(j,k) - u1(j,k+1)
432 tem2 = v1(j,k) - v1(j,k+1)
433 dw2 = tem1*tem1 + tem2*tem2
434 shr2 = max(dw2,dw2min) * rdz * rdz
435 bvf2 = grav2 * rdz * (vtk(i,k+1)-vtk(i,k))/ (vtk(i,k+1)+vtk(i,k))
436
437 bnv2(i,k+1) = max( bvf2, bnv2min )
438 ri_n(i,k+1) = bnv2(i,k)/shr2 ! richardson number consistent with bnv2
439! having ri_n
440! we may place here computation for "ktur" and ogw-dissipation for the spectral ORO-scheme
441!
442 enddo
443 enddo
444 k = 1
445!23456
446 do i = 1, npt
447 bnv2(i,k) = bnv2(i,k+1)
448 enddo
449!
450! level khtop => zmet(j,k) < sigfac * hprime(j) < zmet(j,k+1)
451!23456
452 do i = 1, npt
453 j = ipt(i)
454 k_zlow = izlow(i)
455 if (k_zlow == khtop(i)) k_zlow = 1
456 delks(i) = 1.0 / (prsi(j,k_zlow) - prsi(j,khtop(i)))
457! delks1(i) = 1.0 /(prsl(j,k_zlow) - prsl(j,khtop(i)))
458 ubar(i) = 0.0
459 vbar(i) = 0.0
460 roll(i) = 0.0
461 pe(i) = 0.0
462 ek(i) = 0.0
463 bnv2bar(i) = 0.0
464!
465! computation of the mean flow char zlow < z < ztop =sigfac*hprime
466!23456
467 do k = k_zlow, khtop(i)-1
468 rdelks = del(j,k) * delks(i)
469 ubar(i) = ubar(i) + rdelks * u1(j,k)
470 vbar(i) = vbar(i) + rdelks * v1(j,k)
471 roll(i) = roll(i) + rdelks * ro(i,k)
472 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
473 enddo
474 enddo
475!23456
476 do i = 1, npt
477 j = ipt(i)
478!
479! integrate from ztoph = sigfac*hprime down to zblk if exists
480! find ph_blk, dz_blk as introduced in LM-97 and ifs
481!23456
482 ph_blk =0.
483 do k = khtop(i), 1, -1
484 phiang = atan2(v1(j,k),u1(j,k))
485 phiang = theta(j)*rad_to_deg - phiang
486 if ( phiang > pih ) phiang = phiang - pi
487 if ( phiang < -pih ) phiang = phiang + pi
488 ang(i,k) = phiang
489 uds(i,k) = max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), velmin)
490!23456
491 if (idxzb(i) == 0 ) then
492 dz_blk = zmeti(j,k+1) - zmeti(j,k)
493 pe(i) = pe(i) + bnv2(i,k) *( elvmax(j) - zmet(j,k) ) * dz_blk
494 up(i) = max(uds(i,k) * cos(ang(i,k)), velmin)
495 ek(i) = 0.5 * up(i) * up(i)
496 ph_blk = ph_blk + dz_blk*sqrt(bnv2(i,k))/up(i)
497
498! --- dividing stream lime is found when pe =exceeds ek. oper-l gfs
499! if ( pe(i) >= ek(i) ) then
500! --- LM97/IFS
501 if(ph_blk >= fcrit_v1 ) then
502 idxzb(i) = k
503 zobl(j) = zmet(j, k)
504 rdxzb(j) = real(k, kind=kind_phys)
505 endif
506!23456
507 endif
508 enddo
509!
510! --- the drag for the blocked flow
511!
512 if ( idxzb(i) > 0 ) then
513!
514! (4.16)-ifs description
515!
516 gam2 = gamma(j)*gamma(j)
517 bgam = 1.0 - 0.18*gamma(j) - 0.04*gam2
518 cgam = 0.48*gamma(j) + 0.30*gam2
519 do k = idxzb(i)-1, 1, -1
520!23456
521! empirical height dep-nt "blocking" length from LM-1997/IFS
522!
523 zlen = sqrt( (zobl(j)-zmet(j,k) )/(zmet(j,k ) + hprime(j)) )
524 tem = cos(ang(i,k))
525 cosang2 = tem * tem
526 sinang2 = 1.0 - cosang2
527 rdem = cosang2 + gam2 * sinang2
528 rnom = cosang2*gam2 + sinang2
529!
530! metoffice dec 2010
531! correction of H. Wells & A. Zadra for the
532! aspect ratio of the elliptical hill seen by mean flow
533!
534 rdem = max(rdem, 1.e-6)
535 r = sqrt(rnom/rdem)
536 zr = max( 2. - r, 0. )
537 sigres = max(sigmin, sigma(j))
538 mtbridge = zr * sigres*zlen / hprime(j)
539! dbtmp = cdmb4*mtbridge*max(cos(ang(i,k)), gamma(j)*sin(ang(i,k))) ! (4.15)-ifs
540 dbtmp = cdmb4*mtbridge*(bgam* cosang2 +cgam * sinang2) ! (4.16)-ifs
541!
542! linear damping due to OBL [1/sec]=[U/L_block_orthogonal]
543! more accurate along 2-axes of ellipse, here zr-factor is based on Phillips' analytics
544!
545 db(i,k)= dbtmp * uds(i,k)
546! if (db(i,k) > dbmax) print *, ' db > dbmax ', 1./db(i,k)/3600., uds(i,k)
547 db(i,k)= min(db(i,k), dbmax)
548 enddo
549!23456
550 endif
551 enddo
552!.............................
553!.............................
554! finish the mtn blocking
555!.............................
556!.............................
557!
558!--- OGW section
559!
560! scale cleff between im=384*2 and 192*2 for t126/t170 and t62
561! inside "cires_ugwp_initialize.f90" now
562!
563 kmpbl = km / 2
564 iwk(1:npt) = 2
565!
566! in meto/UK-scheme:
567! k_mtb = max(k_zmtb, k_n*hprime/2] to reduce diurnal variations in taub_ogw
568!23456
569 do k=3,kmpbl
570 do i=1,npt
571 j = ipt(i)
572 tem = (prsi(j,1) - prsi(j,k))
573 if (tem < dpmin) iwk(i) = k ! dpmin=50 mb from the surface
574 enddo
575 enddo
576!
577! in all cires-UGWP-schemes: zogw > zobl
578! in ugwp-v1: we consider option for htop ~ (2-3)*hprime > zmtb
579! the top of hill can be inside the PBL.... if kref = khtop
580!
581
582 kbps = 1
583 kmps = km
584 k_mtb = 1
585!23456
586 do i=1,npt
587 j = ipt(i)
588 k_mtb = max(1, idxzb(i))
589 ! WRF/GSL: kogw = max(kpbl, ktop=2*var)
590 kref(i) = max(iwk(i), kpbl(j)+1 ) ! reference level pbl or smt-else...Zogw= sigfac*Hprime
591 kref(i) = max(kref(i), khtop(i) ) ! khtop => sigfac*hprime
592!zogw > zobl
593 if (kref(i) <= k_mtb) kref(i) = k_mtb + 1 ! OGW-layer above the blocking height
594 kbps = max(kbps, kref(i))
595 kmps = min(kmps, kref(i))
596!
597 delks(i) = 1.0 / (prsi(j,k_mtb) - prsi(j,kref(i)))
598 ubar(i) = 0.0
599 vbar(i) = 0.0
600 roll(i) = 0.0
601 bnv2bar(i)= 0.0
602 enddo
603!23456=====================
604!
605!= we estimate MF-parameters from k= k_mtb to [kref~kpbl] > k_mtb
606!computation of the mean flow for zobl < z < ztop =sigfac*hprime inb GSL ztop =max(hpbl, ztop)
607!23456=====================
608 do i = 1,npt
609 k_mtb = max(1, idxzb(i))
610 do k = k_mtb,kbps !kbps = max(kref) kref = (kpbl+1, khtop)
611 if(k < kref(i)) then
612 j = ipt(i)
613 rdelks = del(j,k) * delks(i)
614 ubar(i) = ubar(i) + rdelks * u1(j,k)
615 vbar(i) = vbar(i) + rdelks * v1(j,k)
616 roll(i) = roll(i) + rdelks * ro(i,k)
617 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
618 endif
619 enddo
620 enddo
621!
622! orographic asymmetry parameters (oa), and (clx) [Kim & Arakawa Kim & Doyle]
623!23456
624 do i = 1,npt
625 j = ipt(i)
626 wdir = atan2(ubar(i),vbar(i)) + pi
627 idir = mod(nint(fdir*wdir),mdir) + 1
628 nwd = nwdir(idir)
629 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
630 clx(i) = clx4(j,mod(nwd-1,4)+1) ! number of "effective" hills (?) in the grid-box KA-95/KD-05
631!
632!GSL-drag ->identical to above
633!
634! wdir = atan2(ubar(i),vbar(i)) + pi
635! idir = mod(nint(fdir*wdir),mdir) + 1
636! nwd = nwdir(idir)
637! oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
638! ol(i) = ol4(i,mod(nwd-1,4)+1)
639!
640 dtfac(i) = 1.0
641 icrilv(i) = .false. ! initialize critical level control Logic
642 ulow(i) = max(sqrt(ubar(i)*ubar(i)+vbar(i)*vbar(i)),velmin)
643 xn(i) = ubar(i) / ulow(i)
644 yn(i) = vbar(i) / ulow(i)
645 enddo
646!23456
647 do k = 1, kmm1
648 do i = 1,npt
649 j = ipt(i)
650 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*xn(i)+ (v1(j,k)+v1(j,k+1))*yn(i))
651 enddo
652 enddo
653
654 do i = 1,npt
655 velco(i,km) = velco(i,kmm1)
656 enddo
657!
658!------------------------------------------------------------------------
659! v0/v1: incorporates modifications for kxridge and heff/hsat
660! and employs taulin for fr <=fcrit_v1
661! concept of "clipped" hill if zmtb > 0. is uded to make
662! the integrated "tau_sso = tau_ogw +tau_mtb" close to reanalysis
663! now it is still used the "single-orowave" along ulow-upwind
664!
665! in contrast ifs/meto/e-canada employ the 2-orthogonal wave (2otw) schemes of
666! it requires "aver angle" and wind projections on axes of ell-hill
667! with 2-stresses: taub_a/b as suggested by analytics of Phillips (1984)
668!------------------------------------------------------------------------
669
670 taub(:) = 0. ; taulin(:)= 0. ;taub_kd05 =0.
671 fcrit2 =fcrit_v1*fcrit_v1
672!
673! taub_oro as in KA-95/KD-05 GSL & EMC includes ALL waves (POGWs, Lee-rotors, etc...)
674! here taub represents mainly OGWs with nonh_fact = 1. -(kx/kz)**2
675! LLWB for Lee-rotors (downslope dynamics and flow-splitting ) is not considered
676!23456
677 do i = 1,npt
678 j = ipt(i)
679 bnv = sqrt( bnv2bar(i) )
680 heff = min(hprime(j),hpmax)
681 if( zobl(j) > 0.) heff = max(sigfac*heff-zobl(j), 0.)/sigfac
682 if (heff <= 0) cycle
683 zw1 = ulow(i)/bnv
684 hsat = fcrit_v1 *zw1
685 heff = min(heff, hsat)
686 fr = heff/zw1
687 fr = min(fr, frmax)
688 fr2 = fr*fr
689 zw2 = fr2 *oc(j) ! oc-values from 0 to 10 (GSL => 100)
690!
691! [Kim & Doyle, 2005]
692!
693 efact = (oa(i) + 2.) ** (ceofrc*fr) ! enhnancement factor due to the resonance ampification/downstream
694 efact = min( max(efact,efmin), efmax )
695 gfobnv = efact* gmax /(zw2 + cg) ! withoutt "multiplication" on <fr2*oc>=zw2
696!
697! ! cleff_max(C768 = 6.28/(12.5 km/5.)) .....
698! xlinv(i) = min(coefm * cleff, cleff_max)
699!
700 mkd05_hills(i) = (1. + clx(i)) ** (oa(i)+1.) ! ex-coefm (1-2) of eff-hills with "some" anizoropy as in KD-2005
701 xlinv(i) = min(cleff * mkd05_hills(i), cleff_max)
702 taub_kd05(i) = roll(i)*xlinv(i) *(gfobnv *zw2)* (zw1 * ulow(i)* ulow(i))
703!
704! old: tem = fr2*oc(j) ; gfobnv = gmax * tem / ((tem + cg)*bnv(i))
705! kx =or max(kxridge, inv_b2eff) ! 6.28/lx ..0.5*sigma(j)/heff = 1./lridge
706!
707 sigres = sigma(j)
708 inv_b2eff = pi*sigres/heff ! pi2/(2b)
709 kxridge = pi /ahdxres(i) ! pi2/(2*dx)
710 xlingfs = max(inv_b2eff, kxridge)
711 nonh_fact = 1. - xlinv(i)*zw1 * xlinv(i)*zw1 ! 1- (kx/kz)^2
712!23456
713 if (nonh_fact <= 0.) cycle ! non-hydrostatic trapping kx = kz = N/U
714!
715 taulin(i) = xlinv(i)*roll(i)*bnv*ulow(i)*heff*heff * nonh_fact
716 tausat(i) = xlinv(i)*roll(i)* zw1*ulow(i)*ulow(i) *fcrit2 * nonh_fact
717!
718! taulin(i) = xlinv(i)*roll(i)*ulow(i)**3/bnv *fr2 =>
719! fr2 = (bnv*heff/Ulow)**2 + non-hydrostatic trapping effects Fr2_nh = Fr2 - kx2*heff^2
720!23456
721 if(fr > fcrit_v1 ) then
722 frnd = fr/fcrit_v1
723 fr_func = frnorm* frnd*frnd/(afr + frnd ** nfr)
724 taub(i) = tausat(i) *max(fr_func, max_frf) ! nonlinear flux tau0...xlinv(i)
725 else
726 taub(i) = taulin(i) ! linear flux for fr <= fcrit_v1
727 endif
728 xlinv(i) = .5*xlinv(i) ! 1/2 rhoint-factor in Ri-solver of PSS-1986
729 k = max(1, kref(i)-1)
730 tem = max(velco(i,k)*velco(i,k), dw2min)
731 scor(i) = bnv2(i,k) / tem ! scorer parameter below kref level Bn2/U2= kz^2
732 zogw(j) = zmeti(j, kref(i) )
733 tau_ogw(j) = taub(i)
734!23456
735 enddo
736!
737!----set up bottom values of stress
738!
739 do i = 1,npt
740 taup(i, 1:kref(i) ) = taub(i)
741 enddo
742!======================================================
743!
744! Having : taub(i)/tau_ogw(j) => solve for OGW-effects
745!
746!======================================================
747 if (strsolver == 'pss-1986') then
748
749!======================================================
750! v0-gfs orogw-solver of Palmer et al 1986 -"pss-1986"
751! modified by KD05 with the emp.expression (11):below k=kref ???
752! tau(k+1) = tau(k)*Scorer(K+1)/Scorer(K)
753! in v1-orogw linsatdis of "wam-2017" for
754! rotational/non-hydrostat ogws; important for
755! highres-fv3gfs with dx < 10 km
756!23456======================================================
757 do k = kmps, kmm1 ! vertical level loop from min(kref)
758 kp1 = k + 1
759 do i = 1, npt
760 if (k >= kref(i)) then
761 icrilv(i) = icrilv(i) .or. ( ri_n(i,k) < ric).or. (velco(i,k) <= 0. )
762 endif
763 enddo
764!
765 do i = 1,npt
766 if (k >= kref(i)) then
767 if (.not.icrilv(i) .and. taup(i,k) > 0.0 ) then
768 zw1 = max(velco(i,k), velmin)
769 temv = 1.0 / zw1
770!===============
771! Condition for levels below kref(i): k+1 < kref(i)) ??? see KD05 expression (11) for LLWB only OA >0
772! k >= kref(i) and .... k+1 <kref(i)
773!
774 if (oa(i) > 0. .and. kp1 < kref(i)) then
775 scork = bnv2(i,k) * temv * temv
776 rscor = min(1.0, scork / scor(i))
777 scor(i) = scork
778 else
779 rscor = 1.
780 endif
781!===============
782 brvf = sqrt(bnv2(i,k)) ! brent-vaisala frequency interface
783 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k)) *brvf* zw1
784 hd = sqrt(taup(i,k) / tem1)
785 fro = brvf * hd * temv
786!
787! rim is the "wave"-richardson number byPalmer,Shutts & Swinbank 1986 , PSS-1986
788!
789 tem2 = sqrt(ri_n(i,k))
790 tem = 1. + tem2 * fro
791 ri_gw = ri_n(i,k) * (1.0-fro) / (tem * tem)
792!
793! check Ri-stability to employ the 'dynamical criterion' of PSS-1986
794! assuming co-existence of simultaneous Dyn-Ins and Conv-Ins
795! cos(GW_phase) =1 and sin(GW_phase)=-1
796!23456
797 if (ri_gw <= ric .and.(oa(i) <= 0. .or. kp1 >= kref(i) )) then
798 temc = 2.0 + 1.0 / tem2
799 hd = zw1 * (2.*sqrt(temc)-temc) / brvf
800 taup(i,kp1) = tem1 * hd * hd
801 else
802 taup(i,kp1) = taup(i,k) * rscor
803 endif
804!
805 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
806 endif ! if (.not.icrilv(i) .and. taup(i,k) > 0.0 )
807 endif ! k >= kref(i))
808 enddo ! oro-points
809 enddo ! do k = kmps, kmm1 vertical level loop
810!23456
811! zero momentum deposition at the top model layer: taup(k+1) = taup(k)
812!
813 taup(1:npt,km+1) = taup(1:npt,km)
814!
815! calculate wave acc-n: - (grav)*d(tau)/d(p) = taud
816!23456
817 do k = 1,km
818 do i = 1,npt
819 zw1 = grav*(taup(i,k+1) - taup(i,k))/del(ipt(i),k)
820!======================================================================================
821! zw1 = zw1 * arhills(i) ! simple scale-awareness nhills=Grid_Area/Hil
822! apply limiters for OGW tendency
823!======================================================================================
824 if (abs(zw1) > max_axyz ) zw1 = sign(max_axyz, zw1)
825 taud(i,k)= zw1
826 enddo
827 enddo
828!23456
829!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
830!------if the gravity wave drag would force a critical line in the
831!------layers below sigma=rlolev during the next deltim timestep,
832!------then only apply drag until that critical line is reached.
833! empirical implementation of the llwb-mechanism: lower level wave breaking
834! by limiting "ax = dtfac*ax" due to possible llwb around kref and 500 mb
835! critical line [v - ax*dtp = 0.] is smt like "llwb" for stationary ogws
836!2019: this option limits sensitivity of taux/tauy to variations in "taub"
837!23456~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
838 do k = 1,kmm1
839 do i = 1,npt
840 if (k >= kref(i) .and. prsi(ipt(i),k) >= rlolev .and. taud(i,k) /= 0.) then
841 tem = dtp * taud(i,k)
842 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem)) ! reduce Ax= Ax*(1, or U/dU <=1)
843!default : dtfac(i) = 1.0
844 endif
845 enddo
846 enddo
847!
848!--------- orogw-solver of gfs PSS-1986 is performed
849 else
850!----------orogw-solver of wam2017 out : taup, taud, pkdis
851
852 call oro_spectral_solver(im, km, npt, ipt, kref, kdt, me, master, &
853 dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsl, &
854 del, sigma, hprime, gamma, theta, sinlat, xlatd, taup, taud, pkdis)
855
856 endif ! oro_linsat - linsatdis-solver for stationary OGWs
857!
858!---- above orogw-solver of wam2017------------
859!
860! tofd as in Beljaars-2004 IFS sep-scale ~5km
861! CESM ~ 6km (TMS + OGW/OBL)
862! sgh30 = varss of GSL
863! ----------------------------------------------
864!23456
865 if( do_tofd ) then
866 do i = 1,npt
867 j = ipt(i)
868 zpbl = zmet( j, kpbl(j) )
869 sigflt = min(sgh30(j), 0.3*hprime(j)) ! cannot exceed 30% of ls-sso GSL-2/limits a) 250 m ; b) var_maxfd =150m
870 zsurf = zmeti(j,1)
871 do k=1,km
872 zpm(k) = zmet(j,k)
873 up1(k) = u1(j,k)
874 vp1(k) = v1(j,k)
875 enddo
876
877 call ugwp_tofd1d(km, cpd, dtp, sigflt, zsurf, zpbl, &
878 up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
879
880 do k=1,km
881 dudt_ofd(j,k) = utofd1(k)
882 dvdt_ofd(j,k) = vtofd1(k)
883!
884! add tofd to gw-tendencies
885!
886 pdvdt(j,k) = pdvdt(j,k) + utofd1(k)
887 pdudt(j,k) = pdudt(j,k) + vtofd1(k)
888 pdtdt(j,k) = pdtdt(j,k) + epstofd1(k)
889 enddo
890 du_ofdcol(j) = sum( utofd1(1:km)* del(j,1:km))
891 dv_ofdcol(j) = sum( vtofd1(1:km)* del(j,1:km))
892
893 dusfc(j) = dusfc(j) + du_ofdcol(j)
894 dvsfc(j) = dvsfc(j) + dv_ofdcol(j)
895 enddo
896 endif ! do_tofd
897!23456
898!--------------------------------------------
899! combine oro-drag effects MB +TOFD + OGWs + diag-3d
900!--------------------------------------------
901!234546
902 do k = 1,km
903 do i = 1,npt
904 j = ipt(i)
905 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
906 if ( k < idxzb(i) .and. idxzb(i) /= 0 ) then
907!
908! if blocking layers -- no ogw effects
909!
910 dbim = db(i,k) / (1.+db(i,k)*dtp)
911 dudt_obl(j,k) = -dbim * u1(j,k)
912 dvdt_obl(j,k) = -dbim * v1(j,k)
913
914 pdudt(j,k) = dudt_obl(j,k) +pdudt(j,k)
915 pdvdt(j,k) = dvdt_obl(j,k) +pdvdt(j,k)
916 du_oblcol(j) = du_oblcol(j) + dudt_obl(j,k)* del(j,k)
917 dv_oblcol(j) = dv_oblcol(j) + dvdt_obl(j,k)* del(j,k)
918 dusfc(j) = dusfc(j) + du_oblcol(j)
919 dvsfc(j) = dvsfc(j) + dv_oblcol(j)
920!23456
921 else
922!
923! ogw-s above blocking height
924!
925 taud(i,k) = taud(i,k) * dtfac(i)
926 dtaux = taud(i,k) * xn(i)
927 dtauy = taud(i,k) * yn(i)
928!
929 dudt_ogw(j,k) = dtaux
930 dvdt_ogw(j,k) = dtauy
931!
932 pdvdt(j,k) = dtauy +pdvdt(j,k)
933 pdudt(j,k) = dtaux +pdudt(j,k)
934
935!
936 du_ogwcol(j) = du_ogwcol(j) + dtaux * del(j,k)
937 dv_ogwcol(j) = dv_ogwcol(j) + dtauy * del(j,k)
938!
939 dusfc(j) = dusfc(j) + du_ogwcol(j)
940 dvsfc(j) = dvsfc(j) + dv_ogwcol(j)
941 endif
942!23456
943!============
944! local energy deposition sso-heat due to loss of kinetic energy
945!============
946 unew = u1(j,k) + pdudt(j,k)*dtp ! pdudt(j,k)*dtp
947 vnew = v1(j,k) + pdvdt(j,k)*dtp ! pdvdt(j,k)*dtp
948 eng1 = 0.5*(unew*unew + vnew*vnew)
949 pdtdt(j,k) = max(eng0-eng1,0.)*rcpdt + pdtdt(j,k)
950 enddo
951 enddo
952! dusfc w/o tofd sign as in the era-i, merra and cfsr
953!23456
954 do i = 1,npt
955 j = ipt(i)
956 dusfc(j) = -rgrav * dusfc(j)
957 dvsfc(j) = -rgrav * dvsfc(j)
958 du_ogwcol(j) = -rgrav *du_ogwcol(j)
959 dv_ogwcol(j) = -rgrav *dv_ogwcol(j)
960 du_oblcol(j) = -rgrav *du_oblcol(j)
961 dv_oblcol(j) = -rgrav *dv_oblcol(j)
962 du_ofdcol(j) = -rgrav * du_ofdcol(j)
963 dv_ofdcol(j) = -rgrav * du_ofdcol(j)
964 enddo
965
966 return
967
968
969!============ print/debug after the RETURN statenemt --------------------------------
970 if (kdt <= 2 .and. me == 0) then
971 print *, 'vgw-oro done gwdps_v0 in ugwp-v0 step-proc ', kdt, me
972!
973 print *, maxval(pdudt)*86400., minval(pdudt)*86400, 'vgw_axoro'
974 print *, maxval(pdvdt)*86400., minval(pdvdt)*86400, 'vgw_ayoro'
975! print *, maxval(kdis), minval(kdis), 'vgw_kdispro m2/sec'
976 print *, maxval(pdtdt)*86400., minval(pdtdt)*86400,'vgw_epsoro'
977! print *, maxval(zobl), ' z_mtb ', maxval(tau_mtb), ' tau_mtb '
978 print *, maxval(zogw), ' z_ogw ', maxval(tau_ogw), ' tau_ogw '
979! print *, maxval(tau_tofd), ' tau_tofd '
980! print *, maxval(axtms)*86400., minval(axtms)*86400, 'vgw_axtms'
981! print *,maxval(dudt_mtb)*86400.,minval(dudt_mtb)*86400,'vgw_axmtb'
982 if (maxval(abs(pdudt))*86400. > 100.) then
983
984 print *, maxval(u1), minval(u1), ' u1 gwdps-v1 '
985 print *, maxval(v1), minval(v1), ' v1 gwdps-v1 '
986 print *, maxval(t1), minval(t1), ' t1 gwdps-v1 '
987 print *, maxval(q1), minval(q1), ' q1 gwdps-v1 '
988 print *, maxval(del), minval(del), ' del gwdps-v1 '
989 print *, maxval(zmet),minval(zmet), 'zmet'
990 print *, maxval(zmeti),minval(zmeti), 'zmeti'
991 print *, maxval(prsi), minval(prsi), ' prsi '
992 print *, maxval(prsl), minval(prsl), ' prsl '
993 print *, maxval(ro), minval(ro), ' ro-dens '
994 print *, maxval(bnv2(1:npt,:)), minval(bnv2(1:npt,:)),' bnv2 '
995 print *, maxval(kpbl), minval(kpbl), ' kpbl '
996 print *, maxval(sgh30), maxval(hprime), maxval(elvmax),'oro-d'
997 print *
998 do i =1, npt
999 j= ipt(i)
1000 print *,zogw(j)/hprime(j), zobl(j)/hprime(j), &
1001 zmet(j,1)*1.e-3, nint(hprime(j)/sigma(j))
1002!
1003 enddo
1004 print *
1005 errflg = 1
1006 errmsg = 'ERROR(orogw_v1): '
1007 return
1008 endif
1009 endif
1010
1011 return
1012 end subroutine orogw_v1
1013!
1014!
1015 subroutine ugwp_tofd1d(levs, con_cp, dtp, sigflt, zsurf, zpbl, u, v, &
1016 zmid, utofd, vtofd, epstofd, krf_tofd)
1017
1018 use machine , only : kind_phys
1019 use ugwp_oro_init, only : n_tofd, const_tofd, ze_tofd, a12_tofd, ztop_tofd
1020!
1021! adding the implicit tendency estimate
1022!
1023 implicit none
1024 integer, intent(in) :: levs
1025 real(kind_phys), intent(in) :: con_cp
1026 real(kind_phys), intent(in) :: dtp
1027
1028 real(kind_phys), intent(in), dimension(levs) :: u, v, zmid
1029 real(kind_phys), intent(in) :: sigflt, zpbl, zsurf
1030
1031 real(kind_phys), intent(out), dimension(levs) :: utofd, vtofd, epstofd, krf_tofd
1032!
1033! locals
1034!
1035 integer :: i, k
1036 real(kind=kind_phys) :: rcpd2, tofd_mag, tofd_zdep
1037 real(kind_phys) :: unew, vnew, eknew
1038
1039 real(kind=kind_phys), parameter :: sghmax = 5. ! dz(1)/5= 25/5 m dz-of the first layer
1040 real(kind=kind_phys), parameter :: tend_imp = 1.
1041
1042
1043 real(kind=kind_phys) :: sgh2, ekin, zdec, rzdec, umag, zmet, zarg, ztexp, krf
1044!
1045 utofd =0.0 ; vtofd = 0.0 ; epstofd =0.0 ; krf_tofd =0.0
1046 rcpd2 = 0.5/con_cp
1047!
1048 zdec = max(n_tofd*sigflt, zpbl) ! ntimes*sgh_turb or Zpbl
1049 zdec = min(ze_tofd, zdec) ! cannot exceed ~1.5 km
1050! H_efold = max(2*varss, hpbl)
1051! H_efold = min(H_efold,1500.)
1052 rzdec = 1.0/zdec
1053
1054 sgh2 = max(sigflt*sigflt, sghmax*sghmax) ! dz ~25m of the first layer in FV3GFS-127L
1055 tofd_mag = const_tofd * a12_tofd * sgh2 ! * scale_res
1056
1057! GSL-darg scheme: varmax_fd, beta_fd ,250.
1058! var_temp = MIN(varss,varmax_fd) + MAX(0., 0.1*(varss-varmax_fd))
1059! var_temp = MIN(var_temp, 250.)
1060! var_temp = var_temp * var_temp
1061!
1062! a12=a1* 0.005363 * 0.0759 * 0.00026615161
1063!
1064! rzdec 1./H_efold
1065! do k=1,levs
1066! zmet = zmid(k)-zsurf
1067! wsp=SQRT(u(k)*u(k) + v(k)*v(k)) ! abs(V)
1068! zarg = zmet*rzdec
1069! var_temp = var_temp * a12 * exp(-zarg*sqrt(zarg))*zmet**(-1.2) ! this > 0
1070! krf = var_temp * wsp /(1. + var_temp*dtp*wsp)
1071! utofd(k) = -u(k) *krf
1072! vtofd(k) = -v(k)/(1. + var_temp*krf
1073! enddo
1074
1075 do k=1, levs
1076 zmet = zmid(k)-zsurf
1077 if (zmet > ztop_tofd) cycle
1078
1079 ekin = u(k)*u(k) + v(k)*v(k)
1080
1081 umag = sqrt(ekin)
1082 zarg = zmet*rzdec
1083 ztexp = exp(-zarg * sqrt(zarg))
1084
1085 tofd_zdep = zmet ** (-1.2) *ztexp
1086 krf = umag * tofd_mag * tofd_zdep
1087
1088 if (tend_imp == 1.) then
1089 krf = krf/(1.+krf*dtp)
1090 endif
1091
1092 utofd(k) = -krf*u(k)
1093 vtofd(k) = -krf*v(k)
1094 if (tend_imp == 1.) then
1095 unew =u(k)+ utofd(k)*dtp ; vnew =v(k)+ vtofd(k)*dtp
1096 eknew =unew*unew + vnew*vnew
1097 epstofd(k) = rcpd2*(ekin-eknew)
1098 else
1099 epstofd(k) = rcpd2*krf*ekin
1100 endif
1101 ! more accurate heat/mom form using "implicit tend-solver"
1102 ! to update momentum and temp-re; epstofd(k) can be skipped
1103 krf_tofd(k) = krf ! can be used as addition to the mesoscale blocking
1104 enddo
1105!
1106 end subroutine ugwp_tofd1d
1107
1108end module cires_ugwpv1_oro