diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl index e59a88af..d0383f82 100644 --- a/examples/near_global_ocean_simulation.jl +++ b/examples/near_global_ocean_simulation.jl @@ -111,20 +111,12 @@ radiation = Radiation(arch) atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(41)) -# ### Sea ice model -# -# This simulation includes a simplified representation of ice cover where the -# air-sea fluxes are shut down whenever the sea surface temperature is below -# the freezing point, - -sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce() - # ## The coupled simulation -# Next we assemble the ocean, sea ice, atmosphere, and radiation +# Next we assemble the ocean, atmosphere, and radiation # into a coupled model, -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) # We then create a coupled simulation, starting with a time step of 10 seconds # and running the simulation for 10 days. diff --git a/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl b/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl index 5e262580..97a2144a 100644 --- a/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl +++ b/experiments/mesoscale_resolving_omip/lat_lon_near_global_simulation.jl @@ -115,12 +115,10 @@ ocean.output_writers[:checkpoint] = Checkpointer(ocean.model, ##### The atmosphere ##### -backend = JRA55NetCDFBackend(4) +backend = JRA55NetCDFBackend(4) atmosphere = JRA55_prescribed_atmosphere(arch; backend) -radiation = Radiation(arch, ocean_albedo=LatitudeDependentAlbedo()) - -sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce() -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +radiation = Radiation(arch, ocean_albedo=LatitudeDependentAlbedo()) +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) @info "Set up coupled model:" @info coupled_model diff --git a/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl b/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl index 90833c65..1e6cd9ac 100644 --- a/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl +++ b/experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl @@ -118,13 +118,10 @@ initial_date = dates[1] ##### The atmosphere ##### -backend = JRA55NetCDFBackend(4) +backend = JRA55NetCDFBackend(4) atmosphere = JRA55_prescribed_atmosphere(arch; backend) -radiation = Radiation(arch) - -sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce() - -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +radiation = Radiation(arch) +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) wall_time = [time_ns()] diff --git a/src/ClimaOcean.jl b/src/ClimaOcean.jl index 0fe306f9..ae1bb084 100644 --- a/src/ClimaOcean.jl +++ b/src/ClimaOcean.jl @@ -9,7 +9,7 @@ end ClimaOcean export OceanSeaIceModel, - MinimumTemperatureSeaIce, + FreezingLimitedOceanTemperature, Radiation, LatitudeDependentAlbedo, SimilarityTheoryTurbulentFluxes, diff --git a/src/OceanSeaIceModels/OceanSeaIceModels.jl b/src/OceanSeaIceModels/OceanSeaIceModels.jl index 19ba472c..ea1b43c7 100644 --- a/src/OceanSeaIceModels/OceanSeaIceModels.jl +++ b/src/OceanSeaIceModels/OceanSeaIceModels.jl @@ -47,8 +47,8 @@ include("CrossRealmFluxes/CrossRealmFluxes.jl") using .CrossRealmFluxes -include("minimum_temperature_sea_ice.jl") include("ocean_sea_ice_model.jl") +include("freezing_limited_ocean_temperature.jl") include("time_step_ocean_sea_ice_model.jl") import .CrossRealmFluxes: @@ -60,42 +60,4 @@ const NoAtmosphereModel = OceanSeaIceModel{<:Any, Nothing} compute_atmosphere_ocean_fluxes!(coupled_model::NoAtmosphereModel) = nothing -##### -##### A fairly dumb, but nevertheless effective "sea ice model" -##### - -struct FreezingLimitedOceanTemperature{L} - liquidus :: L -end - -const FreezingLimitedCoupledModel = OceanSeaIceModel{<:FreezingLimitedOceanTemperature} - -sea_ice_concentration(::FreezingLimitedOceanTemperature) = nothing - -function compute_sea_ice_ocean_fluxes!(cm::FreezingLimitedCoupledModel) - ocean = cm.ocean - liquidus = cm.sea_ice.liquidus - grid = ocean.model.grid - arch = architecture(grid) - Sₒ = ocean.model.tracers.S - Tₒ = ocean.model.tracers.T - - launch!(arch, grid, :xyz, above_freezing_ocean_temperature!, Tₒ, Sₒ, liquidus) - - return nothing -end - -@kernel function above_freezing_ocean_temperature!(Tₒ, Sₒ, liquidus) - - i, j, k = @index(Global, NTuple) - - @inbounds begin - Sᵢ = Sₒ[i, j, k] - Tᵢ = Tₒ[i, j, k] - end - - Tₘ = melting_temperature(liquidus, Sᵢ) - @inbounds Tₒ[i, j, k] = ifelse(Tᵢ < Tₘ, Tₘ, Tᵢ) -end - end # module diff --git a/src/OceanSeaIceModels/freezing_limited_ocean_temperature.jl b/src/OceanSeaIceModels/freezing_limited_ocean_temperature.jl new file mode 100644 index 00000000..1926097c --- /dev/null +++ b/src/OceanSeaIceModels/freezing_limited_ocean_temperature.jl @@ -0,0 +1,103 @@ +using ClimaSeaIce.SeaIceThermodynamics: LinearLiquidus + +##### +##### A simple-minded yet effective "sea ice model" +##### + +struct FreezingLimitedOceanTemperature{L} + liquidus :: L +end + +""" + FreezingLimitedOceanTemperature(FT::DataType) = FreezingLimitedOceanTemperature(LinearLiquidus(FT)) + +The minimal possible sea ice representation, providing an "Insulating layer" on the surface and clipping the +temperature below to the freezing point. Not really a ``model'' per se, however, +it is the most simple way to make sure that temperature does not dip below freezing. +All fluxes are shut down when the surface is below the `T < Tₘ` except for heating to allow temperature to increase. + +the melting temperature is a function of salinity and is controlled by the `liquidus`. +""" +FreezingLimitedOceanTemperature(FT::DataType=Float64) = FreezingLimitedOceanTemperature(LinearLiquidus(FT)) + +const FreezingLimitedCoupledModel = OceanSeaIceModel{<:FreezingLimitedOceanTemperature} + +sea_ice_concentration(::FreezingLimitedOceanTemperature) = nothing + +function compute_sea_ice_ocean_fluxes!(cm::FreezingLimitedCoupledModel) + ocean = cm.ocean + liquidus = cm.sea_ice.liquidus + grid = ocean.model.grid + arch = architecture(grid) + Sₒ = ocean.model.tracers.S + Tₒ = ocean.model.tracers.T + + launch!(arch, grid, :xyz, above_freezing_ocean_temperature!, Tₒ, Sₒ, liquidus) + + return nothing +end + +@kernel function above_freezing_ocean_temperature!(Tₒ, Sₒ, liquidus) + + i, j, k = @index(Global, NTuple) + + @inbounds begin + Sᵢ = Sₒ[i, j, k] + Tᵢ = Tₒ[i, j, k] + end + + Tₘ = melting_temperature(liquidus, Sᵢ) + @inbounds Tₒ[i, j, k] = ifelse(Tᵢ < Tₘ, Tₘ, Tᵢ) +end + +function limit_fluxes_over_sea_ice!(grid, kernel_parameters, + sea_ice::FreezingLimitedOceanTemperature, + centered_velocity_fluxes, + net_tracer_fluxes, + ocean_temperature, + ocean_salinity) + + launch!(architecture(grid), grid, kernel_parameters, _limit_fluxes_over_sea_ice!, + centered_velocity_fluxes, + net_tracer_fluxes, + grid, + sea_ice.liquidus, + ocean_temperature, + ocean_salinity) + + return nothing +end + +@kernel function _limit_fluxes_over_sea_ice!(centered_velocity_fluxes, + net_tracer_fluxes, + grid, + liquidus, + ocean_temperature, + ocean_salinity) + + i, j = @index(Global, NTuple) + kᴺ = size(grid, 3) + + @inbounds begin + Tₒ = ocean_temperature[i, j, kᴺ] + Sₒ = ocean_salinity[i, j, kᴺ] + + Tₘ = melting_temperature(liquidus, Sₒ) + + τx = centered_velocity_fluxes.u + τy = centered_velocity_fluxes.v + Jᵀ = net_tracer_fluxes.T + Jˢ = net_tracer_fluxes.S + + sea_ice = Tₒ < Tₘ + cooling_sea_ice = sea_ice & (Jᵀ[i, j, 1] > 0) + + # Don't allow the ocean to cool below the minimum temperature! (make sure it heats up though!) + Jᵀ[i, j, 1] = ifelse(cooling_sea_ice, zero(grid), Jᵀ[i, j, 1]) + + # If we are in a "sea ice" region we remove all fluxes + Jˢ[i, j, 1] = ifelse(sea_ice, zero(grid), Jˢ[i, j, 1]) + τx[i, j, 1] = ifelse(sea_ice, zero(grid), τx[i, j, 1]) + τy[i, j, 1] = ifelse(sea_ice, zero(grid), τy[i, j, 1]) + end +end diff --git a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl deleted file mode 100644 index 09e2c155..00000000 --- a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl +++ /dev/null @@ -1,66 +0,0 @@ -using Oceananigans.Architectures: architecture - -import ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: limit_fluxes_over_sea_ice! - -""" - struct MinimumTemperatureSeaIce{T} - -The minimal possible sea ice representation, providing an "Insulating layer" on the surface. -Not really a ``model'' per se, however, it is the most simple way to make sure that temperature -does not dip below freezing temperature. -All fluxes are shut down when the surface is below the `minimum_temperature` except for heating. - -# Fields -- `minimum_temperature`: The minimum temperature of water. -""" -struct MinimumTemperatureSeaIce{T} - minimum_temperature :: T -end - -MinimumTemperatureSeaIce() = MinimumTemperatureSeaIce(-1.8) - -function limit_fluxes_over_sea_ice!(grid, kernel_parameters, - sea_ice::MinimumTemperatureSeaIce, - centered_velocity_fluxes, - net_tracer_fluxes, - ocean_temperature, - ocean_salinity) - - launch!(architecture(grid), grid, kernel_parameters, _cap_fluxes_on_sea_ice!, - centered_velocity_fluxes, - net_tracer_fluxes, - grid, - sea_ice.minimum_temperature, - ocean_temperature) - - return nothing -end - -@kernel function _cap_fluxes_on_sea_ice!(centered_velocity_fluxes, - net_tracer_fluxes, - grid, - minimum_temperature, - ocean_temperature) - - i, j = @index(Global, NTuple) - - @inbounds begin - Tₒ = ocean_temperature[i, j, 1] - - τx = centered_velocity_fluxes.u - τy = centered_velocity_fluxes.v - Jᵀ = net_tracer_fluxes.T - Jˢ = net_tracer_fluxes.S - - sea_ice = Tₒ < minimum_temperature - cooling_sea_ice = sea_ice & (Jᵀ[i, j, 1] > 0) - - # Don't allow the ocean to cool below the minimum temperature! (make sure it heats up though!) - Jᵀ[i, j, 1] = ifelse(cooling_sea_ice, zero(grid), Jᵀ[i, j, 1]) - - # If we are in a "sea ice" region we remove all fluxes - Jˢ[i, j, 1] = ifelse(sea_ice, zero(grid), Jˢ[i, j, 1]) - τx[i, j, 1] = ifelse(sea_ice, zero(grid), τx[i, j, 1]) - τy[i, j, 1] = ifelse(sea_ice, zero(grid), τy[i, j, 1]) - end -end diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 71b79844..f8bc3627 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -69,7 +69,7 @@ function heat_capacity(eos::TEOS10EquationOfState{FT}) where FT return convert(FT, cₚ⁰) end -function OceanSeaIceModel(ocean, sea_ice=MinimumTemperatureSeaIce(); +function OceanSeaIceModel(ocean, sea_ice=FreezingLimitedOceanTemperature(); atmosphere = nothing, radiation = nothing, similarity_theory = nothing, diff --git a/test/test_ocean_sea_ice_model_parameter_space.jl b/test/test_ocean_sea_ice_model_parameter_space.jl index d45398ba..ec4c2ccf 100644 --- a/test/test_ocean_sea_ice_model_parameter_space.jl +++ b/test/test_ocean_sea_ice_model_parameter_space.jl @@ -38,8 +38,7 @@ elapsed = 1e-9 * (time_ns() - start_time) # Fluxes are computed when the model is constructed, so we just test that this works. start_time = time_ns() -sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce() -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) elapsed = 1e-9 * (time_ns() - start_time) @info "Coupled model construction time: " * prettytime(elapsed) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 6d7fc302..123199de 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -24,15 +24,13 @@ using OrthogonalSphericalShellGrids free_surface = SplitExplicitFreeSurface(grid; substeps = 20) ocean = ocean_simulation(grid; free_surface) - backend = JRA55NetCDFBackend(4) + backend = JRA55NetCDFBackend(4) atmosphere = JRA55_prescribed_atmosphere(arch; backend) - radiation = Radiation(arch) - - sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce() + radiation = Radiation(arch) # Fluxes are computed when the model is constructed, so we just test that this works. @test begin - coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) true end end diff --git a/test/test_surface_fluxes.jl b/test/test_surface_fluxes.jl index a944b1fa..b1c8cde7 100644 --- a/test/test_surface_fluxes.jl +++ b/test/test_surface_fluxes.jl @@ -137,6 +137,53 @@ end @test turbulent_fluxes.sensible_heat[1, 1, 1] ≈ Qs @test turbulent_fluxes.latent_heat[1, 1, 1] ≈ Ql @test turbulent_fluxes.water_vapor[1, 1, 1] ≈ Mv + + @info " Testing FreezingLimitedOceanTemperature..." + + grid = LatitudeLongitudeGrid(size = (2, 2, 10), + latitude = (-0.5, 0.5), + longitude = (-0.5, 0.5), + z = (-1, 0), + topology = (Periodic, Periodic, Bounded)) + + ocean = ocean_simulation(grid; momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + bottom_drag_coefficient = 0.0) + + atmosphere = JRA55_prescribed_atmosphere(1:2; grid, backend = InMemory()) + + + fill!(ocean.model.tracers.T, -2.0) + + ocean.model.tracers.T[1, 2, 10] = 1.0 + ocean.model.tracers.T[2, 1, 10] = 1.0 + + # Cap all fluxes exept for heating ones where T < 0 + sea_ice = FreezingLimitedOceanTemperature() + + # Always cooling! + fill!(atmosphere.temperature.T, 273.15 - 20) + + coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation=nothing) + + turbulent_fluxes = coupled_model.fluxes.turbulent.fields + + # Make sure that the fluxes are zero when the temperature is below the minimum + # but not zero when it is above + u, v, _ = ocean.model.velocities + T, S = ocean.model.tracers + + for field in (u, v, T, S) + flux = surface_flux(field) + @test flux[1, 2, 10] == 0.0 # below freezing and cooling, no flux + @test flux[2, 1, 10] == 0.0 # below freezing and cooling, no flux + @test flux[1, 1, 10] != 0.0 # above freezing and cooling + @test flux[2, 2, 10] != 0.0 # above freezing and cooling + end + + # Test that the temperature has snapped up to freezing + @test minimum(ocean.model.tracers.T) == 0 end