add snow cover fraction and integrated water to compare in ILAMB (#1087)

Author: kmdeck

experiments/long_runs/land_region.jl

@@ -443,7 +443,7 @@ setup_and_solve_problem(; greet = true);
 # read in diagnostics and make some plots!
 #### ClimaAnalysis ####
 simdir = ClimaAnalysis.SimDir(outdir)
-short_names = ["gpp", "ct", "swc", "si"]
+short_names = ["gpp", "swc", "si"]
 mktempdir(root_path) do tmpdir
     for short_name in short_names
         var = get(simdir; short_name)

experiments/long_runs/snowy_land.jl

@@ -168,10 +168,10 @@ function setup_prob(
     SAI = FT(0.0) # m2/m2
     f_root_to_shoot = FT(3.5)
     RAI = FT(1.0)
-    K_sat_plant = FT(5e-9) # m/s # seems much too small?
-    ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa
-    Weibull_param = FT(4) # unitless, Holtzman's original c param value
-    a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
+    K_sat_plant = FT(7e-8) # m/s 
+    ψ63 = FT(-4 / 0.0098) # / MPa to m
+    Weibull_param = FT(4) # unitless
+    a = FT(0.2 * 0.0098) # 1/m
     conductivity_model =
         Canopy.PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
     retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a)
@@ -494,6 +494,7 @@ short_names = [
     "iwc",
     "trans",
     "msf",
+    "snowc",
 ]
 
 include(

src/diagnostics/Diagnostics.jl

@@ -4,6 +4,8 @@ import Dates: Month, Period
 
 import ClimaComms
 
+import ..Parameters as LP
+
 using ..Bucket: BucketModel
 
 import ..SoilCanopyModel

src/diagnostics/default_diagnostics.jl

@@ -322,11 +322,11 @@ function default_diagnostics(
             "shf",
             "ghf",
             "iwc",
+            "snowc",
         ]
     elseif output_vars == :short
         snowyland_diagnostics = [
             "gpp",
-            "ct",
             "swc",
             "si",
             "sie",
@@ -344,6 +344,7 @@ function default_diagnostics(
             "iwc",
             "swd",
             "lwd",
+            "snowc",
         ]
     end
 

src/diagnostics/define_diagnostics.jl

@@ -820,12 +820,12 @@ function define_diagnostics!(land_model)
 
     add_diagnostic_variable!(
         short_name = "iwc",
-        long_name = "Integrated Soil Water Content in first 1m",
-        standard_name = "soil_1m_water_content",
-        units = "m",
-        comments = "The integrated water content to a depth of 1m",
+        long_name = "Integrated Soil Water Mass in first 10cm",
+        standard_name = "soil_10cm_water_mass",
+        units = "kg/m^2",
+        comments = "The integrated water mass to a depth of 10cm",
         compute! = (out, Y, p, t) ->
-            compute_1m_water_content!(out, Y, p, t, land_model),
+            compute_10cm_water_mass!(out, Y, p, t, land_model),
     )
 
     # Plant water content
@@ -887,4 +887,14 @@ function define_diagnostics!(land_model)
             compute_snow_depth!(out, Y, p, t, land_model),
     )
 
+    # Snow cover fraction
+    add_diagnostic_variable!(
+        short_name = "snowc",
+        long_name = "Snow cover fraction",
+        standard_name = "snow_cover_fraction",
+        units = "",
+        comments = "The snow cover fraction",
+        compute! = (out, Y, p, t) ->
+            compute_snow_cover_fraction!(out, Y, p, t, land_model),
+    )
 end

src/diagnostics/land_compute_methods.jl

@@ -202,7 +202,7 @@ end
 @diagnostic_compute "soil_net_radiation" Union{SoilCanopyModel, LandModel} p.soil.R_n
 @diagnostic_compute "soil_temperature" Union{SoilCanopyModel, LandModel} p.soil.T
 
-function compute_1m_water_content!(
+function compute_10cm_water_mass!(
     out,
     Y,
     p,
@@ -212,13 +212,26 @@ function compute_1m_water_content!(
     ∫Hθdz = p.soil.sfc_scratch
     Hθ = p.soil.sub_sfc_scratch
     z = land_model.soil.domain.fields.z
-    depth = FT(-1)
-    @. Hθ = (Y.soil.ϑ_l + Y.soil.θ_i) * heaviside(z, depth)
+    depth = FT(-0.1)
+    earth_param_set = land_model.soil.parameters.earth_param_set
+    _ρ_liq = LP.ρ_cloud_liq(earth_param_set)
+    _ρ_ice = LP.ρ_cloud_ice(earth_param_set)
+    # Convert from volumetric water content to water mass per unit volume using density
+    @. Hθ = (Y.soil.ϑ_l * _ρ_liq + Y.soil.θ_i * _ρ_ice) * heaviside(z, depth)
     column_integral_definite!(∫Hθdz, Hθ)
+
+    # The layering of the soil model may not coincide with 10 cm exactly, and this could lead
+    # to the integral above not exactly representing 10cm.
+    # To adjust, divide by the ∫heaviside(z, depth) dz, and then multiply by 10cm
+    H = p.subsfc_scratch
+    @. H = heaviside(z, depth)
+    ∫Hdz = p.sfc_scratch
+    column_integral_definite!(∫Hdz, H)
+
     if isnothing(out)
-        return ∫Hθdz
+        return ∫Hθdz ./ ∫Hdz .* FT(0.1)
     else
-        out .= ∫Hθdz
+        @. out = ∫Hθdz / ∫Hdz * FT(0.1)
     end
 end
 function compute_soil_albedo!(
@@ -464,6 +477,7 @@ end
 
 @diagnostic_compute "snow_water_equivalent" LandModel Y.snow.S
 @diagnostic_compute "snow_depth" LandModel p.snow.z_snow
+@diagnostic_compute "snow_cover_fraction" LandModel p.snow.snow_cover_fraction
 
 ### EnergyHydrology ###