Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix axis labels in near-global example + Docstring rendering + Few typos #315

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
4 changes: 1 addition & 3 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -59,7 +57,7 @@ pages = [

makedocs(sitename = "ClimaOcean.jl";
format,
pages,
pages,
modules = [ClimaOcean],
doctest = true,
clean = true,
Expand Down
60 changes: 36 additions & 24 deletions examples/near_global_ocean_simulation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ 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.

Expand All @@ -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)
Expand Down Expand Up @@ -81,7 +81,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),
Expand All @@ -90,8 +90,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
Expand All @@ -101,9 +100,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.

Expand All @@ -116,16 +114,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
Expand Down Expand Up @@ -154,7 +152,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;
Expand All @@ -167,21 +166,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())
Expand Down Expand Up @@ -209,7 +211,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])
Expand All @@ -220,26 +224,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
Expand Down
83 changes: 44 additions & 39 deletions src/Bathymetry.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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``).
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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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))

Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down
Loading
Loading