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