CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
cires_ugwpv1_solv2.F90
1
3
5
6
7contains
8
9
10!---------------------------------------------------
11! Broad spectrum FVS-1993, mkz^nSlope with nSlope = 0, 1,2
12! dissipative solver with NonHyd/ROT-effects
13! reflected GWs treated as waves with "negligible" flux,
14! they are out of given column
15!---------------------------------------------------
16
17 subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, &
18 tau_ngw, tm , um, vm, qm, prsl, prsi, zmet, zmeti, prslk, &
19 xlatd, sinlat, coslat, &
20 pdudt, pdvdt, pdtdt, dked, zngw)
21!
22!--------------------------------------------------------------------------------
23! nov 2015 alternative gw-solver for nggps-wam
24! nov 2017 nh/rotational gw-modes for nh-fv3gfs
25! oct 2019 adding empirical satellite-based
26! source function and *F90 CIRES-style of the code
27! oct 2020 Diagnostics of "tauabs, wrms, trms" is taken out
28! --------------------------------------------------------------------------------
29!
30 use machine, only : kind_phys
31
32 use cires_ugwpv1_module,only : krad, kvg, kion, ktg, ipr_ktgw, pr_kdis, pr_kvkt
33
34 use cires_ugwpv1_module,only : knob_ugwp_doheat, knob_ugwp_dokdis, idebug_gwrms
35
36 use cires_ugwpv1_module,only : psrc => knob_ugwp_palaunch
37
38 use cires_ugwpv1_module,only : maxdudt, maxdtdt, max_eps, dked_min, dked_max
39
40 use ugwp_common , only : rgrav, grav, cpd, rd, rv, rcpdl, grav2cpd, &
41 omega2, rcpd, rcpd2, pi, pi2, fv, &
42 rad_to_deg, deg_to_rad, &
43 rdi, gor, grcp, gocp, &
44 bnv2min, bnv2max, dw2min, velmin, gr2, &
45 hpscale, rhp, rh4, grav2, rgrav2, mkzmin, mkz2min
46!
47 use ugwp_wmsdis_init, only : v_kxw, rv_kxw, v_kxw2, tamp_mpa, tau_min, ucrit, &
48 gw_eff, &
49 nslope, ilaunch, zms, &
50 zci, zdci, zci4, zci3, zci2, &
51 zaz_fct, zcosang, zsinang, nwav, nazd, &
52 zcimin, zcimax, rimin, sc2, sc2u, ric
53!
54 implicit none
55!
56 real(kind=kind_phys), parameter :: zsp_gw = 106.5e3 ! sponge for GWs above the model top
57 real(kind=kind_phys), parameter :: linsat2 = 1.0, dturb_max = 100.0
58 integer, parameter :: ener_norm =0
59 integer, parameter :: ener_lsat=0
60 integer, parameter :: nstdif = 1
61 integer, parameter :: wave_sponge = 1
62
63 integer, intent(in) :: levs ! vertical level
64 integer, intent(in) :: im ! horiz tiles
65 integer, intent(in) :: mpi_id, master, kdt
66
67 real(kind=kind_phys) ,intent(in) :: dtp ! model time step
68 real(kind=kind_phys) ,intent(in) :: tau_ngw(im)
69
70 real(kind=kind_phys) ,intent(in) :: vm(im,levs) ! meridional wind
71 real(kind=kind_phys) ,intent(in) :: um(im,levs) ! zonal wind
72 real(kind=kind_phys) ,intent(in) :: qm(im,levs) ! spec. humidity
73 real(kind=kind_phys) ,intent(in) :: tm(im,levs) ! kinetic temperature
74
75 real(kind=kind_phys) ,intent(in) :: prsl(im,levs) ! mid-layer pressure
76 real(kind=kind_phys) ,intent(in) :: prslk(im,levs) ! mid-layer exner function
77 real(kind=kind_phys) ,intent(in) :: zmet(im,levs) ! meters now !!!!! phil =philg/grav
78 real(kind=kind_phys) ,intent(in) :: prsi(im,levs+1) ! interface pressure
79 real(kind=kind_phys) ,intent(in) :: zmeti(im,levs+1) ! interface geopi/meters
80 real(kind=kind_phys) ,intent(in) :: xlatd(im) ! xlat_d in degrees
81 real(kind=kind_phys) ,intent(in) :: sinlat(im)
82 real(kind=kind_phys) ,intent(in) :: coslat(im)
83!
84! out-gw effects
85!
86 real(kind=kind_phys) ,intent(out) :: pdudt(im,levs) ! zonal momentum tendency
87 real(kind=kind_phys) ,intent(out) :: pdvdt(im,levs) ! meridional momentum tendency
88 real(kind=kind_phys) ,intent(out) :: pdtdt(im,levs) ! gw-heating (u*ax+v*ay)/cp and cooling
89 real(kind=kind_phys) ,intent(out) :: dked(im,levs) ! gw-eddy diffusion
90 real(kind=kind_phys) ,intent(out) :: zngw(im) ! launch height
91!
92!
93!
94! local ===========================================================================================
95
96 real(kind=kind_phys) :: tauabs(im,levs) !
97 real(kind=kind_phys) :: wrms(im,levs) !
98 real(kind=kind_phys) :: trms(im,levs) !
99
100 real(kind=kind_phys) :: zwrms(nwav,nazd), wrk1(levs), wrk2(levs)
101 real(kind=kind_phys) :: atrms(nazd, levs),awrms(nazd, levs), akzw(nwav,nazd, levs+1)
102!
103! local ===========================================================================================
104 real(kind=kind_phys) :: taux(levs+1) ! EW component of vertical momentum flux (pa)
105 real(kind=kind_phys) :: tauy(levs+1) ! NS component of vertical momentum flux (pa)
106 real(kind=kind_phys) :: fpu(nazd, levs+1) ! az-momentum flux
107 real(kind=kind_phys) :: ui(nazd, levs+1) ! azimuthal wind
108
109 real(kind=kind_phys) :: fden_bn(levs+1) ! density/brent
110 real(kind=kind_phys) :: flux(nwav, nazd) , flux_m(nwav, nazd)
111!
112 real(kind=kind_phys) :: bn(levs+1) ! interface BV-frequency
113 real(kind=kind_phys) :: bn2(levs+1) ! interface BV*BV-frequency
114 real(kind=kind_phys) :: rhoint(levs+1) ! interface density
115 real(kind=kind_phys) :: uint(levs+1) ! interface zonal wind
116 real(kind=kind_phys) :: vint(levs+1) ! meridional wind
117 real(kind=kind_phys) :: tint(levs+1) ! temp-re
118
119 real(kind=kind_phys) :: irhodz_mid(levs)
120 real(kind=kind_phys) :: suprf(levs+1) ! RF-super linear dissipation
121 real(kind=kind_phys) :: cstar(levs+1) ,cstar2(levs+1)
122 real(kind=kind_phys) :: v_zmet(levs+1)
123 real(kind=kind_phys) :: vueff(levs+1)
124 real(kind=kind_phys) :: dfdz_v(nazd, levs), dfdz_heat(nazd, levs) ! axj = -df*rho/dz directional Ax
125
126 real(kind=kind_phys), dimension(levs) :: atm , aum, avm, aqm, aprsl, azmet, dz_met
127 real(kind=kind_phys), dimension(levs+1) :: aprsi, azmeti, dz_meti
128
129 real(kind=kind_phys), dimension(levs) :: wrk3
130 real(kind=kind_phys), dimension(levs) :: uold, vold, told, unew, vnew, tnew
131 real(kind=kind_phys), dimension(levs) :: rho, rhomid, adif, cdif, acdif
132 real(kind=kind_phys), dimension(levs) :: qmid, akt
133 real(kind=kind_phys), dimension(levs+1) :: dktur, ktint, kvint
134 real(kind=kind_phys), dimension(levs+1) :: fden_lsat, fden_bnen
135
136 integer, dimension(levs) :: Anstab
137
138 real(kind=kind_phys) :: sig_u2az(nazd), sig_u2az_m(nazd)
139 real(kind=kind_phys) :: wave_dis(nwav, nazd), wave_disaz(nazd)
140 real(kind=kind_phys) :: rdci(nwav), rci(nwav)
141 real(kind=kind_phys) :: wave_act(nwav, nazd) ! active waves at given vert-level
142 real(kind=kind_phys) :: ul(nazd) ! velocity in azimuthal direction at launch level
143!
144! scalars
145!
146 real(kind=kind_phys) :: bvi, bvi2, bvi3, bvi4, rcms ! BV at launch level
147 real(kind=kind_phys) :: c2f2, cf1, wave_distot
148
149
150 real(kind=kind_phys) :: flux_norm ! norm-factor
151 real(kind=kind_phys) :: taub_src, rho_src, zcool, vmdiff
152!
153 real(kind=kind_phys) :: zthm, dtau, cgz, ucrit_maxdc
154 real(kind=kind_phys) :: vm_zflx_mode, vc_zflx_mode
155 real(kind=kind_phys) :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2,v_cdp2
156 real(kind=kind_phys) :: ucrit_max
157 real(kind=kind_phys) :: pwrms, ptrms
158 real(kind=kind_phys) :: zu, zcin, zcin2, zcin3, zcin4, zcinc
159 real(kind=kind_phys) :: zatmp, fluxs, zdep, ze1, ze2
160
161!
162 real(kind=kind_phys) :: zdelp, zdelm, taud_min
163 real(kind=kind_phys) :: tvc, tvm, ptc, ptm
164 real(kind=kind_phys) :: umfp, umfm, umfc, ucrit3
165 real(kind=kind_phys) :: fmode, expdis, fdis
166 real(kind=kind_phys) :: v_kzi, v_kzw, v_cdp, v_wdp, tx1, fcorsat, dzcrit
167 real(kind=kind_phys) :: v_wdi, v_wdpc
168 real(kind=kind_phys) :: ugw, vgw, ek1, ek2, rdtp, rdtp2, rhp_wam
169
170 integer :: j, jj, k, kk, inc, jk, jkp, jl, iaz
171 integer :: ksrc, km2, km1, kp1, ktop
172!
173! Kturb-part
174!
175 real(kind=kind_phys) :: uz, vz, shr2 , ritur, ktur
176
177 real(kind=kind_phys) :: kamp, zmetk, zgrow
178 real(kind=kind_phys) :: stab, stab_dt, dtstab
179 real(kind=kind_phys) :: nslope3
180!
181 integer :: nstab, ist
182 real(kind=kind_phys) :: w1, w2, w3, dtdif
183
184 real(kind=kind_phys) :: dzmetm, dzmetp, dzmetf, bdif, bt_dif, apc, kturp
185 real(kind=kind_phys) :: rstar, rstar2
186
187 real(kind=kind_phys) :: snorm_ener, sigu2, flux_2_sig, ekin_norm
188 real(kind=kind_phys) :: taub_ch, sigu2_ch
189 real(kind=kind_phys) :: pr_kdis_eff, mf_diss_heat, ipr_max
190 real(kind=kind_phys) :: exp_sponge, mi_sponge, gipr
191
192!--------------------------------------------------------------------------
193!
194 nslope3 = nslope + 3.0
195 pr_kdis_eff = gw_eff*pr_kdis
196 ipr_max = max(1.0, ipr_ktgw)
197 gipr = grav* ipr_ktgw
198!
199! test for input fields
200! if (mpi_id == master .and. kdt < -2) then
201! print *, im, levs, dtp, kdt, ' vay-solv2-v1'
202! print *, minval(tm), maxval(tm), ' min-max-tm '
203! print *, minval(vm), maxval(vm), ' min-max-vm '
204! print *, minval(um), maxval(um), ' min-max-um '
205! print *, minval(qm), maxval(qm), ' min-max-qm '
206! print *, minval(prsl), maxval(prsl), ' min-max-Pmid '
207! print *, minval(prsi), maxval(prsi), ' min-max-Pint '
208! print *, minval(zmet), maxval(zmet), ' min-max-Zmid '
209! print *, minval(zmeti), maxval(zmeti), ' min-max-Zint '
210! print *, minval(prslk), maxval(prslk), ' min-max-Exner '
211! print *, minval(tau_ngw), maxval(tau_ngw), ' min-max-taungw '
212! print *, tau_min, ' tau_min ', tamp_mpa, ' tamp_mpa '
213!
214! endif
215
216 if (idebug_gwrms == 1) then
217 tauabs=0.0; wrms =0.0 ; trms =0.0
218 endif
219
220 rci(:) = 1./zci(:)
221 rdci(:) = 1./zdci(:)
222
223 rdtp = 1./dtp
224 rdtp2 = 0.5*rdtp
225
226 ksrc= max(ilaunch, 3)
227 km2 = ksrc - 2
228 km1 = ksrc - 1
229 kp1 = ksrc + 1
230 ktop= levs+1
231
232 suprf(ktop) = kion(levs)
233
234 do k=1,levs
235 suprf(k) = kion(k) ! approximate 1-st order damping with Fast super-RF of FV3
236 pdvdt(:,k) = 0.0
237 pdudt(:,k) = 0.0
238 pdtdt(:,k) = 0.0
239 dked(: ,k) = 0.0
240 enddo
241
242!-----------------------------------------------------------
243! column-based j=1,im pjysics with 1D-arrays
244!-----------------------------------------------------------
245 DO j=1, im
246 jl =j
247 tx1 = omega2 * sinlat(j) *rv_kxw
248 cf1 = abs(tx1)
249 c2f2 = tx1 * tx1
250 ucrit_max = max(ucrit, cf1)
251 ucrit3 = ucrit_max*ucrit_max*ucrit_max
252!
253! ngw-fluxes at all gridpoints (with tau_min at least)
254!
255 aprsl(1:levs) = prsl(jl,1:levs)
256!
257! ksrc-define "aprsi(1:levs+1) redefine "ilaunch"
258!
259 do k=1, levs
260 if (aprsl(k) .lt. psrc ) exit
261 enddo
262 ilaunch = max(k-1, 3)
263 ksrc= max(ilaunch, 3)
264
265 zngw(j) = zmet(j, ksrc)
266
267 km2 = ksrc - 2
268 km1 = ksrc - 1
269 kp1 = ksrc + 1
270
271!=====ksrc
272
273 aum(1:levs) = um(jl,1:levs)
274 avm(1:levs) = vm(jl,1:levs)
275 atm(1:levs) = tm(jl,1:levs)
276 aqm(1:levs) = qm(jl,1:levs)
277 azmet(1:levs) = zmet(jl,1:levs)
278 aprsi(1:levs+1) = prsi(jl,1:levs+1)
279 azmeti(1:levs+1) = zmeti(jl,1:levs+1)
280
281 rho_src = aprsl(ksrc)*rdi/atm(ksrc)
282 taub_ch = max(tau_ngw(jl), tau_min)
283 taub_src = taub_ch
284
285
286 sigu2 = taub_src/rho_src/v_kxw * zms
287 sig_u2az(1:nazd) = sigu2
288!
289! compute diffusion-based arrays km2:levs
290!
291 do jk = km2, levs
292 dz_meti(jk) = azmeti(jk+1)-azmeti(jk)
293 dz_met(jk) = azmet(jk)-azmeti(jk-1)
294 enddo
295! ---------------------------------------------
296! interface mean flow parameters launch -> levs+1
297! ---------------------------------------------
298 do jk= km1,levs
299 tvc = atm(jk)*(1. +fv*aqm(jk))
300 tvm = atm(jk-1)*(1. +fv*aqm(jk-1))
301 ptc = tvc/ prslk(jl, jk)
302 ptm = tvm/prslk(jl,jk-1)
303!
304 zthm = 2.0/(tvc+tvm)
305 rhp_wam = zthm*gor
306!interface
307 uint(jk) = 0.5*(aum(jk-1)+aum(jk))
308 vint(jk) = 0.5*(avm(jk-1)+avm(jk))
309 tint(jk) = 0.5*(tvc+tvm)
310 rhomid(jk) = aprsl(jk)*rdi/atm(jk)
311 rhoint(jk) = aprsi(jk)*rdi*zthm ! rho = p/(RTv)
312 zdelp = dz_meti(jk) ! >0 ...... dz-meters
313 v_zmet(jk) = 2.*zdelp ! 2*kzi*[Z_int(k+1)-Z_int(k)]
314 zdelm = 1./dz_met(jk) ! 1/dz ...... 1/meters
315!
316! bvf2 = grav2*zdelm*(ptc-ptm)/(ptc + ptm) ! N2=[g/PT]*(dPT/dz)
317!
318 bn2(jk) = grav2cpd*zthm*(1.0+rcpdl*(tvc-tvm)*zdelm)
319 bn2(jk) = max(min(bn2(jk), bnv2max), bnv2min)
320 bn(jk) = sqrt(bn2(jk))
321
322
323 wrk3(jk)= 1./zdelp/rhomid(jk) ! 1/rho_mid(k)/[Z_int(k+1)-Z_int(k)]
324 irhodz_mid(jk) = rdtp*zdelp*rhomid(jk)/rho_src
325!
326!
327! diagnostics -Kzz above PBL
328!
329 uz = aum(jk) - aum(jk-1)
330 vz = avm(jk) - avm(jk-1)
331 shr2 = (max(uz*uz+vz*vz, dw2min)) * zdelm *zdelm
332
333 zmetk = azmet(jk)* rh4 ! mid-layer height k_int => k_int+1
334 zgrow = exp(zmetk)
335 ritur = bn2(jk)/shr2
336 kamp = sqrt(shr2)*sc2 *zgrow
337 w1 = 1./(1. + 5*ritur)
338 ktur= min(max(kamp * w1 * w1, dked_min), dked_max)
339 zmetk = azmet(jk)* rhp
340 vueff(jk) = ktur + kvg(jk)
341
342 akt(jk) = gipr/tvc
343 enddo
344
345 if (idebug_gwrms == 1) then
346 do jk= km1,levs
347 wrk1(jk) = rv_kxw/rhoint(jk)
348 wrk2(jk)= rgrav2*zthm*zthm*bn2(jk) ! dimension [K*K]*(c2/m2)
349 enddo
350 endif
351
352!
353! extrapolating values for ktop = levs+1 (lev-interface for prsi(levs+1) =/= 0)
354!
355 jk = levs
356
357 rhoint(ktop) = 0.5*aprsi(levs)*rdi/atm(jk)
358 tint(ktop) = atm(jk)*(1. +fv*aqm(jk))
359 uint(ktop) = aum(jk)
360 vint(ktop) = avm(jk)
361
362 v_zmet(ktop) = v_zmet(jk)
363 vueff(ktop) = vueff(jk)
364 bn2(ktop) = bn2(jk)
365 bn(ktop) = bn(jk)
366!
367! akt_mid *KT = -g*(1/H + 1/T*dT/dz)*KT ... grav/tvc for eddy heat conductivity
368!
369 do jk=km1, levs
370 akt(jk) = -akt(jk)*(gor + (tint(jk+1)-tint(jk))/dz_meti(jk) )
371 enddo
372
373
374 bvi = bn(ksrc); bvi2 = bvi * bvi;
375 bvi3 = bvi2*bvi; bvi4 = bvi2 * bvi2; rcms = zms/bvi
376!
377! project winds at ksrc
378!
379 do iaz=1, nazd
380 ul(iaz) = zcosang(iaz) *uint(ksrc) + zsinang(iaz) *vint(ksrc)
381 enddo
382!
383
384 do jk=ksrc, ktop
385 cstar(jk) = bn(jk)/zms
386 cstar2(jk) = cstar(jk)*cstar(jk)
387
388 fden_lsat(jk) = rhoint(jk)/bn(jk)*v_kxw*linsat2
389
390 do iaz=1, nazd
391 zu = zcosang(iaz)*uint(jk) + zsinang(iaz)*vint(jk)
392 ui(iaz, jk) = zu !- ul(iaz)*0.
393 enddo
394 enddo
395
396 rstar = 1./cstar(ksrc)
397 rstar2 = rstar*rstar
398! -----------------------------------------
399! set launch momentum flux spectral density
400! -----------------------------------------
401
402 fpu(1:nazd, km2:ktop) =0.
403
404 do inc=1,nwav
405
406 zcin = zci(inc)*rstar
407
408!
409! integrate (flux(cin) x dcin ) old tau-flux and normalization
410!
411 flux(inc,1) = rstar*(zcin*zcin)/(1.+ zcin**nslope3)
412!
413! fsat = rstar*(zcin*zcin) * taub_src / SN * [rho/rho_src *N_src/N]
414!
415 fpu(1,ksrc) = fpu(1,ksrc) + flux(inc,1)*zdci(inc) ! dc/cstar = dim-less
416
417 do iaz=1,nazd
418 akzw(inc, iaz, ksrc) = bvi*rci(inc)
419 enddo
420
421 enddo
422!
423! adjust rho/bn vertical factors for saturated fluxes (E(m) ~m^-3)
424
425 flux_norm = taub_src / fpu(1, ksrc) ! [Pa * dc/cstar *dim_less]
426 ze1 = flux_norm * bvi/rhoint(ksrc) *rstar *rstar2
427 do jk=ksrc, ktop
428 fden_bn(jk) = ze1* rhoint(jk) / bn(jk) ! [Pa]/[m/s] * rstar2
429 enddo
430!
431 do inc=1, nwav
432 flux(inc,1) = flux_norm*flux(inc,1)
433 enddo
434
435
436 if (ener_norm == 1) then
437 snorm_ener = 0.
438 do inc=1,nwav
439 zcin = zci(inc)*rstar
440
441 ze2 = zcin /(1.+ zcin**nslope3)
442
443 snorm_ener = snorm_ener + ze2*zdci(inc)*rstar !dim-less
444 flux(inc,1) = ze2 * zcin
445 enddo
446
447 ekin_norm = 1./snorm_ener
448
449! taub_src = sigu2 * rho_src * [v_kxw / zms ]
450! sigu2 = taub_src*zms/(rho_src/v_kxw)
451! ze1 = sigu2*ks*dens/Ns = taub*zms/Ns
452
453 ze1 = taub_src*zms/bvi * ekin_norm
454 taub_src = 0.
455
456 do inc=1,nwav
457 flux(inc,1) = ze1* flux(inc,1)
458 taub_src = taub_src + flux(inc,1)*zdci(inc)
459 enddo
460 ze1 = ekin_norm * v_kxw * rstar2
461 do jk=ksrc, ktop
462 fden_bnen(jk) = rhoint(jk) / bn(jk) *ze1 ! mult on => sigu2(z)*cdf2 => flux_sat
463 enddo
464
465 endif
466!
467 do iaz=1,nazd
468 fpu(iaz, ksrc) = taub_src
469 fpu(iaz, km1) = taub_src
470 enddo
471
472! copy flux-1 into other azimuths
473! --------------------------------
474
475
476 do iaz=2, nazd
477 do inc=1,nwav
478 flux(inc,iaz) = flux(inc,1)
479 enddo
480 enddo
481
482! if (mpi_id == master .and. ener_norm == 1) then
483! print *
484! print *, 'vay_norm: ', taub_src, taub_ch, sigu2, flux_norm, ekin_norm
485! print *
486! endif
487
488 if (idebug_gwrms == 1) then
489 pwrms =0.
490 ptrms =0.
491 tx1 = real(nazd)/rhoint(ksrc)*rv_kxw
492 ze2 = wrk2(ksrc) ! (bvi*atm(ksrc)*rgrav)**2
493 do inc=1, nwav
494 v_kzw = bvi*rci(inc)
495 ze1 = flux(inc,1)*zdci(inc)*tx1*v_kzw
496 pwrms = pwrms + ze1
497 ptrms = ptrms + ze1 * ze2
498 enddo
499 wrms(jl, ksrc) = pwrms
500 trms(jl, ksrc) = ptrms
501 endif
502
503! --------------------------------
504 wave_act(:,:) = 1.0
505! vertical do-loop
506 do jk=ksrc, levs
507
508 jkp = jk+1
509! azimuth do-loop
510 do iaz=1, nazd
511
512 sig_u2az_m(iaz) = sig_u2az(iaz)
513
514 umfp = ui(iaz, jkp)
515 umfm = ui(iaz, jk)
516 umfc = .5*(umfm + umfp)
517! wave-cin loop
518 dfdz_v(iaz, jk) = 0.0
519 dfdz_heat(iaz, jk) = 0.0
520 fpu(iaz, jkp) = 0.0
521 sig_u2az(iaz) =0.0
522!
523! wave_dis(iaz, :) = vueff(jk)
524 do inc=1, nwav
525 flux_m(inc, iaz) = flux(inc, iaz)
526
527 zcin = zci(inc) ! zcin =/0 by definition
528 zcinc = rci(inc)
529
530 if(wave_act(inc,iaz) == 1.0) then
531!=======================================================================
532! discrete mode
533! saturated limit wfit = kzw*kzw*kt; wfdt = wfit/(kxw*cx)*betat
534! & dissipative kzi = 2.*kzw*(wfdm+wfdt)*dzpi(k)
535!=======================================================================
536
537 v_cdp = zcin - umfp
538 v_cdp2=v_cdp*v_cdp
539 cdf2 = v_cdp2 - c2f2
540 if (v_cdp .le. ucrit_max .or. cdf2 .le. 0.0) then
541!
542! between layer [k-1,k or jk-jkp] (Chi - Uk) -> ucrit_max, wave's absorption
543!
544 wave_act(inc,iaz) =0.
545 akzw(inc, iaz, jkp) = pi/dz_meti(jk) ! pi2/dzmet
546 fluxs = 0.0 !max(0., rhobnk(jkp)*ucrit3)*rdci(inc)
547 flux(inc,iaz) = fluxs
548
549 else
550
551 v_wdp = v_kxw*v_cdp
552 wdop2 = v_wdp* v_wdp
553
554!
555! rotational cut-off
556!
557 kzw2 = (bn2(jkp)-wdop2)/cdf2
558!
559!cires_ugwp_initialize.F90: real, parameter :: mkzmin = pi2/80.0e3
560!
561 if ( kzw2 > mkz2min ) then
562 v_kzw = sqrt(kzw2)
563 akzw(inc, iaz, jkp) = v_kzw
564!
565!linsatdis: kzw2, kzw3, kdsat, c2f2, cdf2, cdf1
566!
567!kzw2 = (bn2(k)-wdop2)/Cdf2 - rhp4 - v_kx2w ! full lin DS-NGW (N2-wd2)*k2=(m2+k2+[1/2H]^2)*(wd2-f2)
568! Kds_sat = kxw*Cdf1*rhp2/kzw3
569!krad, kvg, kion, ktg
570 v_cdp = sqrt( cdf2 )
571 v_wdp = v_kxw * v_cdp
572 v_wdi = kzw2*vueff(jk) + kion(jk) ! supRF-diss due for "all" vars
573 v_wdpc = sqrt(v_wdp*v_wdp +v_wdi*v_wdi)
574 v_kzi = v_kzw*v_wdi/v_wdpc
575
576!
577 ze1 = v_kzi*v_zmet(jk)
578
579 if (ze1 .ge. 1.e-2) then
580 expdis = max(exp(-ze1), 0.01)
581 else
582 expdis = 1./(1.+ ze1)
583 endif
584
585!
586 wave_act(inc,iaz) = 1.0
587 fmode = flux(inc,iaz)
588
589 flux_2_sig = v_kzw/v_kxw/rhoint(jkp)
590 w1 = v_wdpc/kzw2/v_kzw/v_zmet(jk)
591 else ! kzw2 <= mkz2min large "Lz"-reflection
592
593 expdis = 1.0
594 v_kzw = mkzmin
595
596 v_cdp = 0. ! no effects of reflected waves
597 wave_act(inc,iaz) = 0.0
598 akzw(inc, iaz, jkp) = v_kzw
599 fmode = 0.
600 w1 =0.
601 endif
602! expdis =1.0
603
604 fdis = fmode*expdis*wave_act(inc,iaz)
605!==============================================================================
606!
607! Saturated Fluxes and Energy: Spectral and Dicrete Modes
608!
609! S2003 fluxs= fden_bn(jk)*(zcin-ui(jk,iaz))**2/zcin
610! WM2001 fluxs= fden_bn(jk)*(zcin-ui(jk,iaz))
611! saturated flux + wave dissipation - Keddy_gwsat in UGWP-V1
612! linsatdis = 1.0 , here: u'^2 ~ linsatdis* [v_cdp*v_cdp]
613!
614! old-sat fluxs= fden_bn(jkp)*cdf2*zcinc*wave_act(inc,iaz)
615! fluxs= fden_bn(jkp)*cdf2*zcinc*wave_act(inc,iaz)
616! new sat fluxs= fden_bn(jkp)*sqrt(cdf2)*wave_act(inc,iaz)
617!
618! fluxs= fden_bn(jkp)*sqrt(cdf2)*wave_act(inc,iaz)
619
620!
621!
622! old spectral sat-limit with "mapping to source-level" sp_tau(cd) = fden_bn(jkp)*sqrt(cdf2)
623! new spectral sat-limit with "mapping to source-level" sp_tau(cd) = fden_bn(jkp)*cdf2*rstar2
624! [fden_bn(jkp)] = Pa/dc
625! fsat = rstar*(zcin*zcin) * [taub_src / SN * [ rstar3*rho/rho_src *N_src/N] = fden_bn ]
626
627 if (ener_norm == 0) fluxs= fden_bn(jkp)*cdf2*wave_act(inc,iaz) ! dim-n: Pa/[m/s]
628!
629! single mode saturation limit: [rho(z)/bn(z)*kx *linsat2* cd^3] /dc
630!
631 if (ener_lsat == 1) fluxs= fden_lsat(jkp)*cdf2*sqrt(cdf2)*rdci(inc)*wave_act(inc,iaz)
632
633 if (ener_norm == 1) then
634
635! spectral saturation limit
636
637 if (ener_lsat == 0) fluxs= fden_bnen(jk)*cdf2*wave_act(inc,iaz)*sig_u2az_m(iaz)
638
639! single mode saturation limit: [rho(z)/bn(z)*kx *linsat2* cd^3] /dc
640
641 if (ener_lsat == 1) fluxs= fden_lsat(jkp)*cdf2*sqrt(cdf2)*rdci(inc)*wave_act(inc,iaz)
642!
643 endif
644!----------------------------------------------------------------------------
645! dicrete mode saturation fden_sat(jkp) = rhoint(jkp)/bn(jkp)*v_kxw
646! fluxs = fden_sat(jkp)*cdf2*sqrt(cdf2)/zdci(inc)*L2sat
647! fluxs_src = fden_sat(ksrc)*cdf2*sqrt(cdf2)/zdci(inc)*L2sat
648!----------------------------------------------------------------------------
649 zdep = fdis-fluxs ! dimension [Pa/dc] *dc = Pa
650 if(zdep > 0.0 ) then
651! subs on sat-limit
652 ze1 = flux(inc,iaz)
653 flux(inc,iaz) = fluxs
654 ze2 = log(ze1/fluxs)*w1 ! Kdsat-compute damping of mode =>df = f-fluxs
655 ! here we can add extra-dissip for the next layer
656 else
657! assign dis-ve flux
658 flux(inc,iaz) = fdis
659 endif
660
661 dtau = flux_m(inc,iaz)-flux(inc,iaz)
662 if (dtau .lt. 0) then
663 flux(inc,iaz) = flux_m(inc,iaz)
664 endif
665!
666! GW-sponge domain: saturate all "GW"-modes above "zsp_gw"
667!
668 if ( azmeti(jkp) .ge. zsp_gw) then
669 mi_sponge = .5/dz_meti(jk)
670 ze2 = v_wdp /v_kzw * mi_sponge ! Ksat*v_kzw2 = [mi_sat*wdp/kzw]
671 v_wdi = ze2 + v_wdi*0.25 ! diss-sat GW-sponge
672 v_wdpc = sqrt(v_wdp*v_wdp +v_wdi*v_wdi)
673 v_kzi = v_kzw*v_wdi/v_wdpc
674!
675 ze1 = v_kzi*v_zmet(jk)
676 exp_sponge = exp(-ze1)
677!
678! additional sponge
679!
680 flux(inc,iaz) = flux(inc,iaz) *exp_sponge
681 endif
682
683 endif ! coriolis or CL condition-checkif => (v_cdp .le. ucrit_max) then
684 endif ! only for waves w/o CL-absorption wave_act=1
685!
686! sum for given (jk, iaz) all active "wave" contributions
687!
688 if (wave_act(inc,iaz) == 1) then
689
690 zcinc =zdci(inc)
691 vc_zflx_mode = flux(inc,iaz)
692 vmdiff = max(0., flux_m(inc,iaz)-vc_zflx_mode)
693 if (vmdiff <= 0. ) vc_zflx_mode = flux_m(inc,iaz)
694 ze1 = vc_zflx_mode*zcinc
695 fpu(iaz, jkp) = fpu(iaz,jkp) + ze1 ! flux (pa) at
696 sig_u2az(iaz) = sig_u2az(iaz) + ze1*flux_2_sig ! ekin(m2/s2) at z+dz
697
698!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
699! (heat deposition integration over spectral mode for each azimuth
700! later sum over selected azimuths as "non-negative" scalars)
701! cdf1 = sqrt( (zci(inc)-umfc)**2-c2f2)
702!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
703! zdelp = wrk3(jk)*cdf1 *zcinc
704
705 zdelp = wrk3(jk)* v_cdp *zcinc * vmdiff
706
707
708! zcool = 1. ! COOL=(-3.5 + Pr)/Pr
709! zcool = [Kv/Pr]*N2*(Pr-Cp/R)/cp
710! edis = (c-u)*ax/cp = Kv_dis*N2/cp
711! cool = -Kt*N2/R
712! add heat-conduction "bulk" impact: 1/Pr*(g*g*rho)* d [rho*Kv(dT/dp- R/Cp *T/p)]
713!
714 dfdz_v(iaz, jk) = dfdz_v(iaz,jk) + zdelp ! +cool !heating & simple cooling < 0
715 dfdz_heat(iaz, jk) = dfdz_heat(iaz,jk) + zdelp ! heating -only > 0
716 endif !wave_act(inc,iaz) == 1)
717!
718 enddo ! wave-inc-loop
719
720 ze1 =fpu(iaz, jk)
721 if (fpu(iaz, jkp) > ze1 ) fpu(iaz, jkp) = ze1
722!
723! compute wind and temp-re rms
724!
725 if (idebug_gwrms == 1) then
726 pwrms =0.
727 ptrms =0.
728 do inc=1, nwav
729 if (wave_act(inc,iaz) > 0.) then
730 v_kzw =akzw(inc, iaz, jk)
731 ze1 = flux(inc,iaz)*v_kzw*zdci(inc)*wrk1(jk)
732 pwrms = pwrms + ze1
733 ptrms = ptrms + ze1*wrk2(jk)
734 endif
735 enddo
736 awrms(iaz, jk) = pwrms
737 atrms(iaz, jk) = ptrms
738 endif
739
740! --------------
741 enddo ! end Azimuth do-loop
742
743!
744! eddy wave dissipation to limit GW-rms
745!
746 tx1 = sum(abs(dfdz_heat(1:nazd, jk)))/bn2(jk)
747 ze1=max(dked_min, tx1)
748 ze2=min(dked_max, ze1)
749 vueff(jkp) = ze2 + vueff(jkp)
750!
751 enddo ! end Vertical do-loop
752!
753! top-layers constant interface-fluxes and zero-heat
754! we allow non-zero momentum fluxes and thermal effects
755! fpu(1:nazd,levs+1) = fpu(1:nazd, levs)
756! dfdz_v(1:nazd, levs) = 0.0
757
758! ---------------------------------------------------------------------
759! sum contribution for total zonal and meridional fluxes +
760! energy dissipation
761! ---------------------------------------------------
762!
763!========================================================================
764! at the source level and below taux = 0 (taux_E=-taux_W by assumption)
765!========================================================================
766
767 do jk=ksrc, levs
768 taux(jk) = 0.0
769 tauy(jk) = 0.0
770 do iaz=1,nazd
771 taux(jk) = taux(jk) + fpu(iaz,jk)*zcosang(iaz)
772 tauy(jk) = tauy(jk) + fpu(iaz,jk)*zsinang(iaz)
773 pdtdt(jl,jk) = pdtdt(jl,jk) + dfdz_v(iaz,jk)
774 dked(jl,jk) = dked(jl,jk) + dfdz_heat(iaz,jk)
775 enddo
776 enddo
777 jk = ktop; taux(jk)=0.; tauy(jk)=0.
778 do iaz=1,nazd
779 taux(jk) = taux(jk) + fpu(iaz,jk)*zcosang(iaz)
780 tauy(jk) = tauy(jk) + fpu(iaz,jk)*zsinang(iaz)
781 enddo
782
783 if (idebug_gwrms == 1) then
784 do jk=kp1, levs
785 do iaz=1,nazd
786 wrms(jl,jk) =wrms(jl,jk) + awrms(iaz,jk)
787 trms(jl,jk) =trms(jl,jk) + atrms(iaz,jk)
788 tauabs(jl,jk)=tauabs(jl,jk) + fpu(iaz,jk)
789 enddo
790 enddo
791 endif
792!
793
794 do jk=ksrc+1,levs
795 jkp = jk + 1
796 zdelp = wrk3(jk)*gw_eff
797 ze1 = (taux(jkp)-taux(jk))* zdelp
798 ze2 = (tauy(jkp)-tauy(jk))* zdelp
799
800 if (abs(ze1) >= maxdudt ) then
801 ze1 = sign(maxdudt, ze1)
802 endif
803 if (abs(ze2) >= maxdudt ) then
804 ze2 = sign(maxdudt, ze2)
805 endif
806
807 pdudt(jl,jk) = -ze1
808 pdvdt(jl,jk) = -ze2
809!
810! Cx =0 based Cx=/= 0. above
811!
812!
813 if (knob_ugwp_doheat == 1) then
814!
815!maxdtdt= dked_max * bnfix2
816!
817 pdtdt(jl,jk) = pdtdt(jl,jk)*gw_eff
818 ze2 = pdtdt(jl,jk)
819 if (abs(ze2) >= max_eps ) pdtdt(jl,jk) = sign(max_eps, ze2)
820
821 dked(jl,jk) = dked(jl,jk)/bn2(jk)
822 ze1 = max(dked_min, dked(jl,jk))
823 dked(jl,jk) = min(dked_max, ze1)
824 qmid(jk) = pdtdt(j,jk)
825 endif
826 enddo
827!----------------------------------------------------------------------------------
828! Update heat = ek_diss/cp and aply 1-2-1 smoother for "dked" => dktur
829! here with "u_new = u +dtp*dudt ; vnew = v + v +dtp*dvdt
830! can check "stability" in the column and "add" ktur-estimation
831! to suppress instability as needed so dked = dked_gw + ktur_ric
832!----------------------------------------------------------------------------------
833
834 dktur(1:levs) = dked(jl,1:levs)
835!
836 do ist= 1, nstdif
837 do jk=ksrc,levs-1
838 adif(jk) =.25*(dktur(jk-1)+ dktur(jk+1)) + .5*dktur(jk)
839 enddo
840 dktur(ksrc:levs-1) = adif(ksrc:levs-1)
841 enddo
842 dktur(levs) = .5*( dked(jl,levs)+ dked(jl,levs-1))
843 dktur(levs+1) = dktur(levs)
844
845 do jk=ksrc,levs+1
846 ze1 = .5*( dktur(jk) +dktur(jk-1) )
847 kvint(jk) = ze1
848 ktint(jk) = ze1*ipr_ktgw
849 enddo
850
851!
852! Thermal budget qmid = qheat + qcool
853!
854 do jk=ksrc+1,levs
855 ze2 = qmid(jk) + dktur(jk)*akt(jk) + grav*(ktint(jk+1)-ktint(jk))/dz_meti(jk)
856 qmid(jk) = ze2
857 if (abs(ze2) >= max_eps ) qmid(jk) = sign(max_eps, ze2)
858 pdtdt(jl,jk) = qmid(jk)*rcpd
859 dked(jl, jk) = dktur(jk)
860 enddo
861!
862! perform explicit eddy "diffusive" 3-point smoothing of "u-v-t"
863! from the surface/launch-gw to the "top"
864!
865!
866! update by source function X(t+dt) = X(t) + dtp * dXdt
867!
868 uold(km2:levs) = aum(km2:levs)+pdudt(jl,km2:levs)*dtp
869 vold(km2:levs) = avm(km2:levs)+pdvdt(jl,km2:levs)*dtp
870 told(km2:levs) = atm(km2:levs)+pdtdt(jl,km2:levs)*dtp
871!
872! diagnose turb-profile using "stability-check" relying on the free-atm diffusion
873! sc2 = 30m x 30m
874!
875 dktur(km2:levs) = dked_min
876
877 do jk=km1,levs
878 uz = uold(jk) - uold(jk-1)
879 vz = vold(jk) - vold(jk-1)
880 ze1 = dz_met(jk)
881 zdelm = 1./ze1
882
883 tvc = told(jk) * (1. +fv*aqm(jk))
884 tvm = told(jk-1) * (1. +fv*aqm(jk-1))
885 zthm = 2.0 / (tvc+tvm)
886 shr2 = (max(uz*uz+vz*vz, dw2min)) * zdelm *zdelm
887
888 bn2(jk) = grav2cpd*zthm * (1.0+rcpdl*(tvc-tvm)*zdelm)
889
890 bn2(jk) = max(min(bn2(jk), bnv2max), bnv2min)
891 zmetk = azmet(jk)* rh4 ! mid-layer height k_int => k_int+1
892 zgrow = exp(zmetk)
893 ritur = bn2(jk)/shr2
894 w1 = 1./(1. + 5*ritur)
895 ze2 = min( sc2 *zgrow, 4.*ze1*ze1)
896!
897! Smag-type of eddy diffusion K_smag = Sqrt(Deformation - N2/Pr)* L2 *const
898!
899 kamp = sqrt(shr2)* ze2 * w1 * w1
900 ktur= min(max(kamp, dked_min), dked_max)
901 dktur(jk) = ktur
902!
903! update of dked = dked_gw + k_turb_mf
904!
905 dked(jl, jk) = dked(jl, jk) +ktur
906
907 enddo
908
909!
910! apply eddy effects due to GWs: explicit scheme Kzz*dt/dz2 < 0.5 stability
911!
912 if (knob_ugwp_dokdis == 2) then
913
914 do jk=ksrc,levs
915 ze1 = min(.5*(dktur(jk) +dktur(jk-1)), dturb_max)
916 kvint(jk) = kvint(jk) + ze1
917! ktint(jk) = ktint(jk) + ze1*iPr_ktgw
918 enddo
919 kvint(km1) = kvint(ksrc)
920 kvint(ktop) = kvint(levs)
921
922 dzmetm = 1./dz_met(km1)
923 adif(km1:levs) = 0.
924 cdif(km1:levs) = 0.
925 do jk=km1,levs-1
926
927 dzmetp = 1./dz_met(jk+1)
928 dzmetf = 1./(dz_meti(jk)*rhomid(jk))
929
930
931 ktur = kvint(jk) *rhoint(jk) * dzmetf
932 kturp =kvint(jk+1)*rhoint(jk+1) * dzmetf
933
934 adif(jk) = ktur * dzmetm
935 cdif(jk) = kturp * dzmetp
936 apc = adif(jk)+cdif(jk)
937 acdif(jk) = apc
938
939 w1 = apc*ipr_max
940 if (rdtp < w1 ) then
941 anstab(jk) = floor(w1*dtp) + 1
942 else
943 anstab(jk) = 1
944 endif
945 dzmetm = dzmetp
946 enddo
947
948 nstab = maxval( anstab(ksrc:levs-1))
949
950! if (nstab .ge. 3) print *, 'nstab ', nstab
951!
952! k instead Jk
953!
954 dtdif = dtp/real(nstab)
955 ze1 = 1./dtdif
956
957 do ist= 1, nstab
958 do k=ksrc,levs-1
959 bdif = ze1 - acdif(k)
960 bt_dif = ze1 - acdif(k)* ipr_ktgw ! ipr_Ktgw = 1./Pr <1
961 unew(k) = uold(k)*bdif + uold(k-1)*adif(k) + uold(k+1)*cdif(k)
962 vnew(k) = vold(k)*bdif + vold(k-1)*adif(k) + vold(k+1)*cdif(k)
963 tnew(k) = told(k)*bt_dif+(told(k-1)*adif(k) + told(k+1)*cdif(k))*ipr_ktgw
964 enddo
965
966 uold(ksrc:levs-1) = unew(ksrc:levs-1)*dtdif ! value du/dtp *dtp = du
967 vold(ksrc:levs-1) = vnew(ksrc:levs-1)*dtdif
968 told(ksrc:levs-1) = tnew(ksrc:levs-1)*dtdif
969!
970! smoothing the boundary points: "k-1" = ksrc-1 and "k+1" = levs
971!
972 uold(levs) = uold(levs-1)
973 vold(levs) = vold(levs-1)
974 told(levs) = told(levs-1)
975 enddo
976!
977! compute "smoothed" tendencies by molecular + GW-eddy diffusions
978!
979 do k=ksrc,levs-1
980!
981! final updates of tendencies and diffusion
982!
983 ze2 = rdtp*(uold(k) - aum(k))
984 ze1 = rdtp*(vold(k) - avm(k))
985 pdtdt(jl,k)= rdtp*( told(k) - atm(k) )
986
987 if (abs(pdtdt(jl,k)) >= maxdtdt ) pdtdt(jl,k) = sign(maxdtdt,pdtdt(jl,k) )
988 if (abs(ze1) >= maxdudt ) then
989 ze1 = sign(maxdudt, ze1)
990 endif
991 if (abs(ze2) >= maxdudt ) then
992 ze2 = sign(maxdudt, ze2)
993 endif
994
995 pdudt(jl, k) = ze2
996 pdvdt(jl, k) = ze1
997 uz = uold(k+1) - uold(k-1)
998 vz = vold(k+1) - vold(k-1)
999 ze2 = 1./(dz_met(k+1)+dz_met(k) )
1000 mf_diss_heat = rcpd*kvint(k)*(uz*uz +vz*vz)*ze2*ze2 ! vert grad heat
1001 pdtdt(jl,k)= pdtdt(jl,k) + mf_diss_heat ! extra heat due to eddy viscosity
1002
1003 enddo
1004
1005
1006 ENDIF ! dissipative IF-loop for vertical eddy difusion u-v-t
1007
1008 enddo ! J-loop
1009!
1010 RETURN
1011
1012!================================= diag print after "return" ======================
1013 if (kdt ==1 .and. mpi_id == master) then
1014!
1015 print *, ' ugwpv1: nazd-nw-ilaunch=', nazd, nwav,ilaunch, maxval(kvg), ' kvg '
1016 print *, 'ugwpv1: zdci(inc)=' , maxval(zdci), minval(zdci)
1017 print *, 'ugwpv1: zcimax=' , maxval(zci) ,' zcimin=' , minval(zci)
1018! print *, 'ugwpv1: tau_ngw=' , maxval(taub_src)*1.e3, minval(taub_src)*1.e3, tau_min
1019
1020 print *
1021
1022 endif
1023
1024 if (kdt == 1 .and. mpi_id == master) then
1025 print *, 'vgw done nstab ', nstab
1026!
1027 print *, maxval(pdudt)*86400., minval(pdudt)*86400, 'vgw ax ugwp'
1028 print *, maxval(pdvdt)*86400., minval(pdvdt)*86400, 'vgw ay ugwp'
1029 print *, maxval(dked)*1., minval(dked)*1, 'vgw keddy m2/sec ugwp'
1030 print *, maxval(pdtdt)*86400., minval(pdtdt)*86400,'vgw eps ugwp'
1031!
1032! print *, ' ugwp -heating rates '
1033 endif
1034!=================================
1035 return
1036 end subroutine cires_ugwpv1_ngw_solv2
1037
1038
1039end module cires_ugwpv1_solv2
subroutine atm(parameters, ep_2, epsm1, sfcprs, sfctmp, q2, prcpconv, prcpnonc, prcpshcv, prcpsnow, prcpgrpl, prcphail, soldn, cosz, thair, qair, eair, rhoair, qprecc, qprecl, solad, solai, swdown, bdfall, rain, snow, fp, fpice, prcp)
re-precess atmospheric forcing.