Merge branch 'main' into ts/checkpointer-for-coupled-model

Author: navidcy

experiments/arctic_simulation.jl

@@ -6,12 +6,13 @@ using Oceananigans.Units
 using Oceananigans.OrthogonalSphericalShellGrids
 using ClimaOcean.OceanSimulations
 using ClimaOcean.ECCO
-using ClimaOcean.ECCO: all_ECCO_dates
+using ClimaOcean.DataWrangling
 using ClimaSeaIce.SeaIceThermodynamics: IceWaterThermalEquilibrium
 using Printf
 
 using CUDA
 CUDA.device!(1)
+arch = GPU()
 
 r_faces = ClimaOcean.exponential_z_faces(; Nz=30, h=10, depth=2000)
 z_faces = MutableVerticalDiscretization(r_faces)
@@ -20,13 +21,13 @@ Nx = 180 # longitudinal direction -> 250 points is about 1.5ᵒ resolution
 Ny = 180 # meridional direction -> same thing, 48 points is about 1.5ᵒ resolution
 Nz = length(r_faces) - 1
 
-grid = RotatedLatitudeLongitudeGrid(GPU(), size = (Nx, Ny, Nz), 
-                                           latitude = (-45, 45),
-                                           longitude = (-45, 45),
-                                           z = r_faces,
-                                           north_pole = (180, 0),
-                                           halo = (5, 5, 4),
-                                           topology = (Bounded, Bounded, Bounded))
+grid = RotatedLatitudeLongitudeGrid(arch, size = (Nx, Ny, Nz), 
+                                          latitude = (-45, 45),
+                                          longitude = (-45, 45),
+                                          z = r_faces,
+                                          north_pole = (180, 0),
+                                          halo = (5, 5, 4),
+                                          topology = (Bounded, Bounded, Bounded))
 
 bottom_height = regrid_bathymetry(grid; minimum_depth=15, major_basins=1)
 
@@ -49,27 +50,47 @@ ocean = ocean_simulation(grid;
                          free_surface,
                          closure)
 
-set!(ocean.model, T=ECCOMetadata(:temperature),
-                  S=ECCOMetadata(:salinity))
+dataset = ECCO4Monthly()
+
+set!(ocean.model, T=Metadata(:temperature; dataset),
+                  S=Metadata(:salinity;    dataset))
 
 #####
 ##### A Prognostic Sea-ice model
 #####
 
+using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.Rheologies
+
 # Remember to pass the SSS as a bottom bc to the sea ice!
 SSS = view(ocean.model.tracers.S, :, :, grid.Nz)
 bottom_heat_boundary_condition = IceWaterThermalEquilibrium(SSS)
 
-sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics=nothing, advection=nothing) 
+SSU = view(ocean.model.velocities.u, :, :, grid.Nz)
+SSV = view(ocean.model.velocities.u, :, :, grid.Nz)
+
+τo  = SemiImplicitStress(uₑ=SSU, vₑ=SSV)
+τua = Field{Face, Center, Nothing}(grid) 
+τva = Field{Center, Face, Nothing}(grid)
+
+dynamics = SeaIceMomentumEquation(grid; 
+                                  coriolis = ocean.model.coriolis,
+                                  top_momentum_stress = (u=τua, v=τva),
+                                  bottom_momentum_stress = τo,
+                                  ocean_velocities = (u=SSU, v=SSV),
+                                  rheology = ElastoViscoPlasticRheology(),
+                                  solver = SplitExplicitSolver(120))
+
+sea_ice = sea_ice_simulation(grid; bottom_heat_boundary_condition, dynamics, advection=WENO(order=7)) 
 
-set!(sea_ice.model, h=ECCOMetadata(:sea_ice_thickness), 
-                    ℵ=ECCOMetadata(:sea_ice_concentration))
+set!(sea_ice.model, h=Metadata(:sea_ice_thickness;     dataset), 
+                    ℵ=Metadata(:sea_ice_concentration; dataset))
 
 ##### 
 ##### A Prescribed Atmosphere model
 #####
 
