clipping sea ice tendencies

experiments/ClimaEarth/components/ocean/prescr_seaice.jl

@@ -14,7 +14,7 @@ include("../slab_utils.jl")
 
 Ice concentration is prescribed, and we solve the following energy equation:
 
-    (h * ρ * c) d T_sfc dt = -(F_turb_energy + F_radiativead) + F_conductive
+    (h * ρ * c) d T_sfc dt = -(F_turb_energy + F_radiative) + F_conductive
 
     with
     F_conductive = k_ice (T_base - T_sfc) / (h)
@@ -24,9 +24,7 @@ Ice concentration is prescribed, and we solve the following energy equation:
     as well as a conductive flux that depends on the temperature difference
     across the ice layer (with `T_base` being prescribed).
 
-In the current version, the sea ice has a prescribed thickness, and we assume that it is not
-sublimating. That contribution has been zeroed out in the atmos fluxes.
-
+In the current version, the sea ice has a prescribed thickness.
 """
 struct PrescribedIceSimulation{P, Y, D, I} <: Interfacer.SeaIceModelSimulation
     params::P
@@ -174,12 +172,11 @@ function ice_rhs!(du, u, p, _)
 
     F_conductive = @. params.k_ice / (params.h) * (params.T_base - Y.T_sfc) # fluxes are defined to be positive when upward
     rhs = @. (-F_turb_energy - F_radiative + F_conductive) / (params.h * params.ρ * params.c)
-
-    # do not count tendencies that lead to temperatures above freezing, and mask out no-ice areas
+    # If tendencies lead to temperature above freezing, set temperature to freezing
+    @. rhs = min(rhs, (T_freeze - Y.T_sfc) / p.dt)
+    # mask out no-ice areas
     area_mask = Regridder.binary_mask.(area_fraction)
-    threshold = zero(FT)
-    unphysical = @. Regridder.binary_mask.(T_freeze - (Y.T_sfc + FT(rhs) * FT(p.dt)), threshold) .* area_mask
-    parent(dY.T_sfc) .= parent(rhs .* unphysical)
+    dY.T_sfc .= rhs .* area_mask
 
     @. p.q_sfc = TD.q_vap_saturation_generic.(p.thermo_params, Y.T_sfc, p.ρ_sfc, TD.Ice())
 end

test/component_model_tests/prescr_seaice_tests.jl

@@ -48,9 +48,10 @@ for FT in (Float32, Float64)
         dT_expected = -1.0 / (p.params.h * p.params.ρ * p.params.c)
         @test sum([i for i in extrema(dY)] .≈ [FT(dT_expected), FT(dT_expected)]) == 2
 
-        # check that tendency not added if T of ice would have done above freezing
+        # check that tendency will not result in above freezing T
         dY, Y, p = test_sea_ice_rhs(F_radiative = 0.0, T_base = 330.0) # Float32 requires a large number here!
-        @test sum([i for i in extrema(dY)] .≈ [FT(0.0), FT(0.0)]) == 2
+        dT_maximum = @. (p.params.T_freeze - Y.T_sfc) / p.dt
+        @test minimum(dT_maximum .- dY.T_sfc) >= FT(0.0)
 
         # check that the correct tendency was added due to basal flux
         dY, Y, p = test_sea_ice_rhs(F_radiative = 0.0, T_base = 269.2, global_mask = 1.0)