consolidate sea ice temperature

Author: simone-silvestri

src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl

@@ -14,6 +14,7 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model)
                       v = ZeroField(),
                       h = sea_ice.model.ice_thickness,
                       ℵ = sea_ice.model.ice_concentration,
+                      hc = sea_ice.model.ice_consolidation_thickness,
                       Tₒ = ocean.model.tracers.T,
                       Sₒ = ocean.model.tracers.S)
 
@@ -97,6 +98,7 @@ end
         uᵢ = zero(FT) # ℑxᶜᵃᵃ(i, j, 1, grid, interior_state.u)
         vᵢ = zero(FT) # ℑyᵃᶜᵃ(i, j, 1, grid, interior_state.v)
         hᵢ = interior_state.h[i, j, 1]
+        hc = interior_state.hc[i, j, 1]
         ℵᵢ = interior_state.ℵ[i, j, 1]
         Tₛ = interface_temperature[i, j, 1]
         Tₛ = convert_to_kelvin(sea_ice_properties.temperature_units, Tₛ)
@@ -117,7 +119,7 @@ end
                               h_bℓ = atmosphere_state.h_bℓ)
 
     downwelling_radiation = (; Qs, Qℓ)
-    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ)
+    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ, hc=hc)
 
     # Estimate initial interface state
     u★ = convert(FT, 1e-4)

src/OceanSeaIceModels/InterfaceComputations/interface_states.jl

@@ -269,9 +269,10 @@ end
 # where Ωc (the sensible heat transfer coefficient) is given by Ωc = Qc / (Tₐ - Tˢ)
 # ⟹  Tₛ = (Tᵢ * k - (Qv + Qu + Qd + Ωc * Tₐ) * h / (k - Ωc * h)
 @inline function flux_balance_temperature(st::SkinTemperature{<:ClimaSeaIce.ConductiveFlux}, Ψₛ, ℙₛ, Qc, Qv, Qu, Qd, Ψᵢ, ℙᵢ, Ψₐ, ℙₐ)
-    F = st.internal_flux
-    k = F.conductivity
-    h = Ψᵢ.h
+    F  = st.internal_flux
+    k  = F.conductivity
+    h  = Ψᵢ.h
+    hc = Ψᵢ.hc # Critical thickness for ice consolidation
 
     # Bottom temperature at the melting temperature
     Tᵢ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S)
@@ -303,6 +304,9 @@ end
     Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ)
     Tₛ⁺ = min(Tₛ⁺, Tₘ)
 
+    # If the ice is not consolidated, use the bottom temperature
+    Tₛ⁺ = ifelse(h ≥ hc, Tₛ⁺, Tᵢ)
+
     return Tₛ⁺
 end
 

src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl

@@ -21,10 +21,10 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_
 
     # TODO after ice time-step:
     #  - Adjust ocean heat flux if the ice completely melts?
-    time_step!(ocean, Δt)
+    !isnothing(ocean) && time_step!(ocean, Δt)
 
     # Time step the atmosphere
-    time_step!(atmosphere, Δt)
+    !isnothing(atmosphere) && time_step!(atmosphere, Δt)
 
     # TODO:
     # - Store fractional ice-free / ice-covered _time_ for more