-atmosphere = JRA55PrescribedAtmosphere(GPU(); backend=JRA55NetCDFBackend(40))
+atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(40))
 radiation  = Radiation()
 
 #####
@@ -80,10 +101,10 @@ arctic = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
 arctic = Simulation(arctic, Δt=5minutes, stop_time=365days)
 
 # Sea-ice variables
-h  = sea_ice.model.ice_thickness
-ℵ  = sea_ice.model.ice_concentration
-Gh = sea_ice.model.timestepper.Gⁿ.h
-Gℵ = sea_ice.model.timestepper.Gⁿ.ℵ
+h = sea_ice.model.ice_thickness
+ℵ = sea_ice.model.ice_concentration
+u = sea_ice.model.velocities.u
+v = sea_ice.model.velocities.v
 
 # Fluxes
 Tu = arctic.model.interfaces.atmosphere_sea_ice_interface.temperature
@@ -93,14 +114,16 @@ 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
+τx = arctic.model.interfaces.net_fluxes.sea_ice_top.u
+τy = arctic.model.interfaces.net_fluxes.sea_ice_top.v
 
 # Output writers
-arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Gh, Gℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ),
+arctic.output_writers[:vars] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, u, v, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, τx, τy),
                                                  filename = "sea_ice_quantities.jld2",
                                                  schedule = IterationInterval(12),
                                                  overwrite_existing=true)
 
-arctic.output_writers[:avrages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ),
+arctic.output_writers[:averages] = JLD2OutputWriter(sea_ice.model, (; h, ℵ, Tu, Qˡ, Qˢ, Qⁱ, Qᶠ, Qᵗ, Qᴮ, u, v, τx, τy),
                                                     filename = "averaged_sea_ice_quantities.jld2",
                                                     schedule = AveragedTimeInterval(1days),
                                                     overwrite_existing=true)
@@ -142,11 +165,11 @@ run!(arctic)
 ##### Comparison to ECCO Climatology
 #####
 
-version = ECCO4Monthly()
-dates   = all_ECCO_dates(version)[1:12]
+dataset = ECCO4Monthly()
+dates   = all_dates(version)[1:12]
 
-h_metadata = ECCOMetadata(:sea_ice_thickness;     version, dates)
-ℵ_metadata = ECCOMetadata(:sea_ice_concentration; version, dates)
+h_metadata = Metadata(:sea_ice_thickness;     dataset, dates)
+ℵ_metadata = Metadata(:sea_ice_concentration; dataset, dates)
 
 # Montly averaged ECCO data
 hE = ECCOFieldTimeSeries(h_metadata, grid; time_indices_in_memory=12)

src/DataWrangling/ECCO/ECCO_metadata.jl

@@ -168,8 +168,6 @@ function download_dataset(metadata::ECCOMetadata; url = urls(metadata))
         # Write down the username and password in a .netrc file
         downloader = netrc_downloader(username, password, "ecco.jpl.nasa.gov", tmp)
         ntasks = Threads.nthreads()
-       
-        missing_files = false
 
         asyncmap(metadata; ntasks) do metadatum # Distribute the download among tasks
 
@@ -194,10 +192,6 @@ function download_dataset(metadata::ECCOMetadata; url = urls(metadata))
                 Downloads.download(fileurl, filepath; downloader, progress=download_progress)
             end
         end
-        
-	    if !missing_files
-            @info "Note: ECCO $(metadata.name) data is in $(metadata.dir)."
-	    end
     end
 
     return nothing

src/OceanSeaIceModels/InterfaceComputations/assemble_net_fluxes.jl

@@ -154,7 +154,7 @@ function compute_net_sea_ice_fluxes!(coupled_model)
     arch  = architecture(grid)
     clock = coupled_model.clock
 
-    top_heat_flux = coupled_model.interfaces.net_fluxes.sea_ice_top.heat
+    top_fluxes = coupled_model.interfaces.net_fluxes.sea_ice_top
     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
