Add consistent environmental buoyancy gradients using mean fields and quadratures.

integration_tests/utils/generate_namelist.jl

@@ -142,6 +142,7 @@ function default_namelist(case_name::String)
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["constant_area"] = false
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["calculate_tke"] = true
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["mixing_length"] = "sbtd_eq"
+    namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["env_buoy_grad"] = "quadratures"
 
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_buoy"] = "normalmode"
     namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_drag"] = "normalmode"

integration_tests/utils/mse_tables.jl

@@ -5,58 +5,58 @@
 all_best_mse = OrderedCollections.OrderedDict()
 #
 all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict()
-all_best_mse["ARM_SGP"]["qt_mean"] = 0.26446718607725306
+all_best_mse["ARM_SGP"]["qt_mean"] = 2.6506823775591898e-01
 all_best_mse["ARM_SGP"]["updraft_area"] = 342.05657568721784
 all_best_mse["ARM_SGP"]["updraft_w"] = 147.98545782669382
 all_best_mse["ARM_SGP"]["updraft_qt"] = 28.294790221283453
 all_best_mse["ARM_SGP"]["updraft_thetal"] = 170.9801640138966
 all_best_mse["ARM_SGP"]["u_mean"] = 5.762470982857257e8
 all_best_mse["ARM_SGP"]["tke_mean"] = 986.3675515161622
 all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00012053009592043601
-all_best_mse["ARM_SGP"]["ql_mean"] = 204.94997013272453
+all_best_mse["ARM_SGP"]["ql_mean"] = 2.0547918525992608e+02
 all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00011122266036009708
-all_best_mse["ARM_SGP"]["Hvar_mean"] = 2633.2221109819784
-all_best_mse["ARM_SGP"]["QTvar_mean"] = 1559.544125199823
+all_best_mse["ARM_SGP"]["Hvar_mean"] = 2.9194184175128917e+03
+all_best_mse["ARM_SGP"]["QTvar_mean"] = 1.7213808112900319e+03
 #
 all_best_mse["Bomex"] = OrderedCollections.OrderedDict()
-all_best_mse["Bomex"]["qt_mean"] = 0.11034670891224195
-all_best_mse["Bomex"]["updraft_area"] = 132.0788865649432
-all_best_mse["Bomex"]["updraft_w"] = 18.392509423445038
-all_best_mse["Bomex"]["updraft_qt"] = 6.288835389580369
-all_best_mse["Bomex"]["updraft_thetal"] = 69.45747637501809
-all_best_mse["Bomex"]["v_mean"] = 66.15195917633564
-all_best_mse["Bomex"]["u_mean"] = 2244.359495472965
-all_best_mse["Bomex"]["tke_mean"] = 55.75185915047925
-all_best_mse["Bomex"]["temperature_mean"] = 4.3858868869668704e-5
-all_best_mse["Bomex"]["ql_mean"] = 8.228711881484596
-all_best_mse["Bomex"]["thetal_mean"] = 4.448304470557014e-5
-all_best_mse["Bomex"]["Hvar_mean"] = 4019.5984832795702
-all_best_mse["Bomex"]["QTvar_mean"] = 1505.300250338699
+all_best_mse["Bomex"]["qt_mean"] = 1.0983320350059607e-01
+all_best_mse["Bomex"]["updraft_area"] = 1.3209194605797390e+02
+all_best_mse["Bomex"]["updraft_w"] = 1.8757855416502004e+01
+all_best_mse["Bomex"]["updraft_qt"] = 6.5576868612678965e+00
+all_best_mse["Bomex"]["updraft_thetal"] = 6.9463116782644903e+01
+all_best_mse["Bomex"]["v_mean"] = 6.6112855769976591e+01
+all_best_mse["Bomex"]["u_mean"] = 2.2443513442723633e+03
+all_best_mse["Bomex"]["tke_mean"] = 5.5835472696934382e+01
+all_best_mse["Bomex"]["temperature_mean"] = 4.3677471984664055e-05
+all_best_mse["Bomex"]["ql_mean"] = 8.8506606580011589e+00
+all_best_mse["Bomex"]["thetal_mean"] = 4.4330671168211216e-05
+all_best_mse["Bomex"]["Hvar_mean"] = 4.2126196703263922e+03
+all_best_mse["Bomex"]["QTvar_mean"] = 1.5794628293458845e+03
 #
 all_best_mse["DryBubble"] = OrderedCollections.OrderedDict()
