diff --git a/docs/make.jl b/docs/make.jl index 317e81e6..0964c0f6 100644 --- 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 index 6db964bc..ae7ce604 100644 --- 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,8 +116,8 @@ 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) @@ -125,7 +125,7 @@ simulation = Simulation(coupled_model; Δt=90, stop_time=10days) 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 index bf9ee264..b9148d06 100644 --- 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 = , + 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 diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl index 3da5b500..46c6e5b7 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/surface_temperature.jl @@ -12,11 +12,12 @@ import Thermodynamics as AtmosphericThermodynamics #### Bulk surface temperature (the easiest case) #### -""" - struct BulkTemperature end +""" + struct BulkTemperature -Represents the surface temperature used in fixed-point iteration for surface fluxes following similarity theory. -The surface temperature is not calculated but provided by either the ocean or the sea ice model. +A type to represent the surface temperature used in fixed-point iteration for surface +fluxes following similarity theory. The surface temperature is not calculated but instead +provided by either the ocean or the sea ice model. """ struct BulkTemperature end @@ -29,17 +30,19 @@ struct BulkTemperature end """ struct SkinTemperature - internal_flux :: I - end A type to represent the surface temperature used in the flux calculation. The surface temperature is calculated from the flux balance at the surface. In particular, the surface temperature ``θₛ`` is the root of: -F(θₛ) - Jᵀ = 0 (all fluxes positive upwards) +```math +F(θₛ) - Jᵀ = 0 +``` + +where ``Jᵀ`` are the fluxes at the top of the surface (turbulent + radiative), and +``F`` is the internal diffusive flux dependent on the surface temperature itself. -where Jᵀ are the fluxes at the top of the surface (turbulent + radiative), and F is the internal diffusive flux -dependent on the surface temperature itself. +Note that all fluxes positive upwards. """ struct SkinTemperature{I} internal_flux :: I @@ -77,19 +80,19 @@ DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert( # that can be explored in the future. @inline flux_balance_temperature(F::DiffusiveFlux, θₒ, Jᵀ) = θₒ - Jᵀ / F.κ * F.δ -# he flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the +# the flaw here is that the ocean emissivity and albedo are fixed, but they might be a function of the # surface temperature, so we might need to pass the radiation and the albedo and emissivity as arguments. -@inline function compute_surface_temperature(st::SkinTemperature, θₛ, ℂ, 𝒬₀, - ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, - prescribed_heat_fluxes, +@inline function compute_surface_temperature(st::SkinTemperature, θₛ, ℂ, 𝒬₀, + ρₐ, cₚ, ℰv, Σ★, ρₒ, cpₒ, g, + prescribed_heat_fluxes, radiation_properties) Rd = prescribed_heat_fluxes # net downwelling radiation (positive out of the ocean) - + # upwelling radiation is calculated explicitly Ru = upwelling_radiation(θₛ, radiation_properties) Rn = Rd + Ru # Net radiation (positive out of the ocean) - + u★ = Σ★.momentum θ★ = Σ★.temperature q★ = Σ★.water_vapor @@ -104,4 +107,4 @@ DiffusiveFlux(FT; κ = 1e-2, δ = 1.0) = DiffusiveFlux(convert(FT, δ), convert( θₛ = flux_balance_temperature(st.internal_flux, θₒ, Jᵀ) return θₛ -end \ No newline at end of file +end