Merge pull request #703 from CliMA/kd/soil_freeze_thaw_timestepping_i…

…mprovements

Increasing phase change timescale for soil

src/shared_utilities/Domains.jl

@@ -857,20 +857,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
 
 """
@@ -912,7 +914,7 @@ 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)
@@ -923,6 +925,7 @@ function get_additional_domain_fields(subsurface_space)
         Δz_bottom = Δz_bottom,
         z_sfc = z_sfc,
         depth = d,
+        Δz = Δz,
     )
     return fields
 end

src/standalone/Soil/energy_hydrology.jl

@@ -120,7 +120,7 @@ end
     ) where {FT, D, PS}
 
 A constructor for a `EnergyHydrology` model, which sets the default value
-of the `lateral_flow` flag to true.
+of the `lateral_flow` flag to false.
 """
 function EnergyHydrology{FT}(;
     parameters::EnergyHydrologyParameters{FT, PSE},
@@ -332,7 +332,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,
                             ),
@@ -354,7 +354,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,
                         ),
@@ -577,7 +577,6 @@ PhaseChange source type.
 """
 struct PhaseChange{FT} <: AbstractSoilSource{FT} end
 
-
 """
      source!(dY::ClimaCore.Fields.FieldVector,
              src::PhaseChange{FT},
@@ -600,7 +599,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,
@@ -613,7 +612,7 @@ 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.κ,
             ),
             ν,
@@ -633,7 +632,7 @@ 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.κ,
             ),
             ν,

test/shared_utilities/domains.jl

@@ -43,7 +43,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
@@ -109,7 +109,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
@@ -260,7 +260,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