-all_best_mse["DryBubble"]["updraft_area"] = 0.0
-all_best_mse["DryBubble"]["updraft_w"] = 0.0
-all_best_mse["DryBubble"]["updraft_thetal"] = 0.0
+all_best_mse["DryBubble"]["updraft_area"] = 8.8530673193084747e-03
+all_best_mse["DryBubble"]["updraft_w"] = 1.2479889474537974e-02
+all_best_mse["DryBubble"]["updraft_thetal"] = 1.2223166344487256e-10
 all_best_mse["DryBubble"]["u_mean"] = 0.0
-all_best_mse["DryBubble"]["tke_mean"] = 0.0
-all_best_mse["DryBubble"]["temperature_mean"] = 0.0
-all_best_mse["DryBubble"]["thetal_mean"] = 0.0
-all_best_mse["DryBubble"]["Hvar_mean"] = 0.0
+all_best_mse["DryBubble"]["tke_mean"] = 2.0135570156946084e-01
+all_best_mse["DryBubble"]["temperature_mean"] = 3.2783493600492222e-07
+all_best_mse["DryBubble"]["thetal_mean"] = 4.5557150098548290e-07
+all_best_mse["DryBubble"]["Hvar_mean"] = 9.6713965877036813e+01
 #
 all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict()
-all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.04290066515065224
-all_best_mse["DYCOMS_RF01"]["ql_mean"] = 0.8976960436098287
-all_best_mse["DYCOMS_RF01"]["updraft_area"] = 31.57502795434024
-all_best_mse["DYCOMS_RF01"]["updraft_w"] = 5.415789637698043
-all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 1.734023675326015
-all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.186016892032455
-all_best_mse["DYCOMS_RF01"]["v_mean"] = 19017.018104720682
-all_best_mse["DYCOMS_RF01"]["u_mean"] = 74833.98050860969
-all_best_mse["DYCOMS_RF01"]["tke_mean"] = 18.837749675831937
-all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 0.00016283538235415605
-all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 0.00016376758990380806
-all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1227.7967628483952
-all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 478.07194167693547
+all_best_mse["DYCOMS_RF01"]["qt_mean"] = 3.0682630810336644e-02
+all_best_mse["DYCOMS_RF01"]["ql_mean"] = 2.7907515016732569e+00
+all_best_mse["DYCOMS_RF01"]["updraft_area"] = 2.9996855083383483e+01
+all_best_mse["DYCOMS_RF01"]["updraft_w"] = 3.6786963992723485e+00
+all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 1.8434661312139911e+00
+all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 4.6186138798587415e+01
+all_best_mse["DYCOMS_RF01"]["v_mean"] = 1.9100567571407046e+04
+all_best_mse["DYCOMS_RF01"]["u_mean"] = 7.5180517589525669e+04
+all_best_mse["DYCOMS_RF01"]["tke_mean"] = 2.4119502493062683e+01
+all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 1.5145561673937695e-04
+all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 1.5230990980765101e-04
+all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1.2554483566036963e+03
+all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 5.1751103217917216e+02
 #
 all_best_mse["GABLS"] = OrderedCollections.OrderedDict()
 all_best_mse["GABLS"]["updraft_thetal"] = 5.584864486190616e-9
@@ -70,80 +70,80 @@ all_best_mse["GABLS"]["QTvar_mean"] = 0.017318287479694217
 all_best_mse["GABLS"]["qt_mean"] = 0.21092607675396896
 #
 all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict()
