(0.4.6) Remove thermodynamics step from ClimaOcean and use new ClimaSeaIce version (#383)

* change

* update to new syntax

* update to new syntax

* reference heat capactity

* bump patch

Author: simone-silvestri

Project.toml

@@ -2,7 +2,7 @@ name = "ClimaOcean"
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
 license = "MIT"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.4.5"
+version = "0.4.6"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -33,7 +33,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 Adapt = "4"
 CFTime = "0.1"
 CUDA = "4, 5"
-ClimaSeaIce = "0.2.1"
+ClimaSeaIce = "0.2.3"
 CubicSplines = "0.2"
 DataDeps = "0.7"
 Downloads = "1.6"

src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

@@ -94,22 +94,19 @@ end
         Sᴺ = Sₒ[i, j, Nz]
     end
 
-    # Adjust temperature 
-    Tₘ = melting_temperature(liquidus, Sᴺ)
-    ΔT = ℵ * (Tₘ - Tᴺ)
-
     # Compute total heat associated with temperature adjustment
-    δE_ice_bath = ρₒ * cₒ * ΔT
+    Tₘ = melting_temperature(liquidus, Sᴺ)
+    δE_ice_bath = ρₒ * cₒ *  (Tₘ - Tᴺ)
 
     # Compute the heat flux from ocean into ice due to sea ice melting.
     # A positive value δQ_melting > 0 corresponds to ocean cooling; ie
     # is fluxing upwards, into the ice. This occurs when applying the
     # ice bath equilibrium condition to cool down a warm ocean (δEₒ < 0).
-    δQ_melting = - δE_ice_bath * uₘ★
+    δQ_melting = - δE_ice_bath * uₘ★  
 
     # Store column-integrated ice-ocean heat flux
     @inbounds Qᶠₒ[i, j, 1] = δQ_frazil
-    @inbounds Qᵢₒ[i, j, 1] = δQ_melting
+    @inbounds Qᵢₒ[i, j, 1] = δQ_melting * ℵ # Melting depends on concentration
 end
 
 function compute_sea_ice_ocean_salinity_flux!(coupled_model)

src/OceanSeaIceModels/OceanSeaIceModels.jl

@@ -78,6 +78,7 @@ compute_atmosphere_ocean_fluxes!(::NoAtmosphereModel) = nothing
 compute_atmosphere_sea_ice_fluxes!(::NoAtmosphereModel) = nothing
 
 const PrescribedAtmosphereModel = OceanSeaIceModel{<:Any, <:PrescribedAtmosphere}
+
 compute_net_atmosphere_fluxes!(::PrescribedAtmosphereModel) = nothing
 
 # "No sea ice" implementation

src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl

@@ -68,124 +68,3 @@ function update_state!(coupled_model::OceanSeaIceModel, callbacks=[]; compute_te
 
     return nothing
 end
-
-function thermodynamic_sea_ice_time_step!(coupled_model)
-    sea_ice = coupled_model.sea_ice
-    model = sea_ice.model
-    Δt = sea_ice.Δt
-    grid = coupled_model.sea_ice.model.grid
-    arch = architecture(grid)
-    clock = model.clock
-    thermodynamics = model.ice_thermodynamics
-    ice_thickness = model.ice_thickness
-    ice_concentration = model.ice_concentration
-    ice_consolidation_thickness = model.ice_consolidation_thickness
-    top_external_heat_flux = model.external_heat_fluxes.top
-    bottom_external_heat_flux = model.external_heat_fluxes.bottom
-    ocean_salinity = coupled_model.ocean.model.tracers.S
-
-    launch!(arch, grid, :xy, update_thickness!,
-            ice_thickness,
-            grid, Δt,
-            ice_concentration,
-            ice_consolidation_thickness,
-            ocean_salinity,
-            thermodynamics,
-            top_external_heat_flux,
-            bottom_external_heat_flux,
-            clock)
-
-    return nothing
-end
-
-@inline function conservative_adjustment(ℵ, h, hᶜ)
-    V = ℵ * h # = ℵ⁺ * (h + dh)
-    dh = max(zero(h), hᶜ - h)
-    ℵ⁺ = V / (h + dh)
-    return ℵ⁺, h + dh
-end
-
-@kernel function update_thickness!(ice_thickness,
-                                   grid, Δt,
-                                   ice_concentration,
-                                   ice_consolidation_thickness,
-                                   ocean_salinity,
-                                   thermodynamics,
-                                   top_external_heat_flux,
-                                   bottom_external_heat_flux,
-                                   clock)
-
-    i, j = @index(Global, NTuple)
-    kᴺ = size(grid, 3)
-
-    phase_transitions = thermodynamics.phase_transitions
-    top_heat_bc = thermodynamics.heat_boundary_conditions.top
-    bottom_heat_bc = thermodynamics.heat_boundary_conditions.bottom
-    liquidus = phase_transitions.liquidus
-
-    Qi = thermodynamics.internal_heat_flux
-    Qu = top_external_heat_flux
-    Qb = bottom_external_heat_flux
-    Tu = thermodynamics.top_surface_temperature
-
-    @inbounds begin
-        hᶜ = ice_consolidation_thickness[i, j, kᴺ]
-        hᵢ = ice_thickness[i, j, kᴺ]
-        ℵᵢ = ice_concentration[i, j, kᴺ]
-    end
-
-    # Volume conserving adjustment to respect minimum thickness
-    ℵᵢ, hᵢ = conservative_adjustment(ℵᵢ, hᵢ, hᶜ)
-
-    # Consolidation criteria
-    @inbounds Tuᵢ = Tu[i, j, 1]
-
-    # Bottom temperature at the melting temperature
-    Sₒ = @inbounds ocean_salinity[i, j, kᴺ]
-    Tbᵢ = SeaIceThermodynamics.melting_temperature(liquidus, Sₒ)
-    ℰb = SeaIceThermodynamics.latent_heat(phase_transitions, Tbᵢ)
-    ℰu = SeaIceThermodynamics.latent_heat(phase_transitions, Tuᵢ)
-
-    # Retrieve fluxes
-    @inbounds begin
-        Quᵢ = Qu[i, j, 1]
-        Qbᵢ = Qb[i, j, 1]
-    end
-
-    # If ice is consolidated, compute tendency for an ice slab; otherwise
-    # just add ocean fluxes from frazil ice formation or melting
-    # wb = - Qbᵢ / ℰb
-
-    𝓀 = Qi.parameters.flux.conductivity
-    Qiᵢ = - 𝓀 * (Tuᵢ - Tbᵢ) / hᵢ * (hᵢ > hᶜ) # getflux(Qi, i, j, grid, Tuᵢ, clock, model_fields)
-
-    # Upper (top) and bottom interface velocities
-    w_top = (Quᵢ - Qiᵢ) / ℰu # < 0 => melting
-    w_bot = +Qiᵢ / ℰb # < 0 => freezing
-    w_frz = -Qbᵢ / ℰb # < 0 => freezing
-
-    Δh_top = w_top * Δt * ℵᵢ
-    Δh_bot = w_bot * Δt * ℵᵢ
-    ΔV_frz = w_frz * Δt # frazil flux contributes from entire cell
-
-    # Compute frazil growth: lateral first, then vertical
-    # dV = dh * ℵ + h * dℵ
-    dℵ = min(1 - ℵᵢ, ΔV_frz / hᵢ)
-    dℵ = max(dℵ, zero(dℵ))
-    ℵ⁺ = ℵᵢ + dℵ
-    Δh_frz = ΔV_frz - hᵢ * ℵ⁺
-
-    Δh = Δh_frz + Δh_bot + Δh_top
-
-    h⁺ = hᵢ + Δh / ℵ⁺ * (Δh > 0)
-    h⁺ = max(zero(h⁺), h⁺)
-
-    # Adjust again to be paranoid?
-    ℵ⁺, h⁺ = conservative_adjustment(ℵ⁺, h⁺, hᶜ)
-
-    @inbounds begin
-        ice_thickness[i, j, kᴺ] = h⁺
-        ice_concentration[i, j, kᴺ] = ℵ⁺
-    end
-end
-