preallocate q_sfc for snow (#978)

Author: kmdeck

src/standalone/Snow/Snow.jl

@@ -246,7 +246,8 @@ density_prog_names(::AbstractDensityModel) = ()
     auxiliary_vars(::SnowModel)
 
 Returns the auxiliary variable names for the snow model. These
-include the mass fraction in liquid water (`q_l`, unitless),
+include the specific humidity at the surface of the snow `(`q_sfc`, unitless),
+the mass fraction in liquid water (`q_l`, unitless),
 the thermal conductivity (`κ`, W/m/K),
 the bulk temperature (`T`, K), the surface temperature (`T_sfc`, K), the bulk snow density (`ρ_snow`, kg/m^3)
 the SHF, LHF, and vapor flux (`turbulent_fluxes.shf`, etc),
@@ -259,6 +260,7 @@ clipped values are what are actually applied as boundary fluxes, and are stored
 `applied_` fluxes.
 """
 auxiliary_vars(::SnowModel) = (
+    :q_sfc,
     :q_l,
     :κ,
     :T,
@@ -281,6 +283,7 @@ auxiliary_types(::SnowModel{FT}) where {FT} = (
     FT,
     FT,
     FT,
+    FT,
     NamedTuple{(:lhf, :shf, :vapor_flux, :r_ae), Tuple{FT, FT, FT, FT}},
     FT,
     FT,
@@ -307,6 +310,7 @@ auxiliary_domain_names(::SnowModel) = (
     :surface,
     :surface,
     :surface,
+    :surface,
 )
 
 
@@ -328,6 +332,13 @@ function ClimaLand.make_update_aux(model::SnowModel{FT}) where {FT}
 
         @. p.snow.T_sfc = snow_surface_temperature(p.snow.T)
 
+        @. p.snow.q_sfc = snow_surface_specific_humidity(
+            p.snow.T_sfc,
+            p.snow.q_l,
+            p.drivers.thermal_state,
+            parameters,
+        )
+
         p.snow.water_runoff .=
             compute_water_runoff.(
                 Y.snow.S,

src/standalone/Snow/snow_parameterizations.jl

@@ -91,39 +91,46 @@ end
 
 
 """
-    ClimaLand.surface_specific_humidity(model::SnowModel, Y, p)
+    ClimaLand.surface_specific_humidity(model::SnowModel, Y, p, _...)
 
-Computes and returns the specific humidity over snow as a weighted
+Returns the precomputed specific humidity over snow as a weighted
 fraction of the saturated specific humidity over liquid and frozen
 water.
 
-This function allocates.
 """
-function ClimaLand.surface_specific_humidity(
-    model::SnowModel,
-    Y,
-    p,
-    T_sfc,
-    ρ_sfc,
-)
-    thermo_params =
-        LP.thermodynamic_parameters(model.parameters.earth_param_set)
-    qsat_over_ice =
-        Thermodynamics.q_vap_saturation_generic.(
-            Ref(thermo_params),
-            T_sfc,
-            ρ_sfc,
-            Ref(Thermodynamics.Ice()),
-        )
-    qsat_over_liq =
-        Thermodynamics.q_vap_saturation_generic.(
-            Ref(thermo_params),
-            T_sfc,
-            ρ_sfc,
-            Ref(Thermodynamics.Liquid()),
-        )
-    q_l = p.snow.q_l
-    return @. qsat_over_ice * (1 - q_l) + q_l * (qsat_over_liq)
+function ClimaLand.surface_specific_humidity(model::SnowModel, Y, p, _...)
+    return p.snow.q_sfc
+end
+
+"""
+    snow_surface_specific_humidity(T_sfc::FT, q_l::FT, atmos_ts, parameters) where {FT}
+
+Computes the snow surface specific humidity at a point, assuming a weighted averaged (by mass fraction)
+of the saturated specific humidity over ice and over liquid, at temperature T_sfc.
+
+This approximates the surface air density using an adiabatic approximation and the current atmospheric state.
+"""
+function snow_surface_specific_humidity(
+    T_sfc::FT,
+    q_l::FT,
+    atmos_ts,
+    parameters,
+) where {FT}
+    thermo_params = LP.thermodynamic_parameters(parameters.earth_param_set)
+    ρ_sfc = compute_ρ_sfc(thermo_params, atmos_ts, T_sfc)
+    qsat_over_ice = Thermodynamics.q_vap_saturation_generic(
+        thermo_params,
+        T_sfc,
+        ρ_sfc,
+        Thermodynamics.Ice(),
+    )
+    qsat_over_liq = Thermodynamics.q_vap_saturation_generic(
+        thermo_params,
+        T_sfc,
+        ρ_sfc,
+        Thermodynamics.Liquid(),
+    )
+    return qsat_over_ice * (1 - q_l) + q_l * (qsat_over_liq)
 end
 
 """

test/standalone/Snow/snow.jl

@@ -53,6 +53,7 @@ import ClimaLand.Parameters as LP
     Y, p, coords = ClimaLand.initialize(model)
     @test (Y.snow |> propertynames) == (:S, :U)
     @test (p.snow |> propertynames) == (
+        :q_sfc,
         :q_l,
         :κ,
         :T,
@@ -124,6 +125,7 @@ import ClimaLand.Parameters as LP
             ρ_sfc,
             Ref(Thermodynamics.Ice()),
         )
+    @test p.snow.q_sfc ≈ q_sfc
     turb_fluxes_copy = copy(p.snow.turbulent_fluxes)
     ClimaLand.turbulent_fluxes!(
         turb_fluxes_copy,