-all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["v_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 0.0
-all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 0.0
+all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 3.9000477871154622e-09
+all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 3.2093573944402762e-06
+all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 1.6350375660088284e-07
+all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 1.5548562332198521e-06
+all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 7.0074885089062560e-10
+all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 2.4326279199247991e-13
+all_best_mse["life_cycle_Tan2018"]["v_mean"] = 5.3858193910463953e-07
+all_best_mse["life_cycle_Tan2018"]["u_mean"] = 3.1825523774843913e-09
+all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 1.9880489190175280e-07
+all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 4.6974217695746801e-12
+all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 4.6149603558343825e-12
+all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 2.7791377140463445e-02
+all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 2.4478780510096252e-02
 #
 all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict()
 all_best_mse["Nieuwstadt"]["updraft_area"] = 110.65702527331422
 all_best_mse["Nieuwstadt"]["updraft_w"] = 12.298314607367429
-all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.3670456600826
+all_best_mse["Nieuwstadt"]["updraft_thetal"] = 1.1736704568882135e+02
 all_best_mse["Nieuwstadt"]["u_mean"] = 59.73642032443851
-all_best_mse["Nieuwstadt"]["tke_mean"] = 222.28508474371915
+all_best_mse["Nieuwstadt"]["tke_mean"] = 2.2229785052178343e+02
 all_best_mse["Nieuwstadt"]["temperature_mean"] = 9.78812700698016e-6
 all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.0097201478188357e-5
-all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1276.761018651613
+all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1.2769135065447590e+03
 #
 all_best_mse["Rico"] = OrderedCollections.OrderedDict()
-all_best_mse["Rico"]["qt_mean"] = 1.3562728436394182
+all_best_mse["Rico"]["qt_mean"] = 1.3576174162823518e+00
 all_best_mse["Rico"]["updraft_area"] = 483.4022083273521
-all_best_mse["Rico"]["updraft_w"] = 105.6746441121623
-all_best_mse["Rico"]["updraft_qt"] = 10.859715229316084
+all_best_mse["Rico"]["updraft_w"] = 1.0605422046690978e+02
+all_best_mse["Rico"]["updraft_qt"] = 1.0914863189252404e+01
 all_best_mse["Rico"]["updraft_thetal"] = 132.79730672922278
-all_best_mse["Rico"]["v_mean"] = 23893.44384637507
-all_best_mse["Rico"]["u_mean"] = 167.0481771976619
-all_best_mse["Rico"]["tke_mean"] = 82.85087674864991
-all_best_mse["Rico"]["temperature_mean"] = 0.0006129402538716675
-all_best_mse["Rico"]["ql_mean"] = 67.02280413541271
+all_best_mse["Rico"]["v_mean"] = 2.3913389154144726e+04
+all_best_mse["Rico"]["u_mean"] = 1.6720325360173803e+02
+all_best_mse["Rico"]["tke_mean"] = 8.2957963768830496e+01
+all_best_mse["Rico"]["temperature_mean"] = 6.1310148877327848e-04
+all_best_mse["Rico"]["ql_mean"] = 6.7286920396505195e+01
 all_best_mse["Rico"]["qr_mean"] = 750.8723638131503
-all_best_mse["Rico"]["thetal_mean"] = 0.0006043563413443778
+all_best_mse["Rico"]["thetal_mean"] = 6.0448834814914290e-04
 all_best_mse["Rico"]["Hvar_mean"] = 18654.16088756812
 all_best_mse["Rico"]["QTvar_mean"] = 4036.079013508425
 #
 all_best_mse["Soares"] = OrderedCollections.OrderedDict()
-all_best_mse["Soares"]["qt_mean"] = 0.12636886303765799
+all_best_mse["Soares"]["qt_mean"] = 1.2659089197243104e-01
 all_best_mse["Soares"]["updraft_area"] = 112.7273316705158
-all_best_mse["Soares"]["updraft_w"] = 11.265484888980165
-all_best_mse["Soares"]["updraft_qt"] = 23.119052296484572
-all_best_mse["Soares"]["updraft_thetal"] = 65.25334973003676
-all_best_mse["Soares"]["u_mean"] = 101.82502556650786
-all_best_mse["Soares"]["tke_mean"] = 184.8671644305064
+all_best_mse["Soares"]["updraft_w"] = 1.1266999168717131e+01
+all_best_mse["Soares"]["updraft_qt"] = 2.3120736199996990e+01
+all_best_mse["Soares"]["updraft_thetal"] = 6.5253349774496741e+01
+all_best_mse["Soares"]["u_mean"] = 1.0182505005403843e+02
+all_best_mse["Soares"]["tke_mean"] = 1.8487029454315532e+02
 all_best_mse["Soares"]["temperature_mean"] = 1.1031965842549474e-5
 all_best_mse["Soares"]["thetal_mean"] = 1.0453842793106344e-5
 all_best_mse["Soares"]["Hvar_mean"] = 1145.5445941618634
 #
 all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict()
 all_best_mse["TRMM_LBA"]["qt_mean"] = 3.696894480356394
-all_best_mse["TRMM_LBA"]["updraft_area"] = 7561.609964448662
-all_best_mse["TRMM_LBA"]["updraft_w"] = 11162.65860693764
+all_best_mse["TRMM_LBA"]["updraft_area"] = 7.5616159743597864e+03
+all_best_mse["TRMM_LBA"]["updraft_w"] = 1.1162662767994134e+04
 all_best_mse["TRMM_LBA"]["updraft_qt"] = 178.62874777200045
 all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1972.5026068596098
 all_best_mse["TRMM_LBA"]["v_mean"] = 816.9725417750085
 all_best_mse["TRMM_LBA"]["u_mean"] = 996.5804518730336
-all_best_mse["TRMM_LBA"]["tke_mean"] = 8195.78791967079
+all_best_mse["TRMM_LBA"]["tke_mean"] = 8.1959925236028794e+03
 all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0011323111668659397
-all_best_mse["TRMM_LBA"]["ql_mean"] = 2623.5510869685795
+all_best_mse["TRMM_LBA"]["ql_mean"] = 2.6235652901417766e+03
 all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.010479008728209362
 all_best_mse["TRMM_LBA"]["Hvar_mean"] = 2511.796313996477
 all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2460.6541168787
 #
 all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict()
-all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.491835618015334
-all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7057083566611517
-all_best_mse["LES_driven_SCM"]["u_mean"] = 1.235432758754188
-all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0028509097874652232
-all_best_mse["LES_driven_SCM"]["ql_mean"] = 252.08833547882853
-all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0031436563856298365
+all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.6765992515646833e+00
+all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7142882513530031e+00
+all_best_mse["LES_driven_SCM"]["u_mean"] = 1.2391794278032586e+00
+all_best_mse["LES_driven_SCM"]["temperature_mean"] = 2.9024394724690302e-03
+all_best_mse["LES_driven_SCM"]["ql_mean"] = 1.1397368402146715e+03
+all_best_mse["LES_driven_SCM"]["thetal_mean"] = 3.2049121067961291e-03
 #
 #################################
 #################################

src/EDMF_Environment.jl

@@ -52,13 +52,15 @@ function update_cloud_dry(en_thermo::EnvironmentThermodynamics, state, k, ts)
     aux_en = center_aux_environment(state)
     if q_liq > 0.0
         aux_en.cloud_fraction[k] = 1.0
-        en_thermo.th_cloudy[k] = TD.liquid_ice_pottemp(ts)
+        en_thermo.th_cloudy[k] = TD.dry_pottemp(ts)
+        en_thermo.thl_cloudy[k] = TD.liquid_ice_pottemp(ts)
         en_thermo.t_cloudy[k] = TD.air_temperature(ts)
         en_thermo.qt_cloudy[k] = TD.total_specific_humidity(ts)
         en_thermo.qv_cloudy[k] = TD.vapor_specific_humidity(ts)
     else
         aux_en.cloud_fraction[k] = 0.0
-        en_thermo.th_dry[k] = TD.liquid_ice_pottemp(ts)
+        en_thermo.th_dry[k] = TD.dry_pottemp(ts)
+        en_thermo.thv_dry[k] = TD.virtual_pottemp(ts)
         en_thermo.qt_dry[k] = TD.total_specific_humidity(ts)
     end
     return
@@ -256,6 +258,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r
             en_thermo.qt_cloudy[k] = outer_env[i_qt_cld]
             ts_cld = TD.PhaseEquil_pTq(param_set, p0_c[k], en_thermo.t_cloudy[k], en_thermo.qt_cloudy[k])
             en_thermo.th_cloudy[k] = TD.dry_pottemp(ts_cld)
+            en_thermo.thl_cloudy[k] = TD.liquid_ice_pottemp(ts_cld)
 
             # update var/covar rain sources
             en_thermo.Hvar_rain_dt[k] = outer_src[i_SH_H] - outer_src[i_SH] * aux_en.θ_liq_ice[k]

src/Operators.jl

@@ -250,6 +250,54 @@ function c∇(f::SA.SVector, grid::Grid, ::BottomBCTag, ::Extrapolate)
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 
+# ∇(center data) for possibly vanishing subdomains.
+
+"""
+    c∇_vanishing_subdomain
+
+Used when the vertical gradient of a field must be computed over a conditional subdomain which may vanish
+below or above the current cell. If the subdomains vanishes, a default user-defined gradient is returned.
+
+Inputs:
+ - f_cut: Slice of field.
+ - sd_cut: Slice of subdomain volume fraction.
+ - default∇: Gradient used in vanishing subdomains.
+"""
+function c∇_vanishing_subdomain(
+    f_cut::SA.SVector,
+    sd_cut::SA.SVector,
+    default∇::FT,
+    grid::Grid,
+    k;
+    bottom = NoBCGivenError(),
+    top = NoBCGivenError(),
+) where {FT <: Real}
+    if is_surface_center(grid, k)
+        return c∇(f_cut, grid, BottomBCTag(), bottom)
+    elseif is_toa_center(grid, k)
+        return c∇(f_cut, grid, TopBCTag(), top)
+    else
+        return c∇_vanishing_subdomain(f_cut, sd_cut, default∇, grid, InteriorTag())
+    end
+end
+function c∇_vanishing_subdomain(
+    f::SA.SVector,
+    sd::SA.SVector,
+    default∇::FT,
+    grid::Grid,
+    ::InteriorTag,
+) where {FT <: Real}
+    @assert length(f) == 3
+    @assert length(sd) == 3
+    if sd[1] * sd[3] ≈ FT(0)
+        return default∇
+    else
+        f_dual⁺ = SA.SVector(f[2], f[3])
+        f_dual⁻ = SA.SVector(f[1], f[2])
+        return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
+    end
+end
+
 #####
 ##### ∇(face data)
 #####

src/closures/buoyancy_gradients.jl

@@ -1,36 +1,71 @@
-function buoyancy_gradients(param_set, bg_model::Tan2018)
+"""
+    buoyancy_gradients(param_set, bg_model::EnvBuoyGrad{FT, EBG}) where {FT <: Real, EBG <: EnvBuoyGradClosure}
+
+Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions.
+The dispatch on EnvBuoyGrad type is performed at the EnvBuoyGrad construction time, and the analytical solutions
+used here are consistent for both mean fields and conditional fields obtained from assumed distributions
+over the conserved thermodynamic variables.
+"""
+function buoyancy_gradients(param_set, bg_model::EnvBuoyGrad{FT, EBG}) where {FT <: Real, EBG <: EnvBuoyGradClosure}
 
     g = CPP.grav(param_set)
     molmass_ratio = CPP.molmass_ratio(param_set)
     R_d = CPP.R_d(param_set)
     R_v = CPP.R_v(param_set)
 
-    # buoyancy_gradients
     phase_part = TD.PhasePartition(0.0, 0.0, 0.0) # assuming R_d = R_m
     Π = TD.exner_given_pressure(param_set, bg_model.p0, phase_part)
 
     ∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p0) * Π
