allocate space in p

src/shared_utilities/models.jl

@@ -164,7 +164,7 @@ end
      make_update_cache(model::AbstractModel)
 
 A helper function which updates all cache variables of a model;
-currently only used in `set_initial_cache` since not all 
+currently only used in `set_initial_cache` since not all
 cache variables are updated at the same time.
 """
 function make_update_cache(model::AbstractModel)
@@ -404,6 +404,15 @@ function get_drivers(model::AbstractModel)
     return (nothing, nothing)
 end
 
+"""
+    add_dtopBCdY_to_cache(p::NamedTuple, model::AbstractModel)
+
+We only need to add this field in the case of a `RichardsModel`, so we
+provide a fallback method for all other models which does nothing.
+"""
+function add_dtopBCdY_to_cache(p::NamedTuple, model::AbstractModel)
+    return p
+end
 
 """
     initialize(model::AbstractModel)
@@ -418,5 +427,6 @@ function initialize(model::AbstractModel{FT}) where {FT}
     Y = initialize_prognostic(model, coords)
     p = initialize_auxiliary(model, coords)
     p = add_drivers_to_cache(p, model, coords)
+    p = add_dtopBCdY_to_cache(p, model)
     return Y, p, coords
 end

src/standalone/Soil/rre.jl

@@ -478,9 +478,25 @@ end
 """
     is_saturated(Y, model::RichardsModel)
 
-A helper function which can be used to indicate whether a layer of soil is 
+A helper function which can be used to indicate whether a layer of soil is
 saturated.
 """
 function is_saturated(Y, model::RichardsModel)
     return @. ClimaLand.heaviside(Y.soil.ϑ_l - model.parameters.ν)
 end
+
+
+"""
+    add_dtopBCdY_to_cache(p::NamedTuple, model::RichardsModel)
+
+Allocate space in `p` for the derivative of the top boundary condition with
+respect to the state variables, and merges it into `p` under the key `dtopBCdY`.
+This is used in the implicit solver for
+`RichardsModel`, where we need to handle the boundary specially in computing
+the Jacobian.
+"""
+function ClimaLand.add_dtopBCdY_to_cache(p::NamedTuple, model::RichardsModel)
+    surface_space = model.domain.space.surface
+    dtopBCdY_field = ClimaCore.Fields.zeros(surface_space)
+    return merge(p, (; dtopBCdtheta = dtopBCdY_field))
+end