A single column sea ice model with prescribed atmosphere and ocean

experiments/arctic_simulation.jl

@@ -81,7 +81,7 @@ dynamics = SeaIceMomentumEquation(grid;
                                   rheology = ElastoViscoPlasticRheology(),
                                   solver = SplitExplicitSolver(120))
 
-sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) 
+sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition) #, dynamics, advection=WENO(order=7)) 
 
 set!(sea_ice.model, h=Metadata(:sea_ice_thickness;     dataset), 
                     ℵ=Metadata(:sea_ice_concentration; dataset))

experiments/flux_climatology/flux_climatology.jl

@@ -115,32 +115,6 @@ end
 ##### A prescribed ocean...
 #####
 
-struct PrescribedOcean{A, G, C, U, T, F} <: AbstractModel{Nothing}
-    architecture :: A       
-    grid :: G        
-    clock :: Clock{C}
-    velocities :: U
-    tracers :: T
-    timeseries :: F
-end
-
-function PrescribedOcean(timeseries; 
-                         grid, 
-                         clock=Clock{Float64}(time = 0)) 
-
-    τx = Field{Face, Center, Nothing}(grid)
-    τy = Field{Center, Face, Nothing}(grid)
-    Jᵀ = Field{Center, Center, Nothing}(grid)
-    Jˢ = Field{Center, Center, Nothing}(grid)
-
-    u = XFaceField(grid,  boundary_conditions=FieldBoundaryConditions(grid, (Face,   Center, Center), top = FluxBoundaryCondition(τx)))
-    v = YFaceField(grid,  boundary_conditions=FieldBoundaryConditions(grid, (Center, Face,   Center), top = FluxBoundaryCondition(τy)))
-    T = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jᵀ)))
-    S = CenterField(grid, boundary_conditions=FieldBoundaryConditions(grid, (Center, Center, Center), top = FluxBoundaryCondition(Jˢ)))
-
-    PrescribedOcean(architecture(grid), grid, clock, (; u, v, w=ZeroField()), (; T, S), timeseries)
-end
-
 # ...with prescribed velocity and tracer fields
 version = ECCO4Monthly()
 dates   = all_ECCO_dates(version)[1:24]
@@ -157,35 +131,6 @@ grid = ECCO.ECCO_immersed_grid(arch)
 ocean_model = PrescribedOcean((; u, v, T, S); grid)
 ocean = Simulation(ocean_model, Δt=3hour, stop_time=365days)
 
-#####
-##### Need to extend a couple of methods
-#####
-
-function time_step!(model::PrescribedOcean, Δt; callbacks=[], euler=true)
-    tick!(model.clock, Δt)
-    time = Time(model.clock.time)
-
-    possible_fts = merge(model.velocities, model.tracers)
-    time_series_tuple = extract_field_time_series(possible_fts)
-
-    for fts in time_series_tuple
-        update_field_time_series!(fts, time)
-    end
-
-    parent(model.velocities.u) .= parent(model.timeseries.u[time])
-    parent(model.velocities.v) .= parent(model.timeseries.v[time])
-    parent(model.tracers.T)    .= parent(model.timeseries.T[time])
-    parent(model.tracers.S)    .= parent(model.timeseries.S[time])
-
-    return nothing
-end
-
-update_state!(::PrescribedOcean) = nothing
-timestepper(::PrescribedOcean) = nothing
-
-reference_density(ocean::Simulation{<:PrescribedOcean}) = 1025.6
-heat_capacity(ocean::Simulation{<:PrescribedOcean}) = 3995.6
-
 #####
 ##### A prescribed atmosphere...
 #####

experiments/single_column_sea_ice.jl