-    ∂b∂θl_dry = ∂b∂θv * (1 + (molmass_ratio - 1) * bg_model.qt_dry)
-    ∂b∂qt_dry = ∂b∂θv * bg_model.th_dry * (molmass_ratio - 1)
 
     if bg_model.en_cld_frac > 0.0
-        ts_cloudy = thermo_state_pθq(param_set, bg_model.p0, bg_model.th_cloudy, bg_model.qt_cloudy)
-        lh = TD.latent_heat_vapor(param_set, bg_model.t_cloudy)
+        ts_cloudy = thermo_state_pθq(param_set, bg_model.p0, bg_model.θ_liq_ice_cloudy, bg_model.qt_cloudy)
+        phase_part = TD.PhasePartition(ts_cloudy)
+        lh = TD.latent_heat_liq_ice(param_set, phase_part)
         cp_m = TD.cp_m(ts_cloudy)
         ∂b∂θl_cld = (
             ∂b∂θv * (1 + molmass_ratio * (1 + lh / R_v / bg_model.t_cloudy) * bg_model.qv_cloudy - bg_model.qt_cloudy) /
             (1 + lh * lh / cp_m / R_v / bg_model.t_cloudy / bg_model.t_cloudy * bg_model.qv_cloudy)
         )
-        ∂b∂qt_cld = (lh / cp_m / bg_model.t_cloudy * ∂b∂θl_cld - ∂b∂θv) * bg_model.th_cloudy
+        ∂b∂qt_cld = (lh / cp_m / bg_model.t_cloudy * ∂b∂θl_cld - ∂b∂θv) * bg_model.θ_cloudy
+    else
+        ∂b∂θl_cld = FT(0)
+        ∂b∂qt_cld = FT(0)
+    end
+
+    ∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy = buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_cld, ∂b∂qt_cld)
+    return GradBuoy(∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy)
+end
+
+"""
+    buoyancy_gradient_chain_rule(
+        bg_model::EnvBuoyGrad{FT, EBG},
+        ∂b∂θv::FT,
+        ∂b∂θl_cld::FT,
+        ∂b∂qt_cld::FT,
+    ) where {FT <: Real, EBG <: EnvBuoyGradClosure}
+
+Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions,
+from the partial derivatives with respect to thermodynamic variables in dry and cloudy volumes.
+"""
+function buoyancy_gradient_chain_rule(
+    bg_model::EnvBuoyGrad{FT, EBG},
+    ∂b∂θv::FT,
+    ∂b∂θl_cld::FT,
+    ∂b∂qt_cld::FT,
+) where {FT <: Real, EBG <: EnvBuoyGradClosure}
+    if bg_model.en_cld_frac > FT(0)
+        ∂b∂z_θl_cld = ∂b∂θl_cld * bg_model.∂θl∂z_cloudy
+        ∂b∂z_qt_cld = ∂b∂qt_cld * bg_model.∂qt∂z_cloudy
     else
-        ∂b∂θl_cld = 0
-        ∂b∂qt_cld = 0
+        ∂b∂z_θl_cld = FT(0)
+        ∂b∂z_qt_cld = FT(0)
     end
-    ∂b∂θl = (bg_model.en_cld_frac * ∂b∂θl_cld + (1 - bg_model.en_cld_frac) * ∂b∂θl_dry)
-    ∂b∂qt = (bg_model.en_cld_frac * ∂b∂qt_cld + (1 - bg_model.en_cld_frac) * ∂b∂qt_dry)
 
-    ∂b∂z_θl = ∂b∂θl * bg_model.∂θl∂z
-    ∂b∂z_qt = ∂b∂qt * bg_model.∂qt∂z
+    ∂b∂z_dry = bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_dry : FT(0)
+
+    ∂b∂z_cloudy = ∂b∂z_θl_cld + ∂b∂z_qt_cld
+    ∂b∂z = (1 - bg_model.en_cld_frac) * ∂b∂z_dry + bg_model.en_cld_frac * ∂b∂z_cloudy
 
-    return GradBuoy(∂b∂θl, ∂b∂qt, ∂b∂z_θl, ∂b∂z_qt)
+    return ∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy
 end