CCPP SciDoc v7.0.0  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!=====================================================================
31
34 SUBROUTINE gwdps_v0(IM, km, imx, do_tofd,
35 & Pdvdt, Pdudt, Pdtdt, Pkdis, U1,V1,T1,Q1,KPBL,
36 & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DTP,KDT,
37 & sgh30, HPRIME,OC,OA4,CLX4,THETA,vSIGMA,vGAMMA,ELVMAXD,
38 & DUSFC, DVSFC, xlatd, sinlat, coslat, sparea,
39 & cdmbgwd, me, master, rdxzb, con_g, con_omega,
40 & zmtb, zogw, tau_mtb, tau_ogw, tau_tofd,
41 & dudt_mtb, dudt_ogw, dudt_tms)
42!----------------------------------------
43! ugwp_v0
44!
45! modified/revised version of gwdps.f (with bug fixes, tofd, appropriate
46! computation of reference level for OGW + COORDE diagnostics
47! all constants/parameters inside cires_ugwp_initialize.F90
48!----------------------------------------
49
50 USE machine , ONLY : kind_phys
51 use ugwp_common_v0,only : rgrav, grav, cpd, rd, rv, rcpd, rcpd2
52 &, pi, rad_to_deg, deg_to_rad, pi2
53 &, rdi, gor, grcp, gocp, fv, gr2
54 &, bnv2min, dw2min, velmin, arad
55
56 use ugwpv0_oro_init, only : rimin, ric, efmin, efmax
57 &, hpmax, hpmin, sigfaci => sigfac
58 &, dpmin, minwnd, hminmt, hncrit
59 &, rlolev, gmax, veleps, factop
60 &, frc, ce, ceofrc, frmax, cg
61 &, fdir, mdir, nwdir
62 &, cdmb, cleff, fcrit_gfs, fcrit_mtb
63 &, n_tofd, ze_tofd, ztop_tofd
64
65 use cires_ugwpv0_module, only : kxw, max_kdis, max_axyz
66
67!----------------------------------------
68 implicit none
69 integer, parameter :: kp = kind_phys
70 character(len=8) :: strsolver='PSS-1986' ! current operational solver or 'WAM-2017'
71 integer, intent(in) :: im, km, imx, kdt
72 integer, intent(in) :: me, master
73 logical, intent(in) :: do_tofd
74 real(kind=kind_phys), parameter :: sigfac = 3, sigfacs = 0.5
75 real(kind=kind_phys) :: ztoph,zlowh,ph_blk, dz_blk
76 integer, intent(in) :: KPBL(IM) ! Index for the PBL top layer!
77 real(kind=kind_phys), intent(in) :: dtp ! time step
78 real(kind=kind_phys), intent(in) :: cdmbgwd(2)
79
80 real(kind=kind_phys), intent(in), dimension(im,km) ::
81 & u1, v1, t1, q1,
82 & del, prsl, prslk, phil
83 real(kind=kind_phys), intent(in),dimension(im,km+1):: prsi, phii
84 real(kind=kind_phys), intent(in) :: xlatd(im),sinlat(im),
85 & coslat(im)
86 real(kind=kind_phys), intent(in) :: sparea(im)
87
88 real(kind=kind_phys), intent(in) :: oc(im), oa4(im,4), clx4(im,4)
89 real(kind=kind_phys), intent(in) :: hprime(im), sgh30(im)
90 real(kind=kind_phys), intent(in) :: elvmaxd(im), theta(im)
91 real(kind=kind_phys), intent(in) :: vsigma(im), vgamma(im)
92 real(kind=kind_phys) :: sigma(im), gamma(im)
93
94 real(kind=kind_phys), intent(in) :: con_g, con_omega
95
96!output -phys-tend
97 real(kind=kind_phys),dimension(im,km),intent(out) ::
98 & pdvdt, pdudt, pkdis, pdtdt
99! output - diag-coorde
100 &, dudt_mtb, dudt_ogw, dudt_tms
101!
102 real(kind=kind_phys),dimension(im) :: rdxzb, zmtb, zogw
103 &, tau_ogw, tau_mtb, tau_tofd
104 &, dusfc, dvsfc
105!
106!---------------------------------------------------------------------
107! # of permissible sub-grid orography hills for "any" resolution < 25
108! correction for "elliptical" hills based on shilmin-area =sgrid/25
109! 4.*gamma*b_ell*b_ell >= shilmin
110! give us limits on [b_ell & gamma *b_ell] > 5 km =sso_min
111! gamma_min = 1/4*shilmin/sso_min/sso_min
112!23.01.2019: cdmb = 4.*192/768_c192=1 x 0.5
113! 192: cdmbgwd = 0.5, 2.5
114! cleff = 2.5*0.5e-5 * sqrt(192./768.) => Lh_eff = 1004. km
115! 6*dx = 240 km 8*dx = 320. ~ 3-5 more effective
116!---------------------------------------------------------------------
117 real(kind=kind_phys) :: gammin = 0.00999999_kp
118 real(kind=kind_phys), parameter :: nhilmax = 25.0_kp
119 real(kind=kind_phys), parameter :: sso_min = 3000.0_kp
120 logical, parameter :: do_adjoro = .true.
121!
122 real(kind=kind_phys) :: shilmin, sgrmax, sgrmin
123 real(kind=kind_phys) :: belpmin, dsmin, dsmax
124! real(kind=kind_phys) :: arhills(im) ! not used why do we need?
125 real(kind=kind_phys) :: xlingfs
126
127!
128! locals
129! mean flow
130 real(kind=kind_phys), dimension(im,km) :: ri_n, bnv2, ro
131 &, vtk, vtj, velco
132!mtb
133 real(kind=kind_phys), dimension(im) :: oa, clx , elvmax, wk
134 &, pe, ek, up
135
136 real(kind=kind_phys), dimension(im,km) :: db, ang, uds
137
138 real(kind=kind_phys) :: zlen, dbtmp, r, phiang, dbim, zr
139 real(kind=kind_phys) :: eng0, eng1, cosang2, sinang2
140 real(kind=kind_phys) :: bgam, cgam, gam2, rnom, rdem
141!
142! TOFD
143! Some constants now in "use ugwp_oro_init" + "use ugwp_common"
144!
145!==================
146 real(kind=kind_phys) :: unew, vnew, zpbl, sigflt, zsurf
147 real(kind=kind_phys), dimension(km) :: utofd1, vtofd1
148 &, epstofd1, krf_tofd1
149 &, up1, vp1, zpm
150 real(kind=kind_phys),dimension(im, km) :: axtms, aytms
151!
152! OGW
153!
154 LOGICAL ICRILV(IM)
155!
156 real(kind=kind_phys), dimension(im) :: xn, yn, ubar, vbar, ulow,
157 & roll, bnv2bar, scor, dtfac, xlinv, delks, delks1
158!
159 real(kind=kind_phys) :: taup(im,km+1), taud(im,km)
160 real(kind=kind_phys) :: taub(im), taulin(im), heff, hsat, hdis
161
162 integer, dimension(im) :: kref, idxzb, ipt, kreflm,
163 & iwklm, iwk, izlow
164!
165!check what we need
166!
167 real(kind=kind_phys) :: bnv, fr, ri_gw
168 &, brvf, tem, tem1, tem2, temc, temv
169 &, ti, rdz, dw2, shr2, bvf2
170 &, rdelks, efact, coefm, gfobnv
171 &, scork, rscor, hd, fro, sira
172 &, dtaux, dtauy, pkp1log, pklog
173 &, grav2, rcpdt, windik, wdir
174 &, sigmin, dxres,sigres,hdxres
175 &, cdmb4, mtbridge
176 &, kxridge, inv_b2eff, zw1, zw2
177 &, belps, aelps, nhills, selps
178
179 integer :: kmm1, kmm2, lcap, lcapp1
180 &, npt, kbps, kbpsp1,kbpsm1
181 &, kmps, idir, nwd, klcap, kp1, kmpbl, kmll
182 &, k_mtb, k_zlow, ktrial, klevm1, i, j, k
183!
184 rcpdt = 1.0 / (cpd*dtp)
185 grav2 = grav + grav
186!
187! mtb-blocking sigma_min and dxres => cires_initialize
188!
189 sgrmax = maxval(sparea) ; sgrmin = minval(sparea)
190 dsmax = sqrt(sgrmax) ; dsmin = sqrt(sgrmin)
191
192 dxres = pi2*arad/float(imx)
193 hdxres = 0.5_kp*dxres
194! shilmin = sgrmin/nhilmax ! not used - Moorthi
195
196! gammin = min(sso_min/dsmax, 1.) ! Moorthi - with this results are not reproducible
197 gammin = min(sso_min/dxres, 1.) ! Moorthi
198
199! sigmin = 2.*hpmin/dsmax !dxres ! Moorthi - this will not reproduce
200 sigmin = 2.*hpmin/dxres !dxres
201
202
203
204 kxridge = float(imx)/arad * cdmbgwd(2)
205
206 do i=1,im
207 idxzb(i) = 0
208 zmtb(i) = 0.0
209 zogw(i) = 0.0
210 rdxzb(i) = 0.0
211 tau_ogw(i) = 0.0
212 tau_mtb(i) = 0.0
213 dusfc(i) = 0.0
214 dvsfc(i) = 0.0
215 tau_tofd(i) = 0.0
216!
217 ipt(i) = 0
218 sigma(i) = max(vsigma(i), sigmin)
219 gamma(i) = max(vgamma(i), gammin)
220 enddo
221
222 do k=1,km
223 do i=1,im
224 pdvdt(i,k) = 0.0
225 pdudt(i,k) = 0.0
226 pdtdt(i,k) = 0.0
227 pkdis(i,k) = 0.0
228 dudt_mtb(i,k) = 0.0
229 dudt_ogw(i,k) = 0.0
230 dudt_tms(i,k) = 0.0
231 enddo
232 enddo
233
234! ---- for lm and gwd calculation points
235
236 npt = 0
237 do i = 1,im
238 if ( elvmaxd(i) >= hminmt .and. hprime(i) >= hpmin ) then
239
240 npt = npt + 1
241 ipt(npt) = i
242! arhills(i) = 1.0
243!
244 sigres = max(sigmin, sigma(i))
245! if (sigma(i) < sigmin) sigma(i)= sigmin
246 dxres = sqrt(sparea(i))
247 if (2.*hprime(i)/sigres > dxres) sigres=2.*hprime(i)/dxres
248 aelps = min(2.*hprime(i)/sigres, 0.5*dxres)
249 if (gamma(i) > 0.0 ) belps = min(aelps/gamma(i),.5*dxres)
250!
251! small-scale "turbulent" oro-scales < sso_min
252!
253 if( aelps < sso_min .and. do_adjoro) then
254
255! a, b > sso_min upscale ellipse a/b > 0.1 a>sso_min & h/b=>new_sigm
256!
257 aelps = sso_min
258 if (belps < sso_min ) then
259 gamma(i) = 1.0
260 belps = aelps*gamma(i)
261 else
262 gamma(i) = min(aelps/belps, 1.0)
263 endif
264 sigma(i) = 2.*hprime(i)/aelps
265 gamma(i) = min(aelps/belps, 1.0)
266 endif
267
268 selps = belps*belps*gamma(i)*4. ! ellipse area of the el-c hill
269 nhills = min(nhilmax, sparea(i)/selps)
270! arhills(i) = max(nhills, 1.0)
271
272!333 format( ' nhil: ', I6, 4(2x, F9.3), 2(2x, E9.3))
273! if (kdt==1 )
274! & write(6,333) nint(nhills)+1,xlatd(i), hprime(i),aelps*1.e-3,
275! & belps*1.e-3, sigma(i),gamma(i)
276
277 endif
278 enddo
279
280 IF (npt == 0) then
281 RETURN ! No gwd/mb calculation done
282 endif
283
284
285 do i=1,npt
286 iwklm(i) = 2
287 idxzb(i) = 0
288 kreflm(i) = 0
289 enddo
290
291 do k=1,km
292 do i=1,im
293 db(i,k) = 0.0
294 ang(i,k) = 0.0
295 uds(i,k) = 0.0
296 enddo
297 enddo
298
299 kmm1 = km - 1 ; kmm2 = km - 2 ; kmll = kmm1
300 lcap = km ; lcapp1 = lcap + 1
301
302 DO i = 1, npt
303 j = ipt(i)
304 elvmax(j) = min(elvmaxd(j)*0. + sigfac * hprime(j), hncrit)
305 izlow(i) = 1 ! surface-level
306 ENDDO
307!
308 DO k = 1, kmm1
309 DO i = 1, npt
310 j = ipt(i)
311 ztoph = sigfac * hprime(j)
312 zlowh = sigfacs* hprime(j)
313 pkp1log = phil(j,k+1) * rgrav
314 pklog = phil(j,k) * rgrav
315! if (( ELVMAX(j) <= pkp1log) .and. (ELVMAX(j).ge.pklog) )
316! & iwklm(I) = MAX(iwklm(I), k+1 )
317 if (( ztoph <= pkp1log) .and. (ztoph >= pklog) )
318 & iwklm(i) = max(iwklm(i), k+1 )
319!
320 if (zlowh <= pkp1log .and. zlowh >= pklog)
321 & izlow(i) = max(izlow(i),k)
322 ENDDO
323 ENDDO
324!
325 DO k = 1,km
326 DO i =1,npt
327 j = ipt(i)
328 vtj(i,k) = t1(j,k) * (1.+fv*q1(j,k))
329 vtk(i,k) = vtj(i,k) / prslk(j,k)
330 ro(i,k) = rdi * prsl(j,k) / vtj(i,k) ! DENSITY mid-levels
331 taup(i,k) = 0.0
332 ENDDO
333 ENDDO
334!
335! check RI_N or RI_MF computation
336!
337 DO k = 1,kmm1
338 DO i =1,npt
339 j = ipt(i)
340 rdz = grav / (phil(j,k+1) - phil(j,k))
341 tem1 = u1(j,k) - u1(j,k+1)
342 tem2 = v1(j,k) - v1(j,k+1)
343 dw2 = tem1*tem1 + tem2*tem2
344 shr2 = max(dw2,dw2min) * rdz * rdz
345! TI = 2.0 / (T1(J,K)+T1(J,K+1))
346! BVF2 = Grav*(GOCP+RDZ*(VTJ(I,K+1)-VTJ(I,K)))* TI
347! RI_N(I,K) = MAX(BVF2/SHR2,RIMIN) ! Richardson number
348!
349 bvf2 = grav2 * rdz * (vtk(i,k+1)-vtk(i,k))
350 & / (vtk(i,k+1)+vtk(i,k))
351 bnv2(i,k+1) = max( bvf2, bnv2min )
352 ri_n(i,k+1) = bnv2(i,k)/shr2 ! Richardson number consistent with BNV2
353!
354! add here computation for Ktur and OGW-dissipation fro VE-GFS
355!
356 ENDDO
357 ENDDO
358 k = 1
359 DO i = 1, npt
360 bnv2(i,k) = bnv2(i,k+1)
361 ENDDO
362!
363! level iwklm =>phil(j,k)/g < sigfac * hprime(j) < phil(j,k+1)/g
364!
365 DO i = 1, npt
366 j = ipt(i)
367 k_zlow = izlow(i)
368 if (k_zlow == iwklm(i)) k_zlow = 1
369 delks(i) = 1.0 / (prsi(j,k_zlow) - prsi(j,iwklm(i)))
370! DELKS1(I) = 1.0 /(PRSL(J,k_zlow) - PRSL(J,iwklm(i)))
371 ubar(i) = 0.0
372 vbar(i) = 0.0
373 roll(i) = 0.0
374 pe(i) = 0.0
375 ek(i) = 0.0
376 bnv2bar(i) = 0.0
377 ENDDO
378!
379 DO i = 1, npt
380 k_zlow = izlow(i)
381 if (k_zlow == iwklm(i)) k_zlow = 1
382 DO k = k_zlow, iwklm(i)-1 ! Kreflm(I)= iwklm(I)-1
383 j = ipt(i) ! laye-aver Rho, U, V
384 rdelks = del(j,k) * delks(i)
385 ubar(i) = ubar(i) + rdelks * u1(j,k) ! trial Mean U below
386 vbar(i) = vbar(i) + rdelks * v1(j,k) ! trial Mean V below
387 roll(i) = roll(i) + rdelks * ro(i,k) ! trial Mean RO below
388!
389 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
390 ENDDO
391 ENDDO
392!
393 DO i = 1, npt
394 j = ipt(i)
395!
396! integrate from Ztoph = sigfac*hprime down to Zblk if exists
397! find ph_blk, dz_blk like in LM-97 and IFS
398!
399 ph_blk =0.
400 DO k = iwklm(i), 1, -1
401 phiang = atan2(v1(j,k),u1(j,k))*rad_to_deg
402 ang(i,k) = ( theta(j) - phiang )
403 if ( ang(i,k) > 90. ) ang(i,k) = ang(i,k) - 180.
404 if ( ang(i,k) < -90. ) ang(i,k) = ang(i,k) + 180.
405 ang(i,k) = ang(i,k) * deg_to_rad
406 uds(i,k) =
407 & max(sqrt(u1(j,k)*u1(j,k) + v1(j,k)*v1(j,k)), velmin)
408!
409 IF (idxzb(i) == 0 ) then
410 dz_blk = ( phii(j,k+1) - phii(j,k) ) *rgrav
411 pe(i) = pe(i) + bnv2(i,k) *
412 & ( elvmax(j) - phil(j,k)*rgrav ) * dz_blk
413
414 up(i) = max(uds(i,k) * cos(ang(i,k)), velmin)
415 ek(i) = 0.5 * up(i) * up(i)
416
417 ph_blk = ph_blk + dz_blk*sqrt(bnv2(i,k))/up(i)
418
419! --- Dividing Stream lime is found when PE =exceeds EK. oper-l GFS
420! IF ( PE(I) >= EK(I) ) THEN
421 IF ( ph_blk >= fcrit_gfs ) THEN
422 idxzb(i) = k
423 zmtb(j) = phil(j, k)*rgrav
424 rdxzb(j) = real(k, kind=kind_phys)
425 ENDIF
426
427 ENDIF
428 ENDDO
429!
430! Alternative expression: ZMTB = max(Heff*(1. -Fcrit_gfs/Fr), 0)
431! fcrit_gfs/fr
432!
433 goto 788
434
435 bnv = sqrt( bnv2bar(i) )
436 heff = 2.*min(hprime(j),hpmax)
437 zw2 = ubar(i)*ubar(i)+vbar(i)*vbar(i)
438 ulow(i) = sqrt(max(zw2,dw2min))
439 fr = heff*bnv/ulow(i)
440 zw1 = max(heff*(1. -fcrit_gfs/fr), 0.0)
441 zw2 = phil(j,2)*rgrav
442 if (fr > fcrit_gfs .and. zw1 > zw2 ) then
443 do k=2, kmm1
444 pkp1log = phil(j,k+1) * rgrav
445 pklog = phil(j,k) * rgrav
446 if (zw1 <= pkp1log .and. zw1 >= pklog) exit
447 enddo
448 idxzb(i) = k
449 zmtb(j) = phil(j, k)*rgrav
450 else
451 zmtb(j) = 0.
452 idxzb(i) = 0
453 endif
454788 continue
455 ENDDO
456
457!
458! --- The drag for mtn blocked flow
459!
460 cdmb4 = 0.25*cdmb
461 DO i = 1, npt
462 j = ipt(i)
463!
464 IF ( idxzb(i) > 0 ) then
465! (4.16)-IFS
466 gam2 = gamma(j)*gamma(j)
467 bgam = 1.0 - 0.18*gamma(j) - 0.04*gam2
468 cgam = 0.48*gamma(j) + 0.30*gam2
469 DO k = idxzb(i)-1, 1, -1
470
471 zlen = sqrt( ( phil(j,idxzb(i)) - phil(j,k) ) /
472 & ( phil(j,k ) + grav * hprime(j) ) )
473
474 tem = cos(ang(i,k))
475 cosang2 = tem * tem
476 sinang2 = 1.0 - cosang2
477!
478! cos =1 sin =0 => 1/R= gam ZR = 2.-gam
479! cos =0 sin =1 => 1/R= 1/gam ZR = 2.- 1/gam
480!
481 rdem = cosang2 + gam2 * sinang2
482 rnom = cosang2*gam2 + sinang2
483!
484! metOffice Dec 2010
485! correction of H. Wells & A. Zadra for the
486! aspect ratio of the hill seen by MF
487! (1/R , R-inverse below: 2-R)
488
489 rdem = max(rdem, 1.e-6)
490 r = sqrt(rnom/rdem)
491 zr = max( 2. - r, 0. )
492
493 sigres = max(sigmin, sigma(j))
494 if (hprime(j)/sigres > dxres) sigres = hprime(j)/dxres
495 mtbridge = zr * sigres*zlen / hprime(j)
496! (4.15)-IFS
497! DBTMP = CDmb4 * mtbridge *
498! & MAX(cos(ANG(I,K)), gamma(J)*sin(ANG(I,K)))
499! (4.16)-IFS
500 dbtmp = cdmb4*mtbridge*(bgam* cosang2 +cgam* sinang2)
501 db(i,k)= dbtmp * uds(i,k)
502 ENDDO
503!
504 endif
505 ENDDO
506!
507!.............................
508!.............................
509! end mtn blocking section
510!.............................
511!.............................
512!
513!--- Orographic Gravity Wave Drag Section
514!
515! Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62
516! inside "cires_ugwp_initialize.F90" now
517!
518 kmpbl = km / 2
519 iwk(1:npt) = 2
520!
521! METO-scheme:
522! k_mtb = max(k_zmtb, k_n*hprime/2] to reduce diurnal variations taub_ogw
523!
524 DO k=3,kmpbl
525 DO i=1,npt
526 j = ipt(i)
527 tem = (prsi(j,1) - prsi(j,k))
528 if (tem < dpmin) iwk(i) = k ! dpmin=50 mb
529
530!===============================================================
531! lev=111 t=311.749 hkm=0.430522 Ps-P(iwk)=52.8958
532! below "Hprime" - source of OGWs and below Zblk !!!
533! 27 2 kpbl ~ 1-2 km < Hprime
534!===============================================================
535 enddo
536 enddo
537!
538! iwk - adhoc GFS-parameter to select OGW-launch level between
539! LEVEL ~0.4-0.5 KM from surface or/and PBL-top
540! in UGWP-V1: options to modify as Htop ~ (2-3)*Hprime > Zmtb
541! in UGWP-V0 we ensured that : Zogw > Zmtb
542!
543
544 kbps = 1
545 kmps = km
546 k_mtb = 1
547 DO i=1,npt
548 j = ipt(i)
549 k_mtb = max(1, idxzb(i))
550
551 kref(i) = max(iwk(i), kpbl(j)+1 ) ! reference level PBL or smt-else ????
552 kref(i) = max(kref(i), iwklm(i) ) ! iwklm => sigfac*hprime
553
554 if (kref(i) <= idxzb(i)) kref(i) = idxzb(i) + 1 ! layer above zmtb
555 kbps = max(kbps, kref(i))
556 kmps = min(kmps, kref(i))
557!
558 delks(i) = 1.0 / (prsi(j,k_mtb) - prsi(j,kref(i)))
559 ubar(i) = 0.0
560 vbar(i) = 0.0
561 roll(i) = 0.0
562 bnv2bar(i)= 0.0
563 ENDDO
564!
565 kbpsp1 = kbps + 1
566 kbpsm1 = kbps - 1
567 k_mtb = 1
568!
569 DO i = 1,npt
570 k_mtb = max(1, idxzb(i))
571 DO k = k_mtb,kbps !KBPS = MAX(kref) ;KMPS= MIN(kref)
572 IF (k < kref(i)) THEN
573 j = ipt(i)
574 rdelks = del(j,k) * delks(i)
575 ubar(i) = ubar(i) + rdelks * u1(j,k) ! Mean U below kref
576 vbar(i) = vbar(i) + rdelks * v1(j,k) ! Mean V below kref
577 roll(i) = roll(i) + rdelks * ro(i,k) ! Mean RO below kref
578 bnv2bar(i) = bnv2bar(i) + .5*(bnv2(i,k)+bnv2(i,k+1))* rdelks
579 ENDIF
580 ENDDO
581 ENDDO
582!
583! orographic asymmetry parameter (OA), and (CLX)
584 DO i = 1,npt
585 j = ipt(i)
586 wdir = atan2(ubar(i),vbar(i)) + pi
587 idir = mod(nint(fdir*wdir),mdir) + 1
588 nwd = nwdir(idir)
589 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(j,mod(nwd-1,4)+1)
590 clx(i) = clx4(j,mod(nwd-1,4)+1)
591 ENDDO
592!
593 DO i = 1,npt
594 dtfac(i) = 1.0
595 icrilv(i) = .false. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR
596 ulow(i) = max(sqrt(ubar(i)*ubar(i)+vbar(i)*vbar(i)),velmin)
597 xn(i) = ubar(i) / ulow(i)
598 yn(i) = vbar(i) / ulow(i)
599 ENDDO
600!
601 DO k = 1, kmm1
602 DO i = 1,npt
603 j = ipt(i)
604 velco(i,k) = 0.5 * ((u1(j,k)+u1(j,k+1))*xn(i)
605 & + (v1(j,k)+v1(j,k+1))*yn(i))
606 ENDDO
607 ENDDO
608!
609!------------------
610! v0: incorporates latest modifications for kxridge and heff/hsat
611! and taulin for Fr <=fcrit_gfs
612! and concept of "clipped" hill if zmtb > 0. to make
613! the integrated "tau_sso = tau_ogw +tau_mtb" close to reanalysis data
614! it is still used the "single-OGWave"-approach along ULOW-upwind
615!
616! in contrast to the 2-orthogonal wave (2OTW) schemes of IFS/METO/E-CANADA
617! 2OTW scheme requires "aver angle" and wind projections on 2 axes of ellipse a-b
618! with 2-stresses: taub_a & taub_b from AS of Phillips et al. (1984)
619!------------------
620 taub(:) = 0. ; taulin(:)= 0.
621 DO i = 1,npt
622 j = ipt(i)
623 bnv = sqrt( bnv2bar(i) )
624 heff = min(hprime(j),hpmax)
625
626 if( zmtb(j) > 0.) heff = max(sigfac*heff-zmtb(j), 0.)/sigfac
627 if (heff <= 0) cycle
628
629 hsat = fcrit_gfs*ulow(i)/bnv
630 heff = min(heff, hsat)
631
632 fr = min(bnv * heff /ulow(i), frmax)
633!
634 efact = (oa(i) + 2.) ** (ceofrc*fr)
635 efact = min( max(efact,efmin), efmax )
636!
637 coefm = (1. + clx(i)) ** (oa(i)+1.)
638!
639 xlinv(i) = coefm * cleff ! effective kxw for Lin-wave
640 xlingfs = coefm * cleff
641!
642 tem = fr * fr * oc(j)
643 gfobnv = gmax * tem / ((tem + cg)*bnv)
644!
645!new specification of XLINV(I) & taulin(i)
646
647 sigres = max(sigmin, sigma(j))
648 if (heff/sigres > hdxres) sigres = heff/hdxres
649 inv_b2eff = 0.5*sigres/heff
650 kxridge = 1.0 / sqrt(sparea(j))
651 xlinv(i) = xlingfs !or max(kxridge, inv_b2eff) ! 6.28/Lx ..0.5*sigma(j)/heff = 1./Lridge
652 taulin(i) = 0.5*roll(i)*xlinv(i)*bnv*ulow(i)*
653 & heff*heff
654
655 if ( fr > fcrit_gfs ) then
656 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
657 & * ulow(i) * gfobnv * efact ! nonlinear FLUX Tau0...XLINV(I)
658!
659 else
660!
661 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i)
662 & * ulow(i) * gfobnv * efact
663!
664! TAUB(I) = taulin(i) ! linear flux for FR <= fcrit_gfs
665!
666 endif
667!
668!
669 k = max(1, kref(i)-1)
670 tem = max(velco(i,k)*velco(i,k), dw2min)
671 scor(i) = bnv2(i,k) / tem ! Scorer parameter below ref level
672!
673! diagnostics for zogw > zmtb
674!
675 zogw(j) = phii(j, kref(i)) *rgrav
676 ENDDO
677!
678!----SET UP BOTTOM VALUES OF STRESS
679!
680 DO k = 1, kbps
681 DO i = 1,npt
682 IF (k <= kref(i)) taup(i,k) = taub(i)
683 ENDDO
684 ENDDO
685
686 if (strsolver == 'PSS-1986') then
687
688!======================================================
689! V0-GFS OROGW-solver of Palmer et al 1986 -"PSS-1986"
690! in V1-OROGW LINSATDIS of "WAM-2017"
691! with LLWB-mechanism for
692! rotational/non-hydrostat OGWs important for
693! HighRES-FV3GFS with dx < 10 km
694!======================================================
695
696 DO k = kmps, kmm1 ! Vertical Level Loop
697 kp1 = k + 1
698 DO i = 1, npt
699!
700 IF (k >= kref(i)) THEN
701 icrilv(i) = icrilv(i) .OR. ( ri_n(i,k) < ric)
702 & .OR. (velco(i,k) <= 0.0)
703 ENDIF
704 ENDDO
705!
706 DO i = 1,npt
707 IF (k >= kref(i)) THEN
708 IF (.NOT.icrilv(i) .AND. taup(i,k) > 0.0 ) THEN
709 temv = 1.0 / max(velco(i,k), velmin)
710!
711 IF (oa(i) > 0. .AND. kp1 < kref(i)) THEN
712 scork = bnv2(i,k) * temv * temv
713 rscor = min(1.0, scork / scor(i))
714 scor(i) = scork
715 ELSE
716 rscor = 1.
717 ENDIF
718!
719 brvf = sqrt(bnv2(i,k)) ! Brent-Vaisala Frequency interface
720! TEM1 = XLINV(I)*(RO(I,KP1)+RO(I,K))*BRVF*VELCO(I,K)*0.5
721
722 tem1 = xlinv(i)*(ro(i,kp1)+ro(i,k))*brvf*0.5
723 & * max(velco(i,k), velmin)
724 hd = sqrt(taup(i,k) / tem1)
725 fro = brvf * hd * temv
726!
727! RIM is the "WAVE"-RICHARDSON NUMBER BY PALMER,Shutts, Swinbank 1986
728!
729
730 tem2 = sqrt(ri_n(i,k))
731 tem = 1. + tem2 * fro
732 ri_gw = ri_n(i,k) * (1.0-fro) / (tem * tem)
733!
734! CHECK STABILITY TO EMPLOY THE 'dynamical SATURATION HYPOTHESIS'
735! OF PALMER,Shutts, Swinbank 1986
736! ----------------------
737 IF (ri_gw <= ric .AND.
738 & (oa(i) <= 0. .OR. kp1 >= kref(i) )) THEN
739 temc = 2.0 + 1.0 / tem2
740 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf
741 taup(i,kp1) = tem1 * hd * hd
742 ELSE
743 taup(i,kp1) = taup(i,k) * rscor
744 ENDIF
745 taup(i,kp1) = min(taup(i,kp1), taup(i,k))
746 ENDIF
747 ENDIF
748 ENDDO
749 ENDDO
750!
751! zero momentum deposition at the top model layer
752!
753 taup(1:npt,km+1) = taup(1:npt,km)
754!
755! Calculate wave acc-n: - (grav)*d(tau)/d(p) = taud
756!
757 DO k = 1,km
758 DO i = 1,npt
759 taud(i,k) = grav*(taup(i,k+1) - taup(i,k))/del(ipt(i),k)
760 ENDDO
761 ENDDO
762!
763!------scale MOMENTUM DEPOSITION AT TOP TO 1/2 VALUE
764! it is zero now
765! DO I = 1,npt
766! TAUD(I, km) = TAUD(I,km) * FACTOP
767! ENDDO
768
769!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
770!------IF THE GRAVITY WAVE DRAG WOULD FORCE A CRITICAL LINE IN THE
771!------LAYERS BELOW SIGMA=RLOLEV DURING THE NEXT DELTIM TIMESTEP,
772!------THEN ONLY APPLY DRAG UNTIL THAT CRITICAL LINE IS REACHED.
773! Empirical implementation of the LLWB-mechanism: Lower Level Wave Breaking
774! by limiting "Ax = Dtfac*Ax" due to possible LLWB around Kref and 500 mb
775! critical line [V - Ax*dtp = 0.] is smt like "LLWB" for stationary OGWs
776!2019: this option limits sensitivity of taux/tauy to increase/decreaseof TAUB
777!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
778 DO k = 1,kmm1
779 DO i = 1,npt
780 IF (k >= kref(i) .and. prsi(ipt(i),k) >= rlolev) THEN
781
782 IF(taud(i,k) /= 0.) THEN
783 tem = dtp * taud(i,k)
784 dtfac(i) = min(dtfac(i),abs(velco(i,k)/tem))
785! DTFAC(I) = 1.0
786 ENDIF
787 ENDIF
788 ENDDO
789 ENDDO
790!
791!--------------------------- OROGW-solver of GFS PSS-1986
792!
793 else
794!
795!--------------------------- OROGW-solver of WAM2017
796!
797! sigres = max(sigmin, sigma(J))
798! if (heff/sigres.gt.dxres) sigres=heff/dxres
799! inv_b2eff = 0.5*sigres/heff
800! XLINV(I) = max(kxridge, inv_b2eff) ! 0.5*sigma(j)/heff = 1./Lridge
801 dtfac(:) = 1.0
802
803 call oro_wam_2017(im, km, npt, ipt, kref, kdt, me, master,
804 & dtp, dxres, taub, u1, v1, t1, xn, yn, bnv2, ro, prsi,prsl,
805 & del, sigma, hprime, gamma, theta,
806 & sinlat, xlatd, taup, taud, pkdis)
807
808 endif ! oro_wam_2017 - LINSATDIS-solver of WAM-2017
809!
810!--------------------------- OROGW-solver of WAM2017
811!
812! TOFD as in BELJAARS-2004
813!
814! ---------------------------
815 IF( do_tofd ) then
816 axtms(:,:) = 0.0 ; aytms(:,:) = 0.0
817
818
819 DO i = 1,npt
820 j = ipt(i)
821 zpbl =rgrav*phil( j, kpbl(j) )
822
823 sigflt = min(sgh30(j), 0.3*hprime(j)) ! cannot exceed 30% of LS-SSO
824
825 zsurf = phii(j,1)*rgrav
826 do k=1,km
827 zpm(k) = phil(j,k)*rgrav
828 up1(k) = u1(j,k)
829 vp1(k) = v1(j,k)
830 enddo
831
832 call ugwpv0_tofd1d(km, sigflt, elvmaxd(j), zsurf, zpbl,
833 & up1, vp1, zpm, utofd1, vtofd1, epstofd1, krf_tofd1)
834
835 do k=1,km
836 axtms(j,k) = utofd1(k)
837 aytms(j,k) = vtofd1(k)
838!
839! add TOFD to GW-tendencies
840!
841 pdvdt(j,k) = pdvdt(j,k) + aytms(j,k)
842 pdudt(j,k) = pdudt(j,k) + axtms(j,k)
843 enddo
844!2018-diag
845 tau_tofd(j) = sum( utofd1(1:km)* del(j,1:km))
846 enddo
847 ENDIF ! do_tofd
848
849!---------------------------
850! combine oro-drag effects
851!---------------------------
852! + diag-3d
853
854 dudt_tms = axtms
855 tau_ogw = 0.
856 tau_mtb = 0.
857
858 DO k = 1,km
859 DO i = 1,npt
860 j = ipt(i)
861!
862 eng0 = 0.5*(u1(j,k)*u1(j,k)+v1(j,k)*v1(j,k))
863!
864 if ( k < idxzb(i) .AND. idxzb(i) /= 0 ) then
865!
866! if blocking layers -- no OGWs
867!
868 dbim = db(i,k) / (1.+db(i,k)*dtp)
869 pdvdt(j,k) = - dbim * v1(j,k) +pdvdt(j,k)
870 pdudt(j,k) = - dbim * u1(j,k) +pdudt(j,k)
871 eng1 = eng0*(1.0-dbim*dtp)*(1.-dbim*dtp)
872
873 dusfc(j) = dusfc(j) - dbim * u1(j,k) * del(j,k)
874 dvsfc(j) = dvsfc(j) - dbim * v1(j,k) * del(j,k)
875!2018-diag
876 dudt_mtb(j,k) = -dbim * u1(j,k)
877 tau_mtb(j) = tau_mtb(j) + dudt_mtb(j,k)* del(j,k)
878
879 else
880!
881! OGW-s above blocking height
882!
883 taud(i,k) = taud(i,k) * dtfac(i)
884 dtaux = taud(i,k) * xn(i)
885 dtauy = taud(i,k) * yn(i)
886
887 pdvdt(j,k) = dtauy +pdvdt(j,k)
888 pdudt(j,k) = dtaux +pdudt(j,k)
889
890 unew = u1(j,k) + dtaux*dtp ! Pdudt(J,K)*DTP
891 vnew = v1(j,k) + dtauy*dtp ! Pdvdt(J,K)*DTP
892 eng1 = 0.5*(unew*unew + vnew*vnew)
893!
894 dusfc(j) = dusfc(j) + dtaux * del(j,k)
895 dvsfc(j) = dvsfc(j) + dtauy * del(j,k)
896!2018-diag
897 dudt_ogw(j,k) = dtaux
898 tau_ogw(j) = tau_ogw(j) +dtaux*del(j,k)
899 endif
900!
901! local energy deposition SSO-heat
902!
903 pdtdt(j,k) = max(eng0-eng1,0.)*rcpdt
904 ENDDO
905 ENDDO
906! dusfc w/o tofd sign as in the ERA-I, MERRA and CFSR
907 DO i = 1,npt
908 j = ipt(i)
909 dusfc(j) = -rgrav * dusfc(j)
910 dvsfc(j) = -rgrav * dvsfc(j)
911 tau_mtb(j) = -rgrav * tau_mtb(j)
912 tau_ogw(j) = -rgrav * tau_ogw(j)
913 tau_tofd(j) = -rgrav * tau_tofd(j)
914 ENDDO
915
916 RETURN
917
918 end subroutine gwdps_v0
919
920
921
922!===============================================================================
923!23456==============================================================================
924
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