Continue sea-ice coupling effort (#365)

* this works best?

* One degree simulation take 4

* add a biharmonic viscosity

* infer maximum delta t

* changes

* add statistics

* infer -> compute

* compute -> infer

* default for distribtued grid

* all_reduce instead of allreduce

* update other packages

* let's see if it works

* try it out

* run simulation

* remove all instances of diffusion

* try a zstar grid

* need a manifest for this

* regrid the bathymetry

* better

* make it work

* Update Project.toml

* fix docs

* fix docs

* update project

* just remove it for now

* try out correct settings

* Update examples/one_degree_simulation.jl

Co-authored-by: Jonathan Lauderdale <[email protected]>

* should work now

* literate one-degree example

* use GPU

* nuke znode

* add CairoMakie 0.13 compat

* explain that there is also e

* add comment for OrthogonalSphericalShellGrids

* add comment for coupled ocean--sea ice

* nuke ExplicitTimeDiscretization and fix typo

* this should work

* add another restoring file

* correct videos

* tamper with the time step

* bugfix

* remove for the moment

* fix examples

* already some bugfix

* this should work?

* go ahead

* some changes

* some more changes

* some fixes

* this should run

* we do not need to remove fluxes as we adjust the ocean interior

* don't need to do that

* also don't need this

* adjust temperature

* remove unused stuff

* remove other stuff

* bugfixes

* changes

* add some tests

* remo this line

* fix tests

* add frazil ice formation

* try it like this

* this should work

* improve model

* this works now?

* this should work

* add default inpainting

* better set! function

* small bug

* bugfix

* bugfix

* correct fluxes

* add the net fluxes

* interface heat

* Update coefficient_based_turbulent_fluxes.jl

* some changes

* add a bottom heat BC

* add a bottom heat BC

* tests should pass

---------

Co-authored-by: Navid C. Constantinou <[email protected]>
Co-authored-by: Jonathan Lauderdale <[email protected]>

Author: 3 people

experiments/coupled_simulation/ocean_sea_ice_coupled_simulation.jl

@@ -88,9 +88,7 @@ set!(ocean.model, T=temperature, S=salinity)
 set!(sea_ice.model.ice_thickness,     ice_thickness,     inpainting=NearestNeighborInpainting(1))
 set!(sea_ice.model.ice_concentration, ice_concentration, inpainting=NearestNeighborInpainting(1))
 
-adjust_ocean_temperature!(ocean, sea_ice)
 earth_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
-
 earth = Simulation(earth_model; Δt=30minutes, stop_iteration=10, stop_time=30days)
 
 u, v, _ = ocean.model.velocities

src/DataWrangling/ECCO/ECCO.jl

@@ -155,7 +155,7 @@ function ECCO_field(metadata::ECCOMetadata;
                     mask = nothing,
                     horizontal_halo = (7, 7),
                     cache_inpainted_data = true)
-
+                    
     field = empty_ECCO_field(metadata; architecture, horizontal_halo)
     inpainted_path = inpainted_metadata_path(metadata)
 

src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl

@@ -4,7 +4,7 @@ using Oceananigans.Operators: ℑxᶠᵃᵃ, ℑyᵃᶠᵃ
 using ClimaOcean.OceanSeaIceModels: sea_ice_concentration
 
 @inline computed_sea_ice_ocean_fluxes(interface) = interface.fluxes
-@inline computed_sea_ice_ocean_fluxes(::Nothing) = (heat = ZeroField(), salt = ZeroField())
+@inline computed_sea_ice_ocean_fluxes(::Nothing) = (interface_heat = ZeroField(), frazil_heat = ZeroField(), salt = ZeroField())
 
 function compute_net_ocean_fluxes!(coupled_model)
     ocean = coupled_model.ocean
@@ -131,7 +131,7 @@ end
 
     @inbounds begin
         ℵᵢ   = ℵ[i, j, 1]
-        Qio  = sea_ice_ocean_fluxes.heat[i, j, 1]
+        Qio  = sea_ice_ocean_fluxes.interface_heat[i, j, 1]
         Jˢio = sea_ice_ocean_fluxes.salt[i, j, 1]
         Jᵀio = Qio * ρₒ⁻¹ / cₒ
 
@@ -142,19 +142,22 @@ end
     end
 end
 
-compute_net_sea_ice_fluxes!(coupled_model) = nothing
-
-#=
 function compute_net_sea_ice_fluxes!(coupled_model)
-    ocean = coupled_model.ocean
     sea_ice = coupled_model.sea_ice
-    grid = ocean.model.grid
-    arch = architecture(grid)
+
+    if !(sea_ice isa SeaIceSimulation)
+        return nothing
+    end
+
+    ocean = coupled_model.ocean
+    grid  = ocean.model.grid
+    arch  = architecture(grid)
     clock = coupled_model.clock
 
-    top_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_top
-    atmos_ocean_fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes
+    top_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_top.heat
+    bottom_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_bottom.heat
     sea_ice_ocean_fluxes = coupled_model.interfaces.sea_ice_ocean_interface.fluxes
+    atmosphere_sea_ice_fluxes = coupled_model.interfaces.atmosphere_sea_ice_interface.fluxes
 
     # Simplify NamedTuple to reduce parameter space consumption.
     # See https://github.com/CliMA/ClimaOcean.jl/issues/116.
@@ -165,69 +168,71 @@ function compute_net_sea_ice_fluxes!(coupled_model)
 
     freshwater_flux = atmosphere_fields.Mp.data
 
-    sea_ice_concentration = sea_ice.model.ice_concentration
-    ocean_salinity = ocean.model.tracers.S
-    atmos_ocean_properties = coupled_model.interfaces.atmosphere_ocean_interface.properties
-    ocean_properties = coupled_model.interfaces.ocean_properties
+    atmos_sea_ice_properties = coupled_model.interfaces.atmosphere_sea_ice_interface.properties
+    sea_ice_properties = coupled_model.interfaces.sea_ice_properties
+
     kernel_parameters = interface_kernel_parameters(grid)
 
-    ocean_surface_temperature = coupled_model.interfaces.atmosphere_ocean_interface.temperature
+    sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_ocean_interface.temperature
 
     launch!(arch, grid, kernel_parameters, 
             _assemble_net_sea_ice_fluxes!,
-            net_ocean_fluxes,
+            top_heat_flux,
+            bottom_heat_flux, 
             grid,
             clock,
-            atmos_ocean_fluxes,
+            atmosphere_sea_ice_fluxes,
             sea_ice_ocean_fluxes,
-            ocean_salinity,
-            ocean_surface_temperature,
-            sea_ice_concentration,
-            downwelling_radiation,
             freshwater_flux,
-            atmos_ocean_properties,
-            ocean_properties)
+            sea_ice_surface_temperature,
+            downwelling_radiation,
+            sea_ice_properties,
+            atmos_sea_ice_properties)
 
     return nothing
 end
 