@@ -173,11 +173,11 @@ function compute_net_sea_ice_fluxes!(coupled_model)
 
     kernel_parameters = interface_kernel_parameters(grid)
 
-    sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_ocean_interface.temperature
+    sea_ice_surface_temperature = coupled_model.interfaces.atmosphere_sea_ice_interface.temperature
 
     launch!(arch, grid, kernel_parameters, 
             _assemble_net_sea_ice_fluxes!,
-            top_heat_flux,
+            top_fluxes,
             bottom_heat_flux, 
             grid,
             clock,
@@ -192,7 +192,7 @@ function compute_net_sea_ice_fluxes!(coupled_model)
     return nothing
 end
 
-@kernel function _assemble_net_sea_ice_fluxes!(top_heat_flux,
+@kernel function _assemble_net_sea_ice_fluxes!(top_fluxes,
                                                bottom_heat_flux, 
                                                grid,
                                                clock,
@@ -220,6 +220,9 @@ end
         Qi = sea_ice_ocean_fluxes.interface_heat[i, j, 1]   # interfacial heat flux
     end
 
+    ρτx = atmosphere_sea_ice_fluxes.x_momentum  # zonal momentum flux                      
+    ρτy = atmosphere_sea_ice_fluxes.y_momentum  # meridional momentum flux
+
     # Compute radiation fluxes
     σ = atmos_sea_ice_properties.radiation.σ
     α = stateindex(atmos_sea_ice_properties.radiation.α, i, j, kᴺ, grid, time)
@@ -233,6 +236,8 @@ end
     # Mask fluxes over land for convenience
     inactive = inactive_node(i, j, kᴺ, grid, c, c, c)
 
-    @inbounds top_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQt)
+    @inbounds top_fluxes.heat[i, j, 1]  = ifelse(inactive, zero(grid), ΣQt)
+    @inbounds top_fluxes.u[i, j, 1]     = ifelse(inactive, zero(grid), ℑxᶠᵃᵃ(i, j, 1, grid, ρτx))
+    @inbounds top_fluxes.v[i, j, 1]     = ifelse(inactive, zero(grid), ℑyᵃᶠᵃ(i, j, 1, grid, ρτy))
     @inbounds bottom_heat_flux[i, j, 1] = ifelse(inactive, zero(grid), ΣQb)
 end

src/OceanSeaIceModels/InterfaceComputations/component_interfaces.jl

@@ -266,9 +266,7 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
                              radiation = Radiation(),
                              freshwater_density = 1000,
                              atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
-                             atmosphere_sea_ice_flux_formulation = CoefficientBasedFluxes(drag_coefficient=2e-3,
-                                                                                          heat_transfer_coefficient=1e-4,
-                                                                                          vapor_flux_coefficient=1e-4),
+                             atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
                              atmosphere_ocean_interface_temperature = BulkTemperature(),
                              atmosphere_ocean_velocity_difference = RelativeVelocity(),
                              atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean),
@@ -319,7 +317,17 @@ function ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
                               liquidus           = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus,
                               temperature_units  = sea_ice_temperature_units)
 
-        net_top_sea_ice_fluxes = (; heat=sea_ice.model.external_heat_fluxes.top)
+        net_momentum_fluxes = if sea_ice.model.dynamics isa Nothing 
+            u = Field{Face, Center, Nothing}(sea_ice.model.grid)
+            v = Field{Center, Face, Nothing}(sea_ice.model.grid)
+            (; u, v) 
+        else
+            u = sea_ice.model.dynamics.external_momentum_stresses.top.u
+            v = sea_ice.model.dynamics.external_momentum_stresses.top.v
+            (; u, v)
+        end
+
+        net_top_sea_ice_fluxes = merge((; heat=sea_ice.model.external_heat_fluxes.top), net_momentum_fluxes)
         net_bottom_sea_ice_fluxes = (; heat=sea_ice.model.external_heat_fluxes.bottom)
     else
         sea_ice_properties = nothing

