Integrate sea ice in the flux computation

Project.toml

@@ -2,7 +2,7 @@ name = "ClimaOcean"
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
 license = "MIT"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.3.2"
+version = "0.4.0"
 
 [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.1.3"
+ClimaSeaIce = "0.1.6"
 CubicSplines = "0.2"
 DataDeps = "0.7"
 Downloads = "1.6"
@@ -42,7 +42,7 @@ JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
 MPI = "0.20"
 NCDatasets = "0.12, 0.13, 0.14"
-Oceananigans = "0.95.4 - 0.99"
+Oceananigans = "0.95.5 - 0.99"
 OffsetArrays = "1.14"
 OrthogonalSphericalShellGrids = "0.2.1"
 Scratch = "1"

examples/idealized_single_column_simulation.jl

@@ -1,21 +1,23 @@
 using ClimaOcean
+using ClimaSeaIce
 using Oceananigans
 using Oceananigans.Units
 using SeawaterPolynomials
 
 # Ocean state parameters
-T₀ = 20   # Surface temperature, ᵒC
+T₀ = 2    # Surface temperature, ᵒC
 S₀ = 35   # Surface salinity
 N² = 1e-4 # Buoyancy gradient due to temperature stratification
 f = 0     # Coriolis parameter
 
 # Atmospheric state parameters
-Tₐ = 273.15 + 20 # Kelvin
+Tₐ = 273.15 - 10 # Kelvin
 u₁₀ = 10 # wind at 10 m, m/s
 qₐ = 0.01 # specific humidity
 Qℓ = 400 # shortwave radiation (W m⁻², positive means heating right now)
 
 # Build the atmosphere
+radiation = Radiation(ocean_albedo=0.1)
 atmosphere_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
 atmosphere_times = range(0, 1day, length=3)
 atmosphere = PrescribedAtmosphere(atmosphere_grid, atmosphere_times)
@@ -29,26 +31,52 @@ parent(atmosphere.downwelling_radiation.shortwave) .= Qℓ # W
 grid = RectilinearGrid(size=100, z=(-400, 0), topology=(Flat, Flat, Bounded))
 ocean = ocean_simulation(grid, coriolis=FPlane(; f))
 
-eos = ocean.model.buoyancy.model.equation_of_state
-g = ocean.model.buoyancy.model.gravitational_acceleration
+eos = ocean.model.buoyancy.formulation.equation_of_state
+g = ocean.model.buoyancy.formulation.gravitational_acceleration
 α = SeawaterPolynomials.thermal_expansion(T₀, S₀, 0, eos)
 dTdz = N² / (α * g)
 Tᵢ(z) = T₀ + dTdz * z
 set!(ocean.model, T=Tᵢ, S=S₀)
 
-radiation = Radiation(ocean_albedo=0.1)
-model = OceanSeaIceModel(ocean; atmosphere, radiation)
-simulation = Simulation(model, Δt=5minutes, stop_time=30days)
-
-Q = model.fluxes.total.ocean.heat
-Ql = model.fluxes.turbulent.fields.latent_heat
-Qs = model.fluxes.turbulent.fields.sensible_heat
-τx = model.fluxes.turbulent.fields.x_momentum
-τy = model.fluxes.turbulent.fields.y_momentum
-Fv = model.fluxes.turbulent.fields.water_vapor
-
-fluxes = (; Q, Ql, Qs, τx, τy, Fv)
-outputs = merge(ocean.model.velocities, ocean.model.tracers, fluxes)
+sea_ice_grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+top_sea_ice_temperature = Field{Center, Center, Nothing}(sea_ice_grid)
+top_heat_boundary_condition = PrescribedTemperature(top_sea_ice_temperature)
+ice_thermodynamics = SlabSeaIceThermodynamics(sea_ice_grid; top_heat_boundary_condition)
+
+top_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid)
+bottom_sea_ice_heat_flux = Field{Center, Center, Nothing}(sea_ice_grid)
+
+sea_ice_model = SeaIceModel(sea_ice_grid;
+                            top_heat_flux = top_sea_ice_heat_flux,
+                            bottom_heat_flux = bottom_sea_ice_heat_flux,
+                            ice_thermodynamics)
+
+set!(sea_ice_model.ice_concentration, 1.0)
+set!(sea_ice_model.ice_thickness, 1.0)
+
+sea_ice = Simulation(sea_ice_model, Δt=10minutes)
+
+model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
+simulation = Simulation(model, Δt=10minutes, stop_time=10days)
+
+Ql  = model.fluxes.turbulent.fields.ocean.latent_heat
+Qs  = model.fluxes.turbulent.fields.ocean.sensible_heat
+τx  = model.fluxes.turbulent.fields.ocean.x_momentum
+τy  = model.fluxes.turbulent.fields.ocean.y_momentum
+Fv  = model.fluxes.turbulent.fields.ocean.water_vapor
+Qib = sea_ice_model.external_heat_fluxes.bottom
+
+# TODO: the total fluxes are defined on _interfaces_ between components:
+# atmopshere_ocean, atmosphere_sea_ice, ocean_sea_ice. They aren't defined wrt to 
+# just one component
+Qo = model.fluxes.total.ocean.heat
+Qi = model.fluxes.total.sea_ice.heat
+
+fluxes = (; Qo, Qi, Qib, Ql, Qs, τx, τy, Fv)
+ocean_outputs = merge(ocean.model.velocities, ocean.model.tracers)
+sea_ice_outputs = (; h=sea_ice.model.ice_thickness)
+outputs = merge(ocean_outputs, sea_ice_outputs, fluxes)
+
 ow = JLD2OutputWriter(ocean.model, outputs,
                       schedule = TimeInterval(1hour),
                       filename = "idealized_atmosphere.jld2",
@@ -58,25 +86,28 @@ simulation.output_writers[:ow] = ow
 
 run!(simulation)
 
-ut = FieldTimeSeries("idealized_atmosphere.jld2", "u")
-vt = FieldTimeSeries("idealized_atmosphere.jld2", "v")
-Tt = FieldTimeSeries("idealized_atmosphere.jld2", "T")
-St = FieldTimeSeries("idealized_atmosphere.jld2", "S")
-τx = FieldTimeSeries("idealized_atmosphere.jld2", "τx")
-Q = FieldTimeSeries("idealized_atmosphere.jld2", "Q")
-Nt = length(St)
+ut  = FieldTimeSeries("idealized_atmosphere.jld2", "u")
+vt  = FieldTimeSeries("idealized_atmosphere.jld2", "v")
+Tt  = FieldTimeSeries("idealized_atmosphere.jld2", "T")
+St  = FieldTimeSeries("idealized_atmosphere.jld2", "S")
+τx  = FieldTimeSeries("idealized_atmosphere.jld2", "τx")
+Qo  = FieldTimeSeries("idealized_atmosphere.jld2", "Qo")
+Qib = FieldTimeSeries("idealized_atmosphere.jld2", "Qib")
+h   = FieldTimeSeries("idealized_atmosphere.jld2", "h")
+Nt  = length(St)
 
 using GLMakie
 
 fig = Figure(size=(1000, 800))
 
 axτ = Axis(fig[1, 1:3], xlabel="Time (days)", ylabel="Zonal momentum \n flux (N m⁻²)", xaxisposition=:top)
 axQ = Axis(fig[2, 1:3], xlabel="Time (days)", ylabel="Heat flux \n (W m⁻²)")
-axT = Axis(fig[3, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)")
-axS = Axis(fig[3, 2], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)")
-axu = Axis(fig[3, 3], xlabel="Velocities (m s⁻¹)", ylabel="z (m)", yaxisposition=:right)
+axh = Axis(fig[3, 1:3], xlabel="Time (days)", ylabel="Sea ice thickness \n (m)")
+axT = Axis(fig[4, 1], xlabel="Temperature (ᵒC)", ylabel="z (m)")
+axS = Axis(fig[4, 2], xlabel="Salinity (g kg⁻¹)", ylabel="z (m)")
+axu = Axis(fig[4, 3], xlabel="Velocities (m s⁻¹)", ylabel="z (m)", yaxisposition=:right)
 
-slider = Slider(fig[4, 1:3], startvalue=1, range=1:Nt)
+slider = Slider(fig[5, 1:3], startvalue=1, range=1:Nt)
 n = slider.value
 
 un = @lift ut[$n]
@@ -96,11 +127,15 @@ tn = @lift t[$n]
 lines!(axτ, t, τx[:], label="τx")
 vlines!(axτ, tn)
 
-lines!(axQ, t, Q[:])
+lines!(axQ, t, Qo[:])
 vlines!(axQ, tn)
 
+lines!(axh, t, h[:])
+vlines!(axh, tn)
+
 rowsize!(fig.layout, 1, Relative(0.2))
 rowsize!(fig.layout, 2, Relative(0.2))
+rowsize!(fig.layout, 3, Relative(0.2))
 
 xlims!(axS, 34.9, 35.1)
 xlims!(axu, -0.1, 4.0)

experiments/mesoscale_resolving_omip/prototype_omip_simulation.jl

@@ -11,7 +11,7 @@ using Oceananigans.Units
 using ClimaOcean
 using ClimaOcean.OceanSimulations
 using ClimaOcean.OceanSeaIceModels
-using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: Radiation, SimilarityTheoryTurbulentFluxes
+using ClimaOcean.OceanSeaIceModels.CrossRealmFluxes: Radiation, SimilarityTheoryFluxes
 using ClimaOcean.VerticalGrids: exponential_z_faces
 using ClimaOcean.JRA55
 using ClimaOcean.ECCO

experiments/one_degree_simulation/one_degree_simulation.jl

@@ -38,7 +38,7 @@ bottom_height = regrid_bathymetry(underlying_grid;
 tampered_bottom_height = deepcopy(bottom_height)
 view(tampered_bottom_height, 102:103, 124, 1) .= -400
 
-grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_height))
+grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_height); active_cells_map=true)
 
 #####
 ##### Closures
