Use a thermodynamic tendency ice-ocean salt flux (#539)

* add a thermodynamic tendency

* Update src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* Update src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

* bugfix

---------

Co-authored-by: Gregory L. Wagner <[email protected]>

Author: simone-silvestri

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -42,11 +42,9 @@ mutable struct AtmosphereInterface{J, F, ST, P}
     properties :: P
 end
 
-mutable struct SeaIceOceanInterface{J, P, H, A}
+mutable struct SeaIceOceanInterface{J, P}
     fluxes :: J
     properties :: P
-    previous_ice_thickness :: H
-    previous_ice_concentration :: A
 end
 
 mutable struct ComponentInterfaces{AO, ASI, SIO, C, AP, OP, SIP, EX, P}
@@ -261,8 +259,6 @@ sea_ice_ocean_interface(sea_ice, ocean) = nothing
 function sea_ice_ocean_interface(sea_ice::SeaIceSimulation, ocean;
                                  characteristic_melting_speed = 1e-5)
 
-    previous_ice_thickness = deepcopy(sea_ice.model.ice_thickness)
-    previous_ice_concentration = deepcopy(sea_ice.model.ice_concentration)
     io_bottom_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid)
     io_frazil_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid)
     io_salt_flux = Field{Center, Center, Nothing}(ocean.model.grid)
@@ -281,10 +277,7 @@ function sea_ice_ocean_interface(sea_ice::SeaIceSimulation, ocean;
 
     io_properties = (; characteristic_melting_speed)
 
-    return SeaIceOceanInterface(io_fluxes,
-                                io_properties,
-                                previous_ice_thickness,
-                                previous_ice_concentration)
+    return SeaIceOceanInterface(io_fluxes, io_properties)
 end
 
 default_ai_temperature(sea_ice) = nothing

src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl

@@ -15,7 +15,7 @@ function compute_sea_ice_ocean_fluxes!(coupled_model)
     Sᵢ = sea_ice.model.tracers.S
     ℵᵢ = sea_ice.model.ice_concentration
     hᵢ = sea_ice.model.ice_thickness
-    h⁻ = coupled_model.interfaces.sea_ice_ocean_interface.previous_ice_thickness
+    Gh = sea_ice.model.ice_thermodynamics.thermodynamic_tendency
 
     ocean_properties = coupled_model.interfaces.ocean_properties
     liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus
@@ -35,7 +35,7 @@ function compute_sea_ice_ocean_fluxes!(coupled_model)
     # What about the latent heat removed from the ocean when ice forms?
     # Is it immediately removed from the ocean? Or is it stored in the ice?
     launch!(arch, grid, :xy, _compute_sea_ice_ocean_fluxes!,
-            sea_ice_ocean_fluxes, grid, clock, hᵢ, h⁻, ℵᵢ, Sᵢ, Tₒ, Sₒ, uᵢ, vᵢ,
+            sea_ice_ocean_fluxes, grid, clock, hᵢ, ℵᵢ, Sᵢ, Gh, Tₒ, Sₒ, uᵢ, vᵢ,
             τs, liquidus, ocean_properties, interface_properties, Δt)
 
     return nothing
@@ -45,9 +45,9 @@ end
                                                 grid,
                                                 clock,
                                                 ice_thickness,
-                                                previous_ice_thickness,
                                                 ice_concentration,
                                                 ice_salinity,
+                                                thermodynamic_tendency,
                                                 ocean_temperature,
                                                 ocean_salinity,
                                                 sea_ice_u_velocity,
@@ -73,7 +73,6 @@ end
     Sᵢ  = ice_salinity
     hᵢ  = ice_thickness
     ℵᵢ  = ice_concentration
-    h⁻  = previous_ice_thickness
     ρₒ  = ocean_properties.reference_density
     cₒ  = ocean_properties.heat_capacity
     uₘ★ = interface_properties.characteristic_melting_speed
@@ -136,16 +135,10 @@ end
     τₒᵢ = sea_ice_ocean_stresses
 
     @inbounds begin
-        # Change in thickness
-        Δh = hᵢ[i, j, 1] - h⁻[i, j, 1]
-
         # Update surface salinity flux.
         # Note: the Δt below is the ocean time-step, eg.
         # ΔS = ⋯ - ∮ Jˢ dt ≈ ⋯ - Δtₒ * Jˢ
-        Jˢ[i, j, 1] = Δh / Δt * (Sᵢ[i, j, 1] - Sₒ[i, j, Nz])
-
-        # Update previous ice thickness
-        h⁻[i, j, 1] = hᵢ[i, j, 1]
+        Jˢ[i, j, 1] = thermodynamic_tendency[i, j, 1] * (Sᵢ[i, j, 1] - Sₒ[i, j, Nz])
 
         # Momentum stresses
         τx[i, j, 1] = x_momentum_stress(i, j, Nz, grid, τₒᵢ, clock, sea_ice_fields)

src/OceanSeaIceModels/freezing_limited_ocean_temperature.jl

@@ -29,6 +29,7 @@ sea_ice_thickness(::FreezingLimitedOceanTemperature) = nothing
 # does not matter
 reference_density(::FreezingLimitedOceanTemperature) = 0
 heat_capacity(::FreezingLimitedOceanTemperature) = 0
+time_step!(::FreezingLimitedOceanTemperature, Δt) = nothing
 
 function compute_sea_ice_ocean_fluxes!(cm::FreezingLimitedCoupledModel)
     ocean = cm.ocean

src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl

@@ -14,24 +14,10 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_
     ocean = coupled_model.ocean
     sea_ice = coupled_model.sea_ice
     atmosphere = coupled_model.atmosphere
-    clock = coupled_model.clock
 
     # Eventually, split out into OceanOnlyModel
-    if sea_ice isa SeaIceSimulation
-        h = sea_ice.model.ice_thickness
-        fill_halo_regions!(h)
-
-        # Initialization
-        if coupled_model.clock.iteration == 0
-            @info "Initializing coupled model ice thickness..."
-            h⁻ = coupled_model.interfaces.sea_ice_ocean_interface.previous_ice_thickness
-            hⁿ = coupled_model.sea_ice.model.ice_thickness
-            parent(h⁻) .= parent(hⁿ)
-        end
-
-        time_step!(sea_ice, Δt)
-    end
-
+    !isnothing(sea_ice) && time_step!(sea_ice, Δt)
+    
     # TODO after ice time-step:
     #  - Adjust ocean heat flux if the ice completely melts?
     time_step!(ocean, Δt)