cleanup and fix division by zero

Author: sajjadazimi

src/prognostic_equations/edmfx_boundary_condition.jl

@@ -50,6 +50,12 @@ function sgs_scalar_flux_bc(
     sfc_local_geometry,
 ) where {FT}
 
+    # when surface-interior difference on the grid mean is zero (negligible), 
+    # return grid mean flux (which is zero or negligible) to avoid division by zero
+    if abs(χ_sfc - ᶜχ_int) < sqrt(floatmin(FT))
+        return ρ_flux_χ
+    end
+
     kinematic_sfc_flux_χ = ρ_flux_χ / ρ_sfc # Convert to kinematic flux [K m/s]
 
     χ_sfc_var = get_scalar_variance(

src/prognostic_equations/surface_flux.jl

@@ -170,29 +170,32 @@ function edmfx_surface_flux_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticE
     thermo_params = CAP.thermodynamics_params(params)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
 
-    ρ_flux_χʲ = p.scratch.sfc_temp_C3
-    ᶜaʲ = p.scratch.ᶜtemp_scalar
-
-    (; ρ_flux_h_tot, ρ_flux_q_tot, ustar, obukhov_length, ts) =
+    (; ρ_flux_h_tot, ustar, obukhov_length, ts) =
         p.precomputed.sfc_conditions
 
-    ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
-    ρ_flux_q_tot_val = Fields.field_values(ρ_flux_q_tot)
+    ᶜaʲ = p.scratch.ᶜtemp_scalar
+    ρ_flux_χʲ = p.scratch.sfc_temp_C3
+    # We need field_values everywhere because we are mixing
+    # information from surface and first interior inside the
+    # sgs_scalar_flux_bc call.
+    ρ_flux_χʲ_val = Fields.field_values(ρ_flux_χʲ)
+
     ustar_val = Fields.field_values(ustar)
     obukhov_length_val = Fields.field_values(obukhov_length)
     sfc_local_geometry_val = Fields.field_values(
         Fields.local_geometry_field(Fields.level(Y.f, Fields.half)),
     )
 
-    ρ_sfc = @. lazy(TD.air_density(thermo_params, ts))
-    h_sfc = @. lazy(TD.specific_enthalpy(thermo_params, ts))
-    q_tot_sfc = @. lazy(TD.total_specific_humidity(thermo_params, ts))
     z_sfc_val =
         Fields.field_values(Fields.level(Fields.coordinate_field(Y.f).z, Fields.half))
+    ρ_sfc = @. lazy(TD.air_density(thermo_params, ts))
     ρ_sfc_val = Fields.field_values(ρ_sfc)
-    h_sfc_val = Fields.field_values(h_sfc)
-    q_tot_sfc_val = Fields.field_values(q_tot_sfc)
 
+    # Based on boundary conditions for updrafts we compute
+    # the tendency due to the surface flux for EDMFX ᶜmseʲ...
+    ρ_flux_h_tot_val = Fields.field_values(ρ_flux_h_tot)
+    h_tot_sfc = @. lazy(TD.specific_enthalpy(thermo_params, ts))
+    h_tot_sfc_val = Fields.field_values(h_tot_sfc)
     ᶜh_tot = @. lazy(
         TD.total_specific_enthalpy(
             thermo_params,
@@ -202,20 +205,16 @@ function edmfx_surface_flux_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticE
     )
     ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1))
     ᶜK_int_val = Fields.field_values(Fields.level(ᶜK, 1))
-    ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
-    ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
-
-    ρ_flux_χʲ_val = Fields.field_values(ρ_flux_χʲ)
 
     for j in 1:n
         @. ᶜaʲ = draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
         ᶜaʲ_int_val = Fields.field_values(Fields.level(ᶜaʲ, 1))
-
         ᶜmseʲ_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).mse, 1))
+
         @. ρ_flux_χʲ_val = sgs_scalar_flux_bc(
             z_sfc_val,
             ρ_sfc_val,
-            h_sfc_val,
+            h_tot_sfc_val,
             ᶜh_tot_int_val - ᶜK_int_val,
             ᶜaʲ_int_val,
             ᶜmseʲ_int_val,
@@ -226,22 +225,36 @@ function edmfx_surface_flux_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticE
         )
         btt = boundary_tendency_scalar(Y.c.sgsʲs.:(1).mse, ρ_flux_χʲ)
         @. Yₜ.c.sgsʲs.:($$j).mse -= btt / p.precomputed.ᶜρʲs.:($$j)
+    end
 
-        ᶜq_totʲ_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).q_tot, 1))
-        @. ρ_flux_χʲ_val = sgs_scalar_flux_bc(
-            z_sfc_val,
-            ρ_sfc_val,
-            q_tot_sfc_val,
-            ᶜq_tot_int_val,
-            ᶜaʲ_int_val,
-            ᶜq_totʲ_int_val,
-            ρ_flux_q_tot_val,
-            ustar_val,
-            obukhov_length_val,
-            sfc_local_geometry_val,
-        )
-        btt = boundary_tendency_scalar(Y.c.sgsʲs.:(1).q_tot, ρ_flux_χʲ)
-        @. Yₜ.c.sgsʲs.:($$j).q_tot -= btt / p.precomputed.ᶜρʲs.:($$j)
+    # ... and the tendency for EDMFX ᶜq_totʲ.
+    if !(p.atmos.moisture_model isa DryModel)
+        ρ_flux_q_tot_val = Fields.field_values(p.precomputed.sfc_conditions.ρ_flux_q_tot)
+        q_tot_sfc = @. lazy(TD.total_specific_humidity(thermo_params, ts))
+        q_tot_sfc_val = Fields.field_values(q_tot_sfc)
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+        ᶜq_tot_int_val = Fields.field_values(Fields.level(ᶜq_tot, 1))
+
+        for j in 1:n
+            @. ᶜaʲ = draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
+            ᶜaʲ_int_val = Fields.field_values(Fields.level(ᶜaʲ, 1))
+            ᶜq_totʲ_int_val = Fields.field_values(Fields.level(Y.c.sgsʲs.:($j).q_tot, 1))
+
+            @. ρ_flux_χʲ_val = sgs_scalar_flux_bc(
+                z_sfc_val,
+                ρ_sfc_val,
+                q_tot_sfc_val,
+                ᶜq_tot_int_val,
+                ᶜaʲ_int_val,
+                ᶜq_totʲ_int_val,
+                ρ_flux_q_tot_val,
+                ustar_val,
+                obukhov_length_val,
+                sfc_local_geometry_val,
+            )
+            btt = boundary_tendency_scalar(Y.c.sgsʲs.:(1).q_tot, ρ_flux_χʲ)
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -= btt / p.precomputed.ᶜρʲs.:($$j)
+        end
     end
 
 end