Merge MinimumTemperatureSeaIce and FreezingLimitedOceanTemperature (

#255) * new sea ice? * done? * add a constructor * docstring * comment * comment * tests * tests * switch includes * no need for FT I guess * Update src/OceanSeaIceModels/freezing_limited_ocean_temperature.jl Co-authored-by: Gregory L. Wagner <[email protected]> * a couple of typos * Remove MinimumTemperatureSeaIce from repo * reorder includes * reorder * Update src/OceanSeaIceModels/ocean_sea_ice_model.jl Co-authored-by: Gregory L. Wagner <[email protected]> * FT default --------- Co-authored-by: Gregory L. Wagner <[email protected]> Co-authored-by: Navid C. Constantinou <[email protected]>
CliMA · Nov 16, 2024 · 50e05f2 · 50e05f2
1 parent b11863f
commit 50e05f2
Show file tree

Hide file tree

Showing 11 changed files with 165 additions and 135 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
diff --git 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
@@ -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
@@ -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
@@ -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