-@kernel function _assemble_atmosphere_sea_ice_fluxes!(net_heat_flux,
-                                                      grid,
-                                                      clock,
-                                                      surface_temperature,
-                                                      surface_temperature_units,
-                                                      turbulent_fluxes,
-                                                      downwelling_radiation,
-                                                      stefan_boltzmann_constant,
-                                                      albedo,
-                                                      emissivity)
+@kernel function _assemble_net_sea_ice_fluxes!(top_heat_flux,
+                                               bottom_heat_flux, 
+                                               grid,
+                                               clock,
+                                               atmosphere_sea_ice_fluxes,
+                                               sea_ice_ocean_fluxes,
+                                               freshwater_flux, # Where do we add this one?
+                                               surface_temperature,
+                                               downwelling_radiation,
+                                               sea_ice_properties,
+                                               atmos_sea_ice_properties)
 
     i, j = @index(Global, NTuple)
     kᴺ = size(grid, 3)
     time = Time(clock.time)
-
+    
     @inbounds begin
         Ts = surface_temperature[i, j, kᴺ]
-        Ts = convert_to_kelvin(surface_temperature_units, Ts)
-
-        Qs  = downwelling_radiation.shortwave[i, j, 1]
-        Qℓ  = downwelling_radiation.longwave[i, j, 1]
-        Qc  = turbulent_fluxes.sensible_heat[i, j, 1] # sensible or "conductive" heat flux
-        Qv  = turbulent_fluxes.latent_heat[i, j, 1]   # latent heat flux
+        Ts = convert_to_kelvin(sea_ice_properties.temperature_units, Ts)
+
+        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
     end
 
     # Compute radiation fluxes
