trying to get soil alone to be NaN free

experiments/long_runs/soil.jl

@@ -353,7 +353,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
 
     soil_args = (domain = domain, parameters = soil_params)
     soil_model_type = Soil.EnergyHydrology{FT}
-    sources = (Soil.PhaseChange{FT}(),)# sublimation and subsurface runoff are added automatically
+    sources = (Soil.PhaseChange{FT}(Δt),)# sublimation and subsurface runoff are added automatically
     top_bc = ClimaLand.Soil.AtmosDrivenFluxBC(atmos, radiation, runoff_model)
     zero_flux = Soil.HeatFluxBC((p, t) -> 0.0)
     boundary_conditions = (;

src/integrated/soil_canopy_model.jl

@@ -88,9 +88,9 @@ function SoilCanopyModel{FT}(;
     MM <: Soil.Biogeochemistry.SoilCO2Model{FT},
 }
 
-    (; atmos, radiation, soil_organic_carbon) = land_args
+    (; atmos, radiation, soil_organic_carbon, Δt) = land_args
     # These should always be set by the constructor.
-    sources = (RootExtraction{FT}(), Soil.PhaseChange{FT}())
+    sources = (RootExtraction{FT}(), Soil.PhaseChange{FT}(Δt))
     if :runoff ∈ propertynames(land_args)
         top_bc = ClimaLand.Soil.AtmosDrivenFluxBC(
             atmos,

src/shared_utilities/Domains.jl

@@ -753,20 +753,22 @@ bottom_center_to_surface(val) = val
 
 A function to return a tuple containing the distance between the top boundary
 and its closest center, and the bottom boundary and its closest center,
-both as Fields.
+both as Fields. It also returns the widths of each layer as a field.
 """
 function get_Δz(z::ClimaCore.Fields.Field)
     # Extract the differences between levels of the face space
     fs = obtain_face_space(axes(z))
     z_face = ClimaCore.Fields.coordinate_field(fs).z
-    Δz = ClimaCore.Fields.Δz_field(z_face)
-
+    Δz_face = ClimaCore.Fields.Δz_field(z_face)
     Δz_top = ClimaCore.Fields.level(
-        Δz,
+        Δz_face,
         ClimaCore.Utilities.PlusHalf(ClimaCore.Spaces.nlevels(fs) - 1),
     )
-    Δz_bottom = ClimaCore.Fields.level(Δz, ClimaCore.Utilities.PlusHalf(0))
-    return Δz_top ./ 2, Δz_bottom ./ 2
+    Δz_bottom = ClimaCore.Fields.level(Δz_face, ClimaCore.Utilities.PlusHalf(0))
+
+    #Layer widths:
+    Δz_center = ClimaCore.Fields.Δz_field(z)
+    return Δz_top ./ 2, Δz_bottom ./ 2, Δz_center
 end
 
 """
@@ -808,11 +810,17 @@ during the simulation.
 function get_additional_domain_fields(subsurface_space)
     surface_space = obtain_surface_space(subsurface_space)
     z = ClimaCore.Fields.coordinate_field(subsurface_space).z
-    Δz_top, Δz_bottom = get_Δz(z)
+    Δz_top, Δz_bottom, Δz = get_Δz(z)
     face_space = obtain_face_space(subsurface_space)
     z_face = ClimaCore.Fields.coordinate_field(face_space).z
     z_sfc = top_face_to_surface(z_face, surface_space)
-    fields = (; z = z, Δz_top = Δz_top, Δz_bottom = Δz_bottom, z_sfc = z_sfc)
+    fields = (;
+        z = z,
+        Δz_top = Δz_top,
+        Δz_bottom = Δz_bottom,
+        z_sfc = z_sfc,
+        Δz = Δz,
+    )
     return fields
 end
 

src/standalone/Soil/energy_hydrology.jl

@@ -331,7 +331,7 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT}
                             ClimaLand.Soil.dψdϑ(
                                 hydrology_cm,
                                 Y.soil.ϑ_l,
-                                ν,
+                                ν - Y.soil.θ_i, #ν_eff
                                 θ_r,
                                 S_s,
                             ),
@@ -353,7 +353,7 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT}
                         ClimaLand.Soil.dψdϑ(
                             hydrology_cm,
                             Y.soil.ϑ_l,
-                            ν,
+                            ν - Y.soil.θ_i, #ν_eff
                             θ_r,
                             S_s,
                         ),
@@ -572,9 +572,11 @@ end
 """
     PhaseChange{FT} <: AbstractSoilSource{FT}
 
-PhaseChange source type.
+PhaseChange source type
 """
