diff --git a/Project.toml b/Project.toml index 5667342b..5732718b 100644 --- a/Project.toml +++ b/Project.toml @@ -44,7 +44,7 @@ ClimaOceanReactantExt = "Reactant" Adapt = "4" CFTime = "0.1, 0.2" CUDA = "4, 5" -ClimaSeaIce = "0.2.6" +ClimaSeaIce = "0.3" CondaPkg = "0.2.28" CubicSplines = "0.2" DataDeps = "0.7" diff --git a/experiments/arctic_simulation.jl b/experiments/arctic_simulation.jl index 31e8e3f0..f93ad70a 100644 --- a/experiments/arctic_simulation.jl +++ b/experiments/arctic_simulation.jl @@ -59,7 +59,7 @@ set!(ocean.model, T=Metadatum(:temperature; dataset), ##### A Prognostic Sea-ice model ##### -using ClimaSeaIce.SeaIceMomentumEquations +using ClimaSeaIce.SeaIceDynamics using ClimaSeaIce.Rheologies # Remember to pass the SSS as a bottom bc to the sea ice! diff --git a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl index 77f103bc..8c068096 100644 --- a/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/InterfaceComputations/sea_ice_ocean_fluxes.jl @@ -1,6 +1,6 @@ using Oceananigans.Operators: Δzᶜᶜᶜ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature -using ClimaSeaIce.SeaIceMomentumEquations: x_momentum_stress, y_momentum_stress +using ClimaSeaIce.SeaIceDynamics: x_momentum_stress, y_momentum_stress function compute_sea_ice_ocean_fluxes!(coupled_model) ocean = coupled_model.ocean diff --git a/src/SeaIceSimulations.jl b/src/SeaIceSimulations.jl index 58766da0..08f3de9e 100644 --- a/src/SeaIceSimulations.jl +++ b/src/SeaIceSimulations.jl @@ -17,19 +17,21 @@ using Oceananigans.Operators using ClimaSeaIce using ClimaSeaIce: SeaIceModel, SlabSeaIceThermodynamics, PhaseTransitions, ConductiveFlux using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium +using ClimaSeaIce.SeaIceDynamics: SplitExplicitSolver, SemiImplicitStress, SeaIceMomentumEquation, StressBalanceFreeDrift +using ClimaSeaIce.Rheologies: IceStrength, ElastoViscoPlasticRheology using ClimaOcean.OceanSimulations: Default -function sea_ice_simulation(grid; +function sea_ice_simulation(grid, ocean=nothing; Δt = 5minutes, - ice_salinity = 0, # psu + ice_salinity = 4, # psu advection = nothing, # for the moment tracers = (), ice_heat_capacity = 2100, # J kg⁻¹ K⁻¹ ice_consolidation_thickness = 0.05, # m ice_density = 900, # kg m⁻³ - dynamics = nothing, - bottom_heat_boundary_condition = IceWaterThermalEquilibrium(), + dynamics = sea_ice_dynamics(grid, ocean), + bottom_heat_boundary_condition = nothing, phase_transitions = PhaseTransitions(; ice_heat_capacity, ice_density), conductivity = 2, # kg m s⁻³ K⁻¹ internal_heat_flux = ConductiveFlux(; conductivity)) @@ -41,6 +43,16 @@ function sea_ice_simulation(grid; top_surface_temperature = Field{Center, Center, Nothing}(grid) top_heat_boundary_condition = PrescribedTemperature(top_surface_temperature) + if isnothing(bottom_heat_boundary_condition) + if isnothing(ocean) + surface_ocean_salinity = 0 + else + kᴺ = size(grid, 3) + surface_ocean_salinity = interior(ocean.model.tracers.S, :, :, kᴺ:kᴺ) + end + bottom_heat_boundary_condition = IceWaterThermalEquilibrium(surface_ocean_salinity) + end + ice_thermodynamics = SlabSeaIceThermodynamics(grid; internal_heat_flux, phase_transitions, @@ -50,16 +62,12 @@ function sea_ice_simulation(grid; bottom_heat_flux = Field{Center, Center, Nothing}(grid) top_heat_flux = Field{Center, Center, Nothing}(grid) - # top_momentum_stress = (u = Field{Face, Center, Nothing}(grid), - # v = Field{Center, Face, Nothing}(grid)) - # Build the sea ice model sea_ice_model = SeaIceModel(grid; ice_salinity, advection, tracers, ice_consolidation_thickness, - # top_momentum_stress, ice_thermodynamics, dynamics, bottom_heat_flux, @@ -73,4 +81,40 @@ function sea_ice_simulation(grid; return sea_ice end -end \ No newline at end of file +function sea_ice_dynamics(grid, ocean=nothing; + sea_ice_ocean_drag_coefficient = 5.5e-3, + rheology = ElastoViscoPlasticRheology(pressure_formulation = IceStrength()), + coriolis = nothing, + free_drift = nothing, + solver = SplitExplicitSolver(120)) + + if isnothing(ocean) + SSU = Oceananigans.Fields.ZeroField() + SSV = Oceananigans.Fields.ZeroField() + else + SSU = view(ocean.model.velocities.u, :, :, grid.Nz) + SSV = view(ocean.model.velocities.v, :, :, grid.Nz) + if isnothing(coriolis) + coriolis = ocean.model.coriolis + end + end + + sea_ice_ocean_drag_coefficient = convert(eltype(grid), sea_ice_ocean_drag_coefficient) + + τo = SemiImplicitStress(uₑ=SSU, vₑ=SSV, Cᴰ=sea_ice_ocean_drag_coefficient) + τua = Field{Face, Center, Nothing}(grid) + τva = Field{Center, Face, Nothing}(grid) + + if isnothing(free_drift) + free_drift = StressBalanceFreeDrift((u=τua, v=τva), τo) + end + + return SeaIceMomentumEquation(grid; + coriolis, + top_momentum_stress = (u=τua, v=τva), + bottom_momentum_stress = τo, + rheology, + solver) +end + +end diff --git a/test/test_ocean_sea_ice_model.jl b/test/test_ocean_sea_ice_model.jl index 71e9e416..6d8edcdb 100644 --- a/test/test_ocean_sea_ice_model.jl +++ b/test/test_ocean_sea_ice_model.jl @@ -4,7 +4,7 @@ using CUDA using Oceananigans.OrthogonalSphericalShellGrids using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature! using ClimaSeaIce.SeaIceThermodynamics: melting_temperature -using ClimaSeaIce.SeaIceMomentumEquations +using ClimaSeaIce.SeaIceDynamics using ClimaSeaIce.Rheologies @inline kernel_melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k]) @@ -81,17 +81,7 @@ using ClimaSeaIce.Rheologies ##### Coupled ocean-sea ice and prescribed atmosphere ##### - # Adding a sea ice model to the coupled model - τua = Field{Face, Center, Nothing}(grid) - τva = Field{Center, Face, Nothing}(grid) - - dynamics = SeaIceMomentumEquation(grid; - coriolis = ocean.model.coriolis, - top_momentum_stress = (u=τua, v=τva), - rheology = ElastoViscoPlasticRheology(), - solver = SplitExplicitSolver(120)) - - sea_ice = sea_ice_simulation(grid; dynamics, advection=WENO(order=7)) + sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=7)) liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus # Set the ocean temperature and salinity diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index e594b513..7f0c3267 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -18,7 +18,7 @@ using Oceananigans.TimeSteppers: update_state! using Oceananigans.Units: hours, days using ClimaOcean.DataWrangling: all_dates -using ClimaSeaIce.SeaIceMomentumEquations +using ClimaSeaIce.SeaIceDynamics using ClimaSeaIce.Rheologies import ClimaOcean.OceanSeaIceModels.InterfaceComputations: surface_specific_humidity