Merge #1364

1364: Fix bottom boundary condition on diffusive fluxes r=costachris a=costachris Fix bottom boundary condition on diffusive fluxes Co-authored-by: costachris <christopouloscosta@gmail.com>
CliMA · Oct 30, 2023 · 2236995 · 2236995
2 parents add5c50 + 6c2c458
commit 2236995
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 7 deletions.
diff --git a/driver/compute_diagnostics.jl b/driver/compute_diagnostics.jl
@@ -124,7 +124,9 @@ function compute_diagnostics!(
     FT = TC.float_type(state)
     N_up = TC.n_updrafts(edmf)
     aux_gm = TC.center_aux_grid_mean(state)
+    prog_gm_f = TC.face_prog_grid_mean(state)
     aux_en = TC.center_aux_environment(state)
+    aux_en_f = TC.face_aux_environment(state)
     aux_up = TC.center_aux_updrafts(state)
     aux_up_f = TC.face_aux_updrafts(state)
     aux_tc_f = TC.face_aux_turbconv(state)
@@ -145,6 +147,8 @@ function compute_diagnostics!(
     diag_tc_f_precip = face_diagnostics_precip(diagnostics)
     ρ_c = prog_gm.ρ
     p_c = aux_gm.p
+    ρ_f = aux_gm_f.ρ
+    w_gm = prog_gm_f.w
 
     diag_tc_svpc = svpc_diagnostics_turbconv(diagnostics)
     diag_svpc = svpc_diagnostics_grid_mean(diagnostics)
@@ -179,15 +183,18 @@ function compute_diagnostics!(
     end
 
     wvec = CC.Geometry.WVector
-    aeKHs_bc = -surf.ρs_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf]
+    aeKHs_bc = -surf.ρs_flux / aux_tc_f.ρ_ae_KH[kf_surf]
 
     If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0)))
     ∇f_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(FT(aeKHs_bc)), top = CCO.SetDivergence(FT(0)))
     massflux_s = aux_gm_f.massflux_s
-    parent(massflux_s) .= 0
     @. aux_gm_f.diffusive_flux_s = -aux_tc_f.ρ_ae_KH * ∇f_en(wvec(aux_en.s))
+
+    @. massflux_s = ρ_f * TC.ᶠinterp_a(a_en) * (aux_en_f.w - TC.toscalar(w_gm)) * (If(aux_en.s) - If(aux_gm.s))
     @inbounds for i in 1:N_up
-        @. massflux_s += aux_up_f[i].massflux * (If(aux_up[i].s) - If(aux_en.s))
+        w_up_i = aux_up_f[i].w
+        a_up = aux_up[i].area
+        @. massflux_s += ρ_f * (TC.ᶠinterp_a(a_up) * (w_up_i - TC.toscalar(w_gm)) * (If(aux_up[i].s) - If(aux_gm.s)))
     end
 
     # Mean water paths for calibration

diff --git a/src/EDMF_functions.jl b/src/EDMF_functions.jl
@@ -54,6 +54,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf
     ∂M∂z = aux_tc.∂M∂z
     ∂lnM∂z = aux_tc.∂lnM∂z
     massflux_c = aux_tc.massflux
+    parent(massflux_c) .= 0
+    parent(∂M∂z) .= 0
+    parent(∂lnM∂z) .= 0
 
     wvec = CC.Geometry.WVector
     ∇c = CCO.DivergenceF2C()
@@ -205,10 +208,10 @@ function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, sur
     @. aux_tc_f.ρ_ae_KH = IfKH(aeKH) * ρ_f
     @. aux_tc_f.ρ_ae_KQ = IfKQ(aeKQ) * ρ_f
 
-    aeKQq_tot_bc = -surf.ρq_tot_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KQ[kf_surf]
-    aeKHθ_liq_ice_bc = -surf.ρθ_liq_ice_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf]
-    aeKMu_bc = -surf.ρu_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf]
-    aeKMv_bc = -surf.ρv_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf]
+    aeKQq_tot_bc = -surf.ρq_tot_flux / aux_tc_f.ρ_ae_KQ[kf_surf]
+    aeKHθ_liq_ice_bc = -surf.ρθ_liq_ice_flux / aux_tc_f.ρ_ae_KH[kf_surf]
+    aeKMu_bc = -surf.ρu_flux / aux_tc_f.ρ_ae_KM[kf_surf]
+    aeKMv_bc = -surf.ρv_flux / aux_tc_f.ρ_ae_KM[kf_surf]
 
     aeKMuₕ_bc = CCG.UVVector(aeKMu_bc, aeKMv_bc)