-struct PhaseChange{FT} <: AbstractSoilSource{FT} end
+struct PhaseChange{FT} <: AbstractSoilSource{FT}
+    Δt::FT
+end
 
 
 """
@@ -599,7 +601,7 @@ function ClimaLand.source!(
     (; ν, ρc_ds, θ_r, hydrology_cm, earth_param_set) = params
     _ρ_l = FT(LP.ρ_cloud_liq(earth_param_set))
     _ρ_i = FT(LP.ρ_cloud_ice(earth_param_set))
-    Δz_top = model.domain.fields.Δz_top # center face distance
+    Δz = model.domain.fields.Δz # center face distance
     @. dY.soil.ϑ_l +=
         -phase_change_source(
             p.soil.θ_l,
@@ -612,8 +614,9 @@ function ClimaLand.source!(
                     ρc_ds,
                     earth_param_set,
                 ),
-                2 * Δz_top, # the factor of 2 appears to get the face-face/layer thickness, Δz_top is center-face distance
+                Δz,
                 p.soil.κ,
+                src.Δt,
             ),
             ν,
             θ_r,
@@ -632,8 +635,9 @@ function ClimaLand.source!(
                     ρc_ds,
                     earth_param_set,
                 ),
-                2 * Δz_top, #the factor of 2 appears to get the face-face/layer thickness, Δz_top is center-face distance
+                Δz,
                 p.soil.κ,
+                src.Δt,
             ),
             ν,
             θ_r,

src/standalone/Soil/soil_heat_parameterizations.jl

@@ -18,9 +18,12 @@ export volumetric_internal_energy,
 
 Returns the thermal timescale for temperature differences across
 a typical thickness Δz to equilibrate.
+
+Clip to 10x the timestep in order to respect CFL condition,
+given that phase change is stepped explicitly.
 """
-function thermal_time(ρc::FT, Δz::FT, κ::FT) where {FT}
-    return ρc * Δz^2 / κ
+function thermal_time(ρc::FT, Δz::FT, κ::FT, Δt::FT) where {FT}
+    return max(ρc * Δz^2 / κ, 10 * Δt)
 end
 
 """
@@ -209,7 +212,7 @@ Compute the expression for relative saturation.
 This is referred to as θ_sat in Balland and Arp's paper.
 """
 function relative_saturation(θ_l::FT, θ_i::FT, ν::FT) where {FT}
-    return (θ_l + θ_i) / ν
+    return max((θ_l + θ_i) / ν, eps(FT))
 end
 
 """

src/standalone/Soil/soil_hydrology_parameterizations.jl

@@ -113,13 +113,14 @@ function pressure_head(
 end
 
 """
-   dψdϑ(cm::vanGenuchten{FT}, ϑ, ν, θ_r, S_s)
+   dψdϑ(cm::vanGenuchten{FT}, ϑ_l, ν, θ_r, S_s)
 
 Computes and returns the derivative of the pressure head 
 with respect to ϑ for the van Genuchten formulation.
 """
-function dψdϑ(cm::vanGenuchten{FT}, ϑ, ν, θ_r, S_s) where {FT}
-    S = effective_saturation(ν, ϑ, θ_r)
+function dψdϑ(cm::vanGenuchten{FT}, ϑ_l, ν, θ_r, S_s) where {FT}
+    ϑ_l_safe = max(ϑ_l, θ_r + eps(FT))
+    S = effective_saturation(ν, ϑ_l_safe, θ_r)
     (; α, m, n) = cm
     if S < 1.0
         return FT(

test/shared_utilities/domains.jl

@@ -42,7 +42,7 @@ FT = Float32
     face_space = obtain_face_space(shell.space.subsurface)
     z_face = ClimaCore.Fields.coordinate_field(face_space).z
     @test shell.fields.z_sfc == top_face_to_surface(z_face, shell.space.surface)
-    Δz_top, Δz_bottom = get_Δz(shell.fields.z)
+    Δz_top, Δz_bottom, Δz = get_Δz(shell.fields.z)
     @test shell.fields.Δz_top == Δz_top
     @test shell.fields.Δz_bottom == Δz_bottom
     @test shell.radius == radius
@@ -107,7 +107,7 @@ FT = Float32
     face_space = obtain_face_space(box.space.subsurface)
     z_face = ClimaCore.Fields.coordinate_field(face_space).z
     @test box.fields.z_sfc == top_face_to_surface(z_face, box.space.surface)
-    Δz_top, Δz_bottom = get_Δz(box.fields.z)
+    Δz_top, Δz_bottom, Δz = get_Δz(box.fields.z)
     @test box.fields.Δz_top == Δz_top
     @test box.fields.Δz_bottom == Δz_bottom
     box_coords = coordinates(box).subsurface
@@ -178,7 +178,9 @@ FT = Float32
     z_face = ClimaCore.Fields.coordinate_field(face_space).z
     @test z_column.fields.z_sfc ==
           top_face_to_surface(z_face, z_column.space.surface)
-    Δz_top, Δz_bottom = get_Δz(z_column.fields.z)
+    Δz_top, Δz_bottom, Δz = get_Δz(z_column.fields.z)
+    z = ClimaCore.Fields.coordinate_field(z_column.space.subsurface).z
+    @test z_column.fields.z == z
     @test z_column.fields.Δz_top == Δz_top
     @test z_column.fields.Δz_bottom == Δz_bottom
     column_coords = coordinates(z_column).subsurface

test/standalone/Soil/soiltest.jl

@@ -663,8 +663,8 @@ for FT in (Float32, Float64)
                 heat = zero_heat_flux_bc,
             ),
         )
-
-        sources = (PhaseChange{FT}(),)
+        Δt = FT(1000)# doesnt matter what this is here, we are checking dY.liquid water and dY.ice are in sync
+        sources = (PhaseChange{FT}(Δt),)
 
         ###
         hyd_off_en_on = Soil.EnergyHydrologyParameters(