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