Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

use total energy as grid mean variable #960

Merged
merged 1 commit into from
May 25, 2022
Merged

use total energy as grid mean variable #960

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions .buildkite/pipeline.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,12 +154,12 @@ steps:
command: "julia --color=yes --project=integration_tests integration_tests/driver.jl --case Rico --thermo_covariance_model prognostic --skip_tests true --suffix _prognostic_covarinaces"
artifact_paths: "Output.Rico.01_prognostic_covarinaces/stats/comparison/*"

- group: "3D Experiments"
steps:
# - group: "3D Experiments"
# steps:

- label: ":package: 3D Bomex"
command: "julia --color=yes --project=integration_tests integration_tests/3dBomex.jl"
artifact_paths: "integration_tests/output/3dBomex/*"
# - label: ":package: 3D Bomex"
# command: "julia --color=yes --project=integration_tests integration_tests/3dBomex.jl"
# artifact_paths: "integration_tests/output/3dBomex/*"

- group: "Performance monitoring"
steps:
Expand Down
6 changes: 3 additions & 3 deletions driver/Forcing.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,19 +9,19 @@ function initialize(self::TC.ForcingBase{TC.ForcingLES}, grid, state, LESDat::TC
zc_les = Array(TC.get_nc_data(data, "zc"))

dTdt_hadv = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dtdt_hadv", imin, imax))
H_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "thetali_mean", imin, imax))
T_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "temperature_mean", imin, imax))
dTdt_fluc = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dtdt_fluc", imin, imax))
dqtdt_hadv = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dqtdt_hadv", imin, imax))
qt_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "qt_mean", imin, imax))
dqtdt_fluc = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dqtdt_fluc", imin, imax))
subsidence = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "ls_subsidence", imin, imax))
u_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "u_mean", imin, imax))
v_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "v_mean", imin, imax))
(; dTdt_hadv, H_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
(; dTdt_hadv, T_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
end
for k in TC.real_center_indices(grid)
aux_gm.dTdt_hadv[k] = nt.dTdt_hadv[k]
aux_gm.H_nudge[k] = nt.H_nudge[k]
aux_gm.T_nudge[k] = nt.T_nudge[k]
aux_gm.dTdt_fluc[k] = nt.dTdt_fluc[k]
aux_gm.dqtdt_hadv[k] = nt.dqtdt_hadv[k]
aux_gm.qt_nudge[k] = nt.qt_nudge[k]
Expand Down
8 changes: 4 additions & 4 deletions driver/Surface.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ function get_surface(
ch = result.Ch,
ρu_flux = surf_params.zero_uv_fluxes ? FT(0) : result.ρτxz,
ρv_flux = surf_params.zero_uv_fluxes ? FT(0) : result.ρτyz,
ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
ρe_tot_flux = shf + lhf,
yairchn marked this conversation as resolved.
Show resolved Hide resolved
ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
wstar = convective_vel,
ρq_liq_flux = FT(0),
Expand Down Expand Up @@ -126,7 +126,7 @@ function get_surface(
ustar = result.ustar,
ρu_flux = result.ρτxz,
ρv_flux = result.ρτyz,
ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
ρe_tot_flux = shf + lhf,
ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
bflux = result.buoy_flux,
wstar = convective_vel,
Expand Down Expand Up @@ -183,7 +183,7 @@ function get_surface(
ustar = result.ustar,
ρu_flux = result.ρτxz,
ρv_flux = result.ρτyz,
ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
ρe_tot_flux = shf + lhf,
ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
bflux = result.buoy_flux,
wstar = convective_vel,
Expand Down Expand Up @@ -245,7 +245,7 @@ function get_surface(
ustar = result.ustar,
ρu_flux = result.ρτxz,
ρv_flux = result.ρτyz,
ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in),
ρe_tot_flux = shf + lhf,
ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in),
bflux = result.buoy_flux,
wstar = convective_vel,
Expand Down
2 changes: 0 additions & 2 deletions driver/compute_diagnostics.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,8 +139,6 @@ function compute_diagnostics!(

@inbounds for k in TC.real_center_indices(grid)
aux_gm.s[k] = TD.specific_entropy(param_set, ts_gm[k])
aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
aux_en.s[k] = TD.specific_entropy(param_set, ts_en[k])
end
@inbounds for k in TC.real_center_indices(grid)
Expand Down
133 changes: 93 additions & 40 deletions driver/dycore.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,15 +62,15 @@ cent_aux_vars_gm(FT, edmf) = (;
subsidence = FT(0), #Large-scale subsidence
dTdt_hadv = FT(0), #Horizontal advection of temperature
dqtdt_hadv = FT(0), #Horizontal advection of moisture
H_nudge = FT(0), #Reference H profile for relaxation tendency
T_nudge = FT(0), #Reference T profile for relaxation tendency
qt_nudge = FT(0), #Reference qt profile for relaxation tendency
dTdt_fluc = FT(0), #Vertical turbulent advection of temperature
dqtdt_fluc = FT(0), #Vertical turbulent advection of moisture
u_nudge = FT(0), #Reference u profile for relaxation tendency
v_nudge = FT(0), #Reference v profile for relaxation tendency
ug = FT(0), #Geostrophic u velocity
vg = FT(0), #Geostrophic v velocity
∇θ_liq_ice_gm = FT(0),
∇MSE_gm = FT(0),
∇q_tot_gm = FT(0),
cent_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
θ_virt = FT(0),
Expand All @@ -89,7 +89,7 @@ face_aux_vars_gm(FT, edmf) = (;
diffusive_flux_s = FT(0),
total_flux_s = FT(0),
f_rad = FT(0),
sgs_flux_θ_liq_ice = FT(0),
sgs_flux_h_tot = FT(0),
sgs_flux_q_tot = FT(0),
face_aux_vars_gm_moisture(FT, edmf.moisture_model)...,
sgs_flux_u = FT(0),
Expand Down Expand Up @@ -139,7 +139,7 @@ cent_prognostic_vars_gm(::Type{FT}, local_geometry, edmf) where {FT} = (;
ρ = FT(0),
u = FT(0),
v = FT(0),
ρθ_liq_ice = FT(0),
ρe_tot = FT(0),
ρq_tot = FT(0),
# TODO: Change to:
# uₕ = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry),
Expand Down Expand Up @@ -221,7 +221,6 @@ function compute_ref_state!(
@info "z_span = $z_span"
prob = ODE.ODEProblem(rhs, logp, z_span)
sol = ODE.solve(prob, ODE.Tsit5(), reltol = 1e-12, abstol = 1e-12)

parent(p_f) .= sol.(vec(grid.zf))
parent(p_c) .= sol.(vec(grid.zc))

Expand All @@ -241,23 +240,43 @@ function compute_ref_state!(
return nothing
end

function set_prog_from_aux!(state, ::TC.ρθVar)

function set_thermo_state_peq!(state, grid, moisture_model, param_set)
Ic = CCO.InterpolateF2C()
FT = eltype(grid)
ts_gm = TC.center_aux_grid_mean(state).ts
prog_gm = TC.center_prog_grid_mean(state)
prog_gm_f = TC.face_prog_grid_mean(state)
aux_gm = TC.center_aux_grid_mean(state)
@. prog_gm.ρθ_liq_ice = prog_gm.ρ * aux_gm.θ_liq_ice
@. prog_gm.ρq_tot .= prog_gm.ρ * aux_gm.q_tot
p_c = aux_gm.p
ρ_c = prog_gm.ρ
w_c = copy(prog_gm.ρe_tot)
charleskawczynski marked this conversation as resolved.
Show resolved Hide resolved
@. w_c = Ic(FT(0) + prog_gm_f.w)
@inbounds for k in TC.real_center_indices(grid)
thermo_args = if moisture_model isa TC.EquilibriumMoisture
()
tapios marked this conversation as resolved.
Show resolved Hide resolved
elseif moisture_model isa TC.NonEquilibriumMoisture
(prog_gm.q_liq[k], prog_gm.q_ice[k])
else
error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium")
end
e_kin = TC.kinetic_energy(prog_gm.u[k], prog_gm.v[k], w_c[k])
e_pot = TC.geopotential(param_set, grid.zc.z[k])
e_int = prog_gm.ρe_tot[k] / ρ_c[k] - e_kin - e_pot
ts_gm[k] = TC.thermo_state_peq(param_set, p_c[k], e_int, aux_gm.q_tot[k], thermo_args...)
aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp(param_set, ts_gm[k])
aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
end
return nothing
end

function set_thermo_state!(state, grid, moisture_model, param_set)
function set_thermo_state_pθq!(state, grid, moisture_model, param_set)
Ic = CCO.InterpolateF2C()
ts_gm = TC.center_aux_grid_mean(state).ts
prog_gm = TC.center_prog_grid_mean(state)
aux_gm = TC.center_aux_grid_mean(state)
p_c = aux_gm.p
ρ_c = prog_gm.ρ
@inbounds for k in TC.real_center_indices(grid)
aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
thermo_args = if moisture_model isa TC.EquilibriumMoisture
()
elseif moisture_model isa TC.NonEquilibriumMoisture
Expand All @@ -270,14 +289,31 @@ function set_thermo_state!(state, grid, moisture_model, param_set)
return nothing
end

function set_grid_mean_from_thermo_state!(param_set, state, grid)
Ic = CCO.InterpolateF2C()
ts_gm = TC.center_aux_grid_mean(state).ts
prog_gm = TC.center_prog_grid_mean(state)
prog_gm_f = TC.face_prog_grid_mean(state)
aux_gm = TC.center_aux_grid_mean(state)
ρ_c = prog_gm.ρ
FT = eltype(grid)
w_c = copy(prog_gm.ρe_tot)
@. w_c = Ic(FT(0) + prog_gm_f.w)
@inbounds for k in TC.real_center_indices(grid)
e_kin = TC.kinetic_energy(prog_gm.u[k], prog_gm.v[k], w_c[k])
e_pot = TC.geopotential(param_set, grid.zc.z[k])
prog_gm.ρe_tot[k] = ρ_c[k] * TD.total_energy(param_set, ts_gm[k], e_kin, e_pot)
prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k]
end
return nothing
end

function assign_thermo_aux!(state, grid, moisture_model, param_set)
aux_gm = TC.center_aux_grid_mean(state)
prog_gm = TC.center_prog_grid_mean(state)
ts_gm = TC.center_aux_grid_mean(state).ts
ρ_c = prog_gm.ρ
@inbounds for k in TC.real_center_indices(grid)
aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k]
aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k]
ts = ts_gm[k]
aux_gm.q_liq[k] = TD.liquid_specific_humidity(param_set, ts)
aux_gm.q_ice[k] = TD.ice_specific_humidity(param_set, ts)
Expand Down Expand Up @@ -311,7 +347,7 @@ function ∑tendencies!(tendencies::FV, prog::FV, params::NT, t::Real) where {NT

state = TC.State(prog, aux, tendencies)

set_thermo_state!(state, grid, edmf.moisture_model, param_set)
set_thermo_state_peq!(state, grid, edmf.moisture_model, param_set)

# TODO: where should this live?
aux_gm = TC.center_aux_grid_mean(state)
Expand Down Expand Up @@ -365,14 +401,19 @@ function compute_gm_tendencies!(
force::TC.ForcingBase,
param_set::APS,
)
Ic = CCO.InterpolateF2C()
R_d = CPP.R_d(param_set)
T_0 = CPP.T_0(param_set)
Lv_0 = CPP.LH_v0(param_set)
tendencies_gm = TC.center_tendencies_grid_mean(state)
kc_toa = TC.kc_top_of_atmos(grid)
kf_surf = TC.kf_surface(grid)
FT = eltype(grid)
prog_gm = TC.center_prog_grid_mean(state)
prog_gm_f = TC.face_prog_grid_mean(state)
aux_gm = TC.center_aux_grid_mean(state)
aux_gm_f = TC.face_aux_grid_mean(state)
∇θ_liq_ice_gm = TC.center_aux_grid_mean(state).∇θ_liq_ice_gm
∇MSE_gm = TC.center_aux_grid_mean(state).∇MSE_gm
∇q_tot_gm = TC.center_aux_grid_mean(state).∇q_tot_gm
aux_en = TC.center_aux_environment(state)
aux_en_f = TC.face_aux_environment(state)
Expand All @@ -384,14 +425,20 @@ function compute_gm_tendencies!(
aux_tc = TC.center_aux_turbconv(state)
ts_gm = TC.center_aux_grid_mean(state).ts

θ_liq_ice_gm_toa = aux_gm.θ_liq_ice[kc_toa]
q_tot_gm_toa = aux_gm.q_tot[kc_toa]
RBθ = CCO.RightBiasedC2F(; top = CCO.SetValue(θ_liq_ice_gm_toa))
e_kin = copy(prog_gm.ρe_tot)
w_c = copy(prog_gm.ρe_tot)
yairchn marked this conversation as resolved.
Show resolved Hide resolved
h_tot_gm = copy(prog_gm.ρe_tot)
@. w_c = Ic(FT(0) + prog_gm_f.w)
@. h_tot_gm = TC.anelastic_total_enthalpy(param_set, prog_gm.ρe_tot / ρ_c, ts_gm)
@. e_kin = TC.kinetic_energy(prog_gm.u, prog_gm.v, w_c)
MSE_gm_toa = h_tot_gm[kc_toa] - e_kin[kc_toa]
q_tot_gm_toa = prog_gm.ρq_tot[kc_toa] / ρ_c[kc_toa]
RBe = CCO.RightBiasedC2F(; top = CCO.SetValue(MSE_gm_toa))
RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa))
wvec = CC.Geometry.WVector
∇c = CCO.DivergenceF2C()
@. ∇θ_liq_ice_gm = ∇c(wvec(RBθ(aux_gm.θ_liq_ice)))
@. ∇q_tot_gm = ∇c(wvec(RBq(aux_gm.q_tot)))
@. ∇MSE_gm = ∇c(wvec(RBe(h_tot_gm - e_kin)))
@. ∇q_tot_gm = ∇c(wvec(RBq(prog_gm.ρq_tot / ρ_c)))

if edmf.moisture_model isa TC.NonEquilibriumMoisture
∇q_liq_gm = TC.center_aux_grid_mean(state).∇q_liq_gm
Expand All @@ -406,28 +453,33 @@ function compute_gm_tendencies!(

# Apply forcing and radiation
@inbounds for k in TC.real_center_indices(grid)
Π = TD.exner(param_set, ts_gm[k])
cp_m = TD.cp_m(param_set, ts_gm[k])
cp_v = CPP.cp_v(param_set)
cv_v = CPP.cv_v(param_set)
R_v = CPP.R_v(param_set)
cv_m = TD.cv_m(param_set, ts_gm[k])
h_v = cv_v * (aux_gm.T[k] - T_0) + Lv_0 - R_v * T_0

# Coriolis
if force.apply_coriolis
tendencies_gm.u[k] -= force.coriolis_param * (aux_gm.vg[k] - prog_gm.v[k])
tendencies_gm.v[k] += force.coriolis_param * (aux_gm.ug[k] - prog_gm.u[k])
end
# LS Subsidence
tendencies_gm.ρθ_liq_ice[k] -= ρ_c[k] * ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k]
tendencies_gm.ρq_tot[k] -= ρ_c[k] * ∇q_tot_gm[k] * aux_gm.subsidence[k]
tendencies_gm.ρe_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇MSE_gm[k]
tendencies_gm.ρq_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇q_tot_gm[k]
if edmf.moisture_model isa TC.NonEquilibriumMoisture
tendencies_gm.q_liq[k] -= ∇q_liq_gm[k] * aux_gm.subsidence[k]
tendencies_gm.q_ice[k] -= ∇q_ice_gm[k] * aux_gm.subsidence[k]
end
# Radiation
if TC.rad_type(radiation) <: Union{TC.RadiationDYCOMS_RF01, TC.RadiationLES, TC.RadiationTRMM_LBA}
tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_rad[k] / Π
tendencies_gm.ρe_tot[k] += ρ_c[k] * cv_m * aux_gm.dTdt_rad[k]
end
# LS advection
yairchn marked this conversation as resolved.
Show resolved Hide resolved
tendencies_gm.ρq_tot[k] += ρ_c[k] * aux_gm.dqtdt_hadv[k]
if !(TC.force_type(force) <: TC.ForcingDYCOMS_RF01)
tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_hadv[k] / Π
tendencies_gm.ρe_tot[k] += ρ_c[k] * (cp_m * aux_gm.dTdt_hadv[k] + h_v * aux_gm.dqtdt_hadv[k])
end
if edmf.moisture_model isa TC.NonEquilibriumMoisture
tendencies_gm.q_liq[k] += aux_gm.dqldt[k]
Expand All @@ -436,21 +488,22 @@ function compute_gm_tendencies!(

# LES specific forcings
if TC.force_type(force) <: TC.ForcingLES
H_fluc = aux_gm.dTdt_fluc[k] / Π
T_fluc = aux_gm.dTdt_fluc[k]

gm_U_nudge_k = (aux_gm.u_nudge[k] - prog_gm.u[k]) / force.wind_nudge_τᵣ
gm_V_nudge_k = (aux_gm.v_nudge[k] - prog_gm.v[k]) / force.wind_nudge_τᵣ

Γᵣ = TC.compute_les_Γᵣ(grid.zc[k], force.scalar_nudge_τᵣ, force.scalar_nudge_zᵢ, force.scalar_nudge_zᵣ)
gm_H_nudge_k = Γᵣ * (aux_gm.H_nudge[k] - aux_gm.θ_liq_ice[k])
gm_T_nudge_k = Γᵣ * (aux_gm.T_nudge[k] - aux_gm.T[k])
gm_q_tot_nudge_k = Γᵣ * (aux_gm.qt_nudge[k] - aux_gm.q_tot[k])
if edmf.moisture_model isa TC.NonEquilibriumMoisture
gm_q_liq_nudge_k = Γᵣ * (aux_gm.ql_nudge[k] - prog_gm.q_liq[k])
gm_q_ice_nudge_k = Γᵣ * (aux_gm.qi_nudge[k] - prog_gm.q_ice[k])
end

tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * (gm_H_nudge_k + H_fluc)
tendencies_gm.ρq_tot[k] += ρ_c[k] * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k])
tendencies_gm.ρe_tot[k] +=
ρ_c[k] * (cv_m * (gm_T_nudge_k + T_fluc) + h_v * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k]))
tendencies_gm.u[k] += gm_U_nudge_k
tendencies_gm.v[k] += gm_V_nudge_k
if edmf.moisture_model isa TC.NonEquilibriumMoisture
Expand All @@ -466,38 +519,38 @@ function compute_gm_tendencies!(
aux_en.qt_tendency_precip_formation[k] +
aux_tc.qt_tendency_precip_sinks[k]
)
tendencies_gm.ρθ_liq_ice[k] +=

tendencies_gm.ρe_tot[k] +=
ρ_c[k] * (
aux_bulk.θ_liq_ice_tendency_precip_formation[k] +
aux_en.θ_liq_ice_tendency_precip_formation[k] +
aux_tc.θ_liq_ice_tendency_precip_sinks[k]
aux_bulk.e_tot_tendency_precip_formation[k] +
aux_en.e_tot_tendency_precip_formation[k] +
aux_tc.e_tot_tendency_precip_sinks[k]
)

if edmf.moisture_model isa TC.NonEquilibriumMoisture
tendencies_gm.q_liq[k] += aux_bulk.ql_tendency_precip_formation[k] + aux_en.ql_tendency_precip_formation[k]
tendencies_gm.q_ice[k] += aux_bulk.qi_tendency_precip_formation[k] + aux_en.qi_tendency_precip_formation[k]
end
end

TC.compute_sgs_flux!(edmf, grid, state, surf)

sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice
TC.compute_sgs_flux!(edmf, grid, state, surf, param_set)
sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot
sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot
sgs_flux_u = aux_gm_f.sgs_flux_u
sgs_flux_v = aux_gm_f.sgs_flux_v
# apply surface BC as SGS flux at lowest level
sgs_flux_θ_liq_ice[kf_surf] = surf.ρθ_liq_ice_flux
sgs_flux_h_tot[kf_surf] = surf.ρe_tot_flux
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You have probably already discussed this, but why do we equate these?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the flux from the surface is the SGS flux at the lowest level (first face), at any other level that flux is given by the EDMF. I named it sgs_flux_h_tot to make sure it is clear that this is the value of our sgs_flux_h_tot at the lowest level and nothing else is missing here

sgs_flux_q_tot[kf_surf] = surf.ρq_tot_flux
sgs_flux_u[kf_surf] = surf.ρu_flux
sgs_flux_v[kf_surf] = surf.ρv_flux

tends_ρθ_liq_ice = tendencies_gm.ρθ_liq_ice
tends_ρe_tot = tendencies_gm.ρe_tot
tends_ρq_tot = tendencies_gm.ρq_tot
tends_u = tendencies_gm.u
tends_v = tendencies_gm.v

∇sgs = CCO.DivergenceF2C()

@. tends_ρθ_liq_ice += -∇sgs(wvec(sgs_flux_θ_liq_ice))
@. tends_ρe_tot += -∇sgs(wvec(sgs_flux_h_tot))
@. tends_ρq_tot += -∇sgs(wvec(sgs_flux_q_tot))
@. tends_u += -∇sgs(wvec(sgs_flux_u)) / ρ_c
@. tends_v += -∇sgs(wvec(sgs_flux_v)) / ρ_c
Expand Down
Loading