Merge remote-tracking branch 'origin/main' into ss/integrate-sea-ice

docs/make.jl

@@ -13,7 +13,6 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
 const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")
 
 to_be_literated = [
-    "ecco_inspect_temperature_salinity.jl",
     "generate_bathymetry.jl",
     "generate_surface_fluxes.jl",
     "single_column_os_papa_simulation.jl",
@@ -41,7 +40,6 @@ pages = [
     "Home" => "index.md",
 
     "Examples" => [
-        "Inspect ECCO2 data" => "literated/ecco_inspect_temperature_salinity.md",
         "Generate bathymetry" => "literated/generate_bathymetry.md",
         "Surface fluxes" => "literated/generate_surface_fluxes.md",
         "Single-column simulation" => "literated/single_column_os_papa_simulation.md",
@@ -59,7 +57,7 @@ pages = [
 
 makedocs(sitename = "ClimaOcean.jl";
          format,
-         pages, 
+         pages,
          modules = [ClimaOcean],
          doctest = true,
          clean = true,

docs/src/library/internals.md

@@ -38,6 +38,13 @@ Modules = [ClimaOcean.ECCO]
 Public = false
 ```
 
+## JRA55
+
+```@autodocs
+Modules = [ClimaOcean.JRA55]
+Public = false
+```
+
 ## Bathymetry
 
 ```@autodocs

docs/src/library/public.md

@@ -39,6 +39,13 @@ Modules = [ClimaOcean.ECCO]
 Private = false
 ```
 
+## JRA55
+
+```@autodocs
+Modules = [ClimaOcean.JRA55]
+Private = false
+```
+
 ## Bathymetry
 
 ```@autodocs

examples/near_global_ocean_simulation.jl

@@ -25,11 +25,11 @@ using Printf
 #
 # We define a global grid with a horizontal resolution of 1/4 degree and 40 vertical levels.
 # The grid is a `LatitudeLongitudeGrid` spanning latitudes from 75°S to 75°N.
-# We use an exponential vertical spacing to better resolve the upper ocean layers.
+# We use an exponential vertical spacing to better resolve the upper-ocean layers.
 # The total depth of the domain is set to 6000 meters.
 # Finally, we specify the architecture for the simulation, which in this case is a GPU.
 
-arch = GPU() 
+arch = GPU()
 
 Nx = 1440
 Ny = 600
@@ -49,10 +49,10 @@ grid = LatitudeLongitudeGrid(arch;
 #
 # We use `regrid_bathymetry` to derive the bottom height from ETOPO1 data.
 # To smooth the interpolated data we use 5 interpolation passes. We also fill in
-# all the minor enclosed basins except the 3 largest `major_basins`, as well as regions
-# that are shallower than `minimum_depth`.
+# (i) all the minor enclosed basins except the 3 largest `major_basins`, as well as
+# (ii) regions that are shallower than `minimum_depth`.
 
-bottom_height = regrid_bathymetry(grid; 
+bottom_height = regrid_bathymetry(grid;
                                   minimum_depth = 10meters,
                                   interpolation_passes = 5,
                                   major_basins = 3)
@@ -64,8 +64,7 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_
 h = grid.immersed_boundary.bottom_height
 
 fig, ax, hm = heatmap(h, colormap=:deep, colorrange=(-depth, 0))
-cb = Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false)
-hidedecorations!(ax)
+Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false)
 save("bathymetry.png", fig)
 nothing #hide
 
@@ -81,7 +80,7 @@ ocean = ocean_simulation(grid)
 
 ocean.model
 
-# We initialize the ocean model to ECCO2 temperature and salinity for January 1, 1993.
+# We initialize the ocean model with ECCO2 temperature and salinity for January 1, 1993.
 
 date = DateTimeProlepticGregorian(1993, 1, 1)
 set!(ocean.model, T=ECCOMetadata(:temperature; dates=date),
@@ -90,8 +89,7 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=date),
 # ### Prescribed atmosphere and radiation
 #
 # Next we build a prescribed atmosphere state and radiation model,
-# which will drive the development of the ocean simulation.
-# We use the default `Radiation` model,
+# which will drive the ocean simulation. We use the default `Radiation` model,
 
 ## The radiation model specifies an ocean albedo emissivity to compute the net radiative
 ## fluxes. The default ocean albedo is based on Payne (1982) and depends on cloud cover
@@ -101,9 +99,8 @@ set!(ocean.model, T=ECCOMetadata(:temperature; dates=date),
 radiation = Radiation(arch)
 
 # The atmospheric data is prescribed using the JRA55 dataset.
-# The JRA55 dataset provides atmospheric
-# data such as temperature, humidity, and wind fields to calculate turbulent fluxes
-# using bulk formulae, see [`CrossRealmFluxes`](@ref).
+# The JRA55 dataset provides atmospheric data such as temperature, humidity, and winds
+# to calculate turbulent fluxes using bulk formulae, see [`CrossRealmFluxes`](@ref).
 # The number of snapshots that are loaded into memory is determined by
 # the `backend`. Here, we load 41 snapshots at a time into memory.
 
@@ -116,16 +113,16 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
 
 coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
 
-# We then create a coupled simulation, starting with a time step of 90 seconds
-# and running the simulation for 10 days.
+# We then create a coupled simulation. We start with a small-ish time step of 90 seconds.
+# We run the simulation for 10 days with this small-ish time step.
 
 simulation = Simulation(coupled_model; Δt=90, stop_time=10days)
 
 # We define a callback function to monitor the simulation's progress,
 
 wall_time = Ref(time_ns())
 
-function progress(sim) 
+function progress(sim)
     ocean = sim.model.ocean
     u, v, w = ocean.model.velocities
     T = ocean.model.tracers.T
@@ -154,7 +151,8 @@ simulation.callbacks[:progress] = Callback(progress, TimeInterval(5days))
 #
 # We define output writers to save the simulation data at regular intervals.
 # In this case, we save the surface fluxes and surface fields at a relatively high frequency (every day).
-# The `indices` keyword argument allows us to save down a slice at the surface, which is located at `k = grid.Nz`
+# The `indices` keyword argument allows us to save only a slice of the three dimensional variable.
+# Below, we use `indices` to save only the values of the variables at the surface, which corresponds to `k = grid.Nz`
 
 outputs = merge(ocean.model.tracers, ocean.model.velocities)
 ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs;
@@ -167,21 +165,24 @@ ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs;
 
 # ### Spinning up the simulation
 #
-# We spin up the simulation with a very short time-step to resolve the "initialization shock"
-# associated with starting from ECCO initial conditions that are both interpolated and also
+# We spin up the simulation with a small-ish time-step to resolve the "initialization shock"
+# associated with starting from ECCO2 initial conditions that are both interpolated and also
 # satisfy a different dynamical balance than our simulation.
 
 run!(simulation)
 
 # ### Running the simulation for real
 
+# After the initial spin up of 10 days, we can increase the time-step and run for longer.
+
 simulation.stop_time = 60days
 simulation.Δt = 10minutes
 run!(simulation)
 
 # ## A pretty movie
 #
-# It's time to make a pretty movie of the simulation. First we plot a snapshot:
+# It's time to make a pretty movie of the simulation. First we load the output we've been saving on
+# disk and plot the final snapshot:
 
 u = FieldTimeSeries("near_global_surface_fields.jld2", "u"; backend = OnDisk())
 v = FieldTimeSeries("near_global_surface_fields.jld2", "v"; backend = OnDisk())
@@ -209,7 +210,9 @@ end
 
 un = Field{Face, Center, Nothing}(u.grid)
 vn = Field{Center, Face, Nothing}(v.grid)
-s = Field(sqrt(un^2 + vn^2))
+
+s = @at (Center, Center, Nothing) sqrt(un^2 + vn^2) # compute √(u²+v²) and interpolate back to Center, Center
+s = Field(s)
 
 sn = @lift begin
     parent(un) .= parent(u[$n])
@@ -220,26 +223,34 @@ sn = @lift begin
     view(sn, :, :, 1)
 end
 
-fig = Figure(size = (800, 1200))
+title = @lift string("Near-global 1/4 degree ocean simulation after ",
+                     prettytime(times[$n] - times[1]))
+
+λ, φ, _ = nodes(T) # T, e, and s all live on the same grid locations
+
+fig = Figure(size = (1000, 1500))
 
 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!(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)
+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)
+hm = heatmap!(axe, λ, φ, en, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray)
 Colorbar(fig[3, 2], hm, label = "Turbulent Kinetic Energy (m² s⁻²)")
+
+Label(fig[0, :], title)
+
 save("snapshot.png", fig)
 nothing #hide
 
 # ![](snapshot.png)
 
-# And now a movie:
+# And now we make a movie:
 
 record(fig, "near_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn
     n[] = nn