Fix axis labels in near-global example + Docstring rendering + Few ty…

…pos (#315) * fix docstring rendering + few typos * add coords in plots + phrasing tweaks * drop inspect ecco * interpolate s -> (Center, Center) * markdown formatting requires 2 spaces for list within a list * there is no default in remove_minor_basins! * fix math rendering * Update near_global_ocean_simulation.jl * fix bug + add title with time * Update make.jl * Update make.jl * Update src/Bathymetry.jl Co-authored-by: Simone Silvestri <[email protected]> * homogenize font label size * homogenize font label size --------- Co-authored-by: Simone Silvestri <[email protected]>
CliMA · Dec 25, 2024 · 0baa94f · 0baa94f
1 parent 3d76db9
commit 0baa94f
Show file tree

Hide file tree

Showing 4 changed files with 104 additions and 84 deletions.
diff --git a/docs/make.jl b/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,

diff --git a/examples/near_global_ocean_simulation.jl b/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)
@@ -63,7 +63,9 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_
 
 h = grid.immersed_boundary.bottom_height
 
-fig, ax, hm = heatmap(h, colormap=:deep, colorrange=(-depth, 0))
+fig = Figure(size = (800, 450))
+ax = Axis(fig[1, 1])
+hm = heatmap!(ax, h, colormap=:deep, colorrange=(-depth, 0))
 cb = Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false)
 hidedecorations!(ax)
 save("bathymetry.png", fig)
@@ -81,7 +83,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 +92,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 +102,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 +116,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 +154,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 +168,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 +213,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 +226,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

diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl
@@ -32,65 +32,70 @@ end
 
 """
     regrid_bathymetry(target_grid;
-                      url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", 
-                      height_above_water = <none>,
+                      height_above_water = nothing,
                       minimum_depth = 0,
-                      dir = download_cache,
-                      filename = "ETOPO_2022_v1_60s_N90W180_surface.nc")
+                      dir = download_bathymetry_cache,
+                      url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf",
+                      filename = "ETOPO_2022_v1_60s_N90W180_surface.nc",
+                      interpolation_passes = 1,
+                      major_basins = 1)
 
-Regrid bathymetry associated with the NetCDF file at `path = joinpath(dir, filename)` to `target_grid`.
+Return bathymetry associated with the NetCDF file at `path = joinpath(dir, filename)` regridded onto `target_grid`.
 If `path` does not exist, then a download is attempted from `joinpath(url, filename)`.
 
 Arguments
 =========
 
-- `target_grid`: grid to interpolate onto
+- `target_grid`: grid to interpolate the bathymetry onto.
 
 Keyword Arguments
 =================
 
-- `height_above_water`: limits the maximum height of above-water topography (where h > 0). If
-                        `nothing` the original topography is retained
+- `height_above_water`: limits the maximum height of above-water topography (where ``h > 0``) before inetrpolating.
+                        Default: `nothing`, which implies that the original topography is retained.
 
 - `minimum_depth`: minimum depth for the shallow regions, defined as a positive value. 
-                   `h > - minimum_depth` will be considered land
+                   `h > - minimum_depth` is considered land. Default: 0.
 
-- `dir`: directory of the bathymetry-containing file
+- `dir`: directory of the bathymetry-containing file. Default: `download_bathymetry_cache`.
 
-- `filename`: file containing bathymetric data. Must be netcdf with fields:
-              (1) `lat` vector of latitude nodes
-              (2) `lon` vector of longitude nodes
-              (3) `z` matrix of depth values
+- `filename`: file containing bathymetric data. Must be netCDF with fields:
+  1. `lat` vector of latitude nodes
+  2. `lon` vector of longitude nodes
+  3. `z` matrix of depth values
 
 - `interpolation_passes`: regridding/interpolation passes. The bathymetry is interpolated in
-                          `interpolation_passes - 1` intermediate steps. With more steps the 
-                          final bathymetry will be smoother.
-                          
-  Example: interpolating from a 400x200 grid to a 100x100 grid in 4 passes involves:
+                          `interpolation_passes - 1` intermediate steps. The more the interpolation
+                          steps the smoother the final bathymetry becomes.
 
-  * 400x200 → 325x175
-  * 325x175 → 250x150
-  * 250x150 → 175x125
-  * 175x125 → 100x100
+  Example
+  =======
+
+  Interpolating from a 400 x 200 grid to a 100 x 100 grid in 4 passes involves:
+
+  * 400 x 200 → 325 x 175
+  * 325 x 175 → 250 x 150
+  * 250 x 150 → 175 x 125
+  * 175 x 125 → 100 x 100
 
   If _coarsening_ the original grid, linear interpolation in passes is equivalent to
   applying a smoothing filter, with more passes increasing the strength of the filter.
-  If _refining_ the original grid, additional passes will not help and no intermediate
-  steps will be performed.
+  If _refining_ the original grid, additional passes do not help and no intermediate
+  steps are performed.
 
 - `major_basins`: Number of "independent major basins", or fluid regions fully encompassed by land,
                   that are retained by [`remove_minor_basins!`](@ref). Basins are removed by order of size:
-                  the smallest basins are removed first. `major_basins=1` will retain only the largest basin.
-                  Default: `Inf`, which does not remove any basins.
+                  the smallest basins are removed first. `major_basins = 1` retains only the largest basin.
+                  If `Inf` then no basins are removed. Default: 1.
 """
 function regrid_bathymetry(target_grid;
                            height_above_water = nothing,
                            minimum_depth = 0,
                            dir = download_bathymetry_cache,
-                           url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf", 
+                           url = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/60s/60s_surface_elev_netcdf",
                            filename = "ETOPO_2022_v1_60s_N90W180_surface.nc",
                            interpolation_passes = 1,
-                           major_basins = 1) # Allow an `Inf` number of ``lakes''
+                           major_basins = 1) # Allow an `Inf` number of "lakes"
 
     filepath = joinpath(dir, filename)
     fileurl  = url * "/" * filename # joinpath on windows creates the wrong url
