A low resolution coupled simulation

docs/Project.toml

@@ -9,9 +9,9 @@ OrthogonalSphericalShellGrids = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0"
 SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
 
 [compat]
-CairoMakie = "0.10.12, 0.11, 0.12"
+CairoMakie = "0.10.12, 0.11, 0.12, 0.13"
 DataDeps = "0.7"
 Documenter = "1"
 Oceananigans = "0.95.7 - 0.99"
 OrthogonalSphericalShellGrids = "0.2.1"
-SeawaterPolynomials = "0.3.4"
+SeawaterPolynomials = "0.3.5"

docs/make.jl

@@ -17,6 +17,7 @@ to_be_literated = [
     "generate_bathymetry.jl",
     "generate_surface_fluxes.jl",
     "single_column_os_papa_simulation.jl",
+    "one_degree_simulation.jl",
     # "mediterranean_simulation_with_ecco_restoring.jl",
     "near_global_ocean_simulation.jl"
 ]
@@ -46,6 +47,7 @@ pages = [
         "Surface fluxes" => "literated/generate_surface_fluxes.md",
         "Single-column simulation" => "literated/single_column_os_papa_simulation.md",
         # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md",
+        "One-degree Ocean simulation" => "literated/one_degree_simulation.md",
         "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md",
         ],
 

examples/ecco_mixed_layer_depth.jl

@@ -11,8 +11,8 @@ using SeawaterPolynomials: TEOS10EquationOfState
 using Oceananigans.BuoyancyFormulations: buoyancy
 
 arch = CPU()
-Nx = 360 ÷ 1
-Ny = 160 ÷ 1
+Nx = 360 
+Ny = 160 
 
 z = ClimaOcean.DataWrangling.ECCO.ECCO_z
 z = z[20:end]

examples/generate_surface_fluxes.jl

@@ -65,10 +65,10 @@ set!(ocean.model; T=T_metadata, S=S_metadata)
 
 coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation=Radiation())
 
-# Now that the surface fluxes are computed, we can extract and visualize them.
-# The turbulent fluxes are stored in `coupled_model.fluxes.turbulent`.
+# # Now that the surface fluxes are computed, we can extract and visualize them.
+# # The turbulent fluxes are stored in `coupled_modelinterfaces.atmosphere_ocean_interface.fluxes`.
 
-fluxes  = coupled_model.fluxes.turbulent.fields
+fluxes  = coupled_modelinterfaces.atmosphere_ocean_interface.fluxes
 λ, φ, z = nodes(fluxes.sensible_heat)
 
 fig = Figure(size = (800, 800), fontsize = 15)