@@ -0,0 +1,114 @@
+using ClimaOcean
+using ClimaSeaIce
+using Oceananigans
+using Oceananigans.Grids
+using Oceananigans.Units
+using ClimaOcean.OceanSimulations
+using ClimaOcean.ECCO
+using ClimaOcean.DataWrangling
+using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium
+using Printf
+
+# A single column grid with one point in the Arctic ocean
+grid = LatitudeLongitudeGrid(size = 1, 
+                             latitude = 87.0,
+                             longitude = 33.0,
+                             z = (-10, 0),
+                             topology = (Flat, Flat, Bounded))
+
+#####
+##### A Prescribed Ocean model
+#####
+
+dataset = ECCO4Monthly()
+dates   = all_dates(dataset)[1:24]
+
+T = ECCOFieldTimeSeries(:temperature, grid; dates, inpainting=nothing, time_indices_in_memory=24)
+S = ECCOFieldTimeSeries(:salinity, grid;    dates, inpainting=nothing, time_indices_in_memory=24)
+
+ocean_model = PrescribedOceanModel((; T, S); grid)
+ocean = Simulation(ocean_model, Δt=30minutes, stop_time=730days)
+
+##### 
+##### A Prescribed Atmosphere model
+#####
+
+atmosphere = JRA55PrescribedAtmosphere(latitude = 87.0,
+                                       longitude = 33.0,
+                                       backend = InMemory())
+
+#####
+##### A Prognostic Sea-ice model with no dynamics
+#####
+
+# Remember to pass the SSS as a bottom bc to the sea ice!
+SSS = view(ocean.model.tracers.S, :, :, 1)
+bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS)
+
+sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition) 
+
+set!(sea_ice.model, h=Metadata(:sea_ice_thickness;     dataset, dates=first(dates)), 
+                    ℵ=Metadata(:sea_ice_concentration; dataset, dates=first(dates)))
+
+#####
+##### Radiation of the Ocean and the Sea ice
+#####
+
+# to mimick snow, as in Semtner (1975) we use a sea ice albedo that is 0.64 for top temperatures 
+# above -0.1 ᵒC and 0.75 for top temperatures below -0.1 ᵒC. We also disable radiation for the ocean.
+radiation = Radiation(sea_ice_albedo=0.7, ocean_albedo=1, ocean_emissivity=0)
+
+#####
+##### Arctic coupled model
+#####
+
+arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
+arctic = Simulation(arctic, Δt=30minutes, stop_time=730days)
+
+# Sea-ice variables
+h = sea_ice.model.ice_thickness
+ℵ = sea_ice.model.ice_concentration
+
+# Fluxes
+Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature
+Qˡ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat
+Qˢ = arctic.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat
+Qⁱ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.interface_heat
+Qᶠ = arctic.model.interfaces.sea_ice_ocean_interface.fluxes.frazil_heat
+Qᵗ = arctic.model.interfaces.net_fluxes.sea_ice_top.heat
+Qᴮ = arctic.model.interfaces.net_fluxes.sea_ice_bottom.heat
+
+# Output writers
+arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ),
+                                                 filename = "sea_ice_quantities.jld2",
+                                                 schedule = IterationInterval(12),
+                                                 overwrite_existing=true)
+
+wall_time = Ref(time_ns())
+
+using Statistics
+
+function progress(sim)
+    sea_ice = sim.model.sea_ice
+    h = first(sea_ice.model.ice_thickness)
+    ℵ = first(sea_ice.model.ice_concentration)
+    T = first(sim.model.interfaces.atmosphere_sea_ice_interface.temperature)
+
+    step_time = 1e-9 * (time_ns() - wall_time[])
+
+    msg1 = @sprintf("time: %s, iteration: %d, Δt: %s, ", prettytime(sim), iteration(sim), prettytime(sim.Δt))
+    msg2 = @sprintf("h: %.2e m, ℵ: %.2e , T: %.2e ᴼC", h, ℵ, T)
+    msg3 = @sprintf("wall time: %s \n", prettytime(step_time))
+
+    @info msg1 * msg2 * msg3
+
+    wall_time[] = time_ns()
+    return nothing
+end
+
+# And add it as a callback to the simulation.
+add_callback!(arctic, progress, IterationInterval(10))
+
+run!(arctic)
+
+h = FieldTimeSeries("sea_ice_quantities.jld2", "h")

src/ClimaOcean.jl

@@ -41,6 +41,16 @@ using DataDeps
 using Oceananigans.OutputReaders: GPUAdaptedFieldTimeSeries, FieldTimeSeries
 using Oceananigans.Grids: node
 