@@ -191,7 +196,7 @@ function regrid_bathymetry(target_grid;
 end
 
 # Here we can either use `regrid!` (three dimensional version) or `interpolate`
-function interpolate_bathymetry_in_passes(native_z, target_grid; 
+function interpolate_bathymetry_in_passes(native_z, target_grid;
                                           passes = 10)
     Nλt, Nφt = Nt = size(target_grid)
     Nλn, Nφn = Nn = size(native_z)
@@ -201,9 +206,9 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
     if Nλt > Nλn || Nφt > Nφn
         target_z = Field{Center, Center, Nothing}(target_grid)
         interpolate!(target_z, native_z)
-        @info string("Skipping passes for interplating bathymetry of size $Nn ", '\n',
+        @info string("Skipping passes for interpolating bathymetry of size $Nn ", '\n',
                      "to target grid of size $Nt. Interpolation passes may only ", '\n',
-                     "be used to refine bathymetryand requires that the bathymetry ", '\n',
+                     "be used to refine bathymetry and require that the bathymetry ", '\n',
                      "is larger than the target grid in both horizontal directions.")
         return target_z
     end
@@ -230,9 +235,9 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
         @debug "Bathymetry interpolation pass $pass with size $new_size"
 
         new_grid = LatitudeLongitudeGrid(architecture(target_grid),
-                                         size = new_size, 
-                                         latitude = (latitude[1],  latitude[2]), 
-                                         longitude = (longitude[1], longitude[2]), 
+                                         size = new_size,
+                                         latitude = (latitude[1],  latitude[2]),
+                                         longitude = (longitude[1], longitude[2]),
                                          z = (0, 1),
                                          topology = (TX, TY, Bounded))
 
@@ -276,7 +281,7 @@ function regrid_bathymetry(target_grid::DistributedGrid; kw...)
     local_bottom_height = Field{Center, Center, Nothing}(target_grid)
     set!(local_bottom_height, bottom_height)
     fill_halo_regions!(local_bottom_height)
-    
+
     return local_bottom_height
 end
 
@@ -291,14 +296,14 @@ Arguments
 =========
 
 - `z_data`: A 2D array representing the bathymetry data.
-- `keep_major_basins`: The maximum number of connected regions to keep. 
-                       Default is `Inf`, which means all connected regions are kept.
+- `keep_major_basins`: The maximum number of connected regions to keep.
+                       If `Inf` is provided then all connected regions are kept.
 
 """
 function remove_minor_basins!(zb::Field, keep_major_basins)
     zb_cpu = on_architecture(CPU(), zb)
     TX     = topology(zb_cpu.grid, 1)
-    
+
     Nx, Ny, _ = size(zb_cpu.grid)
     zb_data   = maybe_extend_longitude(zb_cpu, TX()) # Outputs a 2D AbstractArray
 
@@ -326,7 +331,7 @@ function maybe_extend_longitude(zb_cpu, ::Periodic)
     # Update offsets
     yoffsets = zb_cpu.data.offsets[2]
     xoffsets = - nx
-    
+
     return OffsetArray(zb_parent, xoffsets, yoffsets)
 end