@@ -89,4 +89,4 @@ ax = Axis(fig[3, 1], title = "Water vapor flux (kg m⁻² s⁻¹)", xlabel = "Lo
 heatmap!(ax, λ, φ, interior(fluxes.water_vapor, :, :, 1); colormap = :bwr)
 
 save("fluxes.png", fig)
-# ![](fluxes.png)
+![](fluxes.png)

examples/one_degree_simulation.jl

@@ -0,0 +1,255 @@
+# # One-degree global ocean simulation
+#
+# This example configures a global ocean--sea ice simulation at 1ᵒ horizontal resolution with
+# realistic bathymetry and some closures.
+#
+# For this example, we need Oceananigans, ClimaOcean, OrthogonalSphericalShellGrids, and
+# CairoMakie to visualize the simulation. Also we need CFTime and Dates for date handling.
+
+using ClimaOcean
+using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting
+using Oceananigans
+using Oceananigans.Units
+using OrthogonalSphericalShellGrids
+using CFTime
+using Dates
+using Printf
+
+arch = CPU()
+
+# ### Grid and Bathymetry
+
+Nx = 360
+Ny = 180
+Nz = 100
+
+r_faces = exponential_z_faces(; Nz, depth=5000, h=34)
+z_faces = Oceananigans.MutableVerticalDiscretization(r_faces)
+
+underlying_grid = TripolarGrid(arch;
+                               size = (Nx, Ny, Nz),
+                               z = z_faces,
+                               halo = (5, 5, 4),
+                               first_pole_longitude = 70,
+                               north_poles_latitude = 55)
+
+bottom_height = regrid_bathymetry(underlying_grid;
+                                  minimum_depth = 10,
+                                  interpolation_passes = 75,
+                                  major_basins = 2)
+
+# For this bathymetry at this horizontal resolution we need to manually open the Gibraltar strait.
+tampered_bottom_height = deepcopy(bottom_height)
+view(tampered_bottom_height, 102:103, 124, 1) .= -400
+
+grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_height);
+                            active_cells_map=true)
+
+# ### Restoring
+
+# We include temperature and salinity surface restoring to ECCO2.
+restoring_rate  = 1 / 10days
+z_below_surface = r_faces[end-1]
+
+mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0))
+
+dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1994, 1, 1)
+temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./")
+salinity    = ECCOMetadata(:salinity;    dates, version=ECCO4Monthly(), dir="./")
+
+FT = ECCORestoring(temperature, grid; mask, rate=restoring_rate)
+FS = ECCORestoring(salinity,    grid; mask, rate=restoring_rate)
+forcing = (T=FT, S=FS)
+
+# ### Closures
+
+# We include a Gent-McWilliam isopycnal diffusivity as a parameterization for the mesoscale
+# eddy fluxes. For vertical mixing at the upper-ocean boundary layer we include the CATKE
+# parameterization. We also include some explicit horizontal diffusivity.
+
+using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity,
+                                       DiffusiveFormulation
+
+eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3,
+                                                 skew_flux_formulation=DiffusiveFormulation())
+vertical_mixing = ClimaOcean.OceanSimulations.default_ocean_closure()
+
+closure = (eddy_closure, vertical_mixing)
+
+# ### Ocean simulation
+
+# Now we bring everything together to construct the ocean simulation.
+# We use a split-explicit timestepping with 30 substeps for the barotropic
+# mode.
+
+free_surface = SplitExplicitFreeSurface(grid; substeps=30)
+
+momentum_advection = WENOVectorInvariant(vorticity_order=3)
+tracer_advection   = Centered()
+
+ocean = ocean_simulation(grid;
+                         momentum_advection,
+                         tracer_advection,
+                         closure,
+                         forcing,
+                         free_surface)
+
+# ### Initial condition
+
+# We initialize the ocean from the ECCO2 state estimate.
+
+set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
+                  S=ECCOMetadata(:salinity; dates=first(dates)))
+
+# ### Atmospheric forcing
+
+# We force the simulation with an JRA55-do atmospheric reanalysis.
+radiation  = Radiation(arch)
+atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20))
+
+# ### Coupled simulation
+
+# Now we are ready to build the coupled ocean--sea ice model and bring everything
+# together into a `simulation`.
+
+# We use a relatively short time step initially and only run for a few days to
+# avoid numerical instabilities from the initial "shock" of the adjustment of the
+# flow fields.
+
+coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
+simulation = Simulation(coupled_model; Δt=1minutes, stop_time=10days)
+
+# ### A progress messenger
+#
+# We write a function that prints out a helpful progress message while the simulation runs.
+
+wall_time = Ref(time_ns())
+
+function progress(sim)
+    ocean = sim.model.ocean
+    u, v, w = ocean.model.velocities
+    T = ocean.model.tracers.T
+    Tmax = maximum(interior(T))
+    Tmin = minimum(interior(T))
+    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, iteration: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s \n",
+                   prettytime(sim), iteration(sim), prettytime(sim.Δt),
+                   umax..., Tmax, Tmin, prettytime(step_time))
+
+     wall_time[] = time_ns()
+
+     return nothing
+end
+
+# And add it as a callback to the simulation.
+add_callback!(simulation, progress, IterationInterval(10))
+
+# ### Output
+#
+# We are almost there! We need to save some output. Below we choose to save _only surface_
+# fields using the `indices` keyword argument. We save all velocity and tracer components.
+# Note, that besides temperature and salinity, the CATKE vertical mixing parameterization
+# also uses a prognostic turbulent kinetic energy, `e`, to diagnose the vertical mixing length.
+
+outputs = merge(ocean.model.tracers, ocean.model.velocities)
+ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs;
+                                                  schedule = TimeInterval(5days),
+                                                  filename = "global_surface_fields",
+                                                  indices = (:, :, grid.Nz),
+                                                  with_halos = true,
+                                                  overwrite_existing = true,
+                                                  array_type = Array{Float32})
+
+# ### Ready to run
+
+# We are ready to press the big red button and run the simulation.
+
+# After we run for a short time (here we set up the simulation with `stop_time = 10days`),
+# we increase the timestep and run for longer.
+
+run!(simulation)
+
+simulation.Δt = 20minutes
+simulation.stop_time = 360days
+
+run!(simulation)
+
+# ### A pretty movie
+#
+# We load the saved output and make a pretty movie of the simulation. First we plot a snapshot:
+using CairoMakie
+
+u = FieldTimeSeries("global_surface_fields.jld2", "u"; backend = OnDisk())
+v = FieldTimeSeries("global_surface_fields.jld2", "v"; backend = OnDisk())
+T = FieldTimeSeries("global_surface_fields.jld2", "T"; backend = OnDisk())
+e = FieldTimeSeries("global_surface_fields.jld2", "e"; backend = OnDisk())
+
+times = u.times
+Nt = length(times)
+
+n = Observable(Nt)
+
+# We create a land mask and use it to fill land points with `NaN`s.
+land = interior(T.grid.immersed_boundary.bottom_height) .≥ 0
+
+Tn = @lift begin
+    Tn = interior(T[$n])
+    Tn[land] .= NaN
+    view(Tn, :, :, 1)
+end
+
+en = @lift begin
+    en = interior(e[$n])
+    en[land] .= NaN
+    view(en, :, :, 1)
+end
+
+# We compute the surface speed.
+un = Field{Face, Center, Nothing}(u.grid)
+vn = Field{Center, Face, Nothing}(v.grid)
+s = Field(sqrt(un^2 + vn^2))
+
+sn = @lift begin
+    parent(un) .= parent(u[$n])
+    parent(vn) .= parent(v[$n])
+    compute!(s)
+    sn = interior(s)
+    sn[land] .= NaN
+    view(sn, :, :, 1)
+end
+
+# Finally, we plot a snapshot of the surface speed, temperature, and the turbulent
+# eddy kinetic energy from the CATKE vertical mixing parameterization.
+fig = Figure(size = (800, 1200))
+
+axs = Axis(fig[1, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)")
+axT = Axis(fig[2, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)")
+axe = Axis(fig[3, 1], xlabel="Longitude (deg)", ylabel="Latitude (deg)")
+
+hm = heatmap!(axs, sn, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray)
+Colorbar(fig[1, 2], hm, label = "Surface speed (m s⁻¹)")
+
+hm = heatmap!(axT, Tn, colorrange = (-1, 30), colormap = :magma, nan_color=:lightgray)
+Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)")
+
+hm = heatmap!(axe, en, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray)
+Colorbar(fig[3, 2], hm, label = "Turbulent Kinetic Energy (m² s⁻²)")
+
+save("global_snapshot.png", fig)
+nothing #hide
+
+# ![](global_snapshot.png)
+
+# And now a movie:
+
+record(fig, "one_degree_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn
+    n[] = nn
+end
+nothing #hide
+
+# ![](one_degree_global_ocean_surface.mp4)

examples/single_column_os_papa_simulation.jl

@@ -122,9 +122,9 @@ function progress(sim)
     S = sim.model.ocean.model.tracers.S
     e = sim.model.ocean.model.tracers.e
 
-    τx = first(sim.model.fluxes.total.ocean.momentum.u)
-    τy = first(sim.model.fluxes.total.ocean.momentum.v)
-    Q = first(sim.model.fluxes.total.ocean.heat)
+    τx = first(sim.model.interfaces.net_fluxes.ocean_surface.u)
+    τy = first(sim.model.interfaces.net_fluxes.ocean_surface.v)
+    Q  = first(sim.model.interfaces.net_fluxes.ocean_surface.Q)
 
     u★ = sqrt(sqrt(τx^2 + τy^2))
 
@@ -142,15 +142,15 @@ end
 simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))
 
 # Build flux outputs
-τx = coupled_model.fluxes.total.ocean.momentum.u
-τy = coupled_model.fluxes.total.ocean.momentum.v
-JT = coupled_model.fluxes.total.ocean.tracers.T
-Js = coupled_model.fluxes.total.ocean.tracers.S
-E  = coupled_model.fluxes.turbulent.fields.water_vapor
-Qc = coupled_model.fluxes.turbulent.fields.sensible_heat
-Qv = coupled_model.fluxes.turbulent.fields.latent_heat
-ρₒ = coupled_model.fluxes.ocean_reference_density
-cₚ = coupled_model.fluxes.ocean_heat_capacity
+τx = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.x_momentum
+τy = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.y_momentum
+JT = coupled_model.interfaces.net_fluxes.ocean_surface.T
+Js = coupled_model.interfaces.net_fluxes.ocean_surface.S
+E  = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.water_vapor
+Qc = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.sensible_heat
+Qv = coupled_model.interfaces.atmosphere_ocean_interface.fluxes.latent_heat
+ρₒ = coupled_model.interfaces.fluxes.ocean_reference_density
+cₚ = coupled_model.interfaces.fluxes.ocean_heat_capacity
 
 Q = ρₒ * cₚ * JT
 ρτx = ρₒ * τx