diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index e3fb4e2a9..93e55e34d 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -151,4 +151,3 @@ steps: - wait: ~ continue_on_failure: true - diff --git a/.github/workflows/DocPreviewCleanup.yml b/.github/workflows/DocPreviewCleanup.yml deleted file mode 100644 index 5d38c6eb1..000000000 --- a/.github/workflows/DocPreviewCleanup.yml +++ /dev/null @@ -1,29 +0,0 @@ -name: Doc Preview Cleanup - -on: - pull_request: - types: [closed] - -jobs: - doc-preview-cleanup: - runs-on: ubuntu-latest - steps: - - name: Checkout gh-pages branch - uses: actions/checkout@v4 - with: - repository: CliMA/ClimaOceanDocumentation - ref: gh-pages - - name: Delete preview and history + push changes - run: | - if [ -d "previews/PR$PRNUM" ]; then - git config user.name "Documenter.jl" - git config user.email "documenter@juliadocs.github.io" - git rm -rf "previews/PR$PRNUM" - git commit -m "delete preview" - git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree}) - git push --force origin gh-pages-new:gh-pages - fi - env: - PRNUM: ${{ github.event.number }} - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - SSH_DEPLOY_KEY: ${{ secrets.SSH_DEPLOY_KEY }} # For authentication with SSH deploy key diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml deleted file mode 100644 index 95ea069f7..000000000 --- a/.github/workflows/Documenter.yml +++ /dev/null @@ -1,36 +0,0 @@ -name: Documentation - -on: - push: - branches: - - main - tags: '*' - pull_request: - -jobs: - build: - runs-on: ubuntu-latest - steps: - - name: Cancel Previous Runs - uses: styfle/cancel-workflow-action@0.8.0 - with: - access_token: ${{ github.token }} - - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@latest - with: - version: '1.10' - show-versioninfo: true - - name: Install dependencies - run: | - julia --color=yes --project -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project -e 'using Pkg; Pkg.precompile()' - julia --color=yes --project=docs/ -e 'using Pkg; Pkg.instantiate()' - julia --color=yes --project=docs/ -e 'using Pkg; Pkg.precompile()' - - name: Build and deploy - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key - JULIA_DEBUG: Documenter - ECCO_USERNAME: ${{ secrets.ECCO_USERNAME }} # To download ECCO data from the podaac website - ECCO_PASSWORD: ${{ secrets.ECCO_PASSWORD }} # To download ECCO data from the podaac website - run: julia --color=yes --project=docs/ docs/make_without_examples.jl diff --git a/docs/Project.toml b/docs/Project.toml index c5baf5a70..ec9e24ea8 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,11 +1,14 @@ [deps] +CFTime = "179af706-886a-5703-950a-314cd64e0468" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" OrthogonalSphericalShellGrids = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0" [compat] -CairoMakie = "0.10.12" +CairoMakie = "0.10.12, 0.11, 0.12" DataDeps = "0.7" Documenter = "1" +Oceananigans = "0.93.1" diff --git a/docs/make.jl b/docs/make.jl index b2bee462b..4e1ada332 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,9 +1,7 @@ -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..")) # add ClimaOcean to environment stack - using + ClimaOcean, Documenter, - Literate, - ClimaOcean + Literate ENV["DATADEPS_ALWAYS_ACCEPT"] = "true" @@ -15,11 +13,11 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/literated") to_be_literated = [ - # "inspect_ECCO_data.jl", + "inspect_ecco_data.jl", "generate_bathymetry.jl", - "generate_surface_fluxes.jl", - # "single_column_simulation.jl", - # "mediterranean_simulation_with_ECCO_restoring.jl", + # "generate_surface_fluxes.jl", + "single_column_os_papa_simulation.jl", + # "mediterranean_simulation_with_ecco_restoring.jl", "near_global_ocean_simulation.jl" ] @@ -34,12 +32,10 @@ end ##### Build and deploy docs ##### -format = Documenter.HTML( - collapselevel = 2, - size_threshold = nothing, - prettyurls = get(ENV, "CI", nothing) == "true", - canonical = "https://clima.github.io/ClimaOceanDocumentation/dev/", -) +format = Documenter.HTML(collapselevel = 2, + size_threshold = nothing, + prettyurls = get(ENV, "CI", nothing) == "true", + canonical = "https://clima.github.io/ClimaOceanDocumentation/dev/") pages = [ "Home" => "index.md", @@ -52,25 +48,21 @@ pages = [ ], "Examples" => [ - # "Inspect ECCO2 data" => "literated/inspect_ECCO_data.md", + "Inspect ECCO2 data" => "literated/inspect_ecco_data.md", "Generate bathymetry" => "literated/generate_bathymetry.md", - "Surface fluxes" => "literated/generate_surface_fluxes.md", - # "Single column simulation" => "literated/single_column_simulation.md", - # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ECCO_restoring.md", + # "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", "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md", ] ] -makedocs( - sitename = "ClimaOcean.jl", - modules = [ClimaOcean], - format = format, - pages = pages, - doctest = true, - clean = true, - warnonly = [:cross_references, :missing_docs], - checkdocs = :exports -) +makedocs(sitename = "ClimaOcean.jl"; format, pages, + modules = [ClimaOcean], + doctest = true, + clean = true, + warnonly = [:cross_references, :missing_docs], + checkdocs = :exports) @info "Clean up temporary .jld2 and .nc output created by doctests or literated examples..." @@ -95,8 +87,8 @@ end withenv("GITHUB_REPOSITORY" => "github.com/CliMA/ClimaOceanDocumentation.git") do deploydocs(repo = "github.com/CliMA/ClimaOceanDocumentation.git", - versions = ["stable" => "v^", "v#.#.#", "dev" => "dev"], - forcepush = true, - devbranch = "main", - push_preview = true) + versions = ["stable" => "v^", "dev" => "dev", "v#.#.#"], + forcepush = true, + devbranch = "main", + push_preview = true) end diff --git a/docs/make_without_examples.jl b/docs/make_without_examples.jl deleted file mode 100644 index 8baec2245..000000000 --- a/docs/make_without_examples.jl +++ /dev/null @@ -1,69 +0,0 @@ -pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..")) # add ClimaOcean to environment stack - -using - Documenter, - Literate, - ClimaOcean - -ENV["DATADEPS_ALWAYS_ACCEPT"] = "true" - -##### -##### Build and deploy docs -##### - -format = Documenter.HTML( - collapselevel = 2, - prettyurls = get(ENV, "CI", nothing) == "true", - canonical = "https://clima.github.io/ClimaOceanDocumentation/dev/", -) - -pages = [ - "Home" => "index.md", - - "Library" => [ - "Contents" => "library/outline.md", - "Public" => "library/public.md", - "Private" => "library/internals.md", - "Function index" => "library/function_index.md", - ], -] - -makedocs( - sitename = "ClimaOcean.jl", - modules = [ClimaOcean], - format = format, - pages = pages, - doctest = true, - clean = true, - warnonly = [:cross_references, :missing_docs], - checkdocs = :exports -) - -@info "Clean up temporary .jld2 and .nc output created by doctests or literated examples..." - -""" - recursive_find(directory, pattern) - -Return list of filepaths within `directory` that contains the `pattern::Regex`. -""" -recursive_find(directory, pattern) = - mapreduce(vcat, walkdir(directory)) do (root, dirs, files) - joinpath.(root, filter(contains(pattern), files)) - end - -files = [] -for pattern in [r"\.jld2", r"\.nc"] - global files = vcat(files, recursive_find(@__DIR__, pattern)) -end - -for file in files - rm(file) -end - -withenv("GITHUB_REPOSITORY" => "github.com/CliMA/ClimaOceanDocumentation.git") do - deploydocs(repo = "github.com/CliMA/ClimaOceanDocumentation.git", - versions = ["stable" => "v^", "dev" => "dev", "v#.#.#"], - forcepush = true, - devbranch = "main", - push_preview = true) -end diff --git a/examples/generate_atmos_dataset.jl b/examples/generate_atmos_dataset.jl index 24ef1cc0d..8ba1d3233 100644 --- a/examples/generate_atmos_dataset.jl +++ b/examples/generate_atmos_dataset.jl @@ -15,4 +15,3 @@ T = Array(interior(Tt[1], 1:4:Nx, 1:4:Ny, 1)) p = Array(interior(pt[1], 1:4:Nx, 1:4:Ny, 1)) @save "atmospheric_state.jld2" q T p - diff --git a/examples/generate_bathymetry.jl b/examples/generate_bathymetry.jl index 42713af48..ae063eeeb 100644 --- a/examples/generate_bathymetry.jl +++ b/examples/generate_bathymetry.jl @@ -1,7 +1,7 @@ # # Generate bathymetry data for the Mediterranean Sea # -# This script shows how to configure an Immersed boundary grid with realistic bathymetry using ClimaOcean.jl -# by generating the bathymetry data for the Mediterranean Sea. +# This example shows how to configure an Immersed boundary grid with realistic bathymetry +# using ClimaOcean.jl by generating the bathymetry data for the Mediterranean Sea. # # For this example, we need Oceananigans for the LatitudeLongitudeGrid and Field utilities, # ClimaOcean to donwload and regrid the bathymetry, and CairoMakie to visualize the grid. @@ -25,7 +25,8 @@ Nλ = 25 * (longitude_range[2] - longitude_range[1]) grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1), latitude = latitude_range, longitude = longitude_range, - z = (0, 1)) + z = (0, 1), + halo = (7, 7, 1)) # Next, we generate the bathymetry data for the Mediterranean Sea using the # `regrid_bathymetry` function from ClimaOcean. The function downloads the bathymetry @@ -33,7 +34,7 @@ grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1), # bathymetry field. The three different regidding procedures below demonstrate the effect # of different parameters on the generated bathymetry: # -# - `h_rough` shows the output of the function with default parameters, which means only +# - `h_rough` shows the output of the function with default parameters, which means only # one interpolation passes and no restrictions on connected regions. # - `h_smooth` shows the output of the function with 40 interpolation passes, which results # in a smoother bathymetry. @@ -44,7 +45,7 @@ grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1), h_rough = regrid_bathymetry(grid) h_smooth = regrid_bathymetry(grid; interpolation_passes = 40) h_one_basin = regrid_bathymetry(grid; major_basins = 1) -nothing # hide +nothing #hide # Finally, we visualize the generated bathymetry data for the Mediterranean Sea using CairoMakie. @@ -59,7 +60,7 @@ interior(h_rough)[land_rough] .= NaN land_one_basin = interior(h_one_basin) .>= 0 interior(h_one_basin)[land_one_basin] .= NaN -fig = Figure(resolution=(850, 1150)) +fig = Figure(size=(850, 1150)) ax = Axis(fig[1, 1], title = "Rough bathymetry", xlabel = "Longitude", ylabel = "Latitude") hm = heatmap!(ax, λ, φ, - interior(h_rough, :, :, 1), nan_color=:white, colormap = Reverse(:deep)) diff --git a/examples/generate_surface_fluxes.jl b/examples/generate_surface_fluxes.jl index 2b65d7ac1..f916cbda2 100644 --- a/examples/generate_surface_fluxes.jl +++ b/examples/generate_surface_fluxes.jl @@ -3,7 +3,7 @@ # ClimaOcean uses bulk formulae to estimate the surface exchange of momentum, # heat, and water vapor between the atmosphere and the ocean. # -# This script demonstrates an example of the turbulent surface flux calculations performed in ClimaOcean +# This example demonstrates an example of the turbulent surface flux calculations performed in ClimaOcean # using ECCO2 data for the ocean and JRA55 data for the atmosphere. # # For this example, we need ClimaOcean with its DataWrangling modules: ECCO2 and JRA55. @@ -28,7 +28,7 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(mask)) fig = Figure() ax = Axis(fig[1, 1]) heatmap!(ax, interior(grid.immersed_boundary.mask, :, :, grid.Nz)) -save("ECCO_continents.png", fig) # hide +save("ECCO_continents.png", fig) #hide # ![](ECCO_continents.png) @@ -47,7 +47,7 @@ save("ECCO_continents.png", fig) # hide # We invoke the constructor with only the first two time indices, corresponding to # January 1st (at 00:00 AM and 03:00 AM). -atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory()) +atmosphere = JRA55_prescribed_atmosphere(1:2; backend = InMemory()) ocean = ocean_simulation(grid) # Now that we have an atmosphere and ocean, we `set!` the ocean temperature and salinity @@ -92,4 +92,3 @@ heatmap!(ax, Mv; colormap = :bwr) save("fluxes.png", fig) # ![](fluxes.png) - diff --git a/examples/inspect_ecco_data.jl b/examples/inspect_ecco_data.jl index 2100a73f1..59c53da0e 100644 --- a/examples/inspect_ecco_data.jl +++ b/examples/inspect_ecco_data.jl @@ -3,10 +3,10 @@ # ClimaOcean can download and utilize data from the "ECCO" state estimate, # which stands for "Estimating the Circulation and Climate of the Ocean" --- two! # -# This script shows how to download three-dimensional temperature and salinity fields +# This example shows how to download three-dimensional temperature and salinity fields # from ECCO, and makes a short animation to showcase the fields' content. # -# For this example we need Oceananigans for Field utilities, CairoMakie for plotting +# For this example we need Oceananigans for Field utilities, CairoMakie for plotting, # Printf for nice labeling, and of course ClimaOcean to actually download and construct # the ECCO fields. @@ -16,8 +16,8 @@ using Printf using ClimaOcean.DataWrangling.ECCO: ECCO_field -# The function `ECCO_field` provided by `ClimaOcean.DataWrangling.ECCO` will automatically -# download ECCO data if it doesn't already exist at the default location. +# The function `ECCO_field` provided by `ClimaOcean.DataWrangling.ECCO` automatically +# downloads ECCO data if that doesn't already exist at the default location. T = ECCO_field(:temperature) S = ECCO_field(:salinity) @@ -40,8 +40,6 @@ fig = Figure(size=(1200, 1400)) axT = Axis(fig[1, 1]) axS = Axis(fig[2, 1]) -λ, φ, z = nodes(T) - # To make an animation that scrolls through the 3D temperature # and salinity fields, we make an Observable for the vertical index, # and then construct slices of T, S using the Observable, `k`. diff --git a/examples/mediterranean_simulation_with_ecco_restoring.jl b/examples/mediterranean_simulation_with_ecco_restoring.jl index 6a33fb3ce..9d640820f 100644 --- a/examples/mediterranean_simulation_with_ecco_restoring.jl +++ b/examples/mediterranean_simulation_with_ecco_restoring.jl @@ -1,16 +1,17 @@ # # Mediterranean simulation with restoring to ECCO # -# This Julia script is a comprehensive example of setting up and running a high-resolution ocean simulation for the Mediterranean Sea using -# the Oceananigans and ClimaOcean packages, with a focus on restoring temperature and salinity fields from the -# ECCO (Estimating the Circulation and Climate of the Ocean) dataset. -# -# The script is divided into several sections, each handling a specific part of the simulation setup and execution process. +# This example is a comprehensive example of setting up and running a high-resolution ocean +# simulation for the Mediterranean Sea using the Oceananigans and ClimaOcean packages, with +# a focus on restoring temperature and salinity fields from the ECCO (Estimating the Circulation +# and Climate of the Ocean) dataset. # +# The example is divided into several sections, each handling a specific part of the simulation +# setup and execution process. # ## Initial Setup with Package Imports # -# The script begins by importing necessary Julia packages for visualization (CairoMakie), -# ocean modeling (Oceananigans, ClimaOcean), and handling of dates and times (CFTime, Dates). +# We begin by importing necessary Julia packages for visualization (CairoMakie), ocean modeling +# (Oceananigans, ClimaOcean), and handling of dates and times (CFTime, Dates). # These packages provide the foundational tools for creating the simulation environment, # including grid setup, physical processes modeling, and data visualization. @@ -29,10 +30,10 @@ using Dates # ## Grid Configuration for the Mediterranean Sea # # The script defines a high-resolution grid to represent the Mediterranean Sea, specifying the domain in terms of longitude (λ₁, λ₂), -# latitude (φ₁, φ₂), and a stretched vertical grid to capture the depth variation (z_faces). +# latitude (φ₁, φ₂), and a stretched vertical grid to capture the depth variation (`z_faces`). # The grid resolution is set to approximately 1/15th of a degree, which translates to a spatial resolution of about 7 km. # This section demonstrates the use of the LatitudeLongitudeGrid function to create a grid that matches the -# Mediterranean's geographical and bathymetric features. +# Mediterranean's geographical and bathymetric features. λ₁, λ₂ = ( 0, 42) # domain in longitude φ₁, φ₂ = (30, 45) # domain in latitude @@ -41,8 +42,8 @@ z_faces = stretched_vertical_faces(depth = 5000, stretching = PowerLawStretching(1.070), surface_layer_height = 50) -Nx = 15 * 42 # 1 / 15th of a degree resolution -Ny = 15 * 15 # 1 / 15th of a degree resolution +Nx = 15 * Int(λ₂ - λ₁) # 1/15 th of a degree resolution +Ny = 15 * Int(φ₂ - φ₁) # 1/15 th of a degree resolution Nz = length(z_faces) - 1 grid = LatitudeLongitudeGrid(GPU(); @@ -56,7 +57,7 @@ grid = LatitudeLongitudeGrid(GPU(); # # The script interpolates bathymetric data onto the grid, ensuring that the model accurately represents # the sea floor's topography. Parameters such as `minimum_depth` and `interpolation_passes` -# are adjusted to refine the bathymetry representation. +# are adjusted to refine the bathymetry representation. bottom_height = regrid_bathymetry(grid, height_above_water = 1, @@ -89,9 +90,9 @@ ocean = ocean_simulation(grid; forcing = (T = FT, S = FS)) # Initializing the model # -# the model can be initialized with custom values or with ecco fields. +# The model can be initialized with custom values or with ecco fields. # In this case, our ECCO dataset has access to a temperature and a salinity -# field, so we initialize T and S from ECCO. +# field, so we initialize temperature T and salinity S from ECCO. set!(ocean.model, T = temperature[1], S = salinity[1]) @@ -102,7 +103,7 @@ ax = Axis(fig[1, 2]) heatmap!(ax, interior(model.tracers.S, :, :, Nz), colorrange = (35, 40), colormap = :haline) function progress(sim) - u, v, w = sim.model.velocities + u, v, w = sim.model.velocities T, S = sim.model.tracers @info @sprintf("Time: %s, Iteration %d, Δt %s, max(vel): (%.2e, %.2e, %.2e), max(T, S): %.2f, %.2f\n", @@ -118,10 +119,10 @@ ocean.callbacks[:progress] = Callback(progress, IterationInterval(10)) # ## Simulation warm up! # # We have regridded from the coarse solution of the ECCO dataset (half of a degree) to a -# fine grid (1/15th of a degree). The bathymetry might also have little mismatches -# that might crash the simulation. We warm up the simulation with a little +# fine grid (1/15th of a degree). The bathymetry might also have little mismatches +# that might crash the simulation. We warm up the simulation with a little # time step for few iterations to allow the solution to adjust to the new grid -# bathymetry +# bathymetry. ocean.Δt = 10 ocean.stop_iteration = 1000 @@ -130,7 +131,7 @@ run!(ocean) # ## Run the real simulation # # Now that the solution has adjusted to the bathymetry we can ramp up the time -# step size. We use a `TimeStepWizard` to automatically adapt to a cfl of 0.2 +# step size. We use a `TimeStepWizard` to automatically adapt to a CFL of 0.2. wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 10minutes, max_change = 1.1) diff --git a/examples/near_global_ocean_simulation.jl b/examples/near_global_ocean_simulation.jl index c53e3eb61..641d1fdea 100644 --- a/examples/near_global_ocean_simulation.jl +++ b/examples/near_global_ocean_simulation.jl @@ -1,15 +1,16 @@ # # Near-global ocean simulation # -# This Julia script sets up and runs a near-global ocean simulation using the Oceananigans.jl and ClimaOcean.jl packages. -# The simulation covers latitudes from 75°S to 75°N with a horizontal resolution of 1/4 degree and 40 vertical levels. +# This example sets up and runs a near-global ocean simulation using the Oceananigans.jl and +# ClimaOcean.jl packages. The simulation covers latitudes from 75°S to 75°N with a horizontal +# resolution of 1/4 degree and 40 vertical levels. # -# The simulation runs for one year, and the results are visualized using the CairoMakie.jl package. +# The simulation's results are visualized using the CairoMakie.jl package. # # ## Initial setup with package imports # -# The script begins by importing the necessary Julia packages for visualization (CairoMakie), -# ocean modeling (Oceananigans, ClimaOcean), and handling dates and times (CFTime, Dates). -# These packages provide the foundational tools for setting up the simulation environment, +# We begin by importing the necessary Julia packages for visualization (CairoMakie), +# ocean modeling (Oceananigans, ClimaOcean), and handling dates and times (CFTime, Dates). +# These packages provide the foundational tools for setting up the simulation environment, # including grid setup, physical processes modeling, and data visualization. using Printf @@ -23,9 +24,10 @@ using Dates # ### Grid configuration # -# We define a global grid with a horizontal resolution of 1/4 degree and 40 vertical levels. -# The grid is a `LatitudeLongitudeGrid` capped at 75°S to 75°N. -# We use an exponential vertical spacing to better resolve the upper ocean layers. The total depth of the domain is set to 6000 meters. +# 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. +# 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() @@ -39,130 +41,122 @@ Nz = length(z_faces) - 1 grid = LatitudeLongitudeGrid(arch; size = (Nx, Ny, Nz), halo = (7, 7, 7), - z = z_faces, + z = z_faces, latitude = (-75, 75), longitude = (0, 360)) # ### Bathymetry and immersed boundary # -# We retrieve the bathymetry from the ETOPO1 data, ensuring a minimum depth of 10 meters -# (depths shallower than this are considered land). The `interpolation_passes` parameter -# specifies the number of passes to interpolate the bathymetry data. A larger number -# results in a smoother bathymetry. We also remove all connected regions (such as inland -# lakes) from the bathymetry data by specifying `connected_regions_allowed = 3` (Mediterrean -# sea an North sea in addition to the ocean). +# 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 +# sminor enclosed basins but the 3 largest `major_basins` as well as reasons +# that are shallower than `minimum_depth = 10`. bottom_height = regrid_bathymetry(grid; minimum_depth = 10, interpolation_passes = 5, major_basins = 3) -# For plotting -bottom_height[bottom_height .>= 0] .= NaN +grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) -fig = Figure(size = (1200, 400)) -ax = Axis(fig[1, 1]) -hm = heatmap!(ax, bottom_height, colormap = :deep, colorrange = (-6000, 0)) -cb = Colorbar(fig[0, 1], hm, label = "Bottom height (m)", vertical = false) -hidedecorations!(ax) +# Let's see what the bathymetry looks like: + +h = grid.immersed_boundary.bottom_height +fig, ax, hm = heatmap(h, colormap=:deep, colorrange=(-6000, 0)) +cb = Colorbar(fig[0, 1], hm, label="Bottom height (m)", vertical=false) +hidedecorations!(ax) save("bathymetry.png", fig) nothing #hide # ![](bathymetry.png) -bottom_height[isnan.(bottom_height)] .= 0 -grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) - # ### Ocean model configuration # -# To configure the ocean simulation, we use the `ocean_simulation` function from ClimaOcean.jl. This function allows us to build -# an ocean simulation with default parameters and numerics. The defaults include: -# - CATKE turbulence closure for vertical mixing -# - WENO-based advection scheme for momentum in the vector invariant form -# - WENO-based advection scheme for tracers -# - `SplitExplicitFreeSurfaceSolver` with 75 substeps -# - TEOS-10 equation of state, see [`TEOS10EquationOfState`](https://clima.github.io/SeawaterPolynomials.jl/dev/#The-TEOS-10-standard) -# - Quadratic bottom drag with a drag coefficient of 0.003 -# -# The ocean model is then initialized with the ECCO2 temperature and salinity fields for January 1, 1993. +# We build our ocean model using `ocean_simulation`, ocean = ocean_simulation(grid) -model = ocean.model -date = DateTimeProlepticGregorian(1993, 1, 1) -set!(model, T=ECCOMetadata(:temperature; date), S=ECCOMetadata(:salinity; date)) +# which uses the default `ocean.model`, + +ocean.model + +# We initialize the ocean model to ECCO2 temperature and salinity for January 1, 1993. + +date = DateTimeProlepticGregorian(1993, 1, 1) +set!(ocean.model, T=ECCOMetadata(:temperature; dates=date), S=ECCOMetadata(:salinity; dates=date)) # ### Prescribed atmosphere and radiation # -# The atmospheric data is prescribed using the JRA55 dataset, which is loaded -# into memory in 4 snapshots at a time. The JRA55 dataset provides atmospheric +# 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, + +## 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 +## (calculated from the ratio of maximum possible incident solar radiation to actual +## incident solar radiation) and latitude. The ocean emissivity is set to 0.97. + +radiation = Radiation(arch) + +# The atmospheric data is prescribed using the JRA55 dataset. +# The number of snapshots that are loaded into memory is determined by +# the `backend` + +# into memory in 41 snapshots at a time. The JRA55 dataset provides atmospheric # data such as temperature, humidity, and wind fields to calculate turbulent fluxes # using bulk formulae, see [`CrossRealmFluxes`](@ref). -# -# 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 -# (calculated from the ratio of maximum possible incident solar radiation to actual -# incident solar radiation) and latitude. The ocean emissivity is set to 0.97. - -backend = JRA55NetCDFBackend(41) -atmosphere = JRA55_prescribed_atmosphere(arch; backend) -radiation = Radiation(arch) -nothing #hide + +atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(41)) # ### Sea ice model # # This simulation includes a simplified representation of ice cover where the # air-sea fluxes are shut down whenever the sea surface temperature is below -# the freezing point. Only heating fluxes are allowed. This is not a full -# sea ice model, but it prevents the temperature from dropping excessively -# low by including atmosphere-ocean fluxes. +# the freezing point, sea_ice = ClimaOcean.OceanSeaIceModels.MinimumTemperatureSeaIce() -nothing #hide # ## The coupled simulation -# -# Finally, we define the coupled model, which includes the ocean, atmosphere, -# and radiation parameters. The model is constructed using the `OceanSeaIceModel` -# constructor. -# + +# Next we assemble the ocean, sea ice, atmosphere, and radiation +# into a coupled model, + +coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) + # We then create a coupled simulation, starting with a time step of 10 seconds # and running the simulation for 10 days. -# We will eventually increase the time step size and end time as the simulation -# progresses and initialization shocks dissipate. -# -# We also define a callback function to monitor the simulation's progress. -# This function prints the current time, iteration, time step, -# as well as the maximum velocities and tracers in the domain. The wall time -# is also printed to monitor the time taken for each iteration. -coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -coupled_simulation = Simulation(coupled_model; Δt=10, stop_time = 10days) +simulation = Simulation(coupled_model; Δt=90, stop_time=10days) -wall_time = [time_ns()] +# We define a callback function to monitor the simulation's progress, + +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(T) - Tmin = minimum(T) - umax = maximum(abs, u), maximum(abs, v), maximum(abs, w) - step_time = 1e-9 * (time_ns() - wall_time[1]) + 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(vel): (%.2e, %.2e, %.2e), max(T): %.2f, min(T): %.2f, wtime: %s \n", - prettytime(ocean.model.clock.time), - ocean.model.clock.iteration, - prettytime(ocean.Δt), - umax..., Tmax, Tmin, prettytime(step_time)) + msg = @sprintf("Iter: %d, time: %s, Δt: %s", iteration(sim), prettytime(sim), prettytime(sim.Δt)) + msg *= @sprintf(", max|u|: (%.2e, %.2e, %.2e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s", + umax..., Tmax, Tmin, prettytime(step_time)) - wall_time[1] = time_ns() + @info msg + + wall_time[] = time_ns() end -coupled_simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000)) -nothing #hide +simulation.callbacks[:progress] = Callback(progress, TimeInterval(5days)) # ### Set up output writers # @@ -170,118 +164,80 @@ nothing #hide # 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` -ocean.output_writers[:surface] = JLD2OutputWriter(model, merge(model.tracers, model.velocities); +outputs = merge(ocean.model.tracers, ocean.model.velocities) +ocean.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs; schedule = TimeInterval(1days), - filename = "surface", + filename = "near_global_surface_fields", indices = (:, :, grid.Nz), + with_halos = true, overwrite_existing = true, array_type = Array{Float32}) -nothing #hide # ### Spinning up the simulation # -# As an initial condition, we have interpolated ECCO tracer fields onto our custom grid. -# The bathymetry of the original ECCO data may differ from our grid, so the initialization of the velocity -# field might cause shocks if a large time step is used. -# -# Therefore, we spin up the simulation with a small time step to ensure that the interpolated initial -# conditions adapt to the model numerics and parameterization without causing instability. A 10-day -# integration with a maximum time step of 1.5 minutes should be sufficient to dissipate spurious -# initialization shocks. -# We use an adaptive time step that maintains the [CFL condition](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition) equal to 0.1. -# For this scope, we use the Oceananigans utility `conjure_time_step_wizard!` (see Oceanigans's documentation). - -ocean.stop_time = 10days -conjure_time_step_wizard!(ocean; cfl = 0.1, max_Δt = 90, max_change = 1.1) -run!(coupled_simulation) -nothing #hide +# 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 +# satisfy a different dynamical balance than our simulation. -# ### Running the simulation -# -# Now that the simulation has spun up, we can run it for the full 100 days. -# We increase the maximum time step size to 10 minutes and let the simulation run for 100 days. -# This time, we set the CFL in the time_step_wizard to be 0.25 as this is the maximum recommended CFL to be -# used in conjunction with Oceananigans' hydrostatic time-stepping algorithm ([two step Adams-Bashfort](https://en.wikipedia.org/wiki/Linear_multistep_method)) - -ocean.stop_time = 100days -coupled_simulation.stop_time = 100days -conjure_time_step_wizard!(ocean; cfl = 0.25, max_Δt = 10minutes, max_change = 1.1) -run!(coupled_simulation) -nothing #hide +run!(simulation) + +# ### Running the simulation for real + +simulation.stop_time = 60days +simulation.Δt = 10minutes +run!(simulation) -# ## Visualizing the results -# -# The simulation has finished, let's visualize the results. -# In this section we pull up the saved data and create visualizations using the CairoMakie.jl package. -# In particular, we generate an animation of the evolution of surface fields: -# surface speed (s), surface temperature (T), and turbulent kinetic energy (e). +# ## A pretty movie +# +# It's time to make a pretty movie of the simulation. First we plot a snapshot: -u = FieldTimeSeries("surface.jld2", "u"; backend = OnDisk()) -v = FieldTimeSeries("surface.jld2", "v"; backend = OnDisk()) -T = FieldTimeSeries("surface.jld2", "T"; backend = OnDisk()) -e = FieldTimeSeries("surface.jld2", "e"; backend = OnDisk()) +u = FieldTimeSeries("near_global_surface_fields.jld2", "u"; backend = OnDisk()) +v = FieldTimeSeries("near_global_surface_fields.jld2", "v"; backend = OnDisk()) +T = FieldTimeSeries("near_global_surface_fields.jld2", "T"; backend = OnDisk()) +e = FieldTimeSeries("near_global_surface_fields.jld2", "e"; backend = OnDisk()) -# Set land values to NaN using the hacky method of deducing where T = 0 times = u.times -Nt = length(u) -T1 = interior(T[1], :, :, 1) -land = T1 .== 0 -for n = 1:Nt - ui = interior(u[n], :, :, 1) - vi = interior(v[n], :, :, 1) - Ti = interior(T[n], :, :, 1) - ei = interior(e[n], :, :, 1) - ui[land] .= NaN - vi[land] .= NaN - Ti[land] .= NaN - ei[land] .= NaN -end +Nt = length(times) + +n = Observable(Nt) -iter = Observable(Nt) +Tn = @lift interior(T[$n], :, :, 1) +en = @lift interior(e[$n], :, :, 1) -Ti = @lift T[$iter] -ei = @lift e[$iter] +un = Field{Face, Center, Nothing}(u.grid) +vn = Field{Center, Face, Nothing}(v.grid) +s = Field(sqrt(un^2 + vn^2)) -si = @lift begin - s = Field(sqrt(u[$iter]^2 + v[$iter]^2)) - compute!(s) +sn = @lift begin + parent(un) .= parent(u[$n]) + parent(vn) .= parent(v[$n]) + compute!(s) end -fig = Figure(size = (800, 400)) -ax = Axis(fig[1, 1]) -hm = heatmap!(ax, si, colorrange = (0, 0.5), colormap = :deep) -cb = Colorbar(fig[0, 1], hm, vertical = false, label = "Surface speed (ms⁻¹)") -hidedecorations!(ax) +fig = Figure(size = (800, 1200)) -CairoMakie.record(fig, "near_global_ocean_surface_s.mp4", 1:Nt, framerate = 8) do i - iter[] = i -end -nothing #hide - - # ![](near_global_ocean_surface_s.mp4) - -fig = Figure(size = (800, 400)) -ax = Axis(fig[1, 1]) -hm = heatmap!(ax, Ti, colorrange = (-1, 30), colormap = :magma) -cb = Colorbar(fig[0, 1], hm, vertical = false, label = "Surface Temperature (Cᵒ)") -hidedecorations!(ax) +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)") -CairoMakie.record(fig, "near_global_ocean_surface_T.mp4", 1:Nt, framerate = 8) do i - iter[] = i -end +hm = heatmap!(axs, sn, colorrange = (0, 0.5), colormap = :deep) +Colorbar(fig[1, 2], hm, label = "Surface speed (m s⁻¹)") + +hm = heatmap!(axT, Tn, colorrange = (-1, 30), colormap = :magma) +Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)") + +hm = heatmap!(axe, en, colorrange = (0, 1e-3), colormap = :solar) +Colorbar(fig[3, 2], hm, label = "Turbulent Kinetic Energy (m² s⁻²)") +save("snapshot.png", fig) nothing #hide - -# ![](near_global_ocean_surface_T.mp4) -fig = Figure(size = (800, 400)) -ax = Axis(fig[1, 1]) -hm = heatmap!(ax, ei, colorrange = (0, 1e-3), colormap = :solar) -cb = Colorbar(fig[0, 1], hm, vertical = false, label = "Turbulent Kinetic Energy (m²s⁻²)") -hidedecorations!(ax) +# ![](snapshot.png) + +# And now a movie: -CairoMakie.record(fig, "near_global_ocean_surface_e.mp4", 1:Nt, framerate = 8) do i - iter[] = i +record(fig, "near_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn + n[] = nn end nothing #hide -# ![](near_global_ocean_surface_e.mp4) +# ![](near_global_ocean_surface.mp4) diff --git a/examples/single_column_simulation.jl b/examples/single_column_os_papa_simulation.jl similarity index 76% rename from examples/single_column_simulation.jl rename to examples/single_column_os_papa_simulation.jl index 904c3beee..04775b790 100644 --- a/examples/single_column_simulation.jl +++ b/examples/single_column_os_papa_simulation.jl @@ -14,15 +14,13 @@ # pkg"add Oceananigans, ClimaOcean, CairoMakie" # ``` +using ClimaOcean + using Oceananigans using Oceananigans.Units using Oceananigans.BuoyancyModels: buoyancy_frequency using Oceananigans.Units: Time -using ClimaOcean -using ClimaOcean.ECCO -using ClimaOcean.OceanSimulations - using CairoMakie using Printf @@ -31,78 +29,78 @@ using Printf # First, we construct a single column grid with 2 meter spacing # located at ocean station Papa. -Nz = 200 -H = 400 - # Ocean station papa location location_name = "ocean_station_papa" -λ★, φ★ = 35.1, 50.1 -longitude = λ★ .+ (-0.25, 0.25) -latitude = φ★ .+ (-0.25, 0.25) +λ★, φ★ = 35.1, 50.1 -grid = LatitudeLongitudeGrid(size = (3, 3, Nz); - longitude, latitude, z = (-H, 0), - topology = (Periodic, Periodic, Bounded)) +grid = RectilinearGrid(size = 200, + x = λ★, + y = φ★, + z = (-400, 0), + topology = (Flat, Flat, Bounded)) # # An "ocean simulation" # # Next, we use ClimaOcean's ocean_simulation constructor to build a realistic -# ocean simulation on the single column grid. +# ocean simulation on the single column grid, + +ocean = ocean_simulation(grid; Δt=10minutes, coriolis=FPlane(latitude = φ★)) -ocean = ocean_simulation(grid; - Δt = 10minutes, - coriolis = FPlane(latitude = φ★), - tracer_advection = nothing, - momentum_advection = nothing, - bottom_drag_coefficient = 0) +# which wraps around the ocean model -start_time = time_ns() +ocean.model -# Initial conditions -set!(ocean.model, T = ECCOMetadata(:temperature), - S = ECCOMetadata(:salinity), - e = 1e-6) +# We set initial conditions from ECCO: -elapsed = time_ns() - start_time -@info "Initial condition built. " * prettytime(elapsed * 1e-9) -start_time = time_ns() +set!(ocean.model, T=ECCOMetadata(:temperature), S=ECCOMetadata(:salinity)) # # A prescribed atmosphere based on JRA55 re-analysis # # We build a PrescribedAtmosphere at the same location as the single colunm grid # which is based on the JRA55 reanalysis. -last_time = floor(Int, 31 * 24 / 3) # 31 days in hours divided by JRA55's frequency in hours -backend = InMemory() -atmosphere = JRA55_prescribed_atmosphere(time_indices = 1:last_time; - longitude, latitude, backend, +simulation_days = 31 +snapshots_per_day = 8 # corresponding to JRA55's 3-hour frequency +last_time = simulation_days * snapshots_per_day +atmosphere = JRA55_prescribed_atmosphere(1:last_time; + longitude = λ★, + latitude = φ★, + #longitude = (λ★ - 1/4, λ★ + 1/4), + #latitude = (φ★ - 1/4, φ★ + 1/4), + backend = InMemory(), include_rivers_and_icebergs = false) -ua = atmosphere.velocities.u -va = atmosphere.velocities.v -Ta = atmosphere.tracers.T -qa = atmosphere.tracers.q -times = ua.times +# This builds a representation of the atmosphere on the small grid -fig = Figure(size=(1200, 1800)) -axu = Axis(fig[1, 1]) -axT = Axis(fig[2, 1]) -axq = Axis(fig[3, 1]) +atmosphere.grid -lines!(axu, times ./ days, interior(ua, 1, 1, 1, :)) -lines!(axu, times ./ days, interior(va, 1, 1, 1, :)) -lines!(axT, times ./ days, interior(Ta, 1, 1, 1, :)) -lines!(axq, times ./ days, interior(qa, 1, 1, 1, :)) +# Let's take a look at the atmospheric state + +ua = interior(atmosphere.velocities.u, 1, 1, 1, :) +va = interior(atmosphere.velocities.v, 1, 1, 1, :) +Ta = interior(atmosphere.tracers.T, 1, 1, 1, :) +qa = interior(atmosphere.tracers.q, 1, 1, 1, :) +t_days = atmosphere.times / days + +fig = Figure(size=(800, 600)) +axu = Axis(fig[2, 1], xlabel="Days since Jan 1 1990", ylabel="Atmosphere \n velocity (m s⁻¹)") +axT = Axis(fig[3, 1], xlabel="Days since Jan 1 1990", ylabel="Atmosphere \n temperature (K)") +axq = Axis(fig[4, 1], xlabel="Days since Jan 1 1990", ylabel="Atmosphere \n specific humidity") +Label(fig[1, 1], "Atmospheric state over ocean station Papa", tellwidth=false) + +lines!(axu, t_days, ua, label="Zonal velocity") +lines!(axu, t_days, va, label="Meridional velocity") +ylims!(axu, -6, 6) +axislegend(axu, framevisible=false, nbanks=2, position=:lb) + +lines!(axT, t_days, Ta) +lines!(axq, t_days, qa) display(fig) radiation = Radiation() coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) -coupled_simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) - -elapsed = time_ns() - start_time -@info "Coupled simulation built. " * prettytime(elapsed * 1e-9) -start_time = time_ns() +simulation = Simulation(coupled_model, Δt=ocean.Δt, stop_time=30days) wall_clock = Ref(time_ns()) @@ -129,16 +127,16 @@ function progress(sim) Nz = size(T, 3) msg *= @sprintf(", u★: %.2f m s⁻¹", u★) - msg *= @sprintf(", Q: %.2f W m⁻²", Q) - msg *= @sprintf(", T₀: %.2f ᵒC", first(interior(T, 1, 1, Nz))) + msg *= @sprintf(", Q: %.2f W m⁻²", Q) + msg *= @sprintf(", T₀: %.2f ᵒC", first(interior(T, 1, 1, Nz))) msg *= @sprintf(", extrema(T): (%.2f, %.2f) ᵒC", minimum(T), maximum(T)) - msg *= @sprintf(", S₀: %.2f g/kg", first(interior(S, 1, 1, Nz))) + msg *= @sprintf(", S₀: %.2f g/kg", first(interior(S, 1, 1, Nz))) msg *= @sprintf(", e₀: %.2e m² s⁻²", first(interior(e, 1, 1, Nz))) @info msg end -coupled_simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) +simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) # Build flux outputs τx = coupled_model.fluxes.total.ocean.momentum.u @@ -151,7 +149,7 @@ Qv = coupled_model.fluxes.turbulent.fields.latent_heat ρₒ = coupled_model.fluxes.ocean_reference_density cₚ = coupled_model.fluxes.ocean_heat_capacity -Q = ρₒ * cₚ * JT +Q = ρₒ * cₚ * JT ρτx = ρₒ * τx ρτy = ρₒ * τy N² = buoyancy_frequency(ocean.model) @@ -166,11 +164,11 @@ outputs = merge(fields, fluxes) filename = "single_column_omip_$(location_name)" -coupled_simulation.output_writers[:jld2] = JLD2OutputWriter(ocean.model, outputs; filename, - schedule = TimeInterval(3hours), - overwrite_existing = true) +simulation.output_writers[:jld2] = JLD2OutputWriter(ocean.model, outputs; filename, + schedule = TimeInterval(3hours), + overwrite_existing = true) -run!(coupled_simulation) +run!(simulation) filename *= ".jld2" @@ -329,4 +327,4 @@ record(fig, "single_column_profiles.mp4", 1:8:Nt, framerate=24) do nn end nothing #hide -# ![](single_column_profiles.mp4) \ No newline at end of file +# ![](single_column_profiles.mp4) diff --git a/experiments/omip_components.jl b/experiments/omip_components.jl index 9f354bcb3..616623fe0 100644 --- a/experiments/omip_components.jl +++ b/experiments/omip_components.jl @@ -26,7 +26,7 @@ function regional_omip_grid(arch, ecco_2_temperature_field; Te = ecco_2_temperature_field launch!(architecture(Te), Te.grid, :xyz, nan_land!, Te) - + ΔΛ = last(longitude) - first(longitude) ΔΦ = last(latitude) - first(latitude) Nx = Int(ΔΛ / resolution) @@ -36,7 +36,7 @@ function regional_omip_grid(arch, ecco_2_temperature_field; grid = LatitudeLongitudeGrid(arch; latitude, longitude, halo, size = (Nx, Ny, Nz), z = z) - + Tᵢ = CenterField(grid) interpolate!(Tᵢ, Te) @@ -45,7 +45,7 @@ function regional_omip_grid(arch, ecco_2_temperature_field; # Construct bottom_height depth by analyzing T launch!(arch, grid, :xy, infer_bottom_height!, bottom_height, Tᵢ, grid) - + grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height)) elapsed = 1e-9 * (time_ns() - start_time) @@ -105,7 +105,7 @@ function omip_ocean_component(grid; tracers = tuple(:e, tracers...) end - + ocean_model = HydrostaticFreeSurfaceModel(; grid, buoyancy, closure, tracers, tracer_advection, momentum_advection, free_surface = SplitExplicitFreeSurface(cfl=0.7; grid), @@ -202,4 +202,3 @@ end T[i, j, k] = ifelse(land, NaN, Tᵢ) end end - diff --git a/experiments/plot_freely_decaying_simulation.jl b/experiments/plot_freely_decaying_simulation.jl index 935035c8b..49ce5de04 100644 --- a/experiments/plot_freely_decaying_simulation.jl +++ b/experiments/plot_freely_decaying_simulation.jl @@ -18,7 +18,7 @@ for n = 1:Nt hn[land] .= NaN end -fig = Figure(resolution=(1200, 600)) +fig = Figure(size=(1200, 600)) axT = Axis(fig[1, 1], xlabel="Longitude", ylabel="Latitude") axh = Axis(fig[2, 1], xlabel="Longitude", ylabel="Latitude") diff --git a/src/Bathymetry.jl b/src/Bathymetry.jl index 76bbbb101..586f282c0 100644 --- a/src/Bathymetry.jl +++ b/src/Bathymetry.jl @@ -46,36 +46,36 @@ Arguments: Keyword Arguments: ================== -- height_above_water: limits the maximum height of above-water topography (where h > 0). If - `nothing` 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 - -- dir: directory of the bathymetry-containing file - -- 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 will involve - - 400x200 -> 325x175 - - 325x175 -> 250x150 - - 250x150 -> 175x125 - - 175x125 -> 100x100 - 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. - -- 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. +- `height_above_water`: limits the maximum height of above-water topography (where h > 0). If + `nothing` 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 + +- `dir`: directory of the bathymetry-containing file + +- `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 will involve + - 400x200 -> 325x175 + - 325x175 -> 250x150 + - 250x150 -> 175x125 + - 175x125 -> 100x100 + 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. + +- `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. """ function regrid_bathymetry(target_grid; height_above_water = nothing, diff --git a/src/DataWrangling/JRA55.jl b/src/DataWrangling/JRA55.jl index 4ead552ac..978b7755d 100644 --- a/src/DataWrangling/JRA55.jl +++ b/src/DataWrangling/JRA55.jl @@ -134,6 +134,11 @@ urls = Dict( compute_bounding_nodes(::Nothing, ::Nothing, LH, hnodes) = nothing compute_bounding_nodes(bounds, ::Nothing, LH, hnodes) = bounds +function compute_bounding_nodes(x::Number, ::Nothing, LH, hnodes) + ϵ = convert(typeof(x), 0.001) # arbitrary? + return (x - ϵ, x + ϵ) +end + # TODO: remove the allowscalar function compute_bounding_nodes(::Nothing, grid, LH, hnodes) hg = hnodes(grid, LH()) @@ -147,7 +152,7 @@ function compute_bounding_indices(::Nothing, hc) return 1, Nh end -function compute_bounding_indices(bounds, hc) +function compute_bounding_indices(bounds::Tuple, hc) h₁, h₂ = bounds Nh = length(hc) @@ -256,14 +261,16 @@ new_backend(::JRA55NetCDFBackend, start, length) = JRA55NetCDFBackend(start, len """ JRA55_field_time_series(variable_name; architecture = CPU(), + time_indices = nothing, + latitude = nothing, + longitude = nothing, location = nothing, url = nothing, filename = nothing, shortname = nothing, backend = InMemory(), preprocess_chunk_size = 10, - preprocess_architecture = CPU(), - time_indices = nothing) + preprocess_architecture = CPU()) Return a `FieldTimeSeries` containing atmospheric reanalysis data for `variable_name`, which describes one of the variables in the "repeat year forcing" dataset derived @@ -300,6 +307,14 @@ Keyword arguments `FieldTimeSeries` with the first three time snapshots of `variable_name`. + - `latitude`: Guiding latitude bounds for the resulting grid. + Used to slice the data when loading into memory. + Default: nothing, which retains the latitude range of the native grid. + + - `longitude`: Guiding longitude bounds for the resulting grid. + Used to slice the data when loading into memory. + Default: nothing, which retains the longitude range of the native grid. + - `url`: The url accessed to download the data for `variable_name`. Default: `ClimaOcean.JRA55.urls[variable_name]`. diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl index 11479dda9..fde3a09bc 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/CrossRealmFluxes.jl @@ -26,6 +26,7 @@ function surface_flux(f::Field) end include("radiation.jl") +include("latitude_dependent_albedo.jl") include("tabulated_albedo.jl") include("roughness_lengths.jl") include("stability_functions.jl") diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index ef977af2e..46d5b5635 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -1,4 +1,5 @@ using Oceananigans.Operators: extrinsic_vector, intrinsic_vector +using Oceananigans.Grids: _node # Fallback! limit_fluxes_over_sea_ice!(args...) = nothing @@ -49,8 +50,14 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) atmosphere_time_indexing = u.time_indexing # kernel parameters that compute fluxes in 0:Nx+1 and 0:Ny+1 - kernel_size = (size(grid, 1) + 2, size(grid, 2) + 2) - kernel_parameters = KernelParameters(kernel_size, (-1, -1)) + Nx, Ny, Nz = size(grid) + single_column_grid = Nx == 1 && Ny == 1 + + if single_column_grid + kernel_parameters = KernelParameters(1:1, 1:1) + else + kernel_parameters = KernelParameters(0:Nx+1, 0:Ny+1) + end surface_atmosphere_state = coupled_model.fluxes.surface_atmosphere_state @@ -88,7 +95,7 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) if !isnothing(auxiliary_freshwater_flux) # TODO: do not assume that `auxiliary_freshater_flux` is a tuple - auxiliary_data = map(ϕ -> ϕ.data, auxiliary_freshwater_flux) + auxiliary_data = map(ϕ -> ϕ.data, auxiliary_freshwater_flux) first_auxiliary_flux = first(auxiliary_freshwater_flux) auxiliary_grid = first_auxiliary_flux.grid @@ -198,7 +205,7 @@ end # The third index "k" should not matter but we put the correct index to get # a surface node anyways. atmos_args = (atmos_grid, atmos_times, atmos_backend, atmos_time_indexing) - X = node(i, j, kᴺ + 1, grid, c, c, f) + X = _node(i, j, kᴺ + 1, grid, c, c, f) time = Time(clock.time) uₐ = interp_atmos_time_series(atmos_velocities.u, X, time, atmos_args...) @@ -239,7 +246,7 @@ end kᴺ = size(grid, 3) # index of the top ocean cell @inbounds begin - X = node(i, j, kᴺ + 1, grid, c, c, f) + X = _node(i, j, kᴺ + 1, grid, c, c, f) time = Time(clock.time) Mr = interp_atmos_time_series(auxiliary_freshwater_flux, X, time, auxiliary_grid, diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl b/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl new file mode 100644 index 000000000..dbe91e775 --- /dev/null +++ b/src/OceanSeaIceModels/CrossRealmFluxes/latitude_dependent_albedo.jl @@ -0,0 +1,53 @@ +using Oceananigans.Grids: ηnode + +struct LatitudeDependentAlbedo{FT} + direct :: FT + diffuse :: FT +end + +""" + LatitudeDependentAlbedo([FT::DataType=Float64]; diffuse = 0.069, direct = 0.011) + +Constructs a `LatitudeDependentAlbedo` object. The albedo of the ocean surface is assumed to be a function of the latitude, +obeying the following formula (Large and Yeager, 2009): + + α(φ) = α.diffuse - α.direct * cos(2φ) + +where `φ` is the latitude, `α.diffuse` is the diffuse albedo, and `α.direct` is the direct albedo. + +# Arguments +=========== + +- `FT::DataType`: The data type of the albedo values. Default is `Float64`. + +# Keyword Arguments +=================== +- `diffuse`: The diffuse albedo value. Default is `0.069`. +- `direct`: The direct albedo value. Default is `0.011`. +""" +function LatitudeDependentAlbedo(FT::DataType=Float64; + diffuse = 0.069, + direct = 0.011) + + return LatitudeDependentAlbedo(convert(FT, direct), + convert(FT, diffuse)) +end + +Adapt.adapt_structure(to, α::LatitudeDependentAlbedo) = + LatitudeDependentAlbedo(Adapt.adapt(to, α.direct), + Adapt.adapt(to, α.diffuse)) + +function Base.summary(α::LatitudeDependentAlbedo{FT}) where FT + formula = string(prettysummary(α.diffuse), " - ", prettysummary(α.direct), " ⋅ cos(2φ)") + return string("LatitudeDepedendentAlbedo{$FT}: ", formula) +end + +Base.show(io::IO, α::LatitudeDependentAlbedo) = print(io, summary(α)) + +@inline function stateindex(α::LatitudeDependentAlbedo, i, j, k, grid, time) + φ = ηnode(i, j, k, grid, Center(), Center(), Center()) + α₀ = α.diffuse + α₁ = α.direct + return α₀ - α₁ * hack_cosd(2φ) +end + diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl index e0998b517..9cb303126 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/radiation.jl @@ -18,7 +18,7 @@ Adapt.adapt_structure(to, r :: Radiation) = Radiation([arch = CPU(), FT=Float64]; ocean_emissivity = 0.97, sea_ice_emissivity = 1.0, - ocean_albedo = TabulatedAlbedo(arch, FT), + ocean_albedo = LatitudeDependentAlbedo(FT), sea_ice_albedo = 0.7, stefan_boltzmann_constant = 5.67e-8) @@ -42,7 +42,7 @@ Constructs a `Radiation` object that represents the radiation properties of the function Radiation(arch = CPU(), FT=Float64; ocean_emissivity = 0.97, sea_ice_emissivity = 1.0, - ocean_albedo = TabulatedAlbedo(arch, FT), + ocean_albedo = LatitudeDependentAlbedo(FT), sea_ice_albedo = 0.7, stefan_boltzmann_constant = 5.67e-8) @@ -59,52 +59,19 @@ function Radiation(arch = CPU(), FT=Float64; convert(FT, stefan_boltzmann_constant)) end -Base.summary(r::Radiation) = "Radiation" -Base.show(io::IO, r::Radiation) = print(io, summary(r)) +Base.summary(r::Radiation{FT}) where FT = "Radiation{$FT}" -struct LatitudeDependentAlbedo{FT} - direct :: FT - diffuse :: FT -end - -""" - LatitudeDependentAlbedo([FT::DataType=Float64]; diffuse = 0.069, direct = 0.011) - -Constructs a `LatitudeDependentAlbedo` object. The albedo of the ocean surface is assumed to be a function of the latitude, -obeying the following formula (Large and Yeager, 2009): - - α(φ) = α.diffuse - α.direct * cos(2φ) - -where `φ` is the latitude, `α.diffuse` is the diffuse albedo, and `α_.irect` is the direct albedo. +function Base.show(io::IO, r::Radiation) + σ = r.stefan_boltzmann_constant -# Arguments -=========== - -- `FT::DataType`: The data type of the albedo values. Default is `Float64`. - -# Keyword Arguments -=================== -- `diffuse`: The diffuse albedo value. Default is `0.069`. -- `direct`: The direct albedo value. Default is `0.011`. -""" -function LatitudeDependentAlbedo(FT::DataType=Float64; - diffuse = 0.069, - direct = 0.011) - - return LatitudeDependentAlbedo(convert(FT, direct), - convert(FT, diffuse)) -end - -Adapt.adapt_structure(to, α::LatitudeDependentAlbedo) = - LatitudeDependentAlbedo(Adapt.adapt(to, α.direct), - Adapt.adapt(to, α.diffuse)) - -@inline function stateindex(α::LatitudeDependentAlbedo, i, j, k, grid, time) - φ = φnode(i, j, k, grid, Center(), Center(), Center()) - α_diffuse = α.diffuse - direct_correction = α.direct * hack_cosd(2φ) - - return α_diffuse - direct_correction + print(io, summary(r), ":", '\n') + print(io, "├── stefan_boltzmann_constant: ", prettysummary(σ), '\n') + print(io, "├── emission: ", summary(r.emission), '\n') + print(io, "│ ├── ocean: ", prettysummary(r.emission.ocean), '\n') + print(io, "│ └── sea_ice: ", prettysummary(r.emission.ocean), '\n') + print(io, "└── reflection: ", summary(r.reflection), '\n') + print(io, " ├── ocean: ", prettysummary(r.reflection.ocean), '\n') + print(io, " └── sea_ice: ", prettysummary(r.reflection.sea_ice)) end struct SurfaceProperties{O, I} @@ -115,3 +82,12 @@ end Adapt.adapt_structure(to, s :: SurfaceProperties) = SurfaceProperties(Adapt.adapt(to, s.ocean), Adapt.adapt(to, s.sea_ice)) + +Base.summary(properties::SurfaceProperties) = "SurfaceProperties" + +function Base.show(io::IO, properties::SurfaceProperties) + print(io, "SurfaceProperties:", '\n') + print(io, "├── ocean: ", summary(properties.ocean), '\n') + print(io, "└── sea_ice: ", summary(properties.sea_ice)) +end + diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/tabulated_albedo.jl b/src/OceanSeaIceModels/CrossRealmFluxes/tabulated_albedo.jl index 59379da5e..2f42e0d2c 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/tabulated_albedo.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/tabulated_albedo.jl @@ -9,7 +9,7 @@ using ClimaOcean.OceanSeaIceModels: # Bilinear interpolation of the albedo α in α_table based on a # transmissivity value (𝓉_values) and latitude (φ_values) -struct TabulatedAlbedo{M, P, T, FT} +struct TabulatedAlbedo{FT, M, P, T} α_table :: M φ_values :: P 𝓉_values :: T @@ -18,7 +18,7 @@ struct TabulatedAlbedo{M, P, T, FT} noon_in_seconds :: Int end -Adapt.adapt_structure(to, α :: TabulatedAlbedo) = +Adapt.adapt_structure(to, α::TabulatedAlbedo) = TabulatedAlbedo(Adapt.adapt(to, α.α_table), Adapt.adapt(to, α.φ_values), Adapt.adapt(to, α.𝓉_values), @@ -83,18 +83,24 @@ function TabulatedAlbedo(arch = CPU(), FT = Float64; φ_values = (0:2:90) ./ 180 * π, 𝓉_values = 0:0.05:1, day_to_radians = convert(FT, 2π / 86400), - noon_in_seconds = 86400 ÷ 2 # assumes that midnight is at t = 0 seconds - ) + noon_in_seconds = 86400 ÷ 2) # assumes that midnight is at t = 0 seconds # Make everything GPU - ready α_table = on_architecture(arch, convert.(FT, α_table)) φ_values = on_architecture(arch, convert.(FT, φ_values)) 𝓉_values = on_architecture(arch, convert.(FT, 𝓉_values)) - return TabulatedAlbedo(α_table, φ_values, 𝓉_values, convert(FT, S₀), convert(FT, day_to_radians), noon_in_seconds) + return TabulatedAlbedo(α_table, + φ_values, + 𝓉_values, + convert(FT, S₀), + convert(FT, day_to_radians), + noon_in_seconds) end -Base.eltype(α::TabulatedAlbedo) = Base.eltype(α.S₀) +Base.eltype(::TabulatedAlbedo{FT}) where FT = FT +Base.summary(::TabulatedAlbedo{FT}) where FT = "TabulatedAlbedo{$FT}" +Base.show(io::IO, α::TabulatedAlbedo) = print(io, summary(α)) @inline ϕ₁(ξ, η) = (1 - ξ) * (1 - η) @inline ϕ₂(ξ, η) = (1 - ξ) * η @@ -109,10 +115,8 @@ Base.eltype(α::TabulatedAlbedo) = Base.eltype(α.S₀) @inline function net_downwelling_radiation(i, j, grid, time, radiation::Radiation{<:Any, <:Any, <:SurfaceProperties{<:TabulatedAlbedo}}, Qs, Qℓ) α = radiation.reflection.ocean - FT = eltype(α) - - λ, φ, z = node(i, j, 1, grid, Center(), Center(), Center()) + λ, φ, z = _node(i, j, 1, grid, Center(), Center(), Center()) φ = deg2rad(φ) λ = deg2rad(λ) @@ -162,3 +166,4 @@ Base.eltype(α::TabulatedAlbedo) = Base.eltype(α.S₀) return - (1 - α) * Qs - ϵ * Qℓ end + diff --git a/src/OceanSeaIceModels/PrescribedAtmospheres.jl b/src/OceanSeaIceModels/PrescribedAtmospheres.jl index a288c1480..e99ca44bc 100644 --- a/src/OceanSeaIceModels/PrescribedAtmospheres.jl +++ b/src/OceanSeaIceModels/PrescribedAtmospheres.jl @@ -1,5 +1,6 @@ module PrescribedAtmospheres +using Oceananigans.Grids: grid_name using Oceananigans.Utils: prettysummary using Oceananigans.OutputReaders: update_field_time_series!, extract_field_time_series @@ -285,7 +286,7 @@ const PATP = PrescribedAtmosphereThermodynamicsParameters ##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic) ##### -struct PrescribedAtmosphere{G, U, P, C, F, I, R, TP, TI, FT} +struct PrescribedAtmosphere{FT, G, U, P, C, F, I, R, TP, TI} grid :: G velocities :: U pressure :: P @@ -299,8 +300,19 @@ struct PrescribedAtmosphere{G, U, P, C, F, I, R, TP, TI, FT} boundary_layer_height :: FT end -Base.summary(::PrescribedAtmosphere) = "PrescribedAtmosphere" -Base.show(io::IO, pa::PrescribedAtmosphere) = print(io, summary(pa)) +function Base.summary(pa::PrescribedAtmosphere{FT}) where FT + Nx, Ny, Nz = size(pa.grid) + Nt = length(pa.times) + sz_str = string(Nx, "×", Ny, "×", Nz, "×", Nt) + return string(sz_str, " PrescribedAtmosphere{$FT}") +end + +function Base.show(io::IO, pa::PrescribedAtmosphere) + print(io, summary(pa), " on ", grid_name(pa.grid), ":", '\n') + print(io, "├── times: ", prettysummary(pa.times), '\n') + print(io, "├── reference_height: ", prettysummary(pa.reference_height), '\n') + print(io, "└── boundary_layer_height: ", prettysummary(pa.reference_height)) +end """ PrescribedAtmosphere(times; diff --git a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl index e3ef045f0..954f7b96a 100644 --- a/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl +++ b/src/OceanSeaIceModels/minimum_temperature_sea_ice.jl @@ -63,3 +63,4 @@ end τy[i, j, 1] = ifelse(sea_ice, zero(grid), τy[i, j, 1]) end end + diff --git a/src/OceanSeaIceModels/ocean_sea_ice_model.jl b/src/OceanSeaIceModels/ocean_sea_ice_model.jl index 5e1257742..09b97a797 100644 --- a/src/OceanSeaIceModels/ocean_sea_ice_model.jl +++ b/src/OceanSeaIceModels/ocean_sea_ice_model.jl @@ -34,10 +34,10 @@ function Base.summary(model::OSIM) end function Base.show(io::IO, cm::OSIM) - print(io, summary(cm)) - print(io, "├── ocean: ", summary(cm.ocean.model)) - print(io, "├── atmosphere: ", summary(cm.atmosphere)) - print(io, "└── sea_ice: ", summary(cm.sea_ice)) + print(io, summary(cm), "\n") + print(io, "├── ocean: ", summary(cm.ocean.model), "\n") + print(io, "├── atmosphere: ", summary(cm.atmosphere), "\n") + print(io, "└── sea_ice: ", summary(cm.sea_ice), "\n") return nothing end @@ -79,10 +79,10 @@ function OceanSeaIceModel(ocean, sea_ice=MinimumTemperatureSeaIce(); # Remove some potentially irksome callbacks from the ocean simulation # TODO: also remove these from sea ice simulations - pop!(ocean.callbacks, :stop_time_exceeded) - pop!(ocean.callbacks, :stop_iteration_exceeded) - pop!(ocean.callbacks, :wall_time_limit_exceeded) - pop!(ocean.callbacks, :nan_checker) + pop!(ocean.callbacks, :stop_time_exceeded, nothing) + pop!(ocean.callbacks, :stop_iteration_exceeded, nothing) + pop!(ocean.callbacks, :wall_time_limit_exceeded, nothing) + pop!(ocean.callbacks, :nan_checker, nothing) # In case there was any doubt these are meaningless. ocean.stop_time = Inf diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index dc04c53bb..46578fb16 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -19,6 +19,13 @@ using Oceananigans.BuoyancyModels: g_Earth using Oceananigans.Coriolis: Ω_Earth using Oceananigans.Operators +struct Default{V} + value :: V +end + +Default() = Default(nothing) +default_or_override(default::Default, value=default.value) = value + # Some defaults default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) @@ -58,13 +65,30 @@ function ocean_simulation(grid; Δt = 5minutes, reference_density = 1020, rotation_rate = Ω_Earth, gravitational_acceleration = g_Earth, - bottom_drag_coefficient = 0.003, + bottom_drag_coefficient = Default(0.003), forcing = NamedTuple(), coriolis = HydrostaticSphericalCoriolis(; rotation_rate), momentum_advection = default_momentum_advection(), tracer_advection = default_tracer_advection(), verbose = false) + FT = eltype(grid) + + # Detect whether we are on a single column grid + Nx, Ny, _ = size(grid) + single_column_simulation = Nx == 1 && Ny == 1 + + if single_column_simulation + # Let users put a bottom drag if they want + bottom_drag_coefficient = default_or_override(bottom_drag_coefficient, zero(grid)) + + # Don't let users use advection in a single column model + tracer_advection = nothing + momentum_advection = nothing + else + bottom_drag_coefficient = default_or_override(bottom_drag_coefficient) + end + # Set up boundary conditions using Field top_zonal_momentum_flux = τx = Field{Face, Center, Nothing}(grid) top_meridional_momentum_flux = τy = Field{Center, Face, Nothing}(grid) @@ -89,13 +113,6 @@ function ocean_simulation(grid; Δt = 5minutes, teos10 = TEOS10EquationOfState(; reference_density) buoyancy = SeawaterBuoyancy(; gravitational_acceleration, equation_of_state=teos10) - # Minor simplifications for single column grids - Nx, Ny, _ = size(grid) - if Nx == Ny == 1 # single column grid - tracer_advection = nothing - momentum_advection = nothing - end - tracers = (:T, :S) if closure isa CATKEVerticalDiffusivity tracers = tuple(tracers..., :e) diff --git a/test/runtests_setup.jl b/test/runtests_setup.jl index 304b916c1..5f3d3bd63 100644 --- a/test/runtests_setup.jl +++ b/test/runtests_setup.jl @@ -11,8 +11,5 @@ using ClimaOcean.JRA55: JRA55_field_time_series using Oceananigans.Architectures: architecture, on_architecture using Oceananigans.OutputReaders: interpolate! -using ClimaOcean - gpu_test = parse(Bool, get(ENV, "GPU_TEST", "false")) test_architectures = gpu_test ? [GPU()] : [CPU()] - diff --git a/test/test_bathymetry.jl b/test/test_bathymetry.jl index 9071f6c0a..c85910649 100644 --- a/test/test_bathymetry.jl +++ b/test/test_bathymetry.jl @@ -1,25 +1,26 @@ using Oceananigans +using Statistics using ClimaOcean + using ClimaOcean.Bathymetry: remove_minor_basins! -using Statistics @testset "Availability of Bathymetry" begin @info "Testing Bathymetry utils..." for arch in test_architectures grid = LatitudeLongitudeGrid(arch; - size = (100, 100, 10), - longitude = (0, 100), + size = (100, 100, 10), + longitude = (0, 100), latitude = (0, 50), z = (-6000, 0)) # Test that remove_minor_basins!(Z, Inf) does nothing - control_bottom_height = regrid_bathymetry(grid) + control_bottom_height = regrid_bathymetry(grid) bottom_height = deepcopy(control_bottom_height) @test_throws ArgumentError remove_minor_basins!(bottom_height, Inf) # A fictitiously large number which should presumably keep all the basins - remove_minor_basins!(bottom_height, 10000000) + remove_minor_basins!(bottom_height, 10000000) @test parent(bottom_height) == parent(control_bottom_height) # Test that remove_minor_basins!(Z, 2) remove the correct number of Basins @@ -27,7 +28,7 @@ using Statistics control_bottom_height = Field{Center, Center, Nothing}(grid) # A two basins bathymetry - bottom(x, y) = - 1000 * Int((x < 10) | (x > 50)) + bottom(x, y) = - 1000 * Int((x < 10) | (x > 50)) set!(bottom_height, bottom) set!(control_bottom_height, bottom) @@ -48,8 +49,8 @@ using Statistics @test mean(view(bottom_height, 51:100, :, 1)) == -1000 grid = LatitudeLongitudeGrid(arch; - size = (200, 200, 10), - longitude = (0, 2), + size = (200, 200, 10), + longitude = (0, 2), latitude = (-10, 50), z = (-6000, 0)) @@ -59,4 +60,4 @@ using Statistics # Testing that multiple passes do not change the solution when refining the grid @test parent(control_bottom_height) == parent(interpolated_bottom_height) end -end \ No newline at end of file +end diff --git a/test/test_ecco.jl b/test/test_ecco.jl index 24ed418f6..4c1b73748 100644 --- a/test/test_ecco.jl +++ b/test/test_ecco.jl @@ -1,13 +1,13 @@ include("runtests_setup.jl") +using CFTime +using Dates using ClimaOcean + using ClimaOcean.ECCO using ClimaOcean.ECCO: ECCO_field, metadata_path using Oceananigans.Grids: topology -using CFTime -using Dates - @testset "ECCO fields utilities" begin for arch in test_architectures A = typeof(arch) diff --git a/test/test_jra55.jl b/test/test_jra55.jl index 2c6ec39af..2aa9a1a95 100644 --- a/test/test_jra55.jl +++ b/test/test_jra55.jl @@ -108,6 +108,5 @@ using ClimaOcean.OceanSeaIceModels: PrescribedAtmosphere atmosphere = JRA55_prescribed_atmosphere(arch; backend, include_rivers_and_icebergs=true) @test haskey(atmosphere.auxiliary_freshwater_flux, :rivers) @test haskey(atmosphere.auxiliary_freshwater_flux, :icebergs) - end + end end - diff --git a/test/test_ocean_sea_ice_model_parameter_space.jl b/test/test_ocean_sea_ice_model_parameter_space.jl index 909510096..d45398ba2 100644 --- a/test/test_ocean_sea_ice_model_parameter_space.jl +++ b/test/test_ocean_sea_ice_model_parameter_space.jl @@ -4,9 +4,9 @@ using OrthogonalSphericalShellGrids start_time = time_ns() arch = GPU() -grid = TripolarGrid(arch; - size = (50, 50, 10), - halo = (7, 7, 7), +grid = TripolarGrid(arch; + size = (50, 50, 10), + halo = (7, 7, 7), z = (-6000, 0), first_pole_longitude = 75, north_poles_latitude = 55) diff --git a/test/test_simulations.jl b/test/test_simulations.jl index 161ca53e5..dbe021a95 100644 --- a/test/test_simulations.jl +++ b/test/test_simulations.jl @@ -5,9 +5,9 @@ using OrthogonalSphericalShellGrids @testset "Parameter space test" begin for arch in test_architectures - grid = TripolarGrid(arch; - size = (50, 50, 10), - halo = (7, 7, 7), + grid = TripolarGrid(arch; + size = (50, 50, 10), + halo = (7, 7, 7), z = (-6000, 0), first_pole_longitude = 75, north_poles_latitude = 55) @@ -18,12 +18,12 @@ using OrthogonalSphericalShellGrids interpolation_passes = 20, major_basins = 1) - grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true) + grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map = true) free_surface = SplitExplicitFreeSurface(grid; substeps = 20) ocean = ocean_simulation(grid; free_surface) - backend = JRA55NetCDFBackend(4) + backend = JRA55NetCDFBackend(4) atmosphere = JRA55_prescribed_atmosphere(arch; backend) radiation = Radiation(arch) @@ -36,5 +36,3 @@ using OrthogonalSphericalShellGrids end end end - -