!>\file cires_ugwpv1_solv2.F90 !! module cires_ugwpv1_solv2 contains !--------------------------------------------------- ! Broad spectrum FVS-1993, mkz^nSlope with nSlope = 0, 1,2 ! dissipative solver with NonHyd/ROT-effects ! reflected GWs treated as waves with "negligible" flux, ! they are out of given column !--------------------------------------------------- subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, & tau_ngw, tm , um, vm, qm, prsl, prsi, zmet, zmeti, prslk, & xlatd, sinlat, coslat, & pdudt, pdvdt, pdtdt, dked, zngw) ! !-------------------------------------------------------------------------------- ! nov 2015 alternative gw-solver for nggps-wam ! nov 2017 nh/rotational gw-modes for nh-fv3gfs ! oct 2019 adding empirical satellite-based ! source function and *F90 CIRES-style of the code ! oct 2020 Diagnostics of "tauabs, wrms, trms" is taken out ! -------------------------------------------------------------------------------- ! use machine, only : kind_phys use cires_ugwpv1_module,only : krad, kvg, kion, ktg, iPr_ktgw, Pr_kdis, Pr_kvkt use cires_ugwpv1_module,only : knob_ugwp_doheat, knob_ugwp_dokdis, idebug_gwrms use cires_ugwpv1_module,only : psrc => knob_ugwp_palaunch use cires_ugwpv1_module,only : maxdudt, maxdtdt, max_eps, dked_min, dked_max use ugwp_common , only : rgrav, grav, cpd, rd, rv, rcpdl, grav2cpd, & omega2, rcpd, rcpd2, pi, pi2, fv, & rad_to_deg, deg_to_rad, & rdi, gor, grcp, gocp, & bnv2min, bnv2max, dw2min, velmin, gr2, & hpscale, rhp, rh4, grav2, rgrav2, mkzmin, mkz2min ! use ugwp_wmsdis_init, only : v_kxw, rv_kxw, v_kxw2, tamp_mpa, tau_min, ucrit, & gw_eff, & nslope, ilaunch, zms, & zci, zdci, zci4, zci3, zci2, & zaz_fct, zcosang, zsinang, nwav, nazd, & zcimin, zcimax, rimin, sc2, sc2u, ric ! implicit none ! real(kind=kind_phys), parameter :: zsp_gw = 106.5e3 ! sponge for GWs above the model top real(kind=kind_phys), parameter :: linsat2 = 1.0, dturb_max = 100.0 integer, parameter :: ener_norm =0 integer, parameter :: ener_lsat=0 integer, parameter :: nstdif = 1 integer, parameter :: wave_sponge = 1 integer, intent(in) :: levs ! vertical level integer, intent(in) :: im ! horiz tiles integer, intent(in) :: mpi_id, master, kdt real(kind=kind_phys) ,intent(in) :: dtp ! model time step real(kind=kind_phys) ,intent(in) :: tau_ngw(im) real(kind=kind_phys) ,intent(in) :: vm(im,levs) ! meridional wind real(kind=kind_phys) ,intent(in) :: um(im,levs) ! zonal wind real(kind=kind_phys) ,intent(in) :: qm(im,levs) ! spec. humidity real(kind=kind_phys) ,intent(in) :: tm(im,levs) ! kinetic temperature real(kind=kind_phys) ,intent(in) :: prsl(im,levs) ! mid-layer pressure real(kind=kind_phys) ,intent(in) :: prslk(im,levs) ! mid-layer exner function real(kind=kind_phys) ,intent(in) :: zmet(im,levs) ! meters now !!!!! phil =philg/grav real(kind=kind_phys) ,intent(in) :: prsi(im,levs+1) ! interface pressure real(kind=kind_phys) ,intent(in) :: zmeti(im,levs+1) ! interface geopi/meters real(kind=kind_phys) ,intent(in) :: xlatd(im) ! xlat_d in degrees real(kind=kind_phys) ,intent(in) :: sinlat(im) real(kind=kind_phys) ,intent(in) :: coslat(im) ! ! out-gw effects ! real(kind=kind_phys) ,intent(out) :: pdudt(im,levs) ! zonal momentum tendency real(kind=kind_phys) ,intent(out) :: pdvdt(im,levs) ! meridional momentum tendency real(kind=kind_phys) ,intent(out) :: pdtdt(im,levs) ! gw-heating (u*ax+v*ay)/cp and cooling real(kind=kind_phys) ,intent(out) :: dked(im,levs) ! gw-eddy diffusion real(kind=kind_phys) ,intent(out) :: zngw(im) ! launch height ! ! ! ! local =========================================================================================== real(kind=kind_phys) :: tauabs(im,levs) ! real(kind=kind_phys) :: wrms(im,levs) ! real(kind=kind_phys) :: trms(im,levs) ! real(kind=kind_phys) :: zwrms(nwav,nazd), wrk1(levs), wrk2(levs) real(kind=kind_phys) :: atrms(nazd, levs),awrms(nazd, levs), akzw(nwav,nazd, levs+1) ! ! local =========================================================================================== real(kind=kind_phys) :: taux(levs+1) ! EW component of vertical momentum flux (pa) real(kind=kind_phys) :: tauy(levs+1) ! NS component of vertical momentum flux (pa) real(kind=kind_phys) :: fpu(nazd, levs+1) ! az-momentum flux real(kind=kind_phys) :: ui(nazd, levs+1) ! azimuthal wind real(kind=kind_phys) :: fden_bn(levs+1) ! density/brent real(kind=kind_phys) :: flux (nwav, nazd) , flux_m (nwav, nazd) ! real(kind=kind_phys) :: bn(levs+1) ! interface BV-frequency real(kind=kind_phys) :: bn2(levs+1) ! interface BV*BV-frequency real(kind=kind_phys) :: rhoint(levs+1) ! interface density real(kind=kind_phys) :: uint(levs+1) ! interface zonal wind real(kind=kind_phys) :: vint(levs+1) ! meridional wind real(kind=kind_phys) :: tint(levs+1) ! temp-re real(kind=kind_phys) :: irhodz_mid(levs) real(kind=kind_phys) :: suprf(levs+1) ! RF-super linear dissipation real(kind=kind_phys) :: cstar(levs+1) ,cstar2(levs+1) real(kind=kind_phys) :: v_zmet(levs+1) real(kind=kind_phys) :: vueff(levs+1) real(kind=kind_phys) :: dfdz_v(nazd, levs), dfdz_heat(nazd, levs) ! axj = -df*rho/dz directional Ax real(kind=kind_phys), dimension(levs) :: atm , aum, avm, aqm, aprsl, azmet, dz_met real(kind=kind_phys), dimension(levs+1) :: aprsi, azmeti, dz_meti real(kind=kind_phys), dimension(levs) :: wrk3 real(kind=kind_phys), dimension(levs) :: uold, vold, told, unew, vnew, tnew real(kind=kind_phys), dimension(levs) :: rho, rhomid, adif, cdif, acdif real(kind=kind_phys), dimension(levs) :: Qmid, AKT real(kind=kind_phys), dimension(levs+1) :: dktur, Ktint, Kvint real(kind=kind_phys), dimension(levs+1) :: fden_lsat, fden_bnen integer, dimension(levs) :: Anstab real(kind=kind_phys) :: sig_u2az(nazd), sig_u2az_m(nazd) real(kind=kind_phys) :: wave_dis(nwav, nazd), wave_disaz(nazd) real(kind=kind_phys) :: rdci(nwav), rci(nwav) real(kind=kind_phys) :: wave_act(nwav, nazd) ! active waves at given vert-level real(kind=kind_phys) :: ul(nazd) ! velocity in azimuthal direction at launch level ! ! scalars ! real(kind=kind_phys) :: bvi, bvi2, bvi3, bvi4, rcms ! BV at launch level real(kind=kind_phys) :: c2f2, cf1, wave_distot real(kind=kind_phys) :: flux_norm ! norm-factor real(kind=kind_phys) :: taub_src, rho_src, zcool, vmdiff ! real(kind=kind_phys) :: zthm, dtau, cgz, ucrit_maxdc real(kind=kind_phys) :: vm_zflx_mode, vc_zflx_mode real(kind=kind_phys) :: kzw2, kzw3, kdsat, cdf2, cdf1, wdop2,v_cdp2 real(kind=kind_phys) :: ucrit_max real(kind=kind_phys) :: pwrms, ptrms real(kind=kind_phys) :: zu, zcin, zcin2, zcin3, zcin4, zcinc real(kind=kind_phys) :: zatmp, fluxs, zdep, ze1, ze2 ! real(kind=kind_phys) :: zdelp, zdelm, taud_min real(kind=kind_phys) :: tvc, tvm, ptc, ptm real(kind=kind_phys) :: umfp, umfm, umfc, ucrit3 real(kind=kind_phys) :: fmode, expdis, fdis real(kind=kind_phys) :: v_kzi, v_kzw, v_cdp, v_wdp, tx1, fcorsat, dzcrit real(kind=kind_phys) :: v_wdi, v_wdpc real(kind=kind_phys) :: ugw, vgw, ek1, ek2, rdtp, rdtp2, rhp_wam integer :: j, jj, k, kk, inc, jk, jkp, jl, iaz integer :: ksrc, km2, km1, kp1, ktop ! ! Kturb-part ! real(kind=kind_phys) :: uz, vz, shr2 , ritur, ktur real(kind=kind_phys) :: kamp, zmetk, zgrow real(kind=kind_phys) :: stab, stab_dt, dtstab real(kind=kind_phys) :: nslope3 ! integer :: nstab, ist real(kind=kind_phys) :: w1, w2, w3, dtdif real(kind=kind_phys) :: dzmetm, dzmetp, dzmetf, bdif, bt_dif, apc, kturp real(kind=kind_phys) :: rstar, rstar2 real(kind=kind_phys) :: snorm_ener, sigu2, flux_2_sig, ekin_norm real(kind=kind_phys) :: taub_ch, sigu2_ch real(kind=kind_phys) :: Pr_kdis_eff, mf_diss_heat, iPr_max real(kind=kind_phys) :: exp_sponge, mi_sponge, gipr !-------------------------------------------------------------------------- ! nslope3 = nslope + 3.0 Pr_kdis_eff = gw_eff*pr_kdis iPr_max = max(1.0, iPr_ktgw) gipr = grav* Ipr_ktgw ! ! test for input fields ! if (mpi_id == master .and. kdt < -2) then ! print *, im, levs, dtp, kdt, ' vay-solv2-v1' ! print *, minval(tm), maxval(tm), ' min-max-tm ' ! print *, minval(vm), maxval(vm), ' min-max-vm ' ! print *, minval(um), maxval(um), ' min-max-um ' ! print *, minval(qm), maxval(qm), ' min-max-qm ' ! print *, minval(prsl), maxval(prsl), ' min-max-Pmid ' ! print *, minval(prsi), maxval(prsi), ' min-max-Pint ' ! print *, minval(zmet), maxval(zmet), ' min-max-Zmid ' ! print *, minval(zmeti), maxval(zmeti), ' min-max-Zint ' ! print *, minval(prslk), maxval(prslk), ' min-max-Exner ' ! print *, minval(tau_ngw), maxval(tau_ngw), ' min-max-taungw ' ! print *, tau_min, ' tau_min ', tamp_mpa, ' tamp_mpa ' ! ! endif if (idebug_gwrms == 1) then tauabs=0.0; wrms =0.0 ; trms =0.0 endif rci(:) = 1./zci(:) rdci(:) = 1./zdci(:) rdtp = 1./dtp rdtp2 = 0.5*rdtp ksrc= max(ilaunch, 3) km2 = ksrc - 2 km1 = ksrc - 1 kp1 = ksrc + 1 ktop= levs+1 suprf(ktop) = kion(levs) do k=1,levs suprf(k) = kion(k) ! approximate 1-st order damping with Fast super-RF of FV3 pdvdt(:,k) = 0.0 pdudt(:,k) = 0.0 pdtdt(:,k) = 0.0 dked(: ,k) = 0.0 enddo !----------------------------------------------------------- ! column-based j=1,im pjysics with 1D-arrays !----------------------------------------------------------- DO j=1, im jl =j tx1 = omega2 * sinlat(j) *rv_kxw cf1 = abs(tx1) c2f2 = tx1 * tx1 ucrit_max = max(ucrit, cf1) ucrit3 = ucrit_max*ucrit_max*ucrit_max ! ! ngw-fluxes at all gridpoints (with tau_min at least) ! aprsl(1:levs) = prsl(jl,1:levs) ! ! ksrc-define "aprsi(1:levs+1) redefine "ilaunch" ! do k=1, levs if (aprsl(k) .lt. psrc ) exit enddo ilaunch = max(k-1, 3) ksrc= max(ilaunch, 3) zngw(j) = zmet(j, ksrc) km2 = ksrc - 2 km1 = ksrc - 1 kp1 = ksrc + 1 !=====ksrc aum(1:levs) = um(jl,1:levs) avm(1:levs) = vm(jl,1:levs) atm(1:levs) = tm(jl,1:levs) aqm(1:levs) = qm(jl,1:levs) azmet(1:levs) = zmet(jl,1:levs) aprsi(1:levs+1) = prsi(jl,1:levs+1) azmeti(1:levs+1) = zmeti(jl,1:levs+1) rho_src = aprsl(ksrc)*rdi/atm(ksrc) taub_ch = max(tau_ngw(jl), tau_min) taub_src = taub_ch sigu2 = taub_src/rho_src/v_kxw * zms sig_u2az(1:nazd) = sigu2 ! ! compute diffusion-based arrays km2:levs ! do jk = km2, levs dz_meti(jk) = azmeti(jk+1)-azmeti(jk) dz_met(jk) = azmet(jk)-azmeti(jk-1) enddo ! --------------------------------------------- ! interface mean flow parameters launch -> levs+1 ! --------------------------------------------- do jk= km1,levs tvc = atm(jk)*(1. +fv*aqm(jk)) tvm = atm(jk-1)*(1. +fv*aqm(jk-1)) ptc = tvc/ prslk(jl, jk) ptm = tvm/prslk(jl,jk-1) ! zthm = 2.0/(tvc+tvm) rhp_wam = zthm*gor !interface uint(jk) = 0.5*(aum(jk-1)+aum(jk)) vint(jk) = 0.5*(avm(jk-1)+avm(jk)) tint(jk) = 0.5*(tvc+tvm) rhomid(jk) = aprsl(jk)*rdi/atm(jk) rhoint(jk) = aprsi(jk)*rdi*zthm ! rho = p/(RTv) zdelp = dz_meti(jk) ! >0 ...... dz-meters v_zmet(jk) = 2.*zdelp ! 2*kzi*[Z_int(k+1)-Z_int(k)] zdelm = 1./dz_met(jk) ! 1/dz ...... 1/meters ! ! bvf2 = grav2*zdelm*(ptc-ptm)/(ptc + ptm) ! N2=[g/PT]*(dPT/dz) ! bn2(jk) = grav2cpd*zthm*(1.0+rcpdl*(tvc-tvm)*zdelm) bn2(jk) = max(min(bn2(jk), bnv2max), bnv2min) bn(jk) = sqrt(bn2(jk)) wrk3(jk)= 1./zdelp/rhomid(jk) ! 1/rho_mid(k)/[Z_int(k+1)-Z_int(k)] irhodz_mid(jk) = rdtp*zdelp*rhomid(jk)/rho_src ! ! ! diagnostics -Kzz above PBL ! uz = aum(jk) - aum(jk-1) vz = avm(jk) - avm(jk-1) shr2 = (max(uz*uz+vz*vz, dw2min)) * zdelm *zdelm zmetk = azmet(jk)* rh4 ! mid-layer height k_int => k_int+1 zgrow = exp(zmetk) ritur = bn2(jk)/shr2 kamp = sqrt(shr2)*sc2 *zgrow w1 = 1./(1. + 5*ritur) ktur= min(max(kamp * w1 * w1, dked_min), dked_max) zmetk = azmet(jk)* rhp vueff(jk) = ktur + kvg(jk) akt(jk) = gipr/tvc enddo if (idebug_gwrms == 1) then do jk= km1,levs wrk1(jk) = rv_kxw/rhoint(jk) wrk2(jk)= rgrav2*zthm*zthm*bn2(jk) ! dimension [K*K]*(c2/m2) enddo endif ! ! extrapolating values for ktop = levs+1 (lev-interface for prsi(levs+1) =/= 0) ! jk = levs rhoint(ktop) = 0.5*aprsi(levs)*rdi/atm(jk) tint(ktop) = atm(jk)*(1. +fv*aqm(jk)) uint(ktop) = aum(jk) vint(ktop) = avm(jk) v_zmet(ktop) = v_zmet(jk) vueff(ktop) = vueff(jk) bn2(ktop) = bn2(jk) bn(ktop) = bn(jk) ! ! akt_mid *KT = -g*(1/H + 1/T*dT/dz)*KT ... grav/tvc for eddy heat conductivity ! do jk=km1, levs akt(jk) = -akt(jk)*(gor + (tint(jk+1)-tint(jk))/dz_meti(jk) ) enddo bvi = bn(ksrc); bvi2 = bvi * bvi; bvi3 = bvi2*bvi; bvi4 = bvi2 * bvi2; rcms = zms/bvi ! ! project winds at ksrc ! do iaz=1, nazd ul(iaz) = zcosang(iaz) *uint(ksrc) + zsinang(iaz) *vint(ksrc) enddo ! do jk=ksrc, ktop cstar(jk) = bn(jk)/zms cstar2(jk) = cstar(jk)*cstar(jk) fden_lsat(jk) = rhoint(jk)/bn(jk)*v_kxw*Linsat2 do iaz=1, nazd zu = zcosang(iaz)*uint(jk) + zsinang(iaz)*vint(jk) ui(iaz, jk) = zu !- ul(iaz)*0. enddo enddo rstar = 1./cstar(ksrc) rstar2 = rstar*rstar ! ----------------------------------------- ! set launch momentum flux spectral density ! ----------------------------------------- fpu(1:nazd, km2:ktop) =0. do inc=1,nwav zcin = zci(inc)*rstar ! ! integrate (flux(cin) x dcin ) old tau-flux and normalization ! flux(inc,1) = rstar*(zcin*zcin)/(1.+ zcin**nslope3) ! ! fsat = rstar*(zcin*zcin) * taub_src / SN * [rho/rho_src *N_src/N] ! fpu(1,ksrc) = fpu(1,ksrc) + flux(inc,1)*zdci(inc) ! dc/cstar = dim-less do iaz=1,nazd akzw(inc, iaz, ksrc) = bvi*rci(inc) enddo enddo ! ! adjust rho/bn vertical factors for saturated fluxes (E(m) ~m^-3) flux_norm = taub_src / fpu(1, ksrc) ! [Pa * dc/cstar *dim_less] ze1 = flux_norm * bvi/rhoint(ksrc) *rstar *rstar2 do jk=ksrc, ktop fden_bn(jk) = ze1* rhoint(jk) / bn(jk) ! [Pa]/[m/s] * rstar2 enddo ! do inc=1, nwav flux(inc,1) = flux_norm*flux(inc,1) enddo if (ener_norm == 1) then snorm_ener = 0. do inc=1,nwav zcin = zci(inc)*rstar ze2 = zcin /(1.+ zcin**nslope3) snorm_ener = snorm_ener + ze2*zdci(inc)*rstar !dim-less flux(inc,1) = ze2 * zcin enddo ekin_norm = 1./snorm_ener ! taub_src = sigu2 * rho_src * [v_kxw / zms ] ! sigu2 = taub_src*zms/(rho_src/v_kxw) ! ze1 = sigu2*ks*dens/Ns = taub*zms/Ns ze1 = taub_src*zms/bvi * ekin_norm taub_src = 0. do inc=1,nwav flux(inc,1) = ze1* flux(inc,1) taub_src = taub_src + flux(inc,1)*zdci(inc) enddo ze1 = ekin_norm * v_kxw * rstar2 do jk=ksrc, ktop fden_bnen(jk) = rhoint(jk) / bn(jk) *ze1 ! mult on => sigu2(z)*cdf2 => flux_sat enddo endif ! do iaz=1,nazd fpu(iaz, ksrc) = taub_src fpu(iaz, km1) = taub_src enddo ! copy flux-1 into other azimuths ! -------------------------------- do iaz=2, nazd do inc=1,nwav flux(inc,iaz) = flux(inc,1) enddo enddo ! if (mpi_id == master .and. ener_norm == 1) then ! print * ! print *, 'vay_norm: ', taub_src, taub_ch, sigu2, flux_norm, ekin_norm ! print * ! endif if (idebug_gwrms == 1) then pwrms =0. ptrms =0. tx1 = real(nazd)/rhoint(ksrc)*rv_kxw ze2 = wrk2(ksrc) ! (bvi*atm(ksrc)*rgrav)**2 do inc=1, nwav v_kzw = bvi*rci(inc) ze1 = flux(inc,1)*zdci(inc)*tx1*v_kzw pwrms = pwrms + ze1 ptrms = ptrms + ze1 * ze2 enddo wrms(jl, ksrc) = pwrms trms(jl, ksrc) = ptrms endif ! -------------------------------- wave_act(:,:) = 1.0 ! vertical do-loop do jk=ksrc, levs jkp = jk+1 ! azimuth do-loop do iaz=1, nazd sig_u2az_m(iaz) = sig_u2az(iaz) umfp = ui(iaz, jkp) umfm = ui(iaz, jk) umfc = .5*(umfm + umfp) ! wave-cin loop dfdz_v(iaz, jk) = 0.0 dfdz_heat(iaz, jk) = 0.0 fpu(iaz, jkp) = 0.0 sig_u2az(iaz) =0.0 ! ! wave_dis(iaz, :) = vueff(jk) do inc=1, nwav flux_m(inc, iaz) = flux(inc, iaz) zcin = zci(inc) ! zcin =/0 by definition zcinc = rci(inc) if(wave_act(inc,iaz) == 1.0) then !======================================================================= ! discrete mode ! saturated limit wfit = kzw*kzw*kt; wfdt = wfit/(kxw*cx)*betat ! & dissipative kzi = 2.*kzw*(wfdm+wfdt)*dzpi(k) !======================================================================= v_cdp = zcin - umfp v_cdp2=v_cdp*v_cdp cdf2 = v_cdp2 - c2f2 if (v_cdp .le. ucrit_max .or. cdf2 .le. 0.0) then ! ! between layer [k-1,k or jk-jkp] (Chi - Uk) -> ucrit_max, wave's absorption ! wave_act(inc,iaz) =0. akzw(inc, iaz, jkp) = pi/dz_meti(jk) ! pi2/dzmet fluxs = 0.0 !max(0., rhobnk(jkp)*ucrit3)*rdci(inc) flux(inc,iaz) = fluxs else v_wdp = v_kxw*v_cdp wdop2 = v_wdp* v_wdp ! ! rotational cut-off ! kzw2 = (bn2(jkp)-wdop2)/Cdf2 ! !cires_ugwp_initialize.F90: real, parameter :: mkzmin = pi2/80.0e3 ! if ( kzw2 > mkz2min ) then v_kzw = sqrt(kzw2) akzw(inc, iaz, jkp) = v_kzw ! !linsatdis: kzw2, kzw3, kdsat, c2f2, cdf2, cdf1 ! !kzw2 = (bn2(k)-wdop2)/Cdf2 - rhp4 - v_kx2w ! full lin DS-NGW (N2-wd2)*k2=(m2+k2+[1/2H]^2)*(wd2-f2) ! Kds_sat = kxw*Cdf1*rhp2/kzw3 !krad, kvg, kion, ktg v_cdp = sqrt( cdf2 ) v_wdp = v_kxw * v_cdp v_wdi = kzw2*vueff(jk) + kion(jk) ! supRF-diss due for "all" vars v_wdpc = sqrt(v_wdp*v_wdp +v_wdi*v_wdi) v_kzi = v_kzw*v_wdi/v_wdpc ! ze1 = v_kzi*v_zmet(jk) if (ze1 .ge. 1.e-2) then expdis = max(exp(-ze1), 0.01) else expdis = 1./(1.+ ze1) endif ! wave_act(inc,iaz) = 1.0 fmode = flux(inc,iaz) flux_2_sig = v_kzw/v_kxw/rhoint(jkp) w1 = v_wdpc/kzw2/v_kzw/v_zmet(jk) else ! kzw2 <= mkz2min large "Lz"-reflection expdis = 1.0 v_kzw = mkzmin v_cdp = 0. ! no effects of reflected waves wave_act(inc,iaz) = 0.0 akzw(inc, iaz, jkp) = v_kzw fmode = 0. w1 =0. endif ! expdis =1.0 fdis = fmode*expdis*wave_act(inc,iaz) !============================================================================== ! ! Saturated Fluxes and Energy: Spectral and Dicrete Modes ! ! S2003 fluxs= fden_bn(jk)*(zcin-ui(jk,iaz))**2/zcin ! WM2001 fluxs= fden_bn(jk)*(zcin-ui(jk,iaz)) ! saturated flux + wave dissipation - Keddy_gwsat in UGWP-V1 ! linsatdis = 1.0 , here: u'^2 ~ linsatdis* [v_cdp*v_cdp] ! ! old-sat fluxs= fden_bn(jkp)*cdf2*zcinc*wave_act(inc,iaz) ! fluxs= fden_bn(jkp)*cdf2*zcinc*wave_act(inc,iaz) ! new sat fluxs= fden_bn(jkp)*sqrt(cdf2)*wave_act(inc,iaz) ! ! fluxs= fden_bn(jkp)*sqrt(cdf2)*wave_act(inc,iaz) ! ! ! old spectral sat-limit with "mapping to source-level" sp_tau(cd) = fden_bn(jkp)*sqrt(cdf2) ! new spectral sat-limit with "mapping to source-level" sp_tau(cd) = fden_bn(jkp)*cdf2*rstar2 ! [fden_bn(jkp)] = Pa/dc ! fsat = rstar*(zcin*zcin) * [taub_src / SN * [ rstar3*rho/rho_src *N_src/N] = fden_bn ] if (ener_norm == 0) fluxs= fden_bn(jkp)*cdf2*wave_act(inc,iaz) ! dim-n: Pa/[m/s] ! ! single mode saturation limit: [rho(z)/bn(z)*kx *linsat2* cd^3] /dc ! if (ener_lsat == 1) fluxs= fden_Lsat(jkp)*cdf2*sqrt(cdf2)*rdci(inc)*wave_act(inc,iaz) if (ener_norm == 1) then ! spectral saturation limit if (ener_lsat == 0) fluxs= fden_bnen(jk)*cdf2*wave_act(inc,iaz)*sig_u2az_m(iaz) ! single mode saturation limit: [rho(z)/bn(z)*kx *linsat2* cd^3] /dc if (ener_lsat == 1) fluxs= fden_Lsat(jkp)*cdf2*sqrt(cdf2)*rdci(inc)*wave_act(inc,iaz) ! endif !---------------------------------------------------------------------------- ! dicrete mode saturation fden_sat(jkp) = rhoint(jkp)/bn(jkp)*v_kxw ! fluxs = fden_sat(jkp)*cdf2*sqrt(cdf2)/zdci(inc)*L2sat ! fluxs_src = fden_sat(ksrc)*cdf2*sqrt(cdf2)/zdci(inc)*L2sat !---------------------------------------------------------------------------- zdep = fdis-fluxs ! dimension [Pa/dc] *dc = Pa if(zdep > 0.0 ) then ! subs on sat-limit ze1 = flux(inc,iaz) flux(inc,iaz) = fluxs ze2 = log(ze1/fluxs)*w1 ! Kdsat-compute damping of mode =>df = f-fluxs ! here we can add extra-dissip for the next layer else ! assign dis-ve flux flux(inc,iaz) = fdis endif dtau = flux_m(inc,iaz)-flux(inc,iaz) if (dtau .lt. 0) then flux(inc,iaz) = flux_m(inc,iaz) endif ! ! GW-sponge domain: saturate all "GW"-modes above "zsp_gw" ! if ( azmeti(jkp) .ge. zsp_gw) then mi_sponge = .5/dz_meti(jk) ze2 = v_wdp /v_kzw * mi_sponge ! Ksat*v_kzw2 = [mi_sat*wdp/kzw] v_wdi = ze2 + v_wdi*0.25 ! diss-sat GW-sponge v_wdpc = sqrt(v_wdp*v_wdp +v_wdi*v_wdi) v_kzi = v_kzw*v_wdi/v_wdpc ! ze1 = v_kzi*v_zmet(jk) exp_sponge = exp(-ze1) ! ! additional sponge ! flux(inc,iaz) = flux(inc,iaz) *exp_sponge endif endif ! coriolis or CL condition-checkif => (v_cdp .le. ucrit_max) then endif ! only for waves w/o CL-absorption wave_act=1 ! ! sum for given (jk, iaz) all active "wave" contributions ! if (wave_act(inc,iaz) == 1) then zcinc =zdci(inc) vc_zflx_mode = flux(inc,iaz) vmdiff = max(0., flux_m(inc,iaz)-vc_zflx_mode) if (vmdiff <= 0. ) vc_zflx_mode = flux_m(inc,iaz) ze1 = vc_zflx_mode*zcinc fpu(iaz, jkp) = fpu(iaz,jkp) + ze1 ! flux (pa) at sig_u2az(iaz) = sig_u2az(iaz) + ze1*flux_2_sig ! ekin(m2/s2) at z+dz !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! (heat deposition integration over spectral mode for each azimuth ! later sum over selected azimuths as "non-negative" scalars) ! cdf1 = sqrt( (zci(inc)-umfc)**2-c2f2) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! zdelp = wrk3(jk)*cdf1 *zcinc zdelp = wrk3(jk)* v_cdp *zcinc * vmdiff ! zcool = 1. ! COOL=(-3.5 + Pr)/Pr ! zcool = [Kv/Pr]*N2*(Pr-Cp/R)/cp ! edis = (c-u)*ax/cp = Kv_dis*N2/cp ! cool = -Kt*N2/R ! add heat-conduction "bulk" impact: 1/Pr*(g*g*rho)* d [rho*Kv(dT/dp- R/Cp *T/p)] ! dfdz_v(iaz, jk) = dfdz_v(iaz,jk) + zdelp ! +cool !heating & simple cooling < 0 dfdz_heat(iaz, jk) = dfdz_heat(iaz,jk) + zdelp ! heating -only > 0 endif !wave_act(inc,iaz) == 1) ! enddo ! wave-inc-loop ze1 =fpu(iaz, jk) if (fpu(iaz, jkp) > ze1 ) fpu(iaz, jkp) = ze1 ! ! compute wind and temp-re rms ! if (idebug_gwrms == 1) then pwrms =0. ptrms =0. do inc=1, nwav if (wave_act(inc,iaz) > 0.) then v_kzw =akzw(inc, iaz, jk) ze1 = flux(inc,iaz)*v_kzw*zdci(inc)*wrk1(jk) pwrms = pwrms + ze1 ptrms = ptrms + ze1*wrk2(jk) endif enddo Awrms(iaz, jk) = pwrms Atrms(iaz, jk) = ptrms endif ! -------------- enddo ! end Azimuth do-loop ! ! eddy wave dissipation to limit GW-rms ! tx1 = sum(abs(dfdz_heat(1:nazd, jk)))/bn2(jk) ze1=max(dked_min, tx1) ze2=min(dked_max, ze1) vueff(jkp) = ze2 + vueff(jkp) ! enddo ! end Vertical do-loop ! ! top-layers constant interface-fluxes and zero-heat ! we allow non-zero momentum fluxes and thermal effects ! fpu(1:nazd,levs+1) = fpu(1:nazd, levs) ! dfdz_v(1:nazd, levs) = 0.0 ! --------------------------------------------------------------------- ! sum contribution for total zonal and meridional fluxes + ! energy dissipation ! --------------------------------------------------- ! !======================================================================== ! at the source level and below taux = 0 (taux_E=-taux_W by assumption) !======================================================================== do jk=ksrc, levs taux(jk) = 0.0 tauy(jk) = 0.0 do iaz=1,nazd taux(jk) = taux(jk) + fpu(iaz,jk)*zcosang(iaz) tauy(jk) = tauy(jk) + fpu(iaz,jk)*zsinang(iaz) pdtdt(jl,jk) = pdtdt(jl,jk) + dfdz_v(iaz,jk) dked(jl,jk) = dked(jl,jk) + dfdz_heat(iaz,jk) enddo enddo jk = ktop; taux(jk)=0.; tauy(jk)=0. do iaz=1,nazd taux(jk) = taux(jk) + fpu(iaz,jk)*zcosang(iaz) tauy(jk) = tauy(jk) + fpu(iaz,jk)*zsinang(iaz) enddo if (idebug_gwrms == 1) then do jk=kp1, levs do iaz=1,nazd wrms(jl,jk) =wrms(jl,jk) + Awrms(iaz,jk) trms(jl,jk) =trms(jl,jk) + Atrms(iaz,jk) tauabs(jl,jk)=tauabs(jl,jk) + fpu(iaz,jk) enddo enddo endif ! do jk=ksrc+1,levs jkp = jk + 1 zdelp = wrk3(jk)*gw_eff ze1 = (taux(jkp)-taux(jk))* zdelp ze2 = (tauy(jkp)-tauy(jk))* zdelp if (abs(ze1) >= maxdudt ) then ze1 = sign(maxdudt, ze1) endif if (abs(ze2) >= maxdudt ) then ze2 = sign(maxdudt, ze2) endif pdudt(jl,jk) = -ze1 pdvdt(jl,jk) = -ze2 ! ! Cx =0 based Cx=/= 0. above ! ! if (knob_ugwp_doheat == 1) then ! !maxdtdt= dked_max * bnfix2 ! pdtdt(jl,jk) = pdtdt(jl,jk)*gw_eff ze2 = pdtdt(jl,jk) if (abs(ze2) >= max_eps ) pdtdt(jl,jk) = sign(max_eps, ze2) dked(jl,jk) = dked(jl,jk)/bn2(jk) ze1 = max(dked_min, dked(jl,jk)) dked(jl,jk) = min(dked_max, ze1) qmid(jk) = pdtdt(j,jk) endif enddo !---------------------------------------------------------------------------------- ! Update heat = ek_diss/cp and aply 1-2-1 smoother for "dked" => dktur ! here with "u_new = u +dtp*dudt ; vnew = v + v +dtp*dvdt ! can check "stability" in the column and "add" ktur-estimation ! to suppress instability as needed so dked = dked_gw + ktur_ric !---------------------------------------------------------------------------------- dktur(1:levs) = dked(jl,1:levs) ! do ist= 1, nstdif do jk=ksrc,levs-1 adif(jk) =.25*(dktur(jk-1)+ dktur(jk+1)) + .5*dktur(jk) enddo dktur(ksrc:levs-1) = adif(ksrc:levs-1) enddo dktur(levs) = .5*( dked(jl,levs)+ dked(jl,levs-1)) dktur(levs+1) = dktur(levs) do jk=ksrc,levs+1 ze1 = .5*( dktur(jk) +dktur(jk-1) ) kvint(jk) = ze1 ktint(jk) = ze1*iPr_ktgw enddo ! ! Thermal budget qmid = qheat + qcool ! do jk=ksrc+1,levs ze2 = qmid(jk) + dktur(jk)*Akt(jk) + grav*(ktint(jk+1)-ktint(jk))/dz_meti(jk) qmid(jk) = ze2 if (abs(ze2) >= max_eps ) qmid(jk) = sign(max_eps, ze2) pdtdt(jl,jk) = qmid(jk)*rcpd dked(jl, jk) = dktur(jk) enddo ! ! perform explicit eddy "diffusive" 3-point smoothing of "u-v-t" ! from the surface/launch-gw to the "top" ! ! ! update by source function X(t+dt) = X(t) + dtp * dXdt ! uold(km2:levs) = aum(km2:levs)+pdudt(jl,km2:levs)*dtp vold(km2:levs) = avm(km2:levs)+pdvdt(jl,km2:levs)*dtp told(km2:levs) = atm(km2:levs)+pdtdt(jl,km2:levs)*dtp ! ! diagnose turb-profile using "stability-check" relying on the free-atm diffusion ! sc2 = 30m x 30m ! dktur(km2:levs) = dked_min do jk=km1,levs uz = uold(jk) - uold(jk-1) vz = vold(jk) - vold(jk-1) ze1 = dz_met(jk) zdelm = 1./ze1 tvc = told(jk) * (1. +fv*aqm(jk)) tvm = told(jk-1) * (1. +fv*aqm(jk-1)) zthm = 2.0 / (tvc+tvm) shr2 = (max(uz*uz+vz*vz, dw2min)) * zdelm *zdelm bn2(jk) = grav2cpd*zthm * (1.0+rcpdl*(tvc-tvm)*zdelm) bn2(jk) = max(min(bn2(jk), bnv2max), bnv2min) zmetk = azmet(jk)* rh4 ! mid-layer height k_int => k_int+1 zgrow = exp(zmetk) ritur = bn2(jk)/shr2 w1 = 1./(1. + 5*ritur) ze2 = min( sc2 *zgrow, 4.*ze1*ze1) ! ! Smag-type of eddy diffusion K_smag = Sqrt(Deformation - N2/Pr)* L2 *const ! kamp = sqrt(shr2)* ze2 * w1 * w1 ktur= min(max(kamp, dked_min), dked_max) dktur(jk) = ktur ! ! update of dked = dked_gw + k_turb_mf ! dked(jl, jk) = dked(jl, jk) +ktur enddo ! ! apply eddy effects due to GWs: explicit scheme Kzz*dt/dz2 < 0.5 stability ! if (knob_ugwp_dokdis == 2) then do jk=ksrc,levs ze1 = min(.5*(dktur(jk) +dktur(jk-1)), dturb_max) kvint(jk) = kvint(jk) + ze1 ! ktint(jk) = ktint(jk) + ze1*iPr_ktgw enddo kvint(km1) = kvint(ksrc) kvint(ktop) = kvint(levs) dzmetm = 1./dz_met(km1) Adif(km1:levs) = 0. Cdif(km1:levs) = 0. do jk=km1,levs-1 dzmetp = 1./dz_met(jk+1) dzmetf = 1./(dz_meti(jk)*rhomid(jk)) ktur = kvint(jk) *rhoint(jk) * dzmetf kturp =Kvint(jk+1)*rhoint(jk+1) * dzmetf Adif(jk) = ktur * dzmetm Cdif(jk) = kturp * dzmetp ApC = adif(jk)+cdif(jk) ACdif(jk) = ApC w1 = ApC*iPr_max if (rdtp < w1 ) then Anstab(jk) = floor(w1*dtp) + 1 else Anstab(jk) = 1 endif dzmetm = dzmetp enddo nstab = maxval( Anstab(ksrc:levs-1)) ! if (nstab .ge. 3) print *, 'nstab ', nstab ! ! k instead Jk ! dtdif = dtp/real(nstab) ze1 = 1./dtdif do ist= 1, nstab do k=ksrc,levs-1 Bdif = ze1 - ACdif(k) Bt_dif = ze1 - ACdif(k)* iPr_ktgw ! ipr_Ktgw = 1./Pr <1 unew(k) = uold(k)*Bdif + uold(k-1)*Adif(k) + uold(k+1)*Cdif(k) vnew(k) = vold(k)*Bdif + vold(k-1)*Adif(k) + vold(k+1)*Cdif(k) tnew(k) = told(k)*Bt_dif+(told(k-1)*Adif(k) + told(k+1)*Cdif(k))*iPr_ktgw enddo uold(ksrc:levs-1) = unew(ksrc:levs-1)*dtdif ! value du/dtp *dtp = du vold(ksrc:levs-1) = vnew(ksrc:levs-1)*dtdif told(ksrc:levs-1) = tnew(ksrc:levs-1)*dtdif ! ! smoothing the boundary points: "k-1" = ksrc-1 and "k+1" = levs ! uold(levs) = uold(levs-1) vold(levs) = vold(levs-1) told(levs) = told(levs-1) enddo ! ! compute "smoothed" tendencies by molecular + GW-eddy diffusions ! do k=ksrc,levs-1 ! ! final updates of tendencies and diffusion ! ze2 = rdtp*(uold(k) - aum(k)) ze1 = rdtp*(vold(k) - avm(k)) pdtdt(jl,k)= rdtp*( told(k) - atm(k) ) if (abs(pdtdt(jl,k)) >= maxdtdt ) pdtdt(jl,k) = sign(maxdtdt,pdtdt(jl,k) ) if (abs(ze1) >= maxdudt ) then ze1 = sign(maxdudt, ze1) endif if (abs(ze2) >= maxdudt ) then ze2 = sign(maxdudt, ze2) endif pdudt(jl, k) = ze2 pdvdt(jl, k) = ze1 uz = uold(k+1) - uold(k-1) vz = vold(k+1) - vold(k-1) ze2 = 1./(dz_met(k+1)+dz_met(k) ) mf_diss_heat = rcpd*kvint(k)*(uz*uz +vz*vz)*ze2*ze2 ! vert grad heat pdtdt(jl,k)= pdtdt(jl,k) + mf_diss_heat ! extra heat due to eddy viscosity enddo ENDIF ! dissipative IF-loop for vertical eddy difusion u-v-t enddo ! J-loop ! RETURN !================================= diag print after "return" ====================== if (kdt ==1 .and. mpi_id == master) then ! print *, ' ugwpv1: nazd-nw-ilaunch=', nazd, nwav,ilaunch, maxval(kvg), ' kvg ' print *, 'ugwpv1: zdci(inc)=' , maxval(zdci), minval(zdci) print *, 'ugwpv1: zcimax=' , maxval(zci) ,' zcimin=' , minval(zci) ! print *, 'ugwpv1: tau_ngw=' , maxval(taub_src)*1.e3, minval(taub_src)*1.e3, tau_min print * endif if (kdt == 1 .and. mpi_id == master) then print *, 'vgw done nstab ', nstab ! print *, maxval(pdudt)*86400., minval(pdudt)*86400, 'vgw ax ugwp' print *, maxval(pdvdt)*86400., minval(pdvdt)*86400, 'vgw ay ugwp' print *, maxval(dked)*1., minval(dked)*1, 'vgw keddy m2/sec ugwp' print *, maxval(pdtdt)*86400., minval(pdtdt)*86400,'vgw eps ugwp' ! ! print *, ' ugwp -heating rates ' endif !================================= return end subroutine cires_ugwpv1_ngw_solv2 end module cires_ugwpv1_solv2