-    σ = stefan_boltzmann_constant
-    α = stateindex(albedo, i, j, 1, grid, time)
-    ϵ = stateindex(emissivity, i, j, 1, grid, time)
+    σ = atmos_sea_ice_properties.radiation.σ
+    α = stateindex(atmos_sea_ice_properties.radiation.α, i, j, kᴺ, grid, time)
+    ϵ = stateindex(atmos_sea_ice_properties.radiation.ϵ, i, j, kᴺ, grid, time)
     Qu = upwelling_radiation(Ts, σ, ϵ)
     Qd = net_downwelling_radiation(i, j, grid, time, α, ϵ, Qs, Qℓ)
 
-    ΣQ = Qd + Qu + Qc + Qv
+    ΣQt = Qd + Qu + Qc + Qv
+    ΣQb = Qf + Qi
 
     # Mask fluxes over land for convenience
     inactive = inactive_node(i, j, kᴺ, grid, c, c, c)
 
-    @inbounds net_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQ)
+    @inbounds top_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQt)
+    @inbounds bottom_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQb)
 end
-=#

src/OceanSeaIceModels/InterfaceComputations/atmosphere_sea_ice_fluxes.jl

@@ -33,7 +33,6 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model)
 
     flux_formulation = coupled_model.interfaces.atmosphere_sea_ice_interface.flux_formulation
     interface_fluxes = coupled_model.interfaces.atmosphere_sea_ice_interface.fluxes
-    net_top_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_top.heat
     interface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature
     interface_properties = coupled_model.interfaces.atmosphere_sea_ice_interface.properties
     sea_ice_properties = coupled_model.interfaces.sea_ice_properties
@@ -48,7 +47,6 @@ function compute_atmosphere_sea_ice_fluxes!(coupled_model)
             _compute_atmosphere_sea_ice_interface_state!,
             interface_fluxes,
             interface_temperature,
-            net_top_heat_flux,
             grid,
             clock,
             flux_formulation,