+# Does not really matter if there is no model
+reference_density(::Nothing) = 0
+heat_capacity(::Nothing) = 0
+
+reference_density(unsupported) =
+    throw(ArgumentError("Cannot extract reference density from $(typeof(unsupported))"))
+
+heat_capacity(unsupported) =
+    throw(ArgumentError("Cannot deduce the heat capacity from $(typeof(unsupported))"))
+
 const SomeKindOfFieldTimeSeries = Union{FieldTimeSeries,
                                         GPUAdaptedFieldTimeSeries}
 

src/DataWrangling/ECCO/ECCO_restoring.jl

@@ -138,15 +138,14 @@ function ECCOFieldTimeSeries(metadata::ECCOMetadata, grid::AbstractGrid;
     return fts	
 end
 
-function ECCOFieldTimeSeries(variable_name::Symbol; 
+function ECCOFieldTimeSeries(variable_name::Symbol, arch_or_grid=CPU();
                              dataset = ECCO4Monthly(),
-                             architecture = CPU(),
                              dates = all_dates(dataset, variable_name),
                              dir = download_ECCO_cache,
                              kw...)
 
     metadata = Metadata(variable_name, dates, dataset, dir)
-    return ECCOFieldTimeSeries(metadata, architecture; kw...)
+    return ECCOFieldTimeSeries(metadata, arch_or_grid; kw...)
 end
 
 # Variable names for restorable data

src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl

@@ -174,6 +174,7 @@ function compute_net_sea_ice_fluxes!(coupled_model)
     kernel_parameters = interface_kernel_parameters(grid)
 
     sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature
+    ice_concentration = sea_ice_concentration(sea_ice)
 
     launch!(arch, grid, kernel_parameters, 
             _assemble_net_sea_ice_fluxes!,
@@ -186,6 +187,7 @@ function compute_net_sea_ice_fluxes!(coupled_model)
             freshwater_flux,
             sea_ice_surface_temperature,
             downwelling_radiation,
+            ice_concentration,
             sea_ice_properties,
             atmos_sea_ice_properties)
 
@@ -201,6 +203,7 @@ end
                                                freshwater_flux, # Where do we add this one?
                                                surface_temperature,
                                                downwelling_radiation,
+                                               ice_concentration,
                                                sea_ice_properties,
                                                atmos_sea_ice_properties)
 
@@ -211,13 +214,13 @@ end
     @inbounds begin
         Ts = surface_temperature[i, j, kᴺ]
         Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts)
-
+        ℵi = ice_concentration[i, j, 1]
         Qs = downwelling_radiation.Qs[i, j, 1]
         Qℓ = downwelling_radiation.Qℓ[i, j, 1]
         Qc = atmosphere_sea_ice_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux
         Qv = atmosphere_sea_ice_fluxes.latent_heat[i, j, 1]   # latent heat flux
         Qf = sea_ice_ocean_fluxes.frazil_heat[i, j, 1]        # frazil heat flux
-        Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1]   # interfacial heat flux
+        Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1]     # interfacial heat flux
     end
 
     ρτx = atmosphere_sea_ice_fluxes.x_momentum  # zonal momentum flux                      
@@ -230,8 +233,8 @@ end
     Qu = upwelling_radiation(Ts, σ, ϵ)
     Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ)
 
-    ΣQt = Qd + Qu + Qc + Qv
-    ΣQb = Qf + Qi
+    ΣQt = (Qd + Qu + Qc + Qv) * ℵi # We need to multiply these times the concentration?
+    ΣQb = Qf + Qi 
 
     # Mask fluxes over land for convenience
     inactive = inactive_node(i, j, kᴺ, grid, c, c, c)

src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl

@@ -1,5 +1,6 @@
 using Oceananigans.Operators: intrinsic_vector
 using Oceananigans.Grids: inactive_node
+using ClimaOcean.OceanSimulations
 using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: thermodynamics_parameters,
                                                           surface_layer_height,
                                                           boundary_layer_height
