subroutine ugwp_lsatdis_naz(levs, ksrc, nw, naz, kxw, taub_spect, ch, xaz, yaz, &
fcor, c2f2, dp, zmid, zint, pmid, pint, rho, ui, vi, ti, &
kvg, ktg, krad, kion, bn2i, bvi, rhoi, ax, ay, eps, ked, tau1)
!
! call ugwp_lsatdis_naz(levs, ksrc, nw, naz, kxw, taub_spect, ch, xaz, yaz, &
! fcor(j), c2f2(j), dp, zmid, zint, pmid, pint, rho, ui, vi, ti, &
! kvg, ktg, krad, kion, bn2i, bvi, rhoi, ax1, ay1, eps1, ked1)
use ugwp_common, only : rcpd, grav, rgrav
implicit none
!
integer :: levs, nw, naz, ksrc
real :: kxw
real, dimension(nw) :: taub_spect, ch
real, dimension(naz) :: xaz, yaz
real, dimension(levs+1) :: ui, vi, ti, bn2i, bvi, rhoi, zint, pint
real, dimension(levs ) :: dp, rho, pmid, zmid
real :: fcor, c2f2
real, dimension(levs+1) :: kvg, ktg, kion, krad, kmol
! output/locals
real, dimension(levs ) :: ax, ay, eps
real, dimension(levs+1) :: ked , tau1
real, dimension(levs+1 ) :: uaz
real, dimension(levs, naz ) :: epsd
real, dimension(levs+1, naz ) :: atau, kedd
real, dimension(levs+1 ) :: taux, tauy
real, dimension(levs ) :: dzirho , dzpi
real :: usrc
!
integer :: iaz, k
!
atau=0.0 ; epsd=0.0 ; kedd=0.0
do k=1,levs
dzpi(k) = -(pint(k+1)-pint(k))/rho(k)*rgrav
dzirho(k) = 1./rho(k)/dzpi(k) ! grav/abs(dp(k)) still hydrostatic "UGWP"
enddo
LOOP_IAZ: do iaz =1, naz
usrc = ui(ksrc)*xaz(iaz) +vi(ksrc)*yaz(iaz)
do k=1,levs+1
uaz(k) =ui(k)*xaz(iaz) +vi(k)*yaz(iaz) -usrc
enddo
!
! if (nw .le. 4) call stochastic ..ugwp_lsatdis_az1 only 4-waves ch_ngw1, fuw_ngw1, eff_ngw1=1
!
! multi-wave scheme
!
if (nw .gt. 4) then
call ugwp_lsatdis_az1(levs, ksrc, nw, kxw, ch, taub_spect, &
fcor, c2f2, zmid, zint, rho, uaz, ti, bn2i, bvi, rhoi, dzirho, dzpi, &
kvg, ktg, krad, kion, kmol, epsd(:, iaz), kedd(:,iaz), atau(:, iaz) )
endif
!
ENDDO LOOP_IAZ ! Azimuth of GW propagation directions
!
! sum over azimuth and project aTau(z, iza) =>(taux and tauy)
! for scalars for "wave-drag vector"
!
eps =0. ; ked =0.
do k=ksrc, levs
eps(k) = sum(epsd(k,:))*rcpd
enddo
do k=ksrc, levs+1
taux(k) = sum( atau(k,:)*xaz(:))
tauy(k) = sum( atau(k,:)*yaz(:))
ked(k) = sum(kedd(k,:))
enddo
tau1(ksrc:levs) = taux(ksrc:levs)
tau1(1:ksrc-1) = tau1(ksrc)
!
! end solver: gw_azimuth_solver_LS81
! sign Ax in rho*dU/dt = -d(rho*tau)/dz
! [(k) - (k+1)]
ax =0. ; ay = 0.
do k=ksrc, levs
ax(k) = dzirho(k)*(taux(k)-taux(k+1))
ay(k) = dzirho(k)*(tauy(k)-tauy(k+1))
enddo
call ugwp_limit_1d(ax, ay, eps, ked, levs)
return
!
print *
print *, ' Ax: ', maxval(Ax(ksrc:levs))*86400., minval(Ax(ksrc:levs))*86400.
print *, ' Ay: ', maxval(Ay(ksrc:levs))*86400., minval(Ay(ksrc:levs))*86400.
print *, 'Eps: ', maxval(Eps(ksrc:levs))*86400., minval(Eps(ksrc:levs))*86400.
print *, 'Ked: ', maxval(Ked(ksrc:levs))*1., minval(Ked(ksrc:levs))*1.
! print *, 'Atau ', maxval(atau(ksrc:levs, 1:Naz))*1.e3, minval(atau(ksrc:levs, 1:Naz))*1.e3
! print *, 'taux_gw: ', maxval(taux( ksrc:levs))*1.e3, minval(taux( ksrc:levs))*1.e3
print *
!-----------------------------------------------------------------------
! Here we can apply "ad-hoc" or/and "stability-based" limiters on
! (axy_gw, ked_gw and eps_gw) and check vert-inegrated conservation laws:
! energy and momentum and after that => final update gw-phys tendencies
!-----------------------------------------------------------------------
end subroutine ugwp_lsatdis_naz
!
subroutine ugwp_lsatdis_az1(levs, ksrc, nw, kxw, ch, taub_sp, &
fcor, c2f2, zm, zi, rho, um, tm, bn2, bn, rhoi, &
dzirho, dzpi, kvg, ktg, krad, kion, kmol, eps, ked, tau )
! call ugwp_lsatdis_az1(levs, ksrc, nw, kxw, ch, taub_spect, &
! fcor, c2f2, zmid, zint, rho, uaz, ti, bn2i, bvi, rhoi, dzirho, dzpi, &
! kvg, ktg, krad, kion, kmol, epsd(:, iaz), kedd(:,iaz), atau(:, iaz) )
use cires_ugwp_module, only : F_coriol, F_nonhyd, F_kds, linsat, linsat2
use cires_ugwp_module, only : iPr_ktgw, iPr_spgw, iPr_turb, iPr_mol
use cires_ugwp_module, only : rhp4, rhp2, rhp1, khp, cd_ulim
!
implicit NONE
!
integer, intent(in) :: nw ! number of GW modes in given direction
integer, intent(in) :: levs ! vertical layers
integer, intent(in) :: ksrc ! level of GW-launch layer
real , intent(in) :: kxw ! horizontal wavelength
real , intent(in) :: ch(nw) ! horizontal phase velocities
real , intent(in) :: taub_sp(nw) ! spectral distribution of the mom-flux
!
real, intent(in) :: fcor, c2f2 ! Corilois factors
real , intent(in) :: um(levs+1)
real , intent(in) :: tm(levs+1)
!in
real, intent(in), dimension(levs) :: rho, zm
real, intent(in), dimension(levs+1) :: rhoi, zi
real, intent(in), dimension(levs+1) :: bn2, bn
real, intent(in), dimension(levs) :: dzpi, dzirho
real, intent(in), dimension(levs+1) :: kvg, ktg, krad, kion, kmol
!========================================================================
!out
real, dimension(levs+1) :: tau, ked
real, dimension(levs) :: eps
!=========================================================================
!local
real :: Fd1, Fd2
real, dimension(levs) :: a_mkz
real, dimension(levs+1,nw) :: sp_tau, sp_ked, sp_kth
real, dimension(levs,nw) :: sp_eps
real, dimension(levs,nw) :: sp_mkz, sp_etot
real, dimension(levs,nw) :: sp_ek, sp_ep
real, dimension(levs) :: swg_ep, swg_ek, swg_et, swg_kz
real, dimension(nw) :: rtaus ! spectral distribution at ksrc
real :: sum_rtaus ! total flux in iaz-azimuth
real :: Chnorm, Cx, Cs, Cxs, Cx2sat
real :: Fdis, Fdisat
real :: Cdf2, Cdf1 ! (Cd*cd-f*f) and sqrt
!
! two-level => upward integration for wave-filtering (dissip + breaking)
!
real :: taus, tauk, tau_lin
real :: etws, etwk, etw_lin
real :: epss, epsk
real :: kds, kdk
real :: kzw, kzw2, kzw3, kzi, kzs
real :: wfd, wfi !
!
! for GW dissipation on the rotational sphere
!
real :: Betadis ! Ep/Ek ratio
real :: BetaM, BetaT ! 0.5 or 1./1+b and 1-1/(1+b)
real :: wfdM, wfdT, wfiM, wfiT, wdop2
real :: dzi, keff, keff_m, keff_t, keffs
real :: sf2k2, cf2
real :: Lzkm, Lzsat
integer :: i, k, igw
integer :: ksat1, ksat2
real :: zsat1, zsat2
real :: kx2_nh
real :: rab1, rab2, rab3, rab4, cd_ulim2
integer :: Ind_out(nw, levs+1)
!
logical, parameter :: dbg_print = .false.
!
!===================================================================
! Nullify arrays
! tau, eps, ked
!====================================================================
tau = 0.0
eps = 0.0
ked = 0.0
Ind_out(1:nw,:) = 0
!
! GW-spectral arrays ..... sp_etot ....sp_tau
!
sp_tau = 0.
sp_eps = 0.
sp_ked = 0.
sp_mkz = -99.
sp_etot = 0.
sp_ek = 0.
sp_ep = 0.
sp_kth = 0.
!
swg_et = 0.
swg_ep = 0.
swg_ek = 0.
swg_kz = 0.
cd_ulim2 = cd_ulim*cd_ulim
cf2 = F_coriol*c2f2
kx2_nh = F_nonhyd*kxw*kxw
if (dbg_print) then
write(6,*) linsat , ' eff-linsat & kx ', kxw
write(6,*) maxval(ch), minval(ch), ' ch '
write(6,*)
write(6,*) maxval(rhoi), minval(rhoi), 'rhoi '
write(6,*) zi(ksrc) , ' zi(ksrc) '
write(6,*) cd_ulim, ' crit-level cd_ulim '
write(6,*) F_coriol, ' F_coriol'
write(6,*) F_nonhyd, ' F_nonhyd '
write(6,*) maxval(Bn), minval(BN), ' BN-BV '
write(6,*) Um(ksrc), ' Um-ksrc ', cd_ulim2 , 'cd_ulim2 ', c2f2, ' c2f2 '
pause
endif
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Loop_GW: over GW-spectra
! of individual non-interactive modes
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
Loop_GW: do i=1, nw
!
Kds = 0.0
!
! src-level
!
Cx = ch(i) - Um(ksrc)
Cdf2 = Cx*Cx - cf2
taus = taub_sp(i) ! momentum flux for i-mode w/o rhoi(ksrc)
kzw = Bn(ksrc) / Ch(i) ! ch(i) > 0. Cx(i) < 0. critica
etws = taus*kzw / kxw
rtaus(i) = taus*rhoi(ksrc)
!
IF( Cx <= cd_ulim .or. Cdf2 <= cd_ulim2) THEN
Ind_out(i, ksrc) =-1 ! -1 - diagnostic index for critical levels
cycle Loop_GW ! got to the next mode of GW-spectra
ELSE
!
kzw2 = Bn2(ksrc)/Cdf2 - rhp4 - kx2_nh
!
if (kzw2 <= 0.) then
Ind_out(i, ksrc) =-2 ! -2 - diagnostic index for reflected waves
cycle Loop_GW ! no wave reflection in GW-LSD scheme
endif
kzw = sqrt(kzw2)
kzw3 = kzw2*kzw
etws = taus*kzw/kxw
!
! Here Linsat == Fr_critical
!
Cx2sat = Linsat2*Cdf2
if (etws >= cx2sat) then
Kds = kxw*Cx*rhp2/kzw3
etws = cx2sat
taus = etws*kxw/kzw
Ind_out(i, ksrc) =-3 ! -3 - dignostic index for saturated waves
endif
!
betadis = cdf2/(Cx*Cx+cf2)
betaM = 1.0 /(1.0+betadis)
betaT = 1.0 - BetaM
!
Cxs = Cx
kzs = kzw
! keffs = (kvg(ksrc)+kds)*iPr_turb*.5*khp
! sp_kth(ksrc, i) = rhoi(ksrc)*keffs*(Tm(ksrc)+Tm(ksrc-1))
rtaus(i) = taus*rhoi(ksrc)
sp_tau(ksrc, i) = rtaus(i)
sp_etot(ksrc, i) = etws
sp_mkz(ksrc, i) = kzw
sp_ek(ksrc, i) = etws*betam
sp_ep(ksrc, i) = etws*betaT ! can be transferred to (T'**2) T-rms
!
ENDIF ! vertical propagation of i-mode to the next upper layer = (ksrc+1)
!
! Loop_Zint .................................. VERTICAL "INTERFACE" LOOP from ksrc => ktop_GW
!
Loop_Zi: do k=ksrc+1, levs
!
Cx = ch(i)-Um(k) ! Um(k) is defined at the interface pressure levels
Cdf2 = Cx*Cx -cf2
if( Cx <= cd_ulim .or. Cdf2 <= 0.) then
Ind_out(i, k) =-1 ! 1 - diagnostic index for critical levels
! print*,'crit level C-U ',int(Cx),int(sqrt(cf2)),' Um ',Um(k)
cycle Loop_GW
endif
cdf1 =sqrt(Cdf2)
wdop2 = (kxw*Cx)* (kxw*Cx)
kzw2 = (Bn2(k)-wdop2)/Cdf2 - rhp4 - kx2_nh ! full lin DS-NIGW (N2-wd2)*k2=(m2+k2+[1/2H]^2)*(wd2-f2)
if (kzw2 < 0.) then
Ind_out(i, k) =-2 ! 2 - diagnostic index for reflected waves
cycle Loop_GW
endif
kzw = sqrt(kzw2)
kzw3 =kzw2*kzw
!
keff_m = kvg(k)*kzw2 + kion(k)
! keff_t = kturb(k)*iPr_turb + kmol(k)*iPr_mol
keff_t = ktg(k)*kzw2 + krad(k)
!
!
betadis = cdf2 / (Cx*Cx+cf2)
betaM = 1.0 / (1.0+betadis)
betaT = 1.0 - BetaM
!
!imaginary frequencies of momentum and heat with "kds at (k-1) level"
!
wfiM = kds*kzw2*F_kds + keff_m
wfiT = kds*iPr_ktgw*F_kds * kzw2 + keff_t
!
wfdM = wfiM/(kxw*Cdf1)*BetaM
wfdT = wfiT/(kxw*Cx)*BetaT
! exp-l: "kzi*dz"
kzi = 2.*kzw*(wfdM+wfdT)*dzpi(k) ! 2-factor energy-momentum (U')^2
!-------------------------------------------------------
! dissipative factor: Fdis
! we can replace WKB-solver by Numerical integration of
! tau_gw == etot_gw/kzw*kxw
! d(rho*tau_gw) = -kdis*rho*tau_gw
! |tau_gw| <= |tau_gwsat|
! linear limit for single mode
! generalization for the "broad" spectra
! or treating single mode breaking
! over finite "vertical"-depth with "efficiency"
! Now: time-step + hor-l scale
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Fdis = exp(-kzi)
!
!
! dissipative "wave rms" by WKB
!
etwk = etws*rhoi(k-1)/rhoi(k)*Fdis*kzw/kzs
!
Cx2sat = Linsat2*Cdf2
!
! Linear saturation
!
if (etwk.ge.cx2sat) then
Ind_out(i, k) =-3 ! 3 - dignostic index for saturated waves
! ! saturate energy and "trigger" keddy <kds>
etw_lin = etwk
etwk = cx2sat
Kds = kxw*Cdf1*rhp2/kzw3
tauk = etwk*kxw/kzw
!===================================================================================
! WAM/case with high Kds tau_lin = (etw_lin-etwk)*kxw/kzw !tau_loss by sat theory
! Lzsat = 6,28/kzw Zsat1 = Zi(k)-.5*Lzsat
! Zsat2 = Zi(k)+.5*Lzsat
! in WAM triggering from "kds = 0 m2/s" => "200 m2/s" for Lzw ~ 10 km
!
! call sat_domain(zi, Zsat1, Zsat2, pver, ksat1, ksat2)
!
! to avoid it do the new diss-n factor with eddy "kds" added to the
! background keff_m and keff_t
!
! can be taken out for the strato-mesosphere in GFS
! wfiM = kds*kzw2 + keff_m
! wfiT = kds*iPr_ktgw * kzw2 +keff_t
! wfdM = wfiM/(kxw*Cdf1)*BetaM
! wfdT = wfiT/(kxw*Cx)*BetaT
! kzi = 2.*kzw*(wfdM+wfdT)*dzpi(k)
! Fdisat = exp(-kzi)
! etwk = etws*rhoi(k-1)/rhoi(k)*Fdis*(kzw/Kzs)
! updated breaking in the Lzsat-domain: zsat1 < zi < zsat2
! =================================================================================
else
kds = 0.0
tauk = etwk*kxw/kzw ! <u'w'> = Ekin*kx/kz
ENDIF
!--------------------------------------
!
! Fill in spectral arrays(levs, nw)
!
!--------------------------------------
sp_ked(k,i) = kds ! defined at interfaces
sp_tau(k, i) = tauk*rhoi(k) ! defined at interfaces
! keff = (kds + kvg(k))*iPr_turb*0.5*KHP
! sp_kth(k, i) = rhoi(k)*keff*(Tm(k)+Tm(k-1)) ! defined at mid-layers
sp_etot(k, i) = etwk ! defined at interfaces
sp_mkz(k, i) = kzw ! defined at interfaces
sp_ek(k, i) = etwk*betam ! defined at interfaces
sp_ep(k, i) = etwk*betaT ! can be transferred to (T'**2)
!
!
if (sp_tau(k,i) > sp_tau(k-1,i)) then
sp_tau(k,i) = sp_tau(k-1,i) ! prevent "possible" numerical "noise"
endif
!
! updates for "eps and keff" from
!
rab1 =.5*(cx+cxs)*dzirho(k)
! heating
! due to wave dissipation
!
sp_eps(k,i) = rab1*(sp_tau(k-1,i)- sp_tau(k,i)) ! defined at mid-layers
!
! cooling term due to eddy heat conduction =0 if Keff_cond =>0,
! usually updated by 1D-heat implict tridiagonal solver
! explicit local solver ---->sp_kth(k,i) = Kt*(dT/dz+ R/Cp*T/Hp~>g/cp)
!
! sp_eps(k,i)=sp_eps(k,i)+dzirho(k)*(sp_kth(k,i)- sp_kth(k-1,i))
!
kzs = kzw
cxs = cX
taus = tauk
etws = etwk
! keffs = keff
enddo Loop_Zi ! ++++++++++++++ vertical layer
!
! ................................! stop ' in solver single-mode'
!
enddo Loop_GW ! i-mode of GW-spectra
!
sum_rtaus =sum(rtaus) ! total momentum flux at k=ksrc
! print *, sum_rtaus, ' tau-src ', nint(zi(ksrc)*1.e-3)
! print *, maxval(ch), minval(ch), ' Ch ', ngwv, ' N-modes '
!
!==============================================================================
! Perform spectral integartion (sum) & apply "efficiency/inremittency" factors
!
! eff_factor: ~ 1./[number of modes in 1-direction of model columns]
!
!==============================================================================
do k=ksrc, levs
ked(k) =0.
Eps(k) = 0.
Tau(k) = 0.
swg_et(k) =0.
swg_ep(k) =0.
swg_ek(k) =0.
do i=1,nw
Ked(k) = Ked(k)+sp_ked(k,i)
Eps(k) = Eps(k)+sp_eps(k,i)
Tau(k) = Tau(k)+sp_tau(k,i)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! GW-energy + GW-en flux ~ Cgz*E, diagnostics-only
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
swg_et(k) = swg_et(k)+sp_etot(k,i) !*eff_fact
swg_ep(k) = swg_ep(k)+sp_ep(k,i) !*eff_fact
swg_ek(k) = swg_ek(k)+sp_ek(k,i) !*eff_fact
enddo
enddo
! fill in below the "source" level ..... [1:ksrc-1]
!
do k=1, ksrc-1
! no loss of the total momentum flux
ked(k) =0.
eps(k) = 0.
tau(k) = tau(ksrc)
! lin-theory diagnostics-only
swg_et(k) =swg_et(ksrc)*rhoi(ksrc)/rhoi(k)
swg_ep(k) =swg_ep(ksrc)*rhoi(ksrc)/rhoi(k)
swg_ek(k) =swg_ek(ksrc)*rhoi(ksrc)/rhoi(k)
enddo
!
RETURN
!
! diagnostics below
!
345 FORMAT(2x, F8.2, 4(2x, F10.3), 2x, F8.2)
if (dbg_print) then
print *
print *, ' Zkm EK m2/s2 Ked m2/s Eps m2/s3 tau-Mpa '
do k=ksrc, levs
! Fd1 = maxval(Fdis_modes(1:nw,k))
! Fd2 = minval(Fdis_modes(1:nw,k))
write(6, 345) Zi(k)*1.e-3, sqrt(swg_ek(k)), Ked(k), Eps(k), Tau(k)*1.e3, Um(k) !, Fd1, Fd2
enddo
print *
write(6,*) nw , ' nwaves-linsat '
write(6,*) maxval(sp_ked), minval(sp_ked), 'ked '
write(6,*) maxval(sp_tau), minval(sp_tau), 'sp_tau '
pause
endif
!
end subroutine ugwp_lsatdis_az1
!
subroutine ugwp_limit_1d(ax, ay,eps, ked, levs)
use cires_ugwp_module, only : max_kdis, max_eps, max_axyz
implicit none
integer :: levs
real, dimension(levs) :: ax, ay,eps
real, dimension(levs+1) :: ked
real, parameter :: xtiny = 1.e-30
where (abs(ax) > max_axyz ) ax = ax/abs(ax+xtiny)*max_axyz
where (abs(ay) > max_axyz ) ay = ay/abs(ay+xtiny)*max_axyz
where (abs(eps) > max_eps ) eps = eps/abs(eps+xtiny)*max_eps
where (ked > max_kdis ) ked = max_kdis
end subroutine ugwp_limit_1d