diff --git a/Project.toml b/Project.toml index ab7f75323..92d71b4ee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,14 @@ name = "NumericalEarth" uuid = "904d977b-046a-4731-8b86-9235c0d1ef02" license = "MIT" -version = "0.3.0" +version = "0.3.1" authors = ["NumericalEarth contributors"] [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CFTime = "179af706-886a-5703-950a-314cd64e0468" ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" +CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" CubedSphere = "7445602f-e544-4518-8976-18f8e8ae6cdb" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -29,6 +30,7 @@ Scratch = "6c6a2e73-6563-6170-7368-637461726353" SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" diff --git a/docs/make.jl b/docs/make.jl index 2e5c289a4..0961faf26 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -32,6 +32,7 @@ examples = [ Example("Single-column ocean simulation", "single_column_os_papa_simulation", true), Example("One-degree ocean--sea ice simulation", "one_degree_simulation", false), Example("Near-global ocean simulation", "near_global_ocean_simulation", false), + Example("Two-degree ocean--sea ice simulation", "orca2", false), Example("Global climate simulation", "global_climate_simulation", false), Example("Veros ocean simulation", "veros_ocean_forced_simulation", false), Example("Breeze over two oceans", "breeze_over_two_oceans", false), diff --git a/examples/orca2.jl b/examples/orca2.jl new file mode 100644 index 000000000..9e79519fb --- /dev/null +++ b/examples/orca2.jl @@ -0,0 +1,339 @@ +# # [Coarse global ocean--sea ice simulation](@id coarse-degree-ocean-seaice) +# +# This example configures a global ocean--sea ice simulation on the ["ORCA2" grid](https://www.nemo-ocean.eu/doc/node108.html#Tab_orca_zgr) +# which is tripolar grid which transitions to a Mercator grid in the tropics to better resolve +# equatorial dynamics (as the sign of the coriolis coefficient switches on the equator). +# The grid nominally has a 2° resolution with meridional refinement to 0.5° in the tropics +# and near Antarctica, and refinement to ~0.5° in the Mediterranean, Red, Black and Caspian Seas. +# The grid has been refined carefully designed by NEMO to carefully capture important physics +# and maintain anisotropy close to 1 in the ocean, especially in strongly eddying regions such +# as the gulf stream. They also provide bathymetry which is refined to represent important straits. +# +# The model is forced with by repeat-year JRA55 atmospheric reanalysis and initialized by +# temperature, salinity, sea ice concentration, and sea ice thickness from the ECCO state estimate. +# It includes a few closures: +# - "Gent-McWilliams" `IsopycnalSkewSymmetricDiffusivity`, +# - `CATKEVerticalDiffusivity` for vertical convective mixing, +# - `HorizontalScalarBiharmonicDiffusivity` to damp grid scale noise, +# - and `VerticalScalarDiffusivity` emulating mixing from internal tides +# +# For this example, we need Oceananigans, NumericalEarth, Dates, CUDA, and +# CairoMakie to visualize the simulation. + +using NumericalEarth +using Oceananigans +using Oceananigans.Units +using Dates +using Printf +using Statistics +using CUDA + +# ### Grid and Bathymetry +arch = GPU() +Nz = 30 +z = ExponentialDiscretization(Nz, -5500, 0; scale = 1240) +grid = ORCATripolarGrid(arch; dataset = ORCA2(), z, Nz, remove_closed_basins = true, halo = (5, 5, 4)) + +# ### Closures +# +# We include a Gent-McWilliams isopycnal diffusivity as a parameterization for the mesoscale +# eddy fluxes. For vertical mixing at the upper-ocean boundary layer we include the CATKE +# parameterization. + +@inline νhb(i, j, k, grid, λ) = Oceananigans.Operators.Az(i, j, k, grid, Center(), Center(), Center())^2 / λ +νh = CenterField(grid) +set!(νh, KernelFunctionOperation{Center, Center, Center}(νhb, grid, 40days)) +horizontal_viscosity = HorizontalScalarBiharmonicDiffusivity(ν=νh) + +@inline henyey_diffusivity(x, y, z) = max(2e-6, 3e-5 * abs(sind(y))) +ν_henyey = CenterField(grid) +set!(ν_henyey, henyey_diffusivity) +vertical_diffusivity = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(); κ=ν_henyey, ν=3e-5) + +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity +eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=500, κ_symmetric=500) + +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: + CATKEVerticalDiffusivity + +vertical_mixing = CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization(); + maximum_tracer_diffusivity = 1, + maximum_tke_diffusivity = 1, + maximum_viscosity = 1) + +free_surface = SplitExplicitFreeSurface(grid; substeps=60) +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=5) + +@inline restoring_mask(x, y, z, t) = z>-50 +rate = 1/50days +dates = DateTime(1993, 1, 1) : Month(1) : DateTime(1993, 11, 1) +salinity = Metadata(:salinity; dates, dataset=ECCO4DarwinMonthly(), dir = data_path) +FS = DatasetRestoring(salinity, grid; mask = restoring_mask, rate) + +ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, + closure=(eddy_closure, vertical_mixing, horizontal_viscosity, vertical_diffusivity), + forcing = (; S = FS)) + +sea_ice = sea_ice_simulation(grid, ocean; dynamics = nothing) + +# ### Initial condition + +# We initialize the ocean and sea ice models with data from the ECCO state estimate. + +date = DateTime(1993, 1, 1) +dataset = ECCO4Monthly() +ecco_temperature = Metadatum(:temperature; date, dataset) +ecco_salinity = Metadatum(:salinity; date, dataset) +ecco_sea_ice_thickness = Metadatum(:sea_ice_thickness; date, dataset) +ecco_sea_ice_concentration = Metadatum(:sea_ice_concentration; date, dataset) + +set!(ocean.model, T=ecco_temperature, S=ecco_salinity) +set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) + +# ### Atmospheric forcing + +# We force the simulation with a JRA55-do atmospheric reanalysis. +radiation = Radiation(arch) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(), + include_rivers_and_icebergs = true) + +# ### Coupled simulation + +# Now we are ready to build the coupled ocean--sea ice model and bring everything +# together into a `simulation`. + +# With Runge-Kutta 3rd order time-stepping we can safely use a timestep of 90 minutes. + +coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +simulation = Simulation(coupled_model; Δt=90minutes, stop_time=2.5*365days) + +# ### A progress messenger +# +# We write a function that prints out a helpful progress message while the simulation runs. + +wall_time = Ref(time_ns()) + +function progress(sim) + ocean = sim.model.ocean + u, v, w = ocean.model.velocities + T = ocean.model.tracers.T + e = ocean.model.tracers.e + emax = maximum(e) + umax = (maximum(abs, u), maximum(abs, v), maximum(abs, w)) + + step_time = 1e-9 * (time_ns() - wall_time[]) + + msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim)) + msg2 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", umax...) + msg4 = @sprintf(", max(e): %.2f m² s⁻²", emax) + msg5 = @sprintf(", wall time: %s", prettytime(step_time)) + msg6 = @sprintf(", Δt = %s, %.1f SYPD\n", prettytime(sim.Δt), Units.day / (step_time * 365 / 10)) + + @info msg1 * msg2 * msg4 * msg5 * msg6 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(simulation, progress, TimeInterval(10days)) + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities) +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature)) + +u, v, w = ocean.model.velocities + +diagnostics = (drake_transport = Field(Integral(view(u, 108, 20:33, :), dims = (2, 3))), + mean_T = Field(Average(ocean.model.tracers.T)), + mean_S = Field(Average(ocean.model.tracers.S)), + mean_SSH = Field(Average(ocean.model.free_surface.displacement)), + total_TKE = Field(Integral(ocean.model.tracers.e))) + +vertical_slices = + (atlantic = view(ocean.model.tracers.T, 125, :, :), + pacific = view(ocean.model.tracers.T, 66, :, :)) + +fname_suffix = "orca2" + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(365days/12), + filename = "surface_"*fname_suffix, + indices = (:, :, grid.Nz)) + +ocean.output_writers[:diagnostics] = JLD2Writer(ocean.model, diagnostics; + schedule = TimeInterval(365days/12), + filename = "diagnostics_"*fname_suffix) + +ocean.output_writers[:vslices] = JLD2Writer(ocean.model, vertical_slices; + schedule = TimeInterval(365days/12), + filename = "basin_slice_"*fname_suffix) + +sea_ice.output_writers[:surface] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = TimeInterval(365days/12), + filename = "sea_ice_"*fname_suffix) + +# ### Ready to run + +# We are ready to press the big red button and run the simulation. +run!(simulation) + +# ### Load and plot results + +using CairoMakie + +fds_surface = FieldDataset("surface_"*fname_suffix*".jld2", backend = OnDisk()) +fds_ice = FieldDataset("sea_ice_"*fname_suffix*".jld2", backend = OnDisk()) +fds_diagnostics = FieldDataset("diagnostics_"*fname_suffix*".jld2") +fds_slice = FieldDataset("basin_slice_"*fname_suffix*".jld2") + +grid = fds_surface["T"].grid + +land = interior(grid.immersed_boundary.bottom_height) .≥ 0 + +Nt = length(fds_surface["u"]) + +uoₙ = Field{Face, Center, Nothing}(grid) +voₙ = Field{Center, Face, Nothing}(grid) + +so = Field(sqrt(uoₙ^2 + voₙ^2)) + +n = Observable(Nt) + +spd_plt = @lift begin + parent(uoₙ) .= parent(fds_surface["u"][$n]) + parent(voₙ) .= parent(fds_surface["v"][$n]) + compute!(so) + soₙ = interior(so) + soₙ[land] .= NaN + [view(soₙ, :, :, 1)...] +end + +T_plt = @lift begin + T = interior(fds_surface["T"][$n]) + T[land] .= NaN + [view(T, :, :, 1)...] +end + +ice_plt = @lift begin + hₙ = interior(fds_ice["h"][$n]) + ℵₙ = interior(fds_ice["ℵ"][$n]) + he = hₙ .* ℵₙ + he[he .< 1e-7] .= NaN + he[land] .= NaN + [view(he, :, :, 1)...] +end + +ice_T_plt = @lift begin + T = interior(fds_ice["T"][$n]) + h = interior(fds_ice["h"][$n]) + ℵ = interior(fds_ice["ℵ"][$n]) + he = h .* ℵ + T[he .< 1e-7] .= NaN + T[land] .= NaN + [view(T, :, :, 1)...] +end + +times = fds_diagnostics["drake_transport"].times + +title = @lift string(floor(times[$n]./365days))*" years, "*string(floor(Int, mod(times[$n]/(365days/12), 12)))*" months" + +λ, φ = nodes(so) + +Nx, Ny = size(grid) +Δx = grid.Δxᶜᶜᵃ[1:Nx, 1:Ny] +Δy = grid.Δyᶜᶜᵃ[1:Nx, 1:Ny] +markersize = [[(8*Δx[i, j]/maximum(Δx)./cosd(φ[i, j]), 8*Δy[i, j]/maximum(Δy)) for i in 1:Nx, j in 1:Ny]...]; +polar_markersize = [[(12*Δx[i, j]/maximum(Δx), 12*Δy[i, j]/maximum(Δy) .* sind(90-φ[i, j])) for i in 1:Nx, j in 1:Ny]...]; + +fig = Figure(size=(1200, 1200)) + +ax = Axis(fig[1:2, 1:4], limits = (-180, 180, -80, 80), backgroundcolor = :lightgray, xgridvisible = false, ygridvisible = false; title) +ax2 = Axis(fig[3:4, 1:4], limits = (-180, 180, -80, 80), backgroundcolor = :lightgray, xgridvisible = false, ygridvisible = false) + +ax3 = PolarAxis(fig[1, 5:6], rgridvisible = false, thetagridvisible = false) +ax4 = PolarAxis(fig[2, 5:6], rgridvisible = false, thetagridvisible = false) + +ax5 = PolarAxis(fig[3, 5:6], rgridvisible = false, thetagridvisible = false) +ax6 = PolarAxis(fig[4, 5:6], rgridvisible = false, thetagridvisible = false) + +ax7 = Axis(fig[5, 1:3], title = "Atlantic", ylabel = "Depth (m)", xlabel = "Latitude (°)") +ax8 = Axis(fig[5, 4:6], title = "Pacific", ylabel = "Depth (m)", xlabel = "Latitude (°)") + +sco = scatter!(ax, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray; markersize, colorrange = (0, 0.5), marker = :rect) +sci = scatter!(ax, [λ...], [φ...], color=ice_plt, colormap=:greys; markersize, colorrange = (0, 4), marker = :rect) + +scatter!(ax3, π/180 .* [λ...], 90 .-[φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = polar_markersize, colorrange = (0, 0.5), marker = :rect, rotation = π/180 .* [λ...]) +scatter!(ax3, π/180 .* [λ...], 90 .-[φ...], color=ice_plt, colormap=:greys; markersize = polar_markersize, colorrange = (0, 4), marker = :rect, rotation = π/180 .* [λ...]) + +scatter!(ax4, π/180 .* [λ...], [φ...] .+ 90, color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = polar_markersize, colorrange = (0, 0.5), marker = :rect, rotation = π/180 .* [λ...]) +scatter!(ax4, π/180 .* [λ...], [φ...] .+ 90, color=ice_plt, colormap=:greys; markersize = polar_markersize, colorrange = (0, 4), marker = :rect, rotation = π/180 .* [λ...]) + +Colorbar(fig[1:2, 0], sco, label = "Ocean Surface Speed (m/s)") +Colorbar(fig[1:2, 7], sci, label = "Sea Ice Effective Thickness (m)") + +lo, hi = -50, 0 +center = -1 +frac = (center - lo) / (hi - lo) +compressed_roma = cgrad( + :roma, + [0.0, frac, 1.0], + categorical = false, + rev = true +) + +scT = scatter!(ax2, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect, nan_color=:lightgray) +scatter!(ax5, π/180 .* [λ...], 90 .-[φ...], color=T_plt, colormap=:magma; markersize = polar_markersize, colorrange = (-1, 32), marker = :rect, nan_color=:lightgray) +scatter!(ax6, π/180 .* [λ...], [φ...] .+ 90, color=T_plt, colormap=:magma; markersize = polar_markersize, colorrange = (-1, 32), marker = :rect, nan_color=:lightgray) +scTi = scatter!(ax2, [λ...], [φ...], color=ice_T_plt, colormap=compressed_roma; markersize, colorrange = (-50, 0), marker = :rect) +scatter!(ax5, π/180 .* [λ...], 90 .-[φ...], color=ice_T_plt, colormap=compressed_roma; markersize = polar_markersize, colorrange = (-50, 0), marker = :rect) +scatter!(ax6, π/180 .* [λ...], [φ...] .+ 90, color=ice_T_plt, colormap=compressed_roma; markersize = polar_markersize, colorrange = (-50, 0), marker = :rect) + +title = @lift string(floor(times[$n]./365days))*" years, "*string(mod(times[$n]/(365days/12), 12))*" months" +Colorbar(fig[3:4, 0], scT, label = "Sea Surface Temperature (°C)") +Colorbar(fig[3:4, 7], scTi, label = "Ice Surface Temperature (°C)") + +λ, φ, z = nodes(grid, Center(), Center(), Center()) + +heatmap!(ax7, φ[125, :], z, (@lift interior(fds_slice["atlantic"][$n], 1, :, :)), colormap=:magma, colorrange = (-1, 32), nan_color=:lightgray) +heatmap!(ax8, φ[66, :], z, (@lift interior(fds_slice["pacific"][$n], 1, :, :)), colormap=:magma, colorrange = (-1, 32), nan_color=:lightgray) + +rlims!(ax3, 0, 40); hiderdecorations!(ax3) +rlims!(ax4, 0, 40); hiderdecorations!(ax4) +rlims!(ax5, 0, 40); hiderdecorations!(ax5) +rlims!(ax6, 0, 40); hiderdecorations!(ax6) + +record(fig, "orca2.mp4", 1:Nt, framerate = 10) do i + n[] = i +end + +nothing #hide + +# ![](orca2.mp4) + +fig = Figure() + +ax = Axis(fig[1, 1], title = "Drake transport (Sverdrops)") +ax2 = Axis(fig[2, 1], ylabel = "Temperature (°C)", title = "Global Mean") +ax3 = Axis(fig[2, 1], yaxisposition = :right, ylabel = "Salinity (psu)") +ax4 = Axis(fig[3, 1], ylabel = "TKE (m²/s)") +ax5 = Axis(fig[3, 1], yaxisposition = :right, ylabel = "SSH (m)", xlabel = "Year") +hidexdecorations!(ax3) +hidexdecorations!(ax5) + +lines!(ax, times./(365days), interior(fds_diagnostics["drake_transport"], 1, 1, 1, :)./10^6) +ylims!(ax, 120, 180) + +lines!(ax2, times./(365days), interior(fds_diagnostics["mean_T"], 1, 1, 1, :)) +lines!(ax3, times./(365days), interior(fds_diagnostics["mean_S"], 1, 1, 1, :), color = Makie.wong_colors()[2]) +lines!(ax4, times./(365days), interior(fds_diagnostics["total_TKE"], 1, 1, 1, :)) +lines!(ax5, times./(365days), interior(fds_diagnostics["mean_SSH"], 1, 1, 1, :), color = Makie.wong_colors()[2]) + +save("orca2_diagnostics.png", fig) + +nothing #hide + +# ![](orca2_diagnostics.png) diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index c290a7f44..3c2a69169 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -1,6 +1,6 @@ module Bathymetry -export regrid_bathymetry, ORCAGrid +export regrid_bathymetry, ORCATripolarGrid using Downloads using ImageMorphology @@ -23,6 +23,6 @@ using ..DataWrangling: Metadatum, native_grid, metadata_path, download_dataset using ..DataWrangling.ETOPO: ETOPO2022 include("regrid_bathymetry.jl") -include("orca_grid.jl") +include("orca_tripolar_grid.jl") end # module diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_tripolar_grid.jl similarity index 94% rename from src/Bathymetry/orca_grid.jl rename to src/Bathymetry/orca_tripolar_grid.jl index ce5cec635..4f5ae0552 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_tripolar_grid.jl @@ -9,7 +9,7 @@ using CubedSphere.SphericalGeometry: lat_lon_to_cartesian, cartesian_to_lat_lon, using Distances: haversine using ..DataWrangling: dataset_variable_name, default_download_directory -using ..DataWrangling.ORCA: ORCA1, ORCA12, default_south_rows_to_remove +using ..DataWrangling.ORCA: ORCA1, ORCA12, ORCA2, default_south_rows_to_remove """ read_2d_nemo_variable(ds, name) @@ -278,7 +278,7 @@ function read_orca_staggered_mesh(ds; radius = Oceananigans.defaults.planet_radi end # Detect periodic overlap columns in NEMO data. -# The eORCA grid has `n` trailing columns that are copies of the first `n` columns +# The eORCA grid has `n` trailing columns that are copies of the first `n` columns # (e.g., columns 361:362 repeat columns 1:2 for eORCA1). function periodic_overlap_index(λCC) Nx = size(λCC, 1) @@ -346,8 +346,12 @@ function halo_fill_stagger(CC, FC, CF, FF, helper_grid, bcs) ) end +fold_topology(::ORCA12) = RightFaceFolded +fold_topology(::ORCA1) = RightFaceFolded +fold_topology(::ORCA2) = RightCenterFolded + """ - ORCAGrid(arch = CPU(), FT::DataType = Float64; + ORCATripolarGrid(arch = CPU(), FT::DataType = Float64; dataset, halo = (4, 4, 4), z = (-6000, 0), @@ -399,16 +403,17 @@ Keyword Arguments - `dir`: Directory to store and look up ORCA files (`mesh_mask` and bathymetry). Defaults to the dataset scratch cache via `default_download_directory(dataset)`. """ -function ORCAGrid(arch = CPU(), FT::DataType = Float64; - dataset = ORCA1(), - halo = (4, 4, 4), - z = (-6000, 0), - Nz = 50, - radius = Oceananigans.defaults.planet_radius, - with_bathymetry = true, - active_cells_map = true, - south_rows_to_remove = default_south_rows_to_remove(dataset), - dir = default_download_directory(dataset)) +function ORCATripolarGrid(arch = CPU(), FT::DataType = Float64; + dataset = ORCA1(), + halo = (4, 4, 4), + z = (-6000, 0), + Nz = 50, + radius = Oceananigans.defaults.planet_radius, + with_bathymetry = true, + active_cells_map = true, + south_rows_to_remove = default_south_rows_to_remove(dataset), + remove_closed_basins = false, + dir = default_download_directory(dataset)) # Download mesh_mask via the metadata interface mesh_meta = Metadatum(:mesh_mask; dataset, dir) @@ -491,7 +496,7 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; # Build the grid to_arch(data) = on_architecture(arch, map(FT, data)) - underlying_grid = OrthogonalSphericalShellGrid{Periodic, RightFaceFolded, Bounded}( + underlying_grid = OrthogonalSphericalShellGrid{Periodic, fold_topology(dataset), Bounded}( arch, Nx, Ny, Nz, Hx, Hy, Hz, @@ -527,7 +532,10 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; # Land (bathymetry == 0) gets mapped to +100 so GridFittedBottom masks it. bottom_height = FT.(coalesce.(bathy_data, FT(0))) bottom_height .= ifelse.(isfinite.(bottom_height) .& (bottom_height .> 0), .-bottom_height, FT(100)) - bottom_height = on_architecture(arch, bottom_height) + remove_closed_basins && remove_minor_basins!(bottom_height, 1, (underlying_grid.Nx, underlying_grid.Ny)) + + bathymetry = Field{Center, Center, Nothing}(underlying_grid) + set!(bathymetry, bottom_height) - return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map) + return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry); active_cells_map) end diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index 4a7f3f623..93f8b399b 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -1,8 +1,10 @@ module ORCA -export ORCA1, ORCA12 +export ORCA1, ORCA2, ORCA12 using Downloads +using Tar +using CodecZlib using Oceananigans using Oceananigans.DistributedComputations: @root using Scratch @@ -124,4 +126,51 @@ default_south_rows_to_remove(::ORCA1) = 35 # Conservative default: keep all rows unless user requests trimming. default_south_rows_to_remove(::ORCA12) = 0 +# ORCA2 +struct ORCA2 <:ORCADataset end + +const ORCA2Metadatum = Metadatum{<:ORCA2} + +ORCA2_variable_names = Dict( + :bottom_height => "bathy_metry", + :mesh_mask => "glamt", +) + +dataset_variable_name(data::ORCA2Metadatum) = ORCA2_variable_names[data.name] + +# Zenodo record 15705144 +const ORCA2_url = "https://zenodo.org/records/15705144/files/ORCA2L75_domaincfg_forcings.tar.gz" + +metadata_url(::ORCA2Metadatum) = ORCA2_url + +function metadata_filename(::ORCA2, name, date, bounding_box) + if name == :mesh_mask + return "ORCA2L75/ORCA2L75/domain_cfg_orca2l75_nemo5.nc" + elseif name == :bottom_height + return "ORCA2L75/ORCA2L75/domain_cfg_orca2l75_nemo5.nc" + else + error("Unknown ORCA2 variable: $name") + end +end + +z_interfaces(::ORCA2Metadatum) = nothing + +function download_dataset(metadatum::ORCA2Metadatum) + fileurl = metadata_url(metadatum) + filepath = metadata_path(metadatum) + fileroute = metadatum.dir + + @root if !isfile(filepath) + @info "Downloading ORCA2 data: $(metadatum.name) to $(metadatum.dir)..." + Downloads.download(fileurl, fileroute*"ORCA2L75.tar.gz"; progress=download_progress) + open(GzipDecompressorStream, fileroute*"ORCA2L75.tar.gz") do io + Tar.extract(io, fileroute*"/ORCA2L75") + end + end + + return filepath +end + +default_south_rows_to_remove(::ORCA2) = 0 + end # module diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index fc5731492..7f253cc25 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -46,7 +46,7 @@ export EN4Monthly, WOAClimatology, WOAAnnual, WOAMonthly, GLORYSDaily, GLORYSMonthly, GLORYSStatic, - ORCA1, ORCA12, + ORCA1, ORCA2, ORCA12, RepeatYearJRA55, MultiYearJRA55, first_date, last_date, @@ -55,7 +55,7 @@ export LinearlyTaperedPolarMask, DatasetRestoring, ocean_simulation, - ORCAGrid, + ORCATripolarGrid, sea_ice_simulation, atmosphere_simulation, sea_ice_dynamics, diff --git a/test/test_orca_grid.jl b/test/test_orca_grid.jl deleted file mode 100644 index 3bf33d03c..000000000 --- a/test/test_orca_grid.jl +++ /dev/null @@ -1,178 +0,0 @@ -include("runtests_setup.jl") -include("download_utils.jl") - -using NumericalEarth -using NumericalEarth.DataWrangling: download_dataset, metadata_path -using NumericalEarth.DataWrangling.ORCA: default_south_rows_to_remove -using Oceananigans -using Oceananigans.OrthogonalSphericalShellGrids: TripolarGrid -using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid -using NCDatasets -using Statistics -using Test - -# Pre-download ORCA1 mesh_mask and bathymetry through the artifacts fallback so -# subsequent ORCAGrid(...) calls find the files locally even when Zenodo is down. -for name in (:mesh_mask, :bottom_height) - md = Metadatum(name; dataset=ORCA1()) - download_dataset_with_fallback(metadata_path(md); dataset_name="ORCA1 $name") do - download_dataset(md) - end -end - -@testset "ORCA1 Metadatum construction" begin - bathy_meta = Metadatum(:bottom_height; dataset=ORCA1()) - @test bathy_meta.name == :bottom_height - @test bathy_meta.dataset isa ORCA1 - - mesh_meta = Metadatum(:mesh_mask; dataset=ORCA1()) - @test mesh_meta.name == :mesh_mask - @test mesh_meta.dataset isa ORCA1 -end - - -@testset "ORCA12 Metadatum construction" begin - bathy_meta = Metadatum(:bottom_height; dataset=ORCA12()) - @test bathy_meta.name == :bottom_height - @test bathy_meta.dataset isa ORCA12 - - mesh_meta = Metadatum(:mesh_mask; dataset=ORCA12()) - @test mesh_meta.name == :mesh_mask - @test mesh_meta.dataset isa ORCA12 - - @test default_south_rows_to_remove(ORCA12()) == 0 - @test occursin("eORCA12", metadata_path(mesh_meta)) - @test occursin("eORCA12", metadata_path(bathy_meta)) -end - -@testset "ORCAGrid with ORCA1 dataset on $(arch)" for arch in test_architectures - south_rows_to_remove = 43 - grid = ORCAGrid(arch; dataset=ORCA1(), Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove) - @test grid.underlying_grid.Ny == 332 - south_rows_to_remove - - grid = ORCAGrid(arch; dataset=ORCA1(), Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove=0) - - # Default returns ImmersedBoundaryGrid with bathymetry - @test grid isa ImmersedBoundaryGrid - underlying = grid.underlying_grid - @test underlying isa Oceananigans.Grids.OrthogonalSphericalShellGrid - @test underlying isa TripolarGrid - @test underlying.Nx == 362 - @test underlying.Ny == 332 - @test underlying.Nz == 5 - - # Coordinates span near-global domain - @test minimum(underlying.λᶜᶜᵃ.parent) < -179 - @test maximum(underlying.λᶜᶜᵃ.parent) > 179 - @test minimum(underlying.φᶜᶜᵃ.parent) < -80 - @test maximum(underlying.φᶜᶜᵃ.parent) > 80 -end - -@testset "ORCAGrid without bathymetry on $(arch)" for arch in test_architectures - grid = ORCAGrid(arch; dataset=ORCA1(), Nz=5, z=(-5000, 0), halo=(4, 4, 4), - with_bathymetry=false) - - @test grid isa Oceananigans.Grids.OrthogonalSphericalShellGrid - @test grid isa TripolarGrid - @test !(grid isa ImmersedBoundaryGrid) - @test grid.Nx == 362 - @test grid.Ny == 332 - default_south_rows_to_remove(ORCA1()) - @test grid.Nz == 5 -end - -@testset "ORCAGrid with south_rows_to_remove on $(arch)" for arch in test_architectures - Nremove = 40 - grid = ORCAGrid(arch; dataset=ORCA1(), Nz=5, z=(-5000, 0), halo=(4, 4, 4), - south_rows_to_remove=Nremove) - - @test grid isa ImmersedBoundaryGrid - underlying = grid.underlying_grid - @test underlying.Nx == 362 - @test underlying.Ny == 332 - Nremove - @test underlying.Nz == 5 -end - -@testset "ORCAGrid metric consistency" begin - grid = ORCAGrid(CPU(); dataset=ORCA1(), Nz=5, z=(-5000, 0), halo=(4, 4, 4), with_bathymetry=false) - - Nx, Ny = grid.Nx, grid.Ny - Hx, Hy = grid.Hx, grid.Hy - - # No NaNs or Infs in any grid data - for name in (:λᶜᶜᵃ, :λᶠᶜᵃ, :λᶜᶠᵃ, :λᶠᶠᵃ, - :φᶜᶜᵃ, :φᶠᶜᵃ, :φᶜᶠᵃ, :φᶠᶠᵃ, - :Δxᶜᶜᵃ, :Δxᶠᶜᵃ, :Δxᶜᶠᵃ, :Δxᶠᶠᵃ, - :Δyᶜᶜᵃ, :Δyᶠᶜᵃ, :Δyᶜᶠᵃ, :Δyᶠᶠᵃ, - :Azᶜᶜᵃ, :Azᶠᶜᵃ, :Azᶜᶠᵃ, :Azᶠᶠᵃ) - data = getproperty(grid, name) - @test all(isfinite, Oceananigans.on_architecture(CPU(), data)) == true - end - - # Metrics strictly positive over the full interior. Face-y fields on - # RightFaceFolded have Ny+1 interior rows; the fold row Ny+1 must be checked. - LYs = Dict(:Δxᶜᶜᵃ => Center, :Δxᶠᶜᵃ => Center, :Δxᶜᶠᵃ => Face, :Δxᶠᶠᵃ => Face, - :Δyᶜᶜᵃ => Center, :Δyᶠᶜᵃ => Center, :Δyᶜᶠᵃ => Face, :Δyᶠᶠᵃ => Face, - :Azᶜᶜᵃ => Center, :Azᶠᶜᵃ => Center, :Azᶜᶠᵃ => Face, :Azᶠᶠᵃ => Face) - for (name, LY) in LYs - data = getproperty(grid, name) - Njf = Base.length(LY(), Oceananigans.Grids.RightFaceFolded(), Ny) - interior = Oceananigans.on_architecture(CPU(), data)[1:Nx, 1:Njf] - @test all(x -> x > 0, interior) == true - end - - for name in (:Δxᶜᶠᵃ, :Δxᶠᶠᵃ, :Δyᶜᶠᵃ, :Δyᶠᶠᵃ, :Azᶜᶠᵃ, :Azᶠᶠᵃ) - data = Oceananigans.on_architecture(CPU(), getproperty(grid, name)) - @test all(x -> x > 0, data[1:Nx, Ny+1]) - end - - # Face-x longitude is west of Center-x longitude (stagger check) - # At mid-latitudes (away from poles), Face[i] should be ≤ Center[i] in longitude. - # Check a mid-latitude row (away from fold and south boundary). - jmid = Ny ÷ 2 - λF = grid.λᶠᶜᵃ[1:Nx, jmid] - λC = grid.λᶜᶜᵃ[1:Nx, jmid] - # Most Face longitudes should be less than the corresponding Center longitude - # (allowing for wraparound near ±180°). Count how many satisfy Face < Center. - n_west = count(i -> begin - δ = λC[i] - λF[i] - # Handle wraparound: if δ < -180, add 360 - δ < -180 && (δ += 360) - δ > 180 && (δ -= 360) - δ > 0 - end, 1:Nx) - - @test n_west / Nx > 0.95 # vast majority should satisfy this - - # Face-y latitude is south of Center-y latitude (stagger check) - # At interior points, Face[j] should be < Center[j] in latitude. - imid = Nx ÷ 2 - φF = grid.φᶜᶠᵃ[imid, 1:Ny] - φC = grid.φᶜᶜᵃ[imid, 1:Ny] - nsouth = count(j -> φF[j] < φC[j], 1:Ny) - @test nsouth / length(φC) > 0.95 - - # Periodic overlap: first and last unique columns should be consistent - # After filling halos, the periodic halo should smoothly wrap. - # Check that Δx at the periodic boundary has no discontinuity. - jmid = Ny ÷ 2 - Δx = grid.Δxᶜᶜᵃ[:, jmid] - # The relative jump from column Nx to column 1 (via periodic halo) - # should be similar to the jump between adjacent interior columns - interior_variation = maximum(abs, diff(Array(Δx[1:Nx]))) / mean(Δx[1:Nx]) - boundary_jump = abs(Δx[Nx] - Δx[1]) / mean(Δx[1:Nx]) - @test boundary_jump < 10 * interior_variation + 1e-10 -end - -@testset "ORCA1 bathymetry retrieval" begin - bathy_md = Metadatum(:bottom_height; dataset=ORCA1()) - download_dataset(bathy_md) - path = metadata_path(bathy_md) - @test isfile(path) - - ds = Dataset(path) - bathy = ds["Bathymetry"][:, :] - close(ds) - - @test size(bathy) == (362, 332) - @test maximum(bathy) > 5000 # Deep ocean -end diff --git a/test/test_orca_tripolar_grid.jl b/test/test_orca_tripolar_grid.jl new file mode 100644 index 000000000..5a59d12fb --- /dev/null +++ b/test/test_orca_tripolar_grid.jl @@ -0,0 +1,186 @@ +include("runtests_setup.jl") +include("download_utils.jl") + +using NumericalEarth +using NumericalEarth.DataWrangling: download_dataset, metadata_path +using NumericalEarth.DataWrangling.ORCA: default_south_rows_to_remove +using Oceananigans +using Oceananigans.OrthogonalSphericalShellGrids: TripolarGrid +using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid +using NCDatasets +using Statistics +using Test + +orca_datasets = (ORCA1, ORCA2)#, ORCA12) +orca_sizes = ((362, 332), (180, 148))#, (4322, 3606)) +bathymetry_names = ("Bathymetry", "bathy_metry")#, "bottom_height") + +# Pre-download ORCA mesh_mask and bathymetry through the artifacts fallback so +# subsequent ORCAGrid(...) calls find the files locally even when Zenodo is down. +for dataset in orca_datasets + for name in (:mesh_mask, :bottom_height) + md = Metadatum(name; dataset=dataset()) + download_dataset_with_fallback(metadata_path(md); dataset_name="$dataset $name") do + download_dataset(md) + end + end +end + +for (n, DATASET) in enumerate(orca_datasets) + dataset = DATASET() + Nx0, Ny0 = orca_sizes[n] + name = nameof(DATASET) + + @testset "$name Metadatum construction" begin + bathy_meta = Metadatum(:bottom_height; dataset) + @test bathy_meta.name == :bottom_height + @test bathy_meta.dataset isa DATASET + + mesh_meta = Metadatum(:mesh_mask; dataset) + @test mesh_meta.name == :mesh_mask + @test mesh_meta.dataset isa DATASET + end + + @testset "ORCATripolarGrid with $name dataset on $(arch)" for arch in test_architectures + south_rows_to_remove = 15 + grid = ORCATripolarGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove) + @test grid.underlying_grid.Ny == Ny0 - south_rows_to_remove + + grid = ORCATripolarGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove=0) + + # Default returns ImmersedBoundaryGrid with bathymetry + @test grid isa ImmersedBoundaryGrid + underlying = grid.underlying_grid + @test underlying isa Oceananigans.Grids.OrthogonalSphericalShellGrid + @test underlying isa TripolarGrid + @test underlying.Nx == Nx0 + @test underlying.Ny == Ny0 + @test underlying.Nz == 5 + + # Coordinates span near-global domain + @test minimum(underlying.λᶜᶜᵃ.parent) < -179 + @test maximum(underlying.λᶜᶜᵃ.parent) > 179 + @test minimum(underlying.φᶜᶜᵃ.parent) < -70 + @test maximum(underlying.φᶜᶜᵃ.parent) > 87 + end + + @testset "ORCATripolarGrid ($name) without bathymetry on $(arch)" for arch in test_architectures + grid = ORCATripolarGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), + with_bathymetry=false) + + @test grid isa Oceananigans.Grids.OrthogonalSphericalShellGrid + @test grid isa TripolarGrid + @test !(grid isa ImmersedBoundaryGrid) + @test grid.Nx == Nx0 + @test grid.Ny == Ny0 - default_south_rows_to_remove(dataset) + @test grid.Nz == 5 + end + + @testset "ORCATripolarGrid ($name) with south_rows_to_remove on $(arch)" for arch in test_architectures + Nremove = 40 + grid = ORCATripolarGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), + south_rows_to_remove=Nremove) + + @test grid isa ImmersedBoundaryGrid + underlying = grid.underlying_grid + @test underlying.Nx == Nx0 + @test underlying.Ny == Ny0 - Nremove + @test underlying.Nz == 5 + end + + @testset "ORCATripolarGrid ($name) metric consistency" begin + grid = ORCATripolarGrid(CPU(); dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), with_bathymetry=false) + + Nx, Ny = grid.Nx, grid.Ny + Hx, Hy = grid.Hx, grid.Hy + + # No NaNs or Infs in any grid data + for name in (:λᶜᶜᵃ, :λᶠᶜᵃ, :λᶜᶠᵃ, :λᶠᶠᵃ, + :φᶜᶜᵃ, :φᶠᶜᵃ, :φᶜᶠᵃ, :φᶠᶠᵃ, + :Δxᶜᶜᵃ, :Δxᶠᶜᵃ, :Δxᶜᶠᵃ, :Δxᶠᶠᵃ, + :Δyᶜᶜᵃ, :Δyᶠᶜᵃ, :Δyᶜᶠᵃ, :Δyᶠᶠᵃ, + :Azᶜᶜᵃ, :Azᶠᶜᵃ, :Azᶜᶠᵃ, :Azᶠᶠᵃ) + data = getproperty(grid, name) + @test all(isfinite, Oceananigans.on_architecture(CPU(), data)) == true + end + + # All interior metrics (Δx, Δy, Az) are strictly positive + # Check only interior points to avoid halo issues + for name in (:Δxᶜᶜᵃ, :Δxᶠᶜᵃ, :Δxᶜᶠᵃ, :Δxᶠᶠᵃ, + :Δyᶜᶜᵃ, :Δyᶠᶜᵃ, :Δyᶜᶠᵃ, :Δyᶠᶠᵃ, + :Azᶜᶜᵃ, :Azᶠᶜᵃ, :Azᶜᶠᵃ, :Azᶠᶠᵃ) + data = getproperty(grid, name) + interior = Oceananigans.on_architecture(CPU(), data)[1:Nx, 1:Ny] + @test all(isfinite, Oceananigans.on_architecture(CPU(), data)) + end + + # Metrics strictly positive over the full interior. Face-y fields on + # RightFaceFolded have Ny+1 interior rows; the fold row Ny+1 must be checked. + LYs = Dict(:Δxᶜᶜᵃ => Center, :Δxᶠᶜᵃ => Center, :Δxᶜᶠᵃ => Face, :Δxᶠᶠᵃ => Face, + :Δyᶜᶜᵃ => Center, :Δyᶠᶜᵃ => Center, :Δyᶜᶠᵃ => Face, :Δyᶠᶠᵃ => Face, + :Azᶜᶜᵃ => Center, :Azᶠᶜᵃ => Center, :Azᶜᶠᵃ => Face, :Azᶠᶠᵃ => Face) + for (name, LY) in LYs + data = getproperty(grid, name) + Njf = Base.length(LY(), Oceananigans.Grids.RightFaceFolded(), Ny) + interior = Oceananigans.on_architecture(CPU(), data)[1:Nx, 1:Njf] + @test all(x -> x > 0, interior) == true + end + + for name in (:Δxᶜᶠᵃ, :Δxᶠᶠᵃ, :Δyᶜᶠᵃ, :Δyᶠᶠᵃ, :Azᶜᶠᵃ, :Azᶠᶠᵃ) + data = Oceananigans.on_architecture(CPU(), getproperty(grid, name)) + @test all(x -> x > 0, data[1:Nx, Ny+1]) + end + + # Face-x longitude is west of Center-x longitude (stagger check) + # At mid-latitudes (away from poles), Face[i] should be ≤ Center[i] in longitude. + # Check a mid-latitude row (away from fold and south boundary). + jmid = Ny ÷ 2 + λF = grid.λᶠᶜᵃ[1:Nx, jmid] + λC = grid.λᶜᶜᵃ[1:Nx, jmid] + # Most Face longitudes should be less than the corresponding Center longitude + # (allowing for wraparound near ±180°). Count how many satisfy Face < Center. + n_west = count(i -> begin + δ = λC[i] - λF[i] + # Handle wraparound: if δ < -180, add 360 + δ < -180 && (δ += 360) + δ > 180 && (δ -= 360) + δ > 0 + end, 1:Nx) + + @test n_west / Nx > 0.95 # vast majority should satisfy this + + # Face-y latitude is south of Center-y latitude (stagger check) + # At interior points, Face[j] should be < Center[j] in latitude. + imid = Nx ÷ 2 + φF = grid.φᶜᶠᵃ[imid, 1:Ny] + φC = grid.φᶜᶜᵃ[imid, 1:Ny] + nsouth = count(j -> φF[j] < φC[j], 1:Ny) + @test nsouth / length(φC) > 0.95 + + # Periodic overlap: first and last unique columns should be consistent + # After filling halos, the periodic halo should smoothly wrap. + # Check that Δx at the periodic boundary has no discontinuity. + jmid = Ny ÷ 2 + Δx = grid.Δxᶜᶜᵃ[:, jmid] + # The relative jump from column Nx to column 1 (via periodic halo) + # should be similar to the jump between adjacent interior columns + interior_variation = maximum(abs, diff(Array(Δx[1:Nx]))) / mean(Δx[1:Nx]) + boundary_jump = abs(Δx[Nx] - Δx[1]) / mean(Δx[1:Nx]) + @test boundary_jump < 10 * interior_variation + 1e-10 + end + + @testset "$name bathymetry retrieval" begin + bathy_md = Metadatum(:bottom_height; dataset) + download_dataset(bathy_md) + path = metadata_path(bathy_md) + @test isfile(path) + + ds = Dataset(path) + name = + bathy = ds[bathymetry_names[n]][:, :] + close(ds) + + @test size(bathy) == (Nx0, Ny0) + @test maximum(bathy) > 5000 # Deep ocean + end +end \ No newline at end of file