@@ -70,7 +71,6 @@ end
 
     i, j = @index(Global, NTuple)
     kᴺ   = size(grid, 3) # index of the top ocean cell
-    time = Time(clock.time)
 
     @inbounds begin
         uₐ = atmosphere_state.u[i, j, 1]

src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl

@@ -116,7 +116,7 @@ end
                               h_bℓ = atmosphere_state.h_bℓ)
 
     downwelling_radiation = (; Qs, Qℓ)
-    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ)
+    local_interior_state = (u=uᵢ, v=vᵢ, T=Tᵢ, S=Sᵢ, h=hᵢ, ℵ=ℵᵢ)
 
     # Estimate initial interface state
     u★ = convert(FT, 1e-4)

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -245,7 +245,7 @@ function default_ao_specific_humidity(ocean)
 end
 
 """
-    ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
+    ComponentInterfaces(atmosphere, ocean, sea_ice;
                         radiation = Radiation(),
                         freshwater_density = 1000,
                         atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(),
@@ -262,11 +262,11 @@ end
                         sea_ice_reference_density = reference_density(sea_ice),
                         sea_ice_heat_capacity = heat_capacity(sea_ice))
 """
-function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
+function ComponentInterfaces(atmosphere, ocean, sea_ice;
                              radiation = Radiation(),
                              freshwater_density = 1000,
                              atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
-                             atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
+                             atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid), solver_maxiter=10000, stability_functions=atmosphere_sea_ice_stability_functions()),
                              atmosphere_ocean_interface_temperature = BulkTemperature(),
                              atmosphere_ocean_velocity_difference = RelativeVelocity(),
                              atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean),

src/OceanSeaIceModels/InterfaceComputations/interface_states.jl

@@ -212,20 +212,21 @@ end
     F = st.internal_flux
     k = F.conductivity
     h = Ψᵢ.h
+    ℵ = Ψᵢ.ℵ
 
     # Bottom temperature at the melting temperature
     Tᵢ  = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S)
     Tᵢ  = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ)
     Tₛ⁻ = Ψₛ.T
-
+    
     #=
     σ = ℙₛ.radiation.σ
     ϵ = ℙₛ.radiation.ϵ
     α = σ * ϵ
     Tₛ = (Tᵢ - h / k * (Qₐ + 4α * Tₛ⁻^4)) / (1 + 4α * h * Tₛ⁻^3 / k)
     =#
 
-    T★ = Tᵢ - Qₐ * h / k
+    T★ = Tᵢ - Qₐ * h / k * ℵ # We need to multiply the flux by the concentration?
 
     # Fix a NaN
     T★ = ifelse(isnan(T★), Tₛ⁻, T★)

src/OceanSeaIceModels/InterfaceComputations/radiation.jl

@@ -67,7 +67,7 @@ function Base.show(io::IO, r::Radiation)
     print(io, "├── stefan_boltzmann_constant: ", prettysummary(σ), '\n')
     print(io, "├── emission: ", summary(r.emission), '\n')
     print(io, "│   ├── ocean: ", prettysummary(r.emission.ocean), '\n')
-    print(io, "│   └── sea_ice: ", prettysummary(r.emission.ocean), '\n')
+    print(io, "│   └── sea_ice: ", prettysummary(r.emission.sea_ice), '\n')
     print(io, "└── reflection: ", summary(r.reflection), '\n')
     print(io, "    ├── ocean: ", prettysummary(r.reflection.ocean), '\n')
     print(io, "    └── sea_ice: ", prettysummary(r.reflection.sea_ice))

src/OceanSeaIceModels/OceanSeaIceModels.jl

@@ -22,14 +22,13 @@ using ClimaSeaIce: SeaIceModel
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
 
 using ClimaOcean: stateindex
-
 using KernelAbstractions: @kernel, @index
 using KernelAbstractions.Extras.LoopInfo: @unroll
 
+import ClimaOcean: reference_density, heat_capacity
+
 function downwelling_radiation end
 function freshwater_flux end
-function reference_density end
-function heat_capacity end
 
 const default_gravitational_acceleration = 9.80665
 const default_freshwater_density = 1000