src/OceanSeaIceModels/InterfaceComputations/interface_states.jl

@@ -214,8 +214,8 @@ end
     h = Ψᵢ.h
 
     # Bottom temperature at the melting temperature
-    Tᵢ = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S)
-    Tᵢ = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ)
+    Tᵢ  = ClimaSeaIce.SeaIceThermodynamics.melting_temperature(ℙᵢ.liquidus, Ψᵢ.S)
+    Tᵢ  = convert_to_kelvin(ℙᵢ.temperature_units, Tᵢ)
     Tₛ⁻ = Ψₛ.T
 
     #=
@@ -227,10 +227,6 @@ end
 
     T★ = Tᵢ - Qₐ * h / k
 
-    # Under heating fluxes, cap surface temperature by melting temperature
-    Tₘ = ℙᵢ.liquidus.freshwater_melting_temperature
-    Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ)
-
     # Fix a NaN
     T★ = ifelse(isnan(T★), Tₛ⁻, T★)
 
@@ -241,6 +237,11 @@ end
     abs_ΔT = min(max_ΔT, abs(ΔT★))
     Tₛ⁺ = Tₛ⁻ + abs_ΔT * sign(ΔT★)
 
+    # Under heating fluxes, cap surface temperature by melting temperature
+    Tₘ = ℙᵢ.liquidus.freshwater_melting_temperature
+    Tₘ = convert_to_kelvin(ℙᵢ.temperature_units, Tₘ)
+    Tₛ⁺ = min(Tₛ⁺, Tₘ)
+
     return Tₛ⁺
 end
 

src/OceanSimulations/barotropic_potential_forcing.jl

@@ -25,11 +25,11 @@ forcing_barotropic_potential(something) = nothing
 forcing_barotropic_potential(f::BarotropicPotentialForcing) = f.potential.data
 
 function forcing_barotropic_potential(mf::MultipleForcings)
-    n = findfirst(f -> f isa BarotropicPotentialForcing, mf.forcing)
+    n = findfirst(f -> f isa BarotropicPotentialForcing, mf.forcings)
     if isnothing(n)
         return nothing
     else
-        return forcing_barotropic_potential(mf.forcing[n])
+        return forcing_barotropic_potential(mf.forcings[n])
     end
 end
 

test/test_ocean_sea_ice_model.jl

@@ -4,6 +4,8 @@ using CUDA
 using Oceananigans.OrthogonalSphericalShellGrids
 using ClimaOcean.OceanSeaIceModels: above_freezing_ocean_temperature!
 using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
+using ClimaSeaIce.SeaIceMomentumEquations
+using ClimaSeaIce.Rheologies
 
 @inline kernel_melting_temperature(i, j, k, grid, liquidus, S) = @inbounds melting_temperature(liquidus, S[i, j, k])
 
@@ -66,10 +68,17 @@ using ClimaSeaIce.SeaIceThermodynamics: melting_temperature
         ##### Coupled ocean-sea ice and prescribed atmosphere
         #####
 
-        sea_ice_grid = TripolarGrid(arch; size=(50, 50, 1), z = (-10, 0))
-        sea_ice_grid = ImmersedBoundaryGrid(sea_ice_grid, GridFittedBottom(bottom_height))
+        # Adding a sea ice model to the coupled model
+        τua = Field{Face, Center, Nothing}(grid) 
+        τva = Field{Center, Face, Nothing}(grid)
 
-        sea_ice  = sea_ice_simulation(sea_ice_grid)
+        dynamics = SeaIceMomentumEquation(grid; 
+                                          coriolis = ocean.model.coriolis,
+                                          top_momentum_stress = (u=τua, v=τva),
+                                          rheology = ElastoViscoPlasticRheology(),
+                                          solver = SplitExplicitSolver(120))
+
+        sea_ice  = sea_ice_simulation(grid; dynamics, advection=WENO(order=7)) 
         liquidus = sea_ice.model.ice_thermodynamics.phase_transitions.liquidus
 
         # Set the ocean temperature and salinity