CCPP SciDoc v7.0.x  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
ugwp_driver_v0.F
1
2
4 module ugwp_driver_v0
5 use cires_orowam2017
6 contains
7!
8!=====================================================================
9!
10!ugwp-v0 subroutines: GWDPS_V0 and fv3_ugwp_solv2_v0
11!
12!=====================================================================
30
33 SUBROUTINE gwdps_v0(IM, km, imx, do_tofd,
34 & Pdvdt, Pdudt, Pdtdt, Pkdis, U1,V1,T1,Q1,KPBL,
35 & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DTP,KDT,
36 & sgh30, HPRIME,OC,OA4,CLX4,THETA,vSIGMA,vGAMMA,ELVMAXD,
37 & DUSFC, DVSFC, xlatd, sinlat, coslat, sparea,
38 & cdmbgwd, me, master, rdxzb, con_g, con_omega,
39 & zmtb, zogw, tau_mtb, tau_ogw, tau_tofd,
40 & dudt_mtb, dudt_ogw, dudt_tms)
41!----------------------------------------
42! ugwp_v0
43!
44! modified/revised version of gwdps.f (with bug fixes, tofd, appropriate
45! computation of reference level for OGW + COORDE diagnostics
46! all constants/parameters inside cires_ugwp_initialize.F90
47!----------------------------------------
48
49 USE machine , ONLY : kind_phys
50 use ugwp_common_v0,only : rgrav, grav, cpd, rd, rv, rcpd, rcpd2
51 &, pi, rad_to_deg, deg_to_rad, pi2
52 &, rdi, gor, grcp, gocp, fv, gr2
53 &, bnv2min, dw2min, velmin, arad
54
55 use ugwpv0_oro_init, only : rimin, ric, efmin, efmax
56 &, hpmax, hpmin, sigfaci => sigfac
57 &, dpmin, minwnd, hminmt, hncrit
58 &, rlolev, gmax, veleps, factop
59 &, frc, ce, ceofrc, frmax, cg
60 &, fdir, mdir, nwdir
61 &, cdmb, cleff, fcrit_gfs, fcrit_mtb
62 &, n_tofd, ze_tofd, ztop_tofd
63
64 use cires_ugwpv0_module, only : kxw, max_kdis, max_axyz
65
66!----------------------------------------
67 implicit none
68 integer, parameter :: kp = kind_phys
69 character(len=8) :: strsolver='PSS-1986' ! current operational solver or 'WAM-2017'
70 integer, intent(in) :: im, km, imx, kdt
71 integer, intent(in) :: me, master
72 logical, intent(in) :: do_tofd
73 real(kind=kind_phys), parameter :: sigfac = 3, sigfacs = 0.5
74 real(kind=kind_phys) :: ztoph,zlowh,ph_blk, dz_blk
75 integer, intent(in) :: KPBL(IM) ! Index for the PBL top layer!
76 real(kind=kind_phys), intent(in) :: dtp ! time step
77 real(kind=kind_phys), intent(in) :: cdmbgwd(2)
78
79 real(kind=kind_phys), intent(in), dimension(im,km) ::
80 & u1, v1, t1, q1,
81 & del, prsl, prslk, phil
82 real(kind=kind_phys), intent(in),dimension(im,km+1):: prsi, phii
83 real(kind=kind_phys), intent(in) :: xlatd(im),sinlat(im),
84 & coslat(im)
85 real(kind=kind_phys), intent(in) :: sparea(im)
86
87 real(kind=kind_phys), intent(in) :: oc(im), oa4(im,4), clx4(im,4)
88 real(kind=kind_phys), intent(in) :: hprime(im), sgh30(im)
89 real(kind=kind_phys), intent(in) :: elvmaxd(im), theta(im)
90 real(kind=kind_phys), intent(in) :: vsigma(im), vgamma(im)
91 real(kind=kind_phys) :: sigma(im), gamma(im)
92
93 real(kind=kind_phys), intent(in) :: con_g, con_omega
94
95!output -phys-tend
96 real(kind=kind_phys),dimension(im,km),intent(out) ::
97 & pdvdt, pdudt, pkdis, pdtdt
98! output - diag-coorde
99 &, dudt_mtb, dudt_ogw, dudt_tms
100!
101 real(kind=kind_phys),dimension(im) :: rdxzb, zmtb, zogw
102 &, tau_ogw, tau_mtb, tau_tofd
103 &, dusfc, dvsfc
104!
105!---------------------------------------------------------------------
106! # of permissible sub-grid orography hills for "any" resolution < 25
107! correction for "elliptical" hills based on shilmin-area =sgrid/25
108! 4.*gamma*b_ell*b_ell >= shilmin
109! give us limits on [b_ell & gamma *b_ell] > 5 km =sso_min
110! gamma_min = 1/4*shilmin/sso_min/sso_min
111!23.01.2019: cdmb = 4.*192/768_c192=1 x 0.5
112! 192: cdmbgwd = 0.5, 2.5
113! cleff = 2.5*0.5e-5 * sqrt(192./768.) => Lh_eff = 1004. km
114! 6*dx = 240 km 8*dx = 320. ~ 3-5 more effective
115!---------------------------------------------------------------------
116 real(kind=kind_phys) :: gammin = 0.00999999_kp
117 real(kind=kind_phys), parameter :: nhilmax = 25.0_kp
118 real(kind=kind_phys), parameter :: sso_min = 3000.0_kp
119 logical, parameter :: do_adjoro = .true.
120!
121 real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
122 real(kind=kind_phys) :: belpmin, dsmin, dsmax
123! real(kind=kind_phys) :: arhills(im) ! not used why do we need?
124 real(kind=kind_phys) :: xlingfs
125
126!
127! locals
128! mean flow
129 real(kind=kind_phys), dimension(im,km) :: ri_n, bnv2, ro
130 &, vtk, vtj, velco
131!mtb
132 real(kind=kind_phys), dimension(im) :: oa, clx , elvmax, wk
133 &, pe, ek, up
134
135 real(kind=kind_phys), dimension(im,km) :: db, ang, uds
136
137 real(kind=kind_phys) :: zlen, dbtmp, r, phiang, dbim, zr
138 real(kind=kind_phys) :: eng0, eng1, cosang2, sinang2
139 real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
140!
141! TOFD
142! Some constants now in "use ugwp_oro_init" + "use ugwp_common"
143!
144!==================
145 real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
146 real(kind=kind_phys), dimension(km) :: utofd1, vtofd1
147 &, epstofd1, krf_tofd1
148 &, up1, vp1, zpm
149 real(kind=kind_phys),dimension(im, km) :: axtms, aytms
150!
151! OGW
152!
153 LOGICAL ICRILV(IM)
154!
155 real(kind=kind_phys), dimension(im) :: xn, yn, ubar, vbar, ulow,
156 & roll, bnv2bar, scor, dtfac, xlinv, delks, delks1
157!
158 real(kind=kind_phys) :: taup(im,km+1), taud(im,km)
159 real(kind=kind_phys) :: taub(im), taulin(im), heff, hsat, hdis
160
161 integer, dimension(im) :: kref, idxzb, ipt, kreflm,
162 & iwklm, iwk, izlow
163!
164!check what we need
165!
166 real(kind=kind_phys) :: bnv, fr, ri_gw
167 &, brvf, tem, tem1, tem2, temc, temv
168 &, ti, rdz, dw2, shr2, bvf2
169 &, rdelks, efact, coefm, gfobnv
170 &, scork, rscor, hd, fro, sira
171 &, dtaux, dtauy, pkp1log, pklog
172 &, grav2, rcpdt, windik, wdir
173 &, sigmin, dxres,sigres,hdxres
174 &, cdmb4, mtbridge
175 &, kxridge, inv_b2eff, zw1, zw2
176 &, belps, aelps, nhills, selps
177
178 integer :: kmm1, kmm2, lcap, lcapp1
179 &, npt, kbps, kbpsp1,kbpsm1
180 &, kmps, idir, nwd, klcap, kp1, kmpbl, kmll
181 &, k_mtb, k_zlow, ktrial, klevm1, i, j, k
182!
183 rcpdt = 1.0 / (cpd*dtp)
184 grav2 = grav + grav
185!
186! mtb-blocking sigma_min and dxres => cires_initialize
187!
188 sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
189 dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
190
191 dxres = pi2*arad/float(imx)
192 hdxres = 0.5_kp*dxres
193! shilmin = sgrmin/nhilmax ! not used - Moorthi
194
195! gammin = min(sso_min/dsmax, 1.) ! Moorthi - with this results are not reproducible
196 gammin = min(sso_min/dxres, 1.) ! Moorthi
197
198! sigmin = 2.*hpmin/dsmax !dxres ! Moorthi - this will not reproduce
199 sigmin = 2.*hpmin/dxres !dxres
200
201
202
203 kxridge = float(imx)/arad * cdmbgwd(2)
204
205 do i=1,im
206 idxzb(i) = 0
207 zmtb(i) = 0.0
208 zogw(i) = 0.0
209 rdxzb(i) = 0.0
210 tau_ogw(i) = 0.0
211 tau_mtb(i) = 0.0
212 dusfc(i) = 0.0
213 dvsfc(i) = 0.0
214 tau_tofd(i) = 0.0
215!
216 ipt(i) = 0
217 sigma(i) = max(vsigma(i), sigmin)
218 gamma(i) = max(vgamma(i), gammin)
219 enddo
220
221 do k=1,km
222 do i=1,im
223 pdvdt(i,k) = 0.0
224 pdudt(i,k) = 0.0
225 pdtdt(i,k) = 0.0
226 pkdis(i,k) = 0.0
227 dudt_mtb(i,k) = 0.0
228 dudt_ogw(i,k) = 0.0
229 dudt_tms(i,k) = 0.0
230 enddo
231 enddo
232
233! ---- for lm and gwd calculation points
234
235 npt = 0
236 do i = 1,im
237 if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin ) then
238
239 npt = npt + 1
240 ipt(npt) = i
241! arhills(i) = 1.0
242!
243 sigres = max(sigmin, sigma(i))
244! if (sigma(i) < sigmin) sigma(i)= sigmin
245 dxres = sqrt(sparea(i))
246 if (2.*hprime(i)/sigres > dxres) sigres=2.*hprime(i)/dxres
247 aelps = min(2.*hprime(i)/sigres, 0.5*dxres)
248 if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i),.5*dxres)
249!
250! small-scale "turbulent" oro-scales < sso_min
251!
252 if( aelps < sso_min .and. do_adjoro) then
253
254! a, b > sso_min upscale ellipse a/b > 0.1 a>sso_min & h/b=>new_sigm
255!
256 aelps = sso_min
257 if (belps < sso_min ) then
258 gamma(i) = 1.0
259 belps = aelps*gamma(i)
260 else
261 gamma(i) = min(aelps/belps, 1.0)
262 endif
263 sigma(i) = 2.*hprime(i)/aelps
264 gamma(i) = min(aelps/belps, 1.0)
265 endif
266
267 selps = belps*belps*gamma(i)*4. ! ellipse area of the el-c hill
268 nhills = min(nhilmax, sparea(i)/selps)
269! arhills(i) = max(nhills, 1.0)
270
271!333 format( ' nhil: ', I6, 4(2x, F9.3), 2(2x, E9.3))
272! if (kdt==1 )
273! & write(6,333) nint(nhills)+1,xlatd(i), hprime(i),aelps*1.e-3,
274! & belps*1.e-3, sigma(i),gamma(i)
275
276 endif
277 enddo
278
279 IF (npt == 0) then
280 RETURN ! No gwd/mb calculation done
281 endif
282
283
284 do i=1,npt
285 iwklm(i) = 2
286 idxzb(i) = 0
287 kreflm(i) = 0
288 enddo
289
290 do k=1,km
291 do i=1,im
292 db(i,k) = 0.0
293 ang(i,k) = 0.0
294 uds(i,k) = 0.0
295 enddo
296 enddo
297
298 kmm1 = km - 1 ; kmm2 = km - 2 ; kmll = kmm1
299 lcap = km ; lcapp1 = lcap + 1
300
301 DO i = 1, npt
302 j = ipt(i)
303 elvmax(j) = min(elvmaxd(j)*0. + sigfac * hprime(j), hncrit)
304 izlow(i) = 1 ! surface-level
305 ENDDO
306!
307 DO k = 1, kmm1
308 DO i = 1, npt
309 j = ipt(i)
310 ztoph = sigfac * hprime(j)
311 zlowh = sigfacs* hprime(j)
312 pkp1log = phil(j,k+1) * rgrav
313 pklog = phil(j,k) * rgrav
314! if (( ELVMAX(j) <= pkp1log) .and. (ELVMAX(j).ge.pklog) )
315! & iwklm(I) = MAX(iwklm(I), k+1 )
316 if (( ztoph <= pkp1log) .and. (ztoph >= pklog) )
317 & iwklm(i) = max(iwklm(i), k+1 )
318!
319 if (zlowh <= pkp1log .and. zlowh >= pklog)
320 & izlow(i) = max(izlow(i),k)
321 ENDDO
322 ENDDO
323!
324 DO k = 1,km
325 DO i =1,npt
326 j = ipt(i)
327 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
328 vtk(i,k) = vtj(i,k) / prslk(j,k)
329 ro(i,k) = rdi * prsl(j,k) / vtj(i,k) ! DENSITY mid-levels
330 taup(i,k) = 0.0
331 ENDDO
332 ENDDO
333!
334! check RI_N or RI_MF computation
335!
336 DO k = 1,kmm1
337 DO i =1,npt
338 j = ipt(i)
339 rdz = grav / (phil(j,k+1) - phil(j,k))
340 tem1 = u1(j,k) - u1(j,k+1)
341 tem2 = v1(j,k) - v1(j,k+1)
342 dw2 = tem1*tem1 + tem2*tem2
343 shr2 = max(dw2,dw2min) * rdz * rdz
344! TI = 2.0 / (T1(J,K)+T1(J,K+1))
345! BVF2 = Grav*(GOCP+RDZ*(VTJ(I,K+1)-VTJ(I,K)))* TI
346! RI_N(I,K) = MAX(BVF2/SHR2,RIMIN) ! Richardson number
347!
348 bvf2 = grav2 * rdz * (vtk(i,k+1)-vtk(i,k))
349 & / (vtk(i,k+1)+vtk(i,k))
350 bnv2(i,k+1) = max( bvf2, bnv2min )
351 ! https://github.com/NCAR/ccpp-physics/issues/1103
352 !RI_N(I,K+1) = Bnv2(i,k)/SHR2 ! Richardson number consistent with BNV2
353 ri_n(i,k+1) = bnv2(i,max(k,2))/shr2 ! Richardson number consistent with BNV2
354!
355! add here computation for Ktur and OGW-dissipation fro VE-GFS
356!
357 ENDDO
358 ENDDO
359 k = 1
360 DO i = 1, npt
361 bnv2(i,k) = bnv2(i,k+1)
362 ENDDO
363!
364! level iwklm =>phil(j,k)/g < sigfac * hprime(j) < phil(j,k+1)/g
365!
366 DO i = 1, npt
367 j = ipt(i)
368 k_zlow = izlow(i)
369 if (k_zlow == iwklm(i)) k_zlow = 1
370 delks(i) = 1.0 / (prsi(j,k_zlow) - prsi(j,iwklm(i)))
371! DELKS1(I) = 1.0 /(PRSL(J,k_zlow) - PRSL(J,iwklm(i)))
372 ubar(i) = 0.0
373 vbar(i) = 0.0
374 roll(i) = 0.0
375 pe(i) = 0.0
376 ek(i) = 0.0
377 bnv2bar(i) = 0.0
378 ENDDO
379!
380 DO i = 1, npt
381 k_zlow = izlow(i)
382 if (k_zlow == iwklm(i)) k_zlow = 1
383 DO k = k_zlow, iwklm(i)-1 ! Kreflm(I)= iwklm(I)-1
384 j = ipt(i) ! laye-aver Rho, U, V
385 rdelks = del(j,k) * delks(i)
386 ubar(i) = ubar(i) + rdelks * u1(j,k) ! trial Mean U below
387 vbar(i) = vbar(i) + rdelks * v1(j,k) ! trial Mean V below
388 roll(i) = roll(i) + rdelks * ro(i,k) ! trial Mean RO below
389!
390 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
391 ENDDO
392 ENDDO
393!
394 DO i = 1, npt
395 j = ipt(i)
396!
397! integrate from Ztoph = sigfac*hprime down to Zblk if exists
398! find ph_blk, dz_blk like in LM-97 and IFS
399!
400 ph_blk =0.
401 DO k = iwklm(i), 1, -1
402 phiang = atan2(v1(j,k),u1(j,k))*rad_to_deg
403 ang(i,k) = ( theta(j) - phiang )
404 if ( ang(i,k) > 90. ) ang(i,k) = ang(i,k) - 180.
405 if ( ang(i,k) < -90. ) ang(i,k) = ang(i,k) + 180.
406 ang(i,k) = ang(i,k) * deg_to_rad
407 uds(i,k) =
408 & max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), velmin)
409!
410 IF (idxzb(i) == 0 ) then
411 dz_blk = ( phii(j,k+1) - phii(j,k) ) *rgrav
412 pe(i) = pe(i) + bnv2(i,k) *
413 & ( elvmax(j) - phil(j,k)*rgrav ) * dz_blk
414
415 up(i) = max(uds(i,k) * cos(ang(i,k)), velmin)
416 ek(i) = 0.5 * up(i) * up(i)
417
418 ph_blk = ph_blk + dz_blk*sqrt(bnv2(i,k))/up(i)
419
420! --- Dividing Stream lime is found when PE =exceeds EK. oper-l GFS
421! IF ( PE(I) >= EK(I) ) THEN
422 IF ( ph_blk >= fcrit_gfs ) THEN
423 idxzb(i) = k
424 zmtb(j) = phil(j, k)*rgrav
425 rdxzb(j) = real(k, kind=kind_phys)
426 ENDIF
427
428 ENDIF
429 ENDDO
430!
431! Alternative expression: ZMTB = max(Heff*(1. -Fcrit_gfs/Fr), 0)
432! fcrit_gfs/fr
433!
434 goto 788
435
436 bnv = sqrt( bnv2bar(i) )
437 heff = 2.*min(hprime(j),hpmax)
438 zw2 = ubar(i)*ubar(i)+vbar(i)*vbar(i)
439 ulow(i) = sqrt(max(zw2,dw2min))
440 fr = heff*bnv/ulow(i)
441 zw1 = max(heff*(1. -fcrit_gfs/fr), 0.0)
442 zw2 = phil(j,2)*rgrav
443 if (fr > fcrit_gfs .and. zw1 > zw2 ) then
444 do k=2, kmm1
445 pkp1log = phil(j,k+1) * rgrav
446 pklog = phil(j,k) * rgrav
447 if (zw1 <= pkp1log .and. zw1 >= pklog) exit
448 enddo
449 idxzb(i) = k
450 zmtb(j) = phil(j, k)*rgrav
451 else
452 zmtb(j) = 0.
453 idxzb(i) = 0
454 endif
455788 continue
456 ENDDO
457
458!
459! --- The drag for mtn blocked flow
460!
461 cdmb4 = 0.25*cdmb
462 DO i = 1, npt
463 j = ipt(i)
464!
465 IF ( idxzb(i) > 0 ) then
466! (4.16)-IFS
467 gam2 = gamma(j)*gamma(j)
468 bgam = 1.0 - 0.18*gamma(j) - 0.04*gam2
469 cgam = 0.48*gamma(j) + 0.30*gam2
470 DO k = idxzb(i)-1, 1, -1
471
472 zlen = sqrt( ( phil(j,idxzb(i)) - phil(j,k) ) /
473 & ( phil(j,k ) + grav * hprime(j) ) )
474
475 tem = cos(ang(i,k))
476 cosang2 = tem * tem
477 sinang2 = 1.0 - cosang2
478!
479! cos =1 sin =0 => 1/R= gam ZR = 2.-gam
480! cos =0 sin =1 => 1/R= 1/gam ZR = 2.- 1/gam
481!
482 rdem = cosang2 + gam2 * sinang2
483 rnom = cosang2*gam2 + sinang2
484!
485! metOffice Dec 2010
486! correction of H. Wells & A. Zadra for the
487! aspect ratio of the hill seen by MF
488! (1/R , R-inverse below: 2-R)
489
490 rdem = max(rdem, 1.e-6)
491 r = sqrt(rnom/rdem)
492 zr = max( 2. - r, 0. )
493
494 sigres = max(sigmin, sigma(j))
495 if (hprime(j)/sigres > dxres) sigres = hprime(j)/dxres
496 mtbridge = zr * sigres*zlen / hprime(j)
497! (4.15)-IFS
498! DBTMP = CDmb4 * mtbridge *
499! & MAX(cos(ANG(I,K)), gamma(J)*sin(ANG(I,K)))
500! (4.16)-IFS
501 dbtmp = cdmb4*mtbridge*(bgam* cosang2 +cgam* sinang2)
502 db(i,k)= dbtmp * uds(i,k)
503 ENDDO
504!
505 endif
506 ENDDO
507!
508!.............................
509!.............................
510! end mtn blocking section
511!.............................
512!.............................
513!
514!--- Orographic Gravity Wave Drag Section
515!
516! Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62
517! inside "cires_ugwp_initialize.F90" now
518!
519 kmpbl = km / 2
520 iwk(1:npt) = 2
521!
522! METO-scheme:
523! k_mtb = max(k_zmtb, k_n*hprime/2] to reduce diurnal variations taub_ogw
524!
525 DO k=3,kmpbl
526 DO i=1,npt
527 j = ipt(i)
528 tem = (prsi(j,1) - prsi(j,k))
529 if (tem < dpmin) iwk(i) = k ! dpmin=50 mb
530
531!===============================================================
532! lev=111 t=311.749 hkm=0.430522 Ps-P(iwk)=52.8958
533! below "Hprime" - source of OGWs and below Zblk !!!
534! 27 2 kpbl ~ 1-2 km < Hprime
535!===============================================================
536 enddo
537 enddo
538!
539! iwk - adhoc GFS-parameter to select OGW-launch level between
540! LEVEL ~0.4-0.5 KM from surface or/and PBL-top
541! in UGWP-V1: options to modify as Htop ~ (2-3)*Hprime > Zmtb
542! in UGWP-V0 we ensured that : Zogw > Zmtb
543!
544
545 kbps = 1
546 kmps = km
547 k_mtb = 1
548 DO i=1,npt
549 j = ipt(i)
550 k_mtb = max(1, idxzb(i))
551
552 kref(i) = max(iwk(i), kpbl(j)+1 ) ! reference level PBL or smt-else ????
553 kref(i) = max(kref(i), iwklm(i) ) ! iwklm => sigfac*hprime
554
555 if (kref(i) <= idxzb(i)) kref(i) = idxzb(i) + 1 ! layer above zmtb
556 kbps = max(kbps, kref(i))
557 kmps = min(kmps, kref(i))
558!
559 delks(i) = 1.0 / (prsi(j,k_mtb) - prsi(j,kref(i)))
560 ubar(i) = 0.0
561 vbar(i) = 0.0
562 roll(i) = 0.0
563 bnv2bar(i)= 0.0
564 ENDDO
565!
566 kbpsp1 = kbps + 1
567 kbpsm1 = kbps - 1
568 k_mtb = 1
569!
570 DO i = 1,npt
571 k_mtb = max(1, idxzb(i))
572 DO k = k_mtb,kbps !KBPS = MAX(kref) ;KMPS= MIN(kref)
573 IF (k < kref(i)) THEN
574 j = ipt(i)
575 rdelks = del(j,k) * delks(i)
576 ubar(i) = ubar(i) + rdelks * u1(j,k) ! Mean U below kref
577 vbar(i) = vbar(i) + rdelks * v1(j,k) ! Mean V below kref
578 roll(i) = roll(i) + rdelks * ro(i,k) ! Mean RO below kref
579 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
580 ENDIF
581 ENDDO
582 ENDDO
583!
584! orographic asymmetry parameter (OA), and (CLX)
585 DO i = 1,npt
586 j = ipt(i)
587 wdir = atan2(ubar(i),vbar(i)) + pi
588 idir = mod(nint(fdir*wdir),mdir) + 1
589 nwd = nwdir(idir)
590 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
591 clx(i) = clx4(j,mod(nwd-1,4)+1)
592 ENDDO
593!
594 DO i = 1,npt
595 dtfac(i) = 1.0
596 icrilv(i) = .false. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
597 ulow(i) = max(sqrt(ubar(i)*ubar(i)+vbar(i)*vbar(i)),velmin)
598 xn(i) = ubar(i) / ulow(i)
599 yn(i) = vbar(i) / ulow(i)
600 ENDDO
601!
602 DO k = 1, kmm1
603 DO i = 1,npt
604 j = ipt(i)
605 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*xn(i)
606 & + (v1(j,k)+v1(j,k+1))*yn(i))
607 ENDDO
608 ENDDO
609!
610!------------------
611! v0: incorporates latest modifications for kxridge and heff/hsat
612! and taulin for Fr <=fcrit_gfs
613! and concept of "clipped" hill if zmtb > 0. to make
614! the integrated "tau_sso = tau_ogw +tau_mtb" close to reanalysis data
615! it is still used the "single-OGWave"-approach along ULOW-upwind
616!
617! in contrast to the 2-orthogonal wave (2OTW) schemes of IFS/METO/E-CANADA
618! 2OTW scheme requires "aver angle" and wind projections on 2 axes of ellipse a-b
619! with 2-stresses: taub_a & taub_b from AS of Phillips et al. (1984)
620!------------------
621 taub(:) = 0. ; taulin(:)= 0.
622 DO i = 1,npt
623 j = ipt(i)
624 bnv = sqrt( bnv2bar(i) )
625 heff = min(hprime(j),hpmax)
626
627 if( zmtb(j) > 0.) heff = max(sigfac*heff-zmtb(j), 0.)/sigfac
628 if (heff <= 0) cycle
629
630 hsat = fcrit_gfs*ulow(i)/bnv
631 heff = min(heff, hsat)
632
633 fr = min(bnv * heff /ulow(i), frmax)
634!
635 efact = (oa(i) + 2.) ** (ceofrc*fr)
636 efact = min( max(efact,efmin), efmax )
637!
638 coefm = (1. + clx(i)) ** (oa(i)+1.)
639!
640 xlinv(i) = coefm * cleff ! effective kxw for Lin-wave
641 xlingfs = coefm * cleff
642!
643 tem = fr * fr * oc(j)
644 gfobnv = gmax * tem / ((tem + cg)*bnv)
645!
646!new specification of XLINV(I) & taulin(i)
647
648 sigres = max(sigmin, sigma(j))
649 if (heff/sigres > hdxres) sigres = heff/hdxres
650 inv_b2eff = 0.5*sigres/heff
651 kxridge = 1.0 / sqrt(sparea(j))
652 xlinv(i) = xlingfs !or max(kxridge, inv_b2eff) ! 6.28/Lx ..0.5*sigma(j)/heff = 1./Lridge
653 taulin(i) = 0.5*roll(i)*xlinv(i)*bnv*ulow(i)*
654 & heff*heff
655
656 if ( fr > fcrit_gfs ) then
657 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
658 & * ulow(i) * gfobnv * efact ! nonlinear FLUX Tau0...XLINV(I)
659!
660 else
661!
662 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
663 & * ulow(i) * gfobnv * efact
664!
665! TAUB(I) = taulin(i) ! linear flux for FR <= fcrit_gfs
666!
667 endif
668!
669!
670 k = max(1, kref(i)-1)
671 tem = max(velco(i,k)*velco(i,k), dw2min)
672 scor(i) = bnv2(i,k) / tem ! Scorer parameter below ref level
673!
674! diagnostics for zogw > zmtb
675!
676 zogw(j) = phii(j, kref(i)) *rgrav
677 ENDDO
678!
679!----SET UP BOTTOM VALUES OF STRESS
680!
681 DO k = 1, kbps
682 DO i = 1,npt
683 IF (k <= kref(i)) taup(i,k) = taub(i)
684 ENDDO
685 ENDDO
686
687 if (strsolver == 'PSS-1986') then
688
689!======================================================
690! V0-GFS OROGW-solver of Palmer et al 1986 -"PSS-1986"
691! in V1-OROGW LINSATDIS of "WAM-2017"
692! with LLWB-mechanism for
693! rotational/non-hydrostat OGWs important for
694! HighRES-FV3GFS with dx < 10 km
695!======================================================
696
697 DO k = kmps, kmm1 ! Vertical Level Loop
698 kp1 = k + 1
699 DO i = 1, npt
700!
701 IF (k >= kref(i)) THEN
702 icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) < ric)
703 & .OR. (velco(i,k) <= 0.0)
704 ENDIF
705 ENDDO
706!
707 DO i = 1,npt
708 IF (k >= kref(i)) THEN
709 IF (.NOT.icrilv(i) .AND. taup(i,k) > 0.0 ) THEN
710 temv = 1.0 / max(velco(i,k), velmin)
711!
712 IF (oa(i) > 0. .AND. kp1 < kref(i)) THEN
713 scork = bnv2(i,k) * temv * temv
714 rscor = min(1.0, scork / scor(i))
715 scor(i) = scork
716 ELSE
717 rscor = 1.
718 ENDIF
719!
720 brvf = sqrt(bnv2(i,k)) ! Brent-Vaisala Frequency interface
721! TEM1 = XLINV(I)*(RO(I,KP1)+RO(I,K))*BRVF*VELCO(I,K)*0.5
722
723 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf*0.5
724 & * max(velco(i,k), velmin)
725 hd = sqrt(taup(i,k) / tem1)
726 fro = brvf * hd * temv
727!
728! RIM is the "WAVE"-RICHARDSON NUMBER BY PALMER,Shutts, Swinbank 1986
729!
730
731 tem2 = sqrt(ri_n(i,k))
732 tem = 1. + tem2 * fro
733 ri_gw = ri_n(i,k) * (1.0-fro) / (tem * tem)
734!
735! CHECK STABILITY TO EMPLOY THE 'dynamical SATURATION HYPOTHESIS'
736! OF PALMER,Shutts, Swinbank 1986
737! ----------------------
738 IF (ri_gw <= ric .AND.
739 & (oa(i) <= 0. .OR. kp1 >= kref(i) )) THEN
740 temc = 2.0 + 1.0 / tem2
741 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf
742 taup(i,kp1) = tem1 * hd * hd
743 ELSE
744 taup(i,kp1) = taup(i,k) * rscor
745 ENDIF
746 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
747 ENDIF
748 ENDIF
749 ENDDO
750 ENDDO
751!
752! zero momentum deposition at the top model layer
753!
754 taup(1:npt,km+1) = taup(1:npt,km)
755!
756! Calculate wave acc-n: - (grav)*d(tau)/d(p) = taud
757!
758 DO k = 1,km
759 DO i = 1,npt
760 taud(i,k) = grav*(taup(i,k+1) - taup(i,k))/del(ipt(i),k)
761 ENDDO
762 ENDDO
763!
764!------scale MOMENTUM DEPOSITION AT TOP TO 1/2 VALUE
765! it is zero now
766! DO I = 1,npt
767! TAUD(I, km) = TAUD(I,km) * FACTOP
768! ENDDO
769
770!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
771!------IF THE GRAVITY WAVE DRAG WOULD FORCE A CRITICAL LINE IN THE
772!------LAYERS BELOW SIGMA=RLOLEV DURING THE NEXT DELTIM TIMESTEP,
773!------THEN ONLY APPLY DRAG UNTIL THAT CRITICAL LINE IS REACHED.
774! Empirical implementation of the LLWB-mechanism: Lower Level Wave Breaking
775! by limiting "Ax = Dtfac*Ax" due to possible LLWB around Kref and 500 mb
776! critical line [V - Ax*dtp = 0.] is smt like "LLWB" for stationary OGWs
777!2019: this option limits sensitivity of taux/tauy to increase/decreaseof TAUB
778!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
779 DO k = 1,kmm1
780 DO i = 1,npt
781 IF (k >= kref(i) .and. prsi(ipt(i),k) >= rlolev) THEN
782
783 IF(taud(i,k) /= 0.) THEN
784 tem = dtp * taud(i,k)
785 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
786! DTFAC(I) = 1.0
787 ENDIF
788 ENDIF
789 ENDDO
790 ENDDO
791!
792!--------------------------- OROGW-solver of GFS PSS-1986
793!
794 else
795!
796!--------------------------- OROGW-solver of WAM2017
797!
798! sigres = max(sigmin, sigma(J))
799! if (heff/sigres.gt.dxres) sigres=heff/dxres
800! inv_b2eff = 0.5*sigres/heff
801! XLINV(I) = max(kxridge, inv_b2eff) ! 0.5*sigma(j)/heff = 1./Lridge
802 dtfac(:) = 1.0
803
804 call oro_wam_2017(im, km, npt, ipt, kref, kdt, me, master,
805 & dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsl,
806 & del, sigma, hprime, gamma, theta,
807 & sinlat, xlatd, taup, taud, pkdis)
808
809 endif ! oro_wam_2017 - LINSATDIS-solver of WAM-2017
810!
811!--------------------------- OROGW-solver of WAM2017
812!
813! TOFD as in BELJAARS-2004
814!
815! ---------------------------
816 IF( do_tofd ) then
817 axtms(:,:) = 0.0 ; aytms(:,:) = 0.0
818
819
820 DO i = 1,npt
821 j = ipt(i)
822 zpbl =rgrav*phil( j, kpbl(j) )
823
824 sigflt = min(sgh30(j), 0.3*hprime(j)) ! cannot exceed 30% of LS-SSO
825
826 zsurf = phii(j,1)*rgrav
827 do k=1,km
828 zpm(k) = phil(j,k)*rgrav
829 up1(k) = u1(j,k)
830 vp1(k) = v1(j,k)
831 enddo
832
833 call ugwpv0_tofd1d(km, sigflt, elvmaxd(j), zsurf, zpbl,
834 & up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
835
836 do k=1,km
837 axtms(j,k) = utofd1(k)
838 aytms(j,k) = vtofd1(k)
839!
840! add TOFD to GW-tendencies
841!
842 pdvdt(j,k) = pdvdt(j,k) + aytms(j,k)
843 pdudt(j,k) = pdudt(j,k) + axtms(j,k)
844 enddo
845!2018-diag
846 tau_tofd(j) = sum( utofd1(1:km)* del(j,1:km))
847 enddo
848 ENDIF ! do_tofd
849
850!---------------------------
851! combine oro-drag effects
852!---------------------------
853! + diag-3d
854
855 dudt_tms = axtms
856 tau_ogw = 0.
857 tau_mtb = 0.
858
859 DO k = 1,km
860 DO i = 1,npt
861 j = ipt(i)
862!
863 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
864!
865 if ( k < idxzb(i) .AND. idxzb(i) /= 0 ) then
866!
867! if blocking layers -- no OGWs
868!
869 dbim = db(i,k) / (1.+db(i,k)*dtp)
870 pdvdt(j,k) = - dbim * v1(j,k) +pdvdt(j,k)
871 pdudt(j,k) = - dbim * u1(j,k) +pdudt(j,k)
872 eng1 = eng0*(1.0-dbim*dtp)*(1.-dbim*dtp)
873
874 dusfc(j) = dusfc(j) - dbim * u1(j,k) * del(j,k)
875 dvsfc(j) = dvsfc(j) - dbim * v1(j,k) * del(j,k)
876!2018-diag
877 dudt_mtb(j,k) = -dbim * u1(j,k)
878 tau_mtb(j) = tau_mtb(j) + dudt_mtb(j,k)* del(j,k)
879
880 else
881!
882! OGW-s above blocking height
883!
884 taud(i,k) = taud(i,k) * dtfac(i)
885 dtaux = taud(i,k) * xn(i)
886 dtauy = taud(i,k) * yn(i)
887
888 pdvdt(j,k) = dtauy +pdvdt(j,k)
889 pdudt(j,k) = dtaux +pdudt(j,k)
890
891 unew = u1(j,k) + dtaux*dtp ! Pdudt(J,K)*DTP
892 vnew = v1(j,k) + dtauy*dtp ! Pdvdt(J,K)*DTP
893 eng1 = 0.5*(unew*unew + vnew*vnew)
894!
895 dusfc(j) = dusfc(j) + dtaux * del(j,k)
896 dvsfc(j) = dvsfc(j) + dtauy * del(j,k)
897!2018-diag
898 dudt_ogw(j,k) = dtaux
899 tau_ogw(j) = tau_ogw(j) +dtaux*del(j,k)
900 endif
901!
902! local energy deposition SSO-heat
903!
904 pdtdt(j,k) = max(eng0-eng1,0.)*rcpdt
905 ENDDO
906 ENDDO
907! dusfc w/o tofd sign as in the ERA-I, MERRA and CFSR
908 DO i = 1,npt
909 j = ipt(i)
910 dusfc(j) = -rgrav * dusfc(j)
911 dvsfc(j) = -rgrav * dvsfc(j)
912 tau_mtb(j) = -rgrav * tau_mtb(j)
913 tau_ogw(j) = -rgrav * tau_ogw(j)
914 tau_tofd(j) = -rgrav * tau_tofd(j)
915 ENDDO
916
917 RETURN
918
919 end subroutine gwdps_v0
920
921
922
923!===============================================================================
924!23456==============================================================================
925
929 subroutine fv3_ugwp_solv2_v0(klon, klev, dtime,
930 & tm1 , um1, vm1, qm1,
931 & prsl, prsi, philg, xlatd, sinlat, coslat,
932 & pdudt, pdvdt, pdtdt, dked, tau_ngw,
933 & mpi_id, master, kdt)
934!
935
936
937!=======================================================
938!
939! nov 2015 alternative gw-solver for nggps-wam
940! nov 2017 nh/rotational gw-modes for nh-fv3gfs
941! ---------------------------------------------------------------------------------
942!
943 use machine, only : kind_phys
944
945 use ugwp_common_v0 , only : rgrav, grav, cpd, rd, rv
946 &, omega2, rcpd2, pi, pi2, fv
947 &, rad_to_deg, deg_to_rad
948 &, rdi, gor, grcp, gocp
949 &, bnv2min, dw2min, velmin, gr2
950!
951 use ugwpv0_wmsdis_init, only : hpscale, rhp2, bv2min, gssec
952 &, v_kxw, v_kxw2, tamp_mpa, zfluxglob
953 &, maxdudt, gw_eff, dked_min
954 &, nslope, ilaunch, zmsi
955 &, zci, zdci, zci4, zci3, zci2
956 &, zaz_fct, zcosang, zsinang
957 &, nwav, nazd, zcimin, zcimax
958!
959 implicit none
960!23456
961
962 integer, parameter :: kp = kind_phys
963 integer, intent(in) :: klev ! vertical level
964 integer, intent(in) :: klon ! horiz tiles
965
966 real, intent(in) :: dtime ! model time step
967 real, intent(in) :: vm1(klon,klev) ! meridional wind
968 real, intent(in) :: um1(klon,klev) ! zonal wind
969 real, intent(in) :: qm1(klon,klev) ! spec. humidity
970 real, intent(in) :: tm1(klon,klev) ! kin temperature
971
972 real, intent(in) :: prsl(klon,klev) ! mid-layer pressure
973 real, intent(in) :: philg(klon,klev) ! m2/s2-phil => meters !!!!! phil =philg/grav
974 real, intent(in) :: prsi(klon,klev+1)! prsi interface pressure
975 real, intent(in) :: xlatd(klon) ! lat was in radians, now with xlat_d in degrees
976 real, intent(in) :: sinlat(klon)
977 real, intent(in) :: coslat(klon)
978 real, intent(in) :: tau_ngw(klon)
979
980 integer, intent(in) :: mpi_id, master, kdt
981!
982!
983! out-gw effects
984!
985 real, intent(out) :: pdudt(klon,klev) ! zonal momentum tendency
986 real, intent(out) :: pdvdt(klon,klev) ! meridional momentum tendency
987 real, intent(out) :: pdtdt(klon,klev) ! gw-heating (u*ax+v*ay)/cp
988 real, intent(out) :: dked(klon,klev) ! gw-eddy diffusion
989 real, parameter :: minvel = 0.5_kp !
990 real, parameter :: epsln = 1.0e-12_kp !
991 real, parameter :: zero = 0.0_kp, one = 1.0_kp, half = 0.5_kp
992
993!vay-2018
994
995 real :: taux(klon,klev+1) ! EW component of vertical momentum flux (pa)
996 real :: tauy(klon,klev+1) ! NS component of vertical momentum flux (pa)
997 real :: phil(klon,klev) ! gphil/grav
998!
999! local ===============================================================================================
1000!
1001
1002! real :: zthm1(klon,klev) ! temperature interface levels
1003 real :: zthm1 ! 1.0 / temperature interface levels
1004 real :: zbvfhm1(klon,ilaunch:klev) ! interface BV-frequency
1005 real :: zbn2(klon,ilaunch:klev) ! interface BV-frequency
1006 real :: zrhohm1(klon,ilaunch:klev) ! interface density
1007 real :: zuhm1(klon,ilaunch:klev) ! interface zonal wind
1008 real :: zvhm1(klon,ilaunch:klev) ! meridional wind
1009 real :: v_zmet(klon,ilaunch:klev)
1010 real :: vueff(klon,ilaunch:klev)
1011 real :: zbvfl(klon) ! BV at launch level
1012 real :: c2f2(klon)
1013
1014!23456
1015 real :: zul(klon,nazd) ! velocity in azimuthal direction at launch level
1016 real :: zci_min(klon,nazd)
1017! real :: zcrt(klon,klev,nazd) ! not used - do we need it? Moorthi
1018 real :: zact(klon, nwav, nazd) ! if =1 then critical level encountered => c-u
1019! real :: zacc(klon, nwav, nazd) ! not used!
1020!
1021 real :: zpu(klon,klev, nazd) ! momentum flux
1022! real :: zdfl(klon,klev, nazd)
1023 real :: zfct(klon,klev)
1024 real :: zfnorm(klon) ! normalisation factor
1025
1026 real :: zfluxlaun(klon)
1027 real :: zui(klon, klev,nazd)
1028!
1029 real :: zdfdz_v(klon,klev, nazd) ! axj = -df*rho/dz directional momentum depositiom
1030 real :: zflux(klon, nwav, nazd) ! momentum flux at each level stored as ( ix, mode, iazdim)
1031
1032 real :: zflux_z (klon, nwav,klev) !momentum flux at each azimuth stored as ( ix, mode, klev)
1033!
1034 real :: vm_zflx_mode, vc_zflx_mode
1035 real :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2
1036
1037! real :: zang, znorm, zang1, ztx
1038 real :: zu, zcin, zcpeak, zcin4, zbvfl4
1039 real :: zcin2, zbvfl2, zcin3, zbvfl3, zcinc
1040 real :: zatmp, zfluxs, zdep, zfluxsq, zulm, zdft, ze1, ze2
1041
1042!
1043 real :: zdelp,zrgpts
1044 real :: zthstd,zrhostd,zbvfstd
1045 real :: tvc1, tvm1, tem1, tem2, tem3
1046 real :: zhook_handle
1047 real :: delpi(klon,ilaunch:klev)
1048
1049
1050! real :: rcpd, grav2cpd
1051 real, parameter :: rcpdl = cpd/grav ! 1/[g/cp] == cp/g
1052 &, grav2cpd = grav/rcpdl ! g*(g/cp)= g^2/cp
1053 &, cpdi = one/cpd
1054
1055 real :: expdis, fdis
1056! real :: fmode, expdis, fdis
1057 real :: v_kzi, v_kzw, v_cdp, v_wdp, sc, tx1
1058
1059 integer :: j, k, inc, jk, jl, iazi
1060!
1061!--------------------------------------------------------------------------
1062!
1063 do k=1,klev
1064 do j=1,klon
1065 pdvdt(j,k) = zero
1066 pdudt(j,k) = zero
1067 pdtdt(j,k) = zero
1068 dked(j,k) = zero
1069 phil(j,k) = philg(j,k) * rgrav
1070 enddo
1071 enddo
1072
1073!=================================================
1074 do iazi=1, nazd
1075 do jk=1,klev
1076 do jl=1,klon
1077 zpu(jl,jk,iazi) = zero
1078! zcrt(jl,jk,iazi) = zero
1079! zdfl(jl,jk,iazi) = zero
1080 enddo
1081 enddo
1082 enddo
1083
1084!
1085! set initial min Cxi for critical level absorption
1086 do iazi=1,nazd
1087 do jl=1,klon
1088 zci_min(jl,iazi) = zcimin
1089 enddo
1090 enddo
1091! define half model level winds and temperature
1092! ---------------------------------------------
1093 do jk=max(ilaunch,2),klev
1094 do jl=1,klon
1095 tvc1 = tm1(jl,jk) * (one +fv*qm1(jl,jk))
1096 tvm1 = tm1(jl,jk-1) * (one +fv*qm1(jl,jk-1))
1097! zthm1(jl,jk) = 0.5 *(tvc1+tvm1)
1098 zthm1 = 2.0_kp / (tvc1+tvm1)
1099 zuhm1(jl,jk) = half *(um1(jl,jk-1)+um1(jl,jk))
1100 zvhm1(jl,jk) = half *(vm1(jl,jk-1)+vm1(jl,jk))
1101! zrhohm1(jl,jk) = prsi(jl,jk)*rdi/zthm1(jl,jk) ! rho = p/(RTv)
1102 zrhohm1(jl,jk) = prsi(jl,jk)*rdi*zthm1 ! rho = p/(RTv)
1103 zdelp = phil(jl,jk)-phil(jl,jk-1)
1104 v_zmet(jl,jk) = zdelp + zdelp
1105 delpi(jl,jk) = grav / (prsi(jl,jk-1) - prsi(jl,jk))
1106 vueff(jl,jk) =
1107 & 2.e-5_kp*exp( (phil(jl,jk)+phil(jl,jk-1))*rhp2)+dked_min
1108!
1109! zbn2(jl,jk) = grav2cpd/zthm1(jl,jk)
1110 zbn2(jl,jk) = grav2cpd*zthm1
1111 & * (one+rcpdl*(tm1(jl,jk)-tm1(jl,jk-1))/zdelp)
1112 zbn2(jl,jk) = max(min(zbn2(jl,jk), gssec), bv2min)
1113 zbvfhm1(jl,jk) = sqrt(zbn2(jl,jk)) ! bn = sqrt(bn2)
1114 enddo
1115 enddo
1116
1117 if (ilaunch == 1) then
1118 jk = 1
1119 do jl=1,klon
1120! zthm1(jl,jk) = tm1(jl,jk) * (1. +fv*qm1(jl,jk)) ! not used
1121 zuhm1(jl,jk) = um1(jl,jk)
1122 zvhm1(jl,jk) = vm1(jl,jk)
1123 zbvfhm1(jl,1) = zbvfhm1(jl,2)
1124 v_zmet(jl,1) = v_zmet(jl,2)
1125 vueff(jl,1) = dked_min
1126 zbn2(jl,1) = zbn2(jl,2)
1127 enddo
1128 endif
1129 do jl=1,klon
1130 tx1 = omega2 * sinlat(jl) / v_kxw
1131 c2f2(jl) = tx1 * tx1
1132 zbvfl(jl) = zbvfhm1(jl,ilaunch)
1133 enddo
1134!
1135! define intrinsic velocity (relative to launch level velocity) u(z)-u(zo), and coefficinets
1136! ------------------------------------------------------------------------------------------
1137 do iazi=1, nazd
1138 do jl=1,klon
1139 zul(jl,iazi) = zcosang(iazi) * zuhm1(jl,ilaunch)
1140 & + zsinang(iazi) * zvhm1(jl,ilaunch)
1141 enddo
1142 enddo
1143!
1144 do jk=ilaunch, klev-1 ! from z-launch up model level from which gw spectrum is launched
1145 do iazi=1, nazd
1146 do jl=1,klon
1147 zu = zcosang(iazi)*zuhm1(jl,jk)
1148 & + zsinang(iazi)*zvhm1(jl,jk)
1149 zui(jl,jk,iazi) = zu - zul(jl,iazi)
1150 enddo
1151 enddo
1152
1153 enddo
1154! define rho(zo)/n(zo)
1155! -------------------
1156 do jk=ilaunch, klev-1
1157 do jl=1,klon
1158 zfct(jl,jk) = zrhohm1(jl,jk) / zbvfhm1(jl,jk)
1159 enddo
1160 enddo
1161
1162! -----------------------------------------
1163! set launch momentum flux spectral density
1164! -----------------------------------------
1165
1166 if(nslope == 1) then ! s=1 case
1167 ! --------
1168 do inc=1,nwav
1169 zcin = zci(inc)
1170 zcin4 = zci4(inc)
1171 do jl=1,klon
1172!n4
1173 zbvfl4 = zbvfl(jl) * zbvfl(jl)
1174 zbvfl4 = zbvfl4 * zbvfl4
1175 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl4*zcin
1176 & / (zbvfl4+zcin4)
1177 enddo
1178 enddo
1179 elseif(nslope == 2) then ! s=2 case
1180 ! --------
1181 do inc=1, nwav
1182 zcin = zci(inc)
1183 zcin4 = zci4(inc)
1184 do jl=1,klon
1185 zbvfl4 = zbvfl(jl) * zbvfl(jl)
1186 zbvfl4 = zbvfl4 * zbvfl4
1187 zcpeak = zbvfl(jl) * zmsi
1188 zflux(jl,inc,1) = zfct(jl,ilaunch)*
1189 & zbvfl4*zcin*zcpeak/(zbvfl4*zcpeak+zcin4*zcin)
1190 enddo
1191 enddo
1192 elseif(nslope == -1) then ! s=-1 case
1193 ! --------
1194 do inc=1,nwav
1195 zcin = zci(inc)
1196 zcin2 = zci2(inc)
1197 do jl=1,klon
1198 zbvfl2 = zbvfl(jl)*zbvfl(jl)
1199 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl2*zcin
1200 & / (zbvfl2+zcin2)
1201 enddo
1202 enddo
1203 elseif(nslope == 0) then ! s=0 case
1204 ! --------
1205
1206 do inc=1, nwav
1207 zcin = zci(inc)
1208 zcin3 = zci3(inc)
1209 do jl=1,klon
1210 zbvfl3 = zbvfl(jl)**3
1211 zflux(jl,inc,1) = zfct(jl,ilaunch)*zbvfl3*zcin
1212 & / (zbvfl3+zcin3)
1213 enddo
1214 enddo
1215
1216 endif ! for slopes
1217!
1218! normalize momentum flux at the src-level
1219! ------------------------------
1220! integrate (zflux x dx)
1221 do inc=1, nwav
1222 zcinc = zdci(inc)
1223 do jl=1,klon
1224 zpu(jl,ilaunch,1) = zpu(jl,ilaunch,1) + zflux(jl,inc,1)*zcinc
1225 enddo
1226 enddo
1227!
1228! normalize and include lat-dep (precip or merra-2)
1229! -----------------------------------------------------------
1230! also other options to alter tropical values
1231!
1232 do jl=1,klon
1233 zfluxlaun(jl) = tau_ngw(jl) !*(.5+.75*coslat(JL)) !zfluxglob/2 on poles
1234 zfnorm(jl) = zfluxlaun(jl) / zpu(jl,ilaunch,1)
1235 enddo
1236!
1237 do iazi=1,nazd
1238 do jl=1,klon
1239 zpu(jl,ilaunch,iazi) = zfluxlaun(jl)
1240 enddo
1241 enddo
1242
1243! adjust constant zfct
1244
1245 do jk=ilaunch, klev-1
1246 do jl=1,klon
1247 zfct(jl,jk) = zfnorm(jl)*zfct(jl,jk)
1248 enddo
1249 enddo
1250! renormalize each spectral mode
1251
1252 do inc=1, nwav
1253 do jl=1,klon
1254 zflux(jl,inc,1) = zfnorm(jl)*zflux(jl,inc,1)
1255 enddo
1256 enddo
1257
1258! copy zflux into all other azimuths
1259! --------------------------------
1260! zact(:,:,:) = one ; zacc(:,:,:) = one
1261 zact(:,:,:) = one
1262 do iazi=2, nazd
1263 do inc=1,nwav
1264 do jl=1,klon
1265 zflux(jl,inc,iazi) = zflux(jl,inc,1)
1266 enddo
1267 enddo
1268 enddo
1269
1270! -------------------------------------------------------------
1271! azimuth do-loop
1272! --------------------
1273 do iazi=1, nazd
1274
1275! write(0,*)' iazi=',iazi,' ilaunch=',ilaunch
1276! vertical do-loop
1277! ----------------
1278 do jk=ilaunch, klev-1
1279! first check for critical levels
1280! ------------------------
1281 do jl=1,klon
1282 zci_min(jl,iazi) = max(zci_min(jl,iazi),zui(jl,jk,iazi))
1283 enddo
1284! set zact to zero if critical level encountered
1285! ----------------------------------------------
1286 do inc=1, nwav
1287! zcin = zci(inc)
1288 do jl=1,klon
1289! zatmp = minvel + sign(minvel,zcin-zci_min(jl,iazi))
1290! zacc(jl,inc,iazi) = zact(jl,inc,iazi)-zatmp
1291! zact(jl,inc,iazi) = zatmp
1292 zact(jl,inc,iazi) = minvel
1293 & + sign(minvel,zci(inc)-zci_min(jl,iazi))
1294 enddo
1295 enddo
1296!
1297! zdfl not used! - do we need it? Moorthi
1298! integrate to get critical-level contribution to mom deposition
1299! ---------------------------------------------------------------
1300! do inc=1, nwav
1301! zcinc = zdci(inc)
1302! do jl=1,klon
1303! zdfl(jl,jk,iazi) = zdfl(jl,jk,iazi) +
1304! & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zcinc
1305! enddo
1306! enddo
1307! --------------------------------------------
1308! get weighted average of phase speed in layer zcrt is not used - do we need it? Moorthi
1309! --------------------------------------------
1310! do jl=1,klon
1311! write(0,*)' jk=',jk,' jl=',jl,' iazi=',iazi, zdfl(jl,jk,iazi)
1312! if(zdfl(jl,jk,iazi) > epsln ) then
1313! zatmp = zcrt(jl,jk,iazi)
1314! do inc=1, nwav
1315! zatmp = zatmp + zci(inc) *
1316! & zacc(jl,inc,iazi)*zflux(jl,inc,iazi)*zdci(inc)
1317! enddo
1318!
1319! zcrt(jl,jk,iazi) = zatmp / zdfl(jl,jk,iazi)
1320! else
1321! zcrt(jl,jk,iazi) = zcrt(jl,jk-1,iazi)
1322! endif
1323! enddo
1324
1325!
1326 do inc=1, nwav
1327 zcin = zci(inc)
1328 if (abs(zcin) > epsln) then
1329 zcinc = one / zcin
1330 else
1331 zcinc = one
1332 endif
1333 do jl=1,klon
1334!=======================================================================
1335! saturated limit wfit = kzw*kzw*kt; wfdt = wfit/(kxw*cx)*betat
1336! & dissipative kzi = 2.*kzw*(wfdm+wfdt)*dzpi(k)
1337! define kxw =
1338!=======================================================================
1339 v_cdp = abs(zcin-zui(jl,jk,iazi))
1340 v_wdp = v_kxw*v_cdp
1341 wdop2 = v_wdp* v_wdp
1342 cdf2 = v_cdp*v_cdp - c2f2(jl)
1343 if (cdf2 > zero) then
1344 kzw2 = (zbn2(jl,jk)-wdop2)/cdf2 - v_kxw2
1345 else
1346 kzw2 = zero
1347 endif
1348 if ( kzw2 > zero .and. cdf2 > zero) then
1349 v_kzw = sqrt(kzw2)
1350!
1351!linsatdis: kzw2, kzw3, kdsat, c2f2, cdf2, cdf1
1352!
1353!kzw2 = (zBn2(k)-wdop2)/Cdf2 - rhp4 - v_kx2w ! full lin DS-NiGW (N2-wd2)*k2=(m2+k2+[1/2H]^2)*(wd2-f2)
1354! Kds = kxw*Cdf1*rhp2/kzw3
1355!
1356 v_cdp = sqrt( cdf2 )
1357 v_wdp = v_kxw * v_cdp
1358 v_kzi = abs(v_kzw*v_kzw*vueff(jl,jk)/v_wdp*v_kzw)
1359 expdis = exp(-v_kzi*v_zmet(jl,jk))
1360 else
1361 v_kzi = zero
1362 expdis = one
1363 v_kzw = zero
1364 v_cdp = zero ! no effects of reflected waves
1365 endif
1366
1367! fmode = zflux(jl,inc,iazi)
1368! fdis = fmode*expdis
1369 fdis = expdis * zflux(jl,inc,iazi)
1370!
1371! saturated flux + wave dissipation - Keddy_gwsat in UGWP-V1
1372! linsatdis = 1.0 , here: u'^2 ~ linsatdis* [v_cdp*v_cdp]
1373!
1374 zfluxs = zfct(jl,jk)*v_cdp*v_cdp*zcinc
1375!
1376! zfluxs= zfct(jl,jk)*(zcin-zui(jl,jk,iazi))**2/zcin
1377! flux_tot - sat.flux
1378!
1379 zdep = zact(jl,inc,iazi)* (fdis-zfluxs)
1380 if(zdep > zero ) then
1381! subs on sat-limit
1382 zflux(jl,inc,iazi) = zfluxs
1383 zflux_z(jl,inc,jk) = zfluxs
1384 else
1385! assign dis-ve flux
1386 zflux(jl,inc,iazi) = fdis
1387 zflux_z(jl,inc,jk) = fdis
1388 endif
1389 enddo
1390 enddo
1391!
1392! integrate over spectral modes zpu(y, z, azimuth) zact(jl,inc,iazi)*zflux(jl,inc,iazi)*[d("zcinc")]
1393!
1394 zdfdz_v(:,jk,iazi) = zero
1395
1396 do inc=1, nwav
1397 zcinc = zdci(inc) ! dc-integration
1398 do jl=1,klon
1399 vc_zflx_mode = zact(jl,inc,iazi)*zflux(jl,inc,iazi)
1400 zpu(jl,jk,iazi) = zpu(jl,jk,iazi) + vc_zflx_mode*zcinc
1401
1402!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1403! check monotonic decrease
1404! (heat deposition integration over spectral mode for each azimuth
1405! later sum over selected azimuths as "non-negative" scalars)
1406!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1407 if (jk > ilaunch)then
1408! zdelp = grav/(prsi(jl,jk-1)-prsi(jl,jk))*
1409! & abs(zcin-zui(jl,jk,iazi)) *zcinc
1410 zdelp = delpi(jl,jk) * abs(zcin-zui(jl,jk,iazi)) *zcinc
1411 vm_zflx_mode = zact(jl,inc,iazi)* zflux_z(jl,inc,jk-1)
1412
1413 if (vc_zflx_mode > vm_zflx_mode)
1414 & vc_zflx_mode = vm_zflx_mode ! no-flux increase
1415 zdfdz_v( jl,jk,iazi) = zdfdz_v( jl,jk,iazi) +
1416 & (vm_zflx_mode-vc_zflx_mode)*zdelp ! heating >0
1417!
1418!
1419 endif
1420 enddo !jl=1,klon
1421 enddo !waves inc=1,nwav
1422
1423! --------------
1424 enddo ! end jk do-loop vertical loop
1425! ---------------
1426 enddo ! end nazd do-loop
1427! ----------------------------------------------------------------------------
1428! sum contribution for total zonal and meridional flux +
1429! energy dissipation
1430! ---------------------------------------------------
1431!
1432 do jk=1,klev+1
1433 do jl=1,klon
1434 taux(jl,jk) = zero
1435 tauy(jl,jk) = zero
1436 enddo
1437 enddo
1438
1439 tem3 = zaz_fct*cpdi
1440 do iazi=1,nazd
1441 tem1 = zaz_fct*zcosang(iazi)
1442 tem2 = zaz_fct*zsinang(iazi)
1443 do jk=ilaunch, klev-1
1444 do jl=1,klon
1445 taux(jl,jk) = taux(jl,jk) + tem1 * zpu(jl,jk,iazi) ! zaz_fct - "azimuth"-norm-n
1446 tauy(jl,jk) = tauy(jl,jk) + tem2 * zpu(jl,jk,iazi)
1447 pdtdt(jl,jk) = pdtdt(jl,jk) + tem3 * zdfdz_v(jl,jk,iazi) ! eps_dis =sum( +d(flux_e)/dz) > 0.
1448 enddo
1449 enddo
1450
1451 enddo
1452!
1453! update du/dt and dv/dt tendencies ..... no contribution to heating => keddy/tracer-mom-heat
1454! ----------------------------
1455!
1456
1457 do jk=ilaunch,klev
1458 do jl=1, klon
1459! zdelp = grav / (prsi(jl,jk-1)-prsi(jl,jk))
1460 zdelp = delpi(jl,jk)
1461 ze1 = (taux(jl,jk)-taux(jl,jk-1))*zdelp
1462 ze2 = (tauy(jl,jk)-tauy(jl,jk-1))*zdelp
1463 if (abs(ze1) >= maxdudt ) then
1464 ze1 = sign(maxdudt, ze1)
1465 endif
1466 if (abs(ze2) >= maxdudt ) then
1467 ze2 = sign(maxdudt, ze2)
1468 endif
1469 pdudt(jl,jk) = -ze1
1470 pdvdt(jl,jk) = -ze2
1471!
1472! Cx =0 based Cx=/= 0. above
1473!
1474 pdtdt(jl,jk) = (ze1*um1(jl,jk) + ze2*vm1(jl,jk)) * cpdi
1475!
1476 dked(jl,jk) = max(dked_min, pdtdt(jl,jk)/zbn2(jl,jk))
1477! if (dked(jl,jk) < 0) dked(jl,jk) = dked_min
1478 enddo
1479 enddo
1480!
1481! add limiters/efficiency for "unbalanced ics" if it is needed
1482!
1483 do jk=ilaunch,klev
1484 do jl=1, klon
1485 pdudt(jl,jk) = gw_eff * pdudt(jl,jk)
1486 pdvdt(jl,jk) = gw_eff * pdvdt(jl,jk)
1487 pdtdt(jl,jk) = gw_eff * pdtdt(jl,jk)
1488 dked(jl,jk) = gw_eff * dked(jl,jk)
1489 enddo
1490 enddo
1491!
1492!---------------------------------------------------------------------------
1493 return
1494 end subroutine fv3_ugwp_solv2_v0
1495
1496 end module ugwp_driver_v0
cires_orowam2017
This module includes the OROGW solver of WAM2017.
Definition cires_orowam2017.f:5
cires_ugwpv0_module
This module contains the UGWPv0 driver.
Definition cires_ugwp_module.F90:5
ugwp_common_v0
This module contains UGWP v0 initialization schemes.
Definition cires_ugwp_initialize.F90:10
ugwp_driver_v0
This module contains the UGWP v0 driver module.
Definition ugwp_driver_v0.F:4
ugwpv0_oro_init
This module contains orographic wave source schemes for UGWP v0.
Definition cires_ugwp_initialize.F90:99
ugwpv0_wmsdis_init
This module contains init-solvers for "broad" non-stationary multi-wave spectra.
Definition cires_ugwp_initialize.F90:276