Merge pull request #770 from CliMA/kd/investigate_NaNs_sahara

LAI lower threshold for nonzero LAI
CliMA · Sep 17, 2024 · 632d83d · 632d83d
2 parents b47ef0f + ffc37e2
commit 632d83d
Show file tree

Hide file tree

Showing 6 changed files with 21 additions and 35 deletions.
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
@@ -473,21 +473,12 @@ function setup_prob(t0, tf, Δt; nelements = (101, 15))
         )
     )
     # Set up plant hydraulics
-    # Note that we clip all values of LAI below 0.05 to zero.
-    # This is because we currently run into issues when LAI is
-    # of order eps(FT) in the SW radiation code.
-    # Please see Issue #644
-    # or PR #645 for details.
-    # For now, this clipping is similar to what CLM does.
     LAIfunction = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
         "lai",
         surface_space;
         reference_date = start_date,
         regridder_type,
-        file_reader_kwargs = (;
-            preprocess_func = (data) -> data > 0.05 ? data : 0.0,
-        ),
     )
     ai_parameterization =
         Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)

diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl
@@ -226,16 +226,12 @@ photosynthesis_args = (;
     parameters = Canopy.FarquharParameters(FT, Canopy.C3(); Vcmax25 = Vcmax25)
 )
 # Set up plant hydraulics
-# Not ideal
 LAIfunction = TimeVaryingInput(
     joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
     "lai",
     surface_space;
     reference_date = start_date,
     regridder_type,
-    file_reader_kwargs = (;
-        preprocess_func = (data) -> data > 0.05 ? data : 0.0,
-    ),
 )
 ai_parameterization = Canopy.PrescribedSiteAreaIndex{FT}(LAIfunction, SAI, RAI)
 

diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
@@ -68,7 +68,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         depth = depth,
         nelements = nelements,
         npolynomial = 1,
-        dz_tuple = FT.((10.0, 0.05)),# top layer should ideally be only a few cm!
+        dz_tuple = FT.((10.0, 0.05)),
     )
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
@@ -483,21 +483,12 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
     )
     # Set up plant hydraulics
 
-    # Note that we clip all values of LAI below 0.05 to zero.
-    # This is because we currently run into issues when LAI is
-    # of order eps(FT) in the SW radiation code.
-    # Please see Issue #644
-    # or PR #645 for details.
-    # For now, this clipping is similar to what CLM does.
     LAIfunction = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
         "lai",
         surface_space;
         reference_date = start_date,
         regridder_type,
-        file_reader_kwargs = (;
-            preprocess_func = (data) -> data > 0.05 ? data : 0.0,
-        ),
         method = time_interpolation_method,
     )
     ai_parameterization =

diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl
@@ -483,21 +483,12 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (10, 10, 15))
     )
     # Set up plant hydraulics
 
-    # Note that we clip all values of LAI below 0.05 to zero.
-    # This is because we currently run into issues when LAI is
-    # of order eps(FT) in the SW radiation code.
-    # Please see Issue #644
-    # or PR #645 for details.
-    # For now, this clipping is similar to what CLM does.
     LAIfunction = TimeVaryingInput(
         joinpath(era5_artifact_path, "era5_lai_2021_0.9x1.25_clima.nc"),
         "lai",
         surface_space;
         reference_date = start_date,
         regridder_type,
-        file_reader_kwargs = (;
-            preprocess_func = (data) -> data > 0.05 ? data : 0.0,
-        ),
         method = time_interpolation_method,
     )
     ai_parameterization =

diff --git a/src/standalone/Vegetation/PlantHydraulics.jl b/src/standalone/Vegetation/PlantHydraulics.jl
@@ -304,7 +304,9 @@ ClimaLand.auxiliary_types(model::PlantHydraulicsModel{FT}) where {FT} = (
 ClimaLand.auxiliary_domain_names(::PlantHydraulicsModel) =
     (:surface, :surface, :surface, :surface, :surface)
 
-
+function clip(x::FT, threshold::FT) where {FT}
+    x > threshold ? x : FT(0)
+end
 """
     set_canopy_prescribed_field!(component::PlantHydraulics{FT},
                                  p,
@@ -314,6 +316,13 @@ ClimaLand.auxiliary_domain_names(::PlantHydraulicsModel) =
 
 Sets the canopy prescribed fields pertaining to the PlantHydraulics
 component (the area indices) with their values at time t.
+
+Note that we clip all values of LAI below 0.05 to zero.
+This is because we currently run into issues when LAI is
+of order eps(FT) in the SW radiation code.
+Please see Issue #644
+or PR #645 for details.
+For now, this clipping is similar to what CLM and NOAH MP do.
 """
 function ClimaLand.Canopy.set_canopy_prescribed_field!(
     component::PlantHydraulicsModel{FT},
@@ -322,6 +331,8 @@ function ClimaLand.Canopy.set_canopy_prescribed_field!(
 ) where {FT}
     (; LAIfunction, SAI, RAI) = component.parameters.ai_parameterization
     evaluate!(p.canopy.hydraulics.area_index.leaf, LAIfunction, floor(t))
+    p.canopy.hydraulics.area_index.leaf .=
+        clip.(p.canopy.hydraulics.area_index.leaf, FT(0.05))
     @. p.canopy.hydraulics.area_index.stem = SAI
     @. p.canopy.hydraulics.area_index.root = RAI
 end

diff --git a/test/standalone/Vegetation/canopy_model.jl b/test/standalone/Vegetation/canopy_model.jl
@@ -498,15 +498,21 @@ end
         set_canopy_prescribed_field!(plant_hydraulics, p, FT(200))
         @test all(
             Array(parent(p.canopy.hydraulics.area_index.leaf)) .==
-            FT(LAI * sin(200 * 2π / 365)),
+            ClimaLand.Canopy.PlantHydraulics.clip(
+                FT(LAI * sin(200 * 2π / 365)),
+                FT(0.05),
+            ),
         )
 
         set_canopy_prescribed_field!(Default{FT}(), p, t0)
         set_canopy_prescribed_field!(Default{FT}(), p, t0)
         # Test that they are unchanged
         @test all(
             parent(p.canopy.hydraulics.area_index.leaf) .==
-            FT(LAI * sin(200 * 2π / 365)),
+            ClimaLand.Canopy.PlantHydraulics.clip(
+                FT(LAI * sin(200 * 2π / 365)),
+                FT(0.05),
+            ),
         )
         @test all(
             Array(parent(p.canopy.hydraulics.area_index.stem)) .== FT(1.0),