@@ -65,7 +63,6 @@ end
 """ Compute turbulent fluxes between an atmosphere and a interface state using similarity theory """
 @kernel function _compute_atmosphere_sea_ice_interface_state!(interface_fluxes,
                                                               interface_temperature,
-                                                              net_top_heat_flux,
                                                               grid,
                                                               clock,
                                                               turbulent_flux_formulation,
@@ -151,13 +148,6 @@ end
     u★ = interface_state.u★
     θ★ = interface_state.θ★
     q★ = interface_state.q★
-
-    #=
-    Pr = similarity_theory.turbulent_prandtl_number
-    θ★ = θ★ / Pr
-    q★ = q★ / Pr
-    =#
-
     Ψₛ = interface_state
     Ψₐ = local_atmosphere_state
     Δu, Δv = velocity_difference(turbulent_flux_formulation.bulk_velocity, Ψₐ, Ψₛ)
@@ -169,26 +159,18 @@ end
     cₚ = AtmosphericThermodynamics.cp_m(ℂₐ, 𝒬ₐ) # moist heat capacity
     ℰs = AtmosphericThermodynamics.latent_heat_sublim(ℂₐ, 𝒬ₐ)
 
-    σ = interface_properties.radiation.σ
-    α = stateindex(interface_properties.radiation.α, i, j, 1, grid, time)
-    ϵ = stateindex(interface_properties.radiation.ϵ, i, j, 1, grid, time)
-    Qu = upwelling_radiation(Ψₛ.T, σ, ϵ)
-    Qd = net_downwelling_radiation(downwelling_radiation, α, ϵ)
-
     # Store fluxes
     Qv = interface_fluxes.latent_heat
     Qc = interface_fluxes.sensible_heat
     Fv = interface_fluxes.water_vapor
     ρτx = interface_fluxes.x_momentum
     ρτy = interface_fluxes.y_momentum
     Ts = interface_temperature
-    ΣQ = net_top_heat_flux
 
     @inbounds begin
         # +0: cooling, -0: heating
-        Qv[i, j, 1]  = _Qv = - ρₐ * u★ * q★ * ℰs 
-        Qc[i, j, 1]  = _Qc = - ρₐ * cₚ * u★ * θ★ 
-        ΣQ[i, j, 1]  = Qu + Qd + _Qv + _Qc
+        Qv[i, j, 1]  = - ρₐ * u★ * q★ * ℰs 
+        Qc[i, j, 1]  = - ρₐ * cₚ * u★ * θ★ 
         Fv[i, j, 1]  = - ρₐ * u★ * q★ 
         ρτx[i, j, 1] = + ρₐ * τx 
         ρτy[i, j, 1] = + ρₐ * τy 

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -17,10 +17,9 @@ using ClimaSeaIce: SeaIceModel
 using Oceananigans: HydrostaticFreeSurfaceModel, architecture
 using Oceananigans.Grids: inactive_node, node
 using Oceananigans.BoundaryConditions: fill_halo_regions!
-using Oceananigans.Fields: ConstantField, interpolate, FractionalIndices
-using Oceananigans.Utils: launch!, Time
 
-# using Oceananigans.OutputReaders: extract_field_time_series, update_field_time_series!
+using Oceananigans.Fields: ConstantField, interpolate, FractionalIndices
+using Oceananigans.Utils: launch!, Time, KernelParameters
 
 using Oceananigans.Operators: ℑxᶜᵃᵃ, ℑyᵃᶜᵃ, ℑxᶠᵃᵃ, ℑyᵃᶠᵃ
 
@@ -179,17 +178,24 @@ end
 
 sea_ice_ocean_interface(sea_ice, ocean) = nothing
 
-function sea_ice_ocean_interface(sea_ice::SeaIceSimulation, ocean)
+function sea_ice_ocean_interface(sea_ice::SeaIceSimulation, ocean;
+                                 characteristic_melting_speed = 1e-5)
+
     previous_ice_thickness = deepcopy(sea_ice.model.ice_thickness)
     previous_ice_concentration = deepcopy(sea_ice.model.ice_concentration)
-    io_heat_flux = sea_ice.model.external_heat_fluxes.bottom
+    io_bottom_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid)
+    io_frazil_heat_flux = Field{Center, Center, Nothing}(ocean.model.grid)
     io_salt_flux = Field{Center, Center, Nothing}(ocean.model.grid)
 
-    @assert io_heat_flux isa Field{Center, Center, Nothing}
+    @assert io_frazil_heat_flux isa Field{Center, Center, Nothing}
+    @assert io_bottom_heat_flux isa Field{Center, Center, Nothing}
     @assert io_salt_flux isa Field{Center, Center, Nothing}
 
-    io_fluxes = (heat=io_heat_flux, salt=io_salt_flux)
-    io_properties = nothing
+    io_fluxes = (interface_heat=io_bottom_heat_flux, 
+                 frazil_heat=io_frazil_heat_flux,
+                 salt=io_salt_flux)
+
+    io_properties = (; characteristic_melting_speed)
 
     return SeaIceOceanInterface(io_fluxes,
                                 io_properties,