@@ -75,14 +75,24 @@ forcing = (T=FT, S=FS)
 momentum_advection = VectorInvariant()
 tracer_advection   = Centered(order=2)
 
+free_surface = SplitExplicitFreeSurface(grid; substeps=30)
+
 # Should we add a side drag since this is at a coarser resolution?
 ocean = ocean_simulation(grid; momentum_advection, tracer_advection,
-                         closure, forcing,
-                         tracers = (:T, :S, :e))
+                         free_surface, forcing)
 
 set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
                   S=ECCOMetadata(:salinity;    dates=first(dates)))
 
+#####
+##### Sea ice simulation
+#####
+
+sea_ice = sea_ice_simulation(grid)
+
+set!(sea_ice.model.ice_thickness,     ECCOMetadata(:sea_ice_thickness;     dates=first(dates)); inpainting=NearestNeighborInpainting(0))
+set!(sea_ice.model.ice_concentration, ECCOMetadata(:sea_ice_concentration; dates=first(dates)); inpainting=NearestNeighborInpainting(0))
+
 #####
 ##### Atmospheric forcing
 #####
@@ -94,8 +104,8 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20))
 ##### Coupled simulation
 #####
 
-coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) 
-simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days)
+coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) 
+simulation = Simulation(coupled_model; Δt=1minutes, stop_time=2*365days)
 
 #####
 ##### Run it!
@@ -105,19 +115,22 @@ wall_time = Ref(time_ns())
 
 function progress(sim)
     ocean = sim.model.ocean
+    sea_ice = sim.model.sea_ice
     u, v, w = ocean.model.velocities
     T = ocean.model.tracers.T
+    h = sea_ice.model.ice_thickness
     Tmax = maximum(interior(T))
     Tmin = minimum(interior(T))
+    hmax = maximum(interior(h))
     umax = (maximum(abs, interior(u)),
             maximum(abs, interior(v)),
             maximum(abs, interior(w)))
 
     step_time = 1e-9 * (time_ns() - wall_time[])
 
-    @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s \n",
+    @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, max(h): %.2f m, wall time: %s \n",
                    prettytime(sim), iteration(sim), prettytime(sim.Δt),
-                   umax..., Tmax, Tmin, prettytime(step_time))
+                   umax..., Tmax, Tmin, hmax, prettytime(step_time))
 
      wall_time[] = time_ns()