From ad68434e6f346631cd02f923b10062be91d16d61 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 21:07:59 +0100 Subject: [PATCH 01/37] Update ORCA.jl --- src/DataWrangling/ORCA/ORCA.jl | 59 +++++++++++++++++++++++++++++++++- 1 file changed, 58 insertions(+), 1 deletion(-) diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index 90db416e0..4ea32a693 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -1,8 +1,10 @@ module ORCA -export ORCA1 +export ORCA1, ORCA2 using Downloads +using Tar +using CodecZlib using Oceananigans using Oceananigans.DistributedComputations: @root using Scratch @@ -87,4 +89,59 @@ end default_south_rows_to_remove(::ORCA1) = 35 +# ORCA2 +struct ORCA2 end + +default_download_directory(::ORCA2) = download_ORCA_cache +reversed_vertical_axis(::ORCA2) = false +longitude_interfaces(::ORCA2) = (-180, 180) +latitude_interfaces(::ORCA2) = (-80, 90) + +all_dates(::ORCA2, args...) = nothing +first_date(::ORCA2, args...) = nothing +last_date(::ORCA2, args...) = nothing + +const ORCA1Metadatum = Metadatum{<:ORCA2} + +ORCA1_variable_names = Dict( + :bottom_height => "Bathymetry", + :mesh_mask => "glamt", +) + +dataset_variable_name(data::ORCA2Metadatum) = ORCA1_variable_names[data.name] + +# Zenodo record 15705144 +const ORCA2_url = "https://zenodo.org/records/15705144/files/ORCA2L75_domaincfg_forcings.tar.gz" + +metadata_url(metadatum::ORCA2Metadatum) = ORCA2_url + +function metadata_filename(::ORCA2, name, date, bounding_box) + if name == :mesh_mask + return "coordinates.nc" + elseif name == :bottom_height + return "bathy_meter.nc" + else + error("Unknown ORCA2 variable: $name") + end +end + +z_interfaces(::ORCA1Metadatum) = nothing + +function download_dataset(metadatum::ORCA2Metadatum) + fileurl = metadata_url(metadatum) + filepath = metadata_path(metadatum) + + @root if !isfile(filepath) + @info "Downloading ORCA2 data: $(metadatum.name) to $(metadatum.dir)..." + Downloads.download(fileurl, filepath; progress=download_progress) + open(GzipDecompressorStream, tarfile) do io + Tar.extract(io, filepath*"/ORCA2L75") + end + end + + return filepath +end + +default_south_rows_to_remove(::ORCA2) = 35 + end # module From b2b9359ef7d30b0fd61500e8361f548ebabcdc78 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 22:00:40 +0100 Subject: [PATCH 02/37] working --- Project.toml | 4 +++- src/Bathymetry/orca_grid.jl | 2 +- src/DataWrangling/ORCA/ORCA.jl | 21 +++++++++++---------- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index 211dec35f..1878f3779 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,14 @@ name = "NumericalEarth" uuid = "904d977b-046a-4731-8b86-9235c0d1ef02" license = "MIT" -version = "0.2.2" authors = ["NumericalEarth contributors"] +version = "0.2.2" [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" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" @@ -27,6 +28,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/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 7a2ef5274..cb4031da7 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -27,7 +27,7 @@ function read_2d_nemo_variable(ds, name) 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) diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index 4ea32a693..cac938df8 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -101,10 +101,10 @@ all_dates(::ORCA2, args...) = nothing first_date(::ORCA2, args...) = nothing last_date(::ORCA2, args...) = nothing -const ORCA1Metadatum = Metadatum{<:ORCA2} +const ORCA2Metadatum = Metadatum{<:ORCA2} -ORCA1_variable_names = Dict( - :bottom_height => "Bathymetry", +ORCA1_variable_names = Dict( + :bottom_height => "bathy_metry", :mesh_mask => "glamt", ) @@ -113,13 +113,13 @@ dataset_variable_name(data::ORCA2Metadatum) = ORCA1_variable_names[data.name] # Zenodo record 15705144 const ORCA2_url = "https://zenodo.org/records/15705144/files/ORCA2L75_domaincfg_forcings.tar.gz" -metadata_url(metadatum::ORCA2Metadatum) = ORCA2_url +metadata_url(::ORCA2Metadatum) = ORCA2_url function metadata_filename(::ORCA2, name, date, bounding_box) if name == :mesh_mask - return "coordinates.nc" + return "ORCA2L75/ORCA2L75/domain_cfg_orca2l75_nemo5.nc" elseif name == :bottom_height - return "bathy_meter.nc" + return "ORCA2L75/ORCA2L75/domain_cfg_orca2l75_nemo5.nc" else error("Unknown ORCA2 variable: $name") end @@ -130,18 +130,19 @@ z_interfaces(::ORCA1Metadatum) = 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, filepath; progress=download_progress) - open(GzipDecompressorStream, tarfile) do io - Tar.extract(io, filepath*"/ORCA2L75") + 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) = 35 +default_south_rows_to_remove(::ORCA2) = 0 end # module From b3136d8f0257e5fe61c6de1a2a89779eaeff775a Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 22:04:46 +0100 Subject: [PATCH 03/37] maybe now --- src/Bathymetry/orca_grid.jl | 5 +++-- src/DataWrangling/ORCA/ORCA.jl | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index cb4031da7..544268da4 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -355,10 +355,11 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; bathy_data = Array(bathy_ds[bathy_varname][:, :]) close(bathy_ds) - # Chop off the same southern rows from bathymetry + # TODO: reinstate + #= # Chop off the same southern rows from bathymetry if jr > 0 bathy_data = chop_south(bathy_data) - end + end=# # NEMO stores bathymetry as positive depth; convert to negative bottom height # (Oceananigans convention: z < 0 below sea level). diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index cac938df8..81e7d1fac 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -143,6 +143,6 @@ function download_dataset(metadatum::ORCA2Metadatum) return filepath end -default_south_rows_to_remove(::ORCA2) = 0 +default_south_rows_to_remove(::ORCA2) = 1 end # module From 10c51bf220dfe25c2d1689fe25df09ed61eab321 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 22:47:18 +0100 Subject: [PATCH 04/37] maybe now --- src/Bathymetry/orca_grid.jl | 14 ++++++++------ src/DataWrangling/ORCA/ORCA.jl | 4 ++-- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 544268da4..3b9dd6b47 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -6,7 +6,7 @@ using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom using Oceananigans.OrthogonalSphericalShellGrids: Tripolar, continue_south! using ..DataWrangling: dataset_variable_name -using ..DataWrangling.ORCA: ORCA1, default_south_rows_to_remove +using ..DataWrangling.ORCA: ORCA1, ORCA2, default_south_rows_to_remove """ read_2d_nemo_variable(ds, name) @@ -70,6 +70,7 @@ function halo_filled_data(data, helper_grid, bcs, LX, LY, overlap) field.data[1:Ni, 1:Nj, 1] .= shifted_data[1:Ni, 1:Nj] else # Face-y: shift +1 in y field.data[1:Ni, 2:Nj+1, 1] .= shifted_data[1:Ni, 1:Nj] + field.data[1:Ni, 1, 1] .= shifted_data[1:Ni, 1] end fill_halo_regions!(field) @@ -133,7 +134,8 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; radius = Oceananigans.defaults.planet_radius, with_bathymetry = true, active_cells_map = true, - south_rows_to_remove = default_south_rows_to_remove(dataset)) + south_rows_to_remove = default_south_rows_to_remove(dataset), + chop_bathymetry = false) # Validate z specification against Nz (mirrors Oceananigans' input_validation.jl) if z isa AbstractVector @@ -355,11 +357,11 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; bathy_data = Array(bathy_ds[bathy_varname][:, :]) close(bathy_ds) - # TODO: reinstate - #= # Chop off the same southern rows from bathymetry - if jr > 0 + + # Chop off the same southern rows from bathymetry + if chop_bathymetry & (jr > 0) bathy_data = chop_south(bathy_data) - end=# + end # NEMO stores bathymetry as positive depth; convert to negative bottom height # (Oceananigans convention: z < 0 below sea level). diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index 81e7d1fac..6eea3ca54 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -103,12 +103,12 @@ last_date(::ORCA2, args...) = nothing const ORCA2Metadatum = Metadatum{<:ORCA2} -ORCA1_variable_names = Dict( +ORCA2_variable_names = Dict( :bottom_height => "bathy_metry", :mesh_mask => "glamt", ) -dataset_variable_name(data::ORCA2Metadatum) = ORCA1_variable_names[data.name] +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" From 998a1f99d4da9aa2c2d9a4e831376e4ef4746d5c Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 22:51:32 +0100 Subject: [PATCH 05/37] oops --- src/DataWrangling/ORCA/ORCA.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index 6eea3ca54..5f89133d4 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -125,7 +125,7 @@ function metadata_filename(::ORCA2, name, date, bounding_box) end end -z_interfaces(::ORCA1Metadatum) = nothing +z_interfaces(::ORCA2Metadatum) = nothing function download_dataset(metadatum::ORCA2Metadatum) fileurl = metadata_url(metadatum) From 0341c291355081c973db7f3c0d82ef721e9b55a9 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 23:17:57 +0100 Subject: [PATCH 06/37] remove closed basins --- src/Bathymetry/orca_grid.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 3b9dd6b47..f2a0b5d72 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -135,7 +135,8 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; with_bathymetry = true, active_cells_map = true, south_rows_to_remove = default_south_rows_to_remove(dataset), - chop_bathymetry = false) + chop_bathymetry = false, + remove_closed_basins = true) # Validate z specification against Nz (mirrors Oceananigans' input_validation.jl) if z isa AbstractVector @@ -371,5 +372,7 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; bottom_height .= ifelse.(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)) + return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map) end From 7dc3e0b980f86678063fb2bfce210e2958db3818 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 23:27:55 +0100 Subject: [PATCH 07/37] remove closed basins --- src/Bathymetry/orca_grid.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index f2a0b5d72..b4470c227 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -370,9 +370,8 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; # (above sea level) so that GridFittedBottom correctly masks them as land. bottom_height = convert.(FT, bathy_data) bottom_height .= ifelse.(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)) - + bottom_height = on_architecture(arch, bottom_height) + return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map) end From 83bb1ffd661c45b2a9942bf42f8b07b9fc4cbaf3 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 30 Mar 2026 23:32:58 +0100 Subject: [PATCH 08/37] weird order required --- src/Bathymetry/orca_grid.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index b4470c227..9ea483fa7 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -358,11 +358,6 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; bathy_data = Array(bathy_ds[bathy_varname][:, :]) close(bathy_ds) - - # Chop off the same southern rows from bathymetry - if chop_bathymetry & (jr > 0) - bathy_data = chop_south(bathy_data) - end # NEMO stores bathymetry as positive depth; convert to negative bottom height # (Oceananigans convention: z < 0 below sea level). @@ -371,7 +366,13 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; bottom_height = convert.(FT, bathy_data) bottom_height .= ifelse.(bottom_height .> 0, .-bottom_height, FT(100)) remove_closed_basins && remove_minor_basins!(bottom_height, 1, (underlying_grid.Nx, underlying_grid.Ny)) - bottom_height = on_architecture(arch, bottom_height) + + # Chop off the same southern rows from bathymetry + if chop_bathymetry & (jr > 0) + bathy_data = chop_south(bathy_data) + end + bottom_height = on_architecture(arch, bottom_height) + return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map) end From b802a28df03f0a71d6cd9e320dfff30496f59fe9 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 31 Mar 2026 12:39:57 +0100 Subject: [PATCH 09/37] still doesnt work --- src/Bathymetry/orca_grid.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 9ea483fa7..24cd948fa 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -369,10 +369,12 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; # Chop off the same southern rows from bathymetry if chop_bathymetry & (jr > 0) - bathy_data = chop_south(bathy_data) + bottom_height = chop_south(bottom_height) end - - bottom_height = on_architecture(arch, bottom_height) + + bottom_height_f = Field{Center, Center, Nothing}(underlying_grid) + + set!(bottom_height_f, bottom_height) return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map) end From b2edac5100417017e01141400ef2fa446941a0bb Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 2 Apr 2026 20:04:30 +0100 Subject: [PATCH 10/37] Removed random attampts to fix things --- src/Bathymetry/orca_grid.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 24cd948fa..1006127a3 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -70,7 +70,6 @@ function halo_filled_data(data, helper_grid, bcs, LX, LY, overlap) field.data[1:Ni, 1:Nj, 1] .= shifted_data[1:Ni, 1:Nj] else # Face-y: shift +1 in y field.data[1:Ni, 2:Nj+1, 1] .= shifted_data[1:Ni, 1:Nj] - field.data[1:Ni, 1, 1] .= shifted_data[1:Ni, 1] end fill_halo_regions!(field) @@ -134,9 +133,7 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; radius = Oceananigans.defaults.planet_radius, with_bathymetry = true, active_cells_map = true, - south_rows_to_remove = default_south_rows_to_remove(dataset), - chop_bathymetry = false, - remove_closed_basins = true) + south_rows_to_remove = default_south_rows_to_remove(dataset)) # Validate z specification against Nz (mirrors Oceananigans' input_validation.jl) if z isa AbstractVector @@ -368,7 +365,7 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; remove_closed_basins && remove_minor_basins!(bottom_height, 1, (underlying_grid.Nx, underlying_grid.Ny)) # Chop off the same southern rows from bathymetry - if chop_bathymetry & (jr > 0) + if jr > 0 bottom_height = chop_south(bottom_height) end From 35c8dc5e196197ac1c98d17c793064754a1c6e68 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 2 Apr 2026 20:12:30 +0100 Subject: [PATCH 11/37] fix bathymetry --- src/Bathymetry/orca_grid.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 1006127a3..af49015bc 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -133,7 +133,8 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; radius = Oceananigans.defaults.planet_radius, with_bathymetry = true, active_cells_map = true, - south_rows_to_remove = default_south_rows_to_remove(dataset)) + south_rows_to_remove = default_south_rows_to_remove(dataset), + remove_closed_basins = false) # Validate z specification against Nz (mirrors Oceananigans' input_validation.jl) if z isa AbstractVector @@ -360,18 +361,16 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; # (Oceananigans convention: z < 0 below sea level). # In NEMO, bathymetry == 0 means land. We map these to bottom_height = 100 # (above sea level) so that GridFittedBottom correctly masks them as land. + bottom_height = convert.(FT, bathy_data) bottom_height .= ifelse.(bottom_height .> 0, .-bottom_height, FT(100)) remove_closed_basins && remove_minor_basins!(bottom_height, 1, (underlying_grid.Nx, underlying_grid.Ny)) # Chop off the same southern rows from bathymetry - if jr > 0 + if (jr > 0) & !(dataset isa ORCA2) bottom_height = chop_south(bottom_height) end - - bottom_height_f = Field{Center, Center, Nothing}(underlying_grid) - - set!(bottom_height_f, bottom_height) + return ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map) end From b694ff802bd454690e0db1312f9202a1078ee75f Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Thu, 2 Apr 2026 20:59:07 +0100 Subject: [PATCH 12/37] export ORCA2 --- src/NumericalEarth.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 8b23f354a..33f21b1d6 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -38,7 +38,7 @@ export EN4Monthly, WOAClimatology, WOAAnnual, WOAMonthly, GLORYSDaily, GLORYSMonthly, GLORYSStatic, - ORCA1, + ORCA1, ORCA2, RepeatYearJRA55, MultiYearJRA55, first_date, last_date, From d60df12a253b3251efadb52a459d6d45d2bd65d6 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 3 Apr 2026 14:51:30 +0100 Subject: [PATCH 13/37] attempt at validation --- examples/coarse_example.jl | 294 +++++++++++++++++++++++++++++++++++++ 1 file changed, 294 insertions(+) create mode 100644 examples/coarse_example.jl diff --git a/examples/coarse_example.jl b/examples/coarse_example.jl new file mode 100644 index 000000000..0604572a7 --- /dev/null +++ b/examples/coarse_example.jl @@ -0,0 +1,294 @@ +# # [One-degree global ocean--sea ice simulation](@id one-degree-ocean-seaice) +# +# This example configures a global ocean--sea ice simulation at 1ᵒ horizontal resolution with +# realistic bathymetry and a few closures including the "Gent-McWilliams" `IsopycnalSkewSymmetricDiffusivity`. +# The simulation is forced by repeat-year JRA55 atmospheric reanalysis +# and initialized by temperature, salinity, sea ice concentration, and sea ice thickness +# from the ECCO state estimate. +# +# 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 + +# We start by constructing an underlying TripolarGrid at ~1 degree resolution, + +arch = GPU() +Nz = 31 + +depth = 5000meters +z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4, mutable = true) + +grid = ORCAGri(arch, z = z.face, Nz, dataset = ORCA2()) + +# ### 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. + +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, AdvectiveFormulation + +eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, skew_flux_formulation=AdvectiveFormulation()) +vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() + +# ### Ocean simulation +# Now we bring everything together to construct the ocean simulation. +# We use a split-explicit timestepping with 70 substeps for the barotropic mode. + +free_surface = SplitExplicitFreeSurface(grid; substeps=70) +momentum_advection = WENOVectorInvariant(order=5) +tracer_advection = WENO(order=5) + +ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, + closure=(eddy_closure, vertical_mixing)) + +@info "We've built an ocean simulation with model:" +@show ocean.model + +# ### Sea Ice simulation +# +# We also build a sea ice simulation. We use the default configuration: +# EVP rheology and a zero-layer thermodynamic model that advances thickness +# and concentration. + +sea_ice = sea_ice_simulation(grid, ocean; advection=tracer_advection) + +# ### 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(80), + include_rivers_and_icebergs = false) + +# ### 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 20 minutes. + +coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) +simulation = Simulation(coupled_model; Δt=20minutes, stop_time=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 + Tmin, Tmax, Tavg = minimum(T), maximum(T), mean(view(T, :, :, ocean.model.grid.Nz)) + 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...) + msg3 = @sprintf(", extrema(To): (%.1f, %.1f) ᵒC, mean(To(z=0)): %.1f ᵒC", Tmin, Tmax, Tavg) + msg4 = @sprintf(", max(e): %.2f m² s⁻²", emax) + msg5 = @sprintf(", wall time: %s \n", prettytime(step_time)) + + @info msg1 * msg2 * msg3 * msg4 * msg5 + + wall_time[] = time_ns() + + return nothing +end + +# And add it as a callback to the simulation. +add_callback!(simulation, progress, TimeInterval(5days)) + +# ### Output +# +# We are almost there! We need to save some output. Below we choose to save _only surface_ +# fields using the `indices` keyword argument. We save all the velocity and tracer components. +# Note, that besides temperature and salinity, the CATKE vertical mixing parameterization +# also uses a prognostic turbulent kinetic energy, `e`, to diagnose the vertical mixing length. + +ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities) +free_surface = ocean.model.free_surface.displacement +sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, + ℵ = sea_ice.model.ice_concentration, + T = sea_ice.model.ice_thermodynamics.top_surface_temperature), + sea_ice.model.velocities) + +ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; + schedule = TimeInterval(1days), + filename = "ocean_one_degree_surface_fields", + indices = (:, :, grid.Nz), + overwrite_existing = true) + +ocean.output_writers[:free_surface] = JLD2Writer(ocean.model, (; η = free_surface); + schedule = TimeInterval(1days), + filename = "ocean_one_degree_free_surface", + overwrite_existing = true) + +sea_ice.output_writers[:surface] = JLD2Writer(sea_ice.model, sea_ice_outputs; + schedule = TimeInterval(1days), + filename = "sea_ice_one_degree_surface_fields", + overwrite_existing = true) + +# ### Ready to run + +# We are ready to press the big red button and run the simulation. +run!(simulation) + +# ### A movie +# +# We load the saved output and make a movie of the simulation. First we plot a snapshot: +using CairoMakie + +# We suffix the ocean fields with "o": +uo = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "u"; backend = OnDisk()) +vo = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "v"; backend = OnDisk()) +To = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "T"; backend = OnDisk()) +eo = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "e"; backend = OnDisk()) +ηo = FieldTimeSeries("ocean_one_degree_free_surface.jld2", "η"; backend = OnDisk()) + +# and sea ice fields with "i": +ui = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "u"; backend = OnDisk()) +vi = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "v"; backend = OnDisk()) +hi = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "h"; backend = OnDisk()) +ℵi = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "ℵ"; backend = OnDisk()) +Ti = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "T"; backend = OnDisk()) + +times = uo.times +Nt = length(times) +n = Observable(Nt) + +# We create a land mask and use it to fill land points with `NaN`s. +land = interior(To.grid.immersed_boundary.bottom_height) .≥ 0 + +Toₙ = @lift begin + Tₙ = interior(To[$n]) + Tₙ[land] .= NaN + view(Tₙ, :, :, 1) +end + +eoₙ = @lift begin + eₙ = interior(eo[$n]) + eₙ[land] .= NaN + view(eₙ, :, :, 1) +end + +ηoₙ = @lift begin + ηₙ = interior(ηo[$n]) + ηₙ[land] .= NaN + view(ηₙ, :, :, 1) +end + +heₙ = @lift begin + hₙ = interior(hi[$n]) + ℵₙ = interior(ℵi[$n]) + hₙ[land] .= NaN + view(hₙ, :, :, 1) .* view(ℵₙ, :, :, 1) +end + +# We compute the surface speeds for the ocean and the sea ice. +uoₙ = Field{Face, Center, Nothing}(uo.grid) +voₙ = Field{Center, Face, Nothing}(vo.grid) + +uiₙ = Field{Face, Center, Nothing}(ui.grid) +viₙ = Field{Center, Face, Nothing}(vi.grid) + +so = Field(sqrt(uoₙ^2 + voₙ^2)) +si = Field(sqrt(uiₙ^2 + viₙ^2)) + +soₙ = @lift begin + parent(uoₙ) .= parent(uo[$n]) + parent(voₙ) .= parent(vo[$n]) + compute!(so) + soₙ = interior(so) + soₙ[land] .= NaN + view(soₙ, :, :, 1) +end + +siₙ = @lift begin + parent(uiₙ) .= parent(ui[$n]) + parent(viₙ) .= parent(vi[$n]) + compute!(si) + siₙ = interior(si) + hₙ = interior(hi[$n]) + ℵₙ = interior(ℵi[$n]) + he = hₙ .* ℵₙ + siₙ[he .< 1e-7] .= 0 + siₙ[land] .= NaN + view(siₙ, :, :, 1) +end + +# Finally, we plot a snapshot of the surface speed, temperature, and the turbulent +# eddy kinetic energy from the CATKE vertical mixing parameterization as well as the +# sea ice speed and the effective sea ice thickness. +fig = Figure(size=(1200, 1000)) + +title = @lift string("Global 1ᵒ ocean simulation after ", prettytime(times[$n] - times[1])) + +axso = Axis(fig[1, 1]) +axηo = Axis(fig[1, 3]) +axTo = Axis(fig[2, 1]) +axeo = Axis(fig[2, 3]) +axsi = Axis(fig[3, 1]) +axhi = Axis(fig[3, 3]) + +hm = heatmap!(axso, soₙ, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) +Colorbar(fig[1, 2], hm, label = "Ocean Surface speed (m s⁻¹)") +hm = heatmap!(axηo, ηoₙ, colorrange = (-1.2, 1.2), colormap = :balance, nan_color=:lightgray) +Colorbar(fig[1, 4], hm, label = "Sea Surface Height (m)") +hm = heatmap!(axTo, Toₙ, colorrange = (-1, 32), colormap = :magma, nan_color=:lightgray) +Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)") +hm = heatmap!(axeo, eoₙ, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray) +Colorbar(fig[2, 4], hm, label = "Turbulent Kinetic Energy (m² s⁻²)") + +hm = heatmap!(axsi, siₙ, colorrange = (0, 0.5), colormap = :greys, nan_color=:lightgray) +Colorbar(fig[1, 4], hm, label = "Sea ice speed (m s⁻¹)") +hm = heatmap!(axhi, heₙ, colorrange = (0, 4), colormap = :blues, nan_color=:lightgray) +Colorbar(fig[2, 4], hm, label = "Effective ice thickness (m)") + + +for ax in (axso, axsi, axTo, axhi, axeo) + hidedecorations!(ax) +end + +Label(fig[0, :], title) + +save("global_snapshot.png", fig) +nothing #hide + +# ![](global_snapshot.png) + +# And now a movie: + +CairoMakie.record(fig, "one_degree_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn + n[] = nn +end +nothing #hide + +# ![](one_degree_global_ocean_surface.mp4) From 51705f415854fc324677466a704a058c2dd606a8 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 3 Apr 2026 15:13:11 +0100 Subject: [PATCH 14/37] Wrap bathymetry in field to make it work on GPU --- src/Bathymetry/orca_grid.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 09fe929f9..84df1400a 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -276,6 +276,8 @@ function ORCAGrid(arch = CPU(), FT::DataType = Float64; bottom_height .= ifelse.(bottom_height .> 0, .-bottom_height, FT(100)) 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 From 7376478a639542a764278a2451253f9efbb62637 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 3 Apr 2026 20:54:14 +0100 Subject: [PATCH 15/37] added tests --- test/test_orca_grid.jl | 271 +++++++++++++++++++++-------------------- 1 file changed, 141 insertions(+), 130 deletions(-) diff --git a/test/test_orca_grid.jl b/test/test_orca_grid.jl index 823e7f125..828854e62 100644 --- a/test/test_orca_grid.jl +++ b/test/test_orca_grid.jl @@ -10,137 +10,148 @@ using NCDatasets using Statistics using Test -@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 "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 +orca_datasets = (ORCA1, ORCA2) +orca_sizes = ((362, 332), (180, 148)) +bathymetry_names = ("Bathymetry", "bathy_metry") + +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 - # 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(x -> x > 0, interior) == true + @testset "ORCAGrid with $name dataset on $(arch)" for arch in test_architectures + south_rows_to_remove = 15 + grid = ORCAGrid(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 = ORCAGrid(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) > 80 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 + @testset "ORCAGrid ($name) without bathymetry on $(arch)" for arch in test_architectures + grid = ORCAGrid(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 "ORCAGrid ($name) with south_rows_to_remove on $(arch)" for arch in test_architectures + Nremove = 40 + grid = ORCAGrid(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 "ORCAGrid ($name) metric consistency" begin + grid = ORCAGrid(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(x -> x > 0, interior) == true + 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 From 676d3e111b2c24cb77fc851591e9ca5c5ce633b8 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 3 Apr 2026 20:54:42 +0100 Subject: [PATCH 16/37] changed halos --- examples/coarse_example.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/coarse_example.jl b/examples/coarse_example.jl index 0604572a7..8c3ed6727 100644 --- a/examples/coarse_example.jl +++ b/examples/coarse_example.jl @@ -27,7 +27,7 @@ Nz = 31 depth = 5000meters z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4, mutable = true) -grid = ORCAGri(arch, z = z.face, Nz, dataset = ORCA2()) +grid = ORCAGrid(arch; z, Nz, dataset = ORCA2(), halo = (5, 5, 4)) # ### Closures # From 61cf02c40c0edb604aca944604ffec1d34602052 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 10 Apr 2026 14:55:35 +0100 Subject: [PATCH 17/37] change name to ORCATripolarGrid --- examples/coarse_example.jl | 22 +++++++++++++--------- src/Bathymetry/Bathymetry.jl | 2 +- src/Bathymetry/orca_grid.jl | 9 ++++++--- src/NumericalEarth.jl | 2 +- test/test_orca_grid.jl | 18 +++++++++--------- 5 files changed, 30 insertions(+), 23 deletions(-) diff --git a/examples/coarse_example.jl b/examples/coarse_example.jl index 8c3ed6727..8b7625a41 100644 --- a/examples/coarse_example.jl +++ b/examples/coarse_example.jl @@ -21,13 +21,14 @@ using CUDA # We start by constructing an underlying TripolarGrid at ~1 degree resolution, +data_path = "/cephfs/home/js2430/store/Global/data" + arch = GPU() Nz = 31 depth = 5000meters z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4, mutable = true) - -grid = ORCAGrid(arch; z, Nz, dataset = ORCA2(), halo = (5, 5, 4)) +grid = ORCATripolarGrid(arch; dataset = ORCA2(), south_rows_to_remove=0) # ### Closures # @@ -39,6 +40,7 @@ using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, Advect eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, skew_flux_formulation=AdvectiveFormulation()) vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() +background_vertical_mixing = VerticalScalarDiffusivity(ν = 3e-5) # ### Ocean simulation # Now we bring everything together to construct the ocean simulation. @@ -49,7 +51,7 @@ momentum_advection = WENOVectorInvariant(order=5) tracer_advection = WENO(order=5) ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, - closure=(eddy_closure, vertical_mixing)) + closure=(eddy_closure, vertical_mixing, background_vertical_mixing)) @info "We've built an ocean simulation with model:" @show ocean.model @@ -68,10 +70,10 @@ sea_ice = sea_ice_simulation(grid, ocean; advection=tracer_advection) 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) +ecco_temperature = Metadatum(:temperature; date, dataset, dir = data_path) +ecco_salinity = Metadatum(:salinity; date, dataset, dir = data_path) +ecco_sea_ice_thickness = Metadatum(:sea_ice_thickness; date, dataset, dir = data_path) +ecco_sea_ice_concentration = Metadatum(:sea_ice_concentration; date, dataset, dir = data_path) set!(ocean.model, T=ecco_temperature, S=ecco_salinity) set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) @@ -81,7 +83,8 @@ set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) # We force the simulation with a JRA55-do atmospheric reanalysis. radiation = Radiation(arch) atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80), - include_rivers_and_icebergs = false) + include_rivers_and_icebergs = false, + dir = data_path) # ### Coupled simulation @@ -157,7 +160,7 @@ sea_ice.output_writers[:surface] = JLD2Writer(sea_ice.model, sea_ice_outputs; overwrite_existing = true) # ### Ready to run - +#= # We are ready to press the big red button and run the simulation. run!(simulation) @@ -292,3 +295,4 @@ end nothing #hide # ![](one_degree_global_ocean_surface.mp4) +=# \ No newline at end of file diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index 40eed858d..40fec7766 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 diff --git a/src/Bathymetry/orca_grid.jl b/src/Bathymetry/orca_grid.jl index 84df1400a..6b2070325 100644 --- a/src/Bathymetry/orca_grid.jl +++ b/src/Bathymetry/orca_grid.jl @@ -85,8 +85,11 @@ function halo_fill_stagger(CC, FC, CF, FF, helper_grid, bcs, overlap) ) end +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), @@ -134,7 +137,7 @@ Keyword Arguments contains degenerate padding rows near Antarctica that are entirely land. Removing them reduces memory usage and computation. """ -function ORCAGrid(arch = CPU(), FT::DataType = Float64; +function ORCATripolarGrid(arch = CPU(), FT::DataType = Float64; dataset = ORCA1(), halo = (4, 4, 4), z = (-6000, 0), @@ -242,7 +245,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, diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 33f21b1d6..453c49418 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -47,7 +47,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 index 828854e62..e015883a1 100644 --- a/test/test_orca_grid.jl +++ b/test/test_orca_grid.jl @@ -29,12 +29,12 @@ for (n, DATASET) in enumerate(orca_datasets) @test mesh_meta.dataset isa DATASET end - @testset "ORCAGrid with $name dataset on $(arch)" for arch in test_architectures + @testset "ORCATripolarGrid with $name dataset on $(arch)" for arch in test_architectures south_rows_to_remove = 15 - grid = ORCAGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove) + 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 = ORCAGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove=0) + 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 @@ -52,8 +52,8 @@ for (n, DATASET) in enumerate(orca_datasets) @test maximum(underlying.φᶜᶜᵃ.parent) > 80 end - @testset "ORCAGrid ($name) without bathymetry on $(arch)" for arch in test_architectures - grid = ORCAGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), + @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 @@ -64,9 +64,9 @@ for (n, DATASET) in enumerate(orca_datasets) @test grid.Nz == 5 end - @testset "ORCAGrid ($name) with south_rows_to_remove on $(arch)" for arch in test_architectures + @testset "ORCATripolarGrid ($name) with south_rows_to_remove on $(arch)" for arch in test_architectures Nremove = 40 - grid = ORCAGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), + grid = ORCATripolarGrid(arch; dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), south_rows_to_remove=Nremove) @test grid isa ImmersedBoundaryGrid @@ -76,8 +76,8 @@ for (n, DATASET) in enumerate(orca_datasets) @test underlying.Nz == 5 end - @testset "ORCAGrid ($name) metric consistency" begin - grid = ORCAGrid(CPU(); dataset, Nz=5, z=(-5000, 0), halo=(4, 4, 4), with_bathymetry=false) + @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 From d84062ec7a975d37973a3f7fdd68185e03e1b8fd Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 10 Apr 2026 14:55:49 +0100 Subject: [PATCH 18/37] correct south chopping - now seems to run --- src/DataWrangling/ORCA/ORCA.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index 5f89133d4..eaec494ac 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -143,6 +143,6 @@ function download_dataset(metadatum::ORCA2Metadatum) return filepath end -default_south_rows_to_remove(::ORCA2) = 1 +default_south_rows_to_remove(::ORCA2) = 0 end # module From cef02b7a6cf80400e7f4e938958171e5ab910602 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 10 Apr 2026 22:02:15 +0100 Subject: [PATCH 19/37] change test file name --- test/{test_orca_grid.jl => test_orca_tripolar_grid.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{test_orca_grid.jl => test_orca_tripolar_grid.jl} (100%) diff --git a/test/test_orca_grid.jl b/test/test_orca_tripolar_grid.jl similarity index 100% rename from test/test_orca_grid.jl rename to test/test_orca_tripolar_grid.jl From dbee27df9579573d0e12ee63c9c68fd386e09410 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 13 Apr 2026 12:05:57 +0100 Subject: [PATCH 20/37] rename file --- src/Bathymetry/Bathymetry.jl | 2 +- src/Bathymetry/{orca_grid.jl => orca_tripolar_grid.jl} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename src/Bathymetry/{orca_grid.jl => orca_tripolar_grid.jl} (100%) diff --git a/src/Bathymetry/Bathymetry.jl b/src/Bathymetry/Bathymetry.jl index 40fec7766..8e53a088e 100644 --- a/src/Bathymetry/Bathymetry.jl +++ b/src/Bathymetry/Bathymetry.jl @@ -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 100% rename from src/Bathymetry/orca_grid.jl rename to src/Bathymetry/orca_tripolar_grid.jl From 8a1a66dbefdae087bc670a854a8421e9e1a66c8a Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 13 Apr 2026 12:08:31 +0100 Subject: [PATCH 21/37] Update test/test_orca_tripolar_grid.jl Co-authored-by: Simone Silvestri --- test/test_orca_tripolar_grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_orca_tripolar_grid.jl b/test/test_orca_tripolar_grid.jl index e015883a1..c88289a88 100644 --- a/test/test_orca_tripolar_grid.jl +++ b/test/test_orca_tripolar_grid.jl @@ -49,7 +49,7 @@ for (n, DATASET) in enumerate(orca_datasets) @test minimum(underlying.λᶜᶜᵃ.parent) < -179 @test maximum(underlying.λᶜᶜᵃ.parent) > 179 @test minimum(underlying.φᶜᶜᵃ.parent) < -70 - @test maximum(underlying.φᶜᶜᵃ.parent) > 80 + @test maximum(underlying.φᶜᶜᵃ.parent) > 87 end @testset "ORCATripolarGrid ($name) without bathymetry on $(arch)" for arch in test_architectures From b108268c3459f7673fa04140bee337e4d1bdde5a Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 14 Apr 2026 15:51:53 +0100 Subject: [PATCH 22/37] added coarse example --- examples/coarse_example.jl | 330 ++++++++++++++++++------------------- 1 file changed, 164 insertions(+), 166 deletions(-) diff --git a/examples/coarse_example.jl b/examples/coarse_example.jl index 8b7625a41..769d8ac53 100644 --- a/examples/coarse_example.jl +++ b/examples/coarse_example.jl @@ -4,7 +4,7 @@ # realistic bathymetry and a few closures including the "Gent-McWilliams" `IsopycnalSkewSymmetricDiffusivity`. # The simulation is forced by repeat-year JRA55 atmospheric reanalysis # and initialized by temperature, salinity, sea ice concentration, and sea ice thickness -# from the ECCO state estimate. +# from the ECCO state estimate. # # For this example, we need Oceananigans, NumericalEarth, Dates, CUDA, and # CairoMakie to visualize the simulation. @@ -21,14 +21,12 @@ using CUDA # We start by constructing an underlying TripolarGrid at ~1 degree resolution, -data_path = "/cephfs/home/js2430/store/Global/data" +const data_path = has_cuda_gpu() ? "/cephfs/home/js2430/store/Global/data" : "data" #"/home/js2430/rds/hpc-work/GlobalOceanBioME/data"# arch = GPU() -Nz = 31 - -depth = 5000meters -z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4, mutable = true) -grid = ORCATripolarGrid(arch; dataset = ORCA2(), south_rows_to_remove=0) +Nz = 30 +z = ExponentialDiscretization(Nz, -5500, 0; scale = 1240) +grid = ORCATripolarGrid(arch; dataset = ORCA2(), z, Nz, remove_closed_basins = true)# ORCATripolarGrid # ### Closures # @@ -36,33 +34,38 @@ grid = ORCATripolarGrid(arch; dataset = ORCA2(), south_rows_to_remove=0) # eddy fluxes. For vertical mixing at the upper-ocean boundary layer we include the CATKE # parameterization. -using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, AdvectiveFormulation +@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) -eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, skew_flux_formulation=AdvectiveFormulation()) -vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() -background_vertical_mixing = VerticalScalarDiffusivity(ν = 3e-5) +@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) -# ### Ocean simulation -# Now we bring everything together to construct the ocean simulation. -# We use a split-explicit timestepping with 70 substeps for the barotropic mode. +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity +eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=500, κ_symmetric=500) -free_surface = SplitExplicitFreeSurface(grid; substeps=70) -momentum_advection = WENOVectorInvariant(order=5) -tracer_advection = WENO(order=5) +vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() -ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, - closure=(eddy_closure, vertical_mixing, background_vertical_mixing)) +free_surface = SplitExplicitFreeSurface(grid; substeps=20*3) +momentum_advection = VectorInvariant() +tracer_advection = UpwindBiased(order = 3) -@info "We've built an ocean simulation with model:" -@show ocean.model +dates = DateTime(1993, 1, 1) : Month(1) : DateTime(1993, 11, 1) +mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(-100, 0)) +salinity = Metadata(:salinity; dates, dataset=ECCO4DarwinMonthly(), dir = data_path) +rate = 1/6days +FS = DatasetRestoring(salinity, grid; mask, rate) -# ### Sea Ice simulation -# -# We also build a sea ice simulation. We use the default configuration: -# EVP rheology and a zero-layer thermodynamic model that advances thickness -# and concentration. +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; advection=tracer_advection) +sea_ice = sea_ice_simulation(grid, ocean; + advection=WENO(order=3, minimum_buffer_upwind_order=1), + dynamics = nothing) # ### Initial condition @@ -82,8 +85,8 @@ set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) # We force the simulation with a JRA55-do atmospheric reanalysis. radiation = Radiation(arch) -atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80), - include_rivers_and_icebergs = false, +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(2920), + include_rivers_and_icebergs = true, dir = data_path) # ### Coupled simulation @@ -94,7 +97,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80), # With Runge-Kutta 3rd order time-stepping we can safely use a timestep of 20 minutes. coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=20minutes, stop_time=365days) +simulation = Simulation(coupled_model; Δt=60minutes, stop_time=10*365days) # ### A progress messenger # @@ -107,7 +110,6 @@ function progress(sim) u, v, w = ocean.model.velocities T = ocean.model.tracers.T e = ocean.model.tracers.e - Tmin, Tmax, Tavg = minimum(T), maximum(T), mean(view(T, :, :, ocean.model.grid.Nz)) emax = maximum(e) umax = (maximum(abs, u), maximum(abs, v), maximum(abs, w)) @@ -115,184 +117,180 @@ function progress(sim) msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim)) msg2 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", umax...) - msg3 = @sprintf(", extrema(To): (%.1f, %.1f) ᵒC, mean(To(z=0)): %.1f ᵒC", Tmin, Tmax, Tavg) msg4 = @sprintf(", max(e): %.2f m² s⁻²", emax) - msg5 = @sprintf(", wall time: %s \n", prettytime(step_time)) + 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 * msg3 * msg4 * msg5 + @info msg1 * msg2 * msg4 * msg5 * msg6 wall_time[] = time_ns() - return nothing + return nothing end # And add it as a callback to the simulation. -add_callback!(simulation, progress, TimeInterval(5days)) - -# ### Output -# -# We are almost there! We need to save some output. Below we choose to save _only surface_ -# fields using the `indices` keyword argument. We save all the velocity and tracer components. -# Note, that besides temperature and salinity, the CATKE vertical mixing parameterization -# also uses a prognostic turbulent kinetic energy, `e`, to diagnose the vertical mixing length. +add_callback!(simulation, progress, TimeInterval(10days)) ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities) -free_surface = ocean.model.free_surface.displacement sea_ice_outputs = merge((h = sea_ice.model.ice_thickness, ℵ = sea_ice.model.ice_concentration, - T = sea_ice.model.ice_thermodynamics.top_surface_temperature), - sea_ice.model.velocities) + 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, :, :), + amoc = view(v, 100:135, 80, :)) + +fname_suffix = "orca2" ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs; - schedule = TimeInterval(1days), - filename = "ocean_one_degree_surface_fields", - indices = (:, :, grid.Nz), - overwrite_existing = true) + 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[:free_surface] = JLD2Writer(ocean.model, (; η = free_surface); - schedule = TimeInterval(1days), - filename = "ocean_one_degree_free_surface", - overwrite_existing = true) +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(1days), - filename = "sea_ice_one_degree_surface_fields", - overwrite_existing = true) + 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) -# ### A movie -# -# We load the saved output and make a movie of the simulation. First we plot a snapshot: -using CairoMakie - -# We suffix the ocean fields with "o": -uo = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "u"; backend = OnDisk()) -vo = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "v"; backend = OnDisk()) -To = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "T"; backend = OnDisk()) -eo = FieldTimeSeries("ocean_one_degree_surface_fields.jld2", "e"; backend = OnDisk()) -ηo = FieldTimeSeries("ocean_one_degree_free_surface.jld2", "η"; backend = OnDisk()) - -# and sea ice fields with "i": -ui = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "u"; backend = OnDisk()) -vi = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "v"; backend = OnDisk()) -hi = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "h"; backend = OnDisk()) -ℵi = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "ℵ"; backend = OnDisk()) -Ti = FieldTimeSeries("sea_ice_one_degree_surface_fields.jld2", "T"; backend = OnDisk()) - -times = uo.times -Nt = length(times) -n = Observable(Nt) - -# We create a land mask and use it to fill land points with `NaN`s. -land = interior(To.grid.immersed_boundary.bottom_height) .≥ 0 +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") -Toₙ = @lift begin - Tₙ = interior(To[$n]) - Tₙ[land] .= NaN - view(Tₙ, :, :, 1) -end +grid = fds_surface["T"].grid -eoₙ = @lift begin - eₙ = interior(eo[$n]) - eₙ[land] .= NaN - view(eₙ, :, :, 1) -end +land = interior(grid.immersed_boundary.bottom_height) .≥ 0 -ηoₙ = @lift begin - ηₙ = interior(ηo[$n]) - ηₙ[land] .= NaN - view(ηₙ, :, :, 1) -end +Nt = length(fds_surface["u"]) -heₙ = @lift begin - hₙ = interior(hi[$n]) - ℵₙ = interior(ℵi[$n]) - hₙ[land] .= NaN - view(hₙ, :, :, 1) .* view(ℵₙ, :, :, 1) -end +uoₙ = Field{Face, Center, Nothing}(grid) +voₙ = Field{Center, Face, Nothing}(grid) -# We compute the surface speeds for the ocean and the sea ice. -uoₙ = Field{Face, Center, Nothing}(uo.grid) -voₙ = Field{Center, Face, Nothing}(vo.grid) - -uiₙ = Field{Face, Center, Nothing}(ui.grid) -viₙ = Field{Center, Face, Nothing}(vi.grid) +uiₙ = Field{Face, Center, Nothing}(grid) +viₙ = Field{Center, Face, Nothing}(grid) so = Field(sqrt(uoₙ^2 + voₙ^2)) si = Field(sqrt(uiₙ^2 + viₙ^2)) -soₙ = @lift begin - parent(uoₙ) .= parent(uo[$n]) - parent(voₙ) .= parent(vo[$n]) +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) + [view(soₙ, :, :, 1)...] +end + +T_plt = @lift begin + T = interior(fds_surface["T"][$n]) + T[land] .= NaN + [view(T, :, :, 1)...] end -siₙ = @lift begin - parent(uiₙ) .= parent(ui[$n]) - parent(viₙ) .= parent(vi[$n]) - compute!(si) - siₙ = interior(si) - hₙ = interior(hi[$n]) - ℵₙ = interior(ℵi[$n]) +ice_plt = @lift begin + hₙ = interior(fds_ice["h"][$n]) + ℵₙ = interior(fds_ice["ℵ"][$n]) he = hₙ .* ℵₙ - siₙ[he .< 1e-7] .= 0 - siₙ[land] .= NaN - view(siₙ, :, :, 1) + he[he .< 1e-7] .= NaN + he[land] .= NaN + [view(he, :, :, 1)...] end -# Finally, we plot a snapshot of the surface speed, temperature, and the turbulent -# eddy kinetic energy from the CATKE vertical mixing parameterization as well as the -# sea ice speed and the effective sea ice thickness. -fig = Figure(size=(1200, 1000)) - -title = @lift string("Global 1ᵒ ocean simulation after ", prettytime(times[$n] - times[1])) - -axso = Axis(fig[1, 1]) -axηo = Axis(fig[1, 3]) -axTo = Axis(fig[2, 1]) -axeo = Axis(fig[2, 3]) -axsi = Axis(fig[3, 1]) -axhi = Axis(fig[3, 3]) - -hm = heatmap!(axso, soₙ, colorrange = (0, 0.5), colormap = :deep, nan_color=:lightgray) -Colorbar(fig[1, 2], hm, label = "Ocean Surface speed (m s⁻¹)") -hm = heatmap!(axηo, ηoₙ, colorrange = (-1.2, 1.2), colormap = :balance, nan_color=:lightgray) -Colorbar(fig[1, 4], hm, label = "Sea Surface Height (m)") -hm = heatmap!(axTo, Toₙ, colorrange = (-1, 32), colormap = :magma, nan_color=:lightgray) -Colorbar(fig[2, 2], hm, label = "Surface Temperature (ᵒC)") -hm = heatmap!(axeo, eoₙ, colorrange = (0, 1e-3), colormap = :solar, nan_color=:lightgray) -Colorbar(fig[2, 4], hm, label = "Turbulent Kinetic Energy (m² s⁻²)") - -hm = heatmap!(axsi, siₙ, colorrange = (0, 0.5), colormap = :greys, nan_color=:lightgray) -Colorbar(fig[1, 4], hm, label = "Sea ice speed (m s⁻¹)") -hm = heatmap!(axhi, heₙ, colorrange = (0, 4), colormap = :blues, nan_color=:lightgray) -Colorbar(fig[2, 4], hm, label = "Effective ice thickness (m)") - - -for ax in (axso, axsi, axTo, axhi, axeo) - hidedecorations!(ax) +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 -Label(fig[0, :], title) +times = fds_diagnostics["drake_transport"].times -save("global_snapshot.png", fig) -nothing #hide +λ, φ = nodes(so) -# ![](global_snapshot.png) +fig = Figure(size=(1200, 1200)) -# And now a movie: +ax = GeoAxis(fig[1:2, 1:2]; dest="+proj=denoy", limits = (-180, 180, -80, 80))#Axis(fig[1, 1]) +ax2 = GeoAxis(fig[1:2, 3:4]; dest="+proj=denoy", limits = (-180, 180, -80, 80))#Axis(fig[1, 3]) -CairoMakie.record(fig, "one_degree_global_ocean_surface.mp4", 1:Nt, framerate = 8) do nn - n[] = nn -end -nothing #hide +ax3 = GeoAxis(fig[3, 1]; dest="+proj=ortho +lat_0=90", limits = (-180, 180, 40, 90)) +ax4 = GeoAxis(fig[3, 2]; dest="+proj=ortho +lat_0=-90", limits = (-180, 180, -90, -40)) + +ax5 = GeoAxis(fig[3, 3]; dest="+proj=ortho +lat_0=90", limits = (-180, 180, 40, 90)) +ax6 = GeoAxis(fig[3, 4]; dest="+proj=ortho +lat_0=-90", limits = (-180, 180, -90, -40)) + +ax7 = Axis(fig[4, 1:2], title = "Atlantic", ylabel = "Depth (m)", xlabel = "Latitude (°)") +ax8 = Axis(fig[4, 3:4], title = "Pacific", ylabel = "Depth (m)", xlabel = "Latitude (°)") + +ax9 = Axis(fig[6, 1:2], title = "Drake transport (Sverdrops)") +ax10 = Axis(fig[6, 3:4], ylabel = "Temperature (°C)", title = "Global Mean", xlabel = "Year") +ax11 = Axis(fig[6, 3:4], yaxisposition = :right, ylabel = "Salinity (psu)", xlabel = "Year") +hidexdecorations!(ax11) + +sco = scatter!(ax, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) +sci = scatter!(ax, [λ...], [φ...], color=ice_plt, colormap=:greys, markersize = 3, colorrange = (0, 4), marker = :rect) + +scatter!(ax3, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) +scatter!(ax3, [λ...], [φ...], color=ice_plt, colormap=:greys, markersize = 3, colorrange = (0, 4), marker = :rect) + +scatter!(ax4, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) +scatter!(ax4, [λ...], [φ...], color=ice_plt, colormap=:greys, markersize = 3, colorrange = (0, 4), marker = :rect) + +Colorbar(fig[0, 1:2], sco, label = "Ocean Surface Speed (m/s)", vertical = false) +Colorbar(fig[5, 1:2], sci, label = "Sea Ice Effective Thickness (m)", vertical = false) + +scT = scatter!(ax2, [λ...], [φ...], color=T_plt, colormap=:magma, markersize = 3, colorrange = (-1, 32), marker = :rect) +scatter!(ax5, [λ...], [φ...], color=T_plt, colormap=:magma, markersize = 3, colorrange = (-1, 32), marker = :rect) +scatter!(ax6, [λ...], [φ...], color=T_plt, colormap=:magma, markersize = 3, colorrange = (-1, 32), marker = :rect) +scTi = scatter!(ax2, [λ...], [φ...], color=ice_T_plt, colormap=:vik, markersize = 3, colorrange = (-2, 0), marker = :rect) +scatter!(ax5, [λ...], [φ...], color=ice_T_plt, colormap=:vik, markersize = 3, colorrange = (-2, 0), marker = :rect) +scatter!(ax6, [λ...], [φ...], color=ice_T_plt, colormap=:vik, markersize = 3, colorrange = (-2, 0), marker = :rect) + +Colorbar(fig[0, 3:4], scT, label = "Sea Surface Temperature (°C)", vertical = false) +Colorbar(fig[5, 3:4], scTi, label = "Ice Surface Temperature (°C)", vertical = false) + +λ, φ, 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) + +lines!(ax9, times./(365days), interior(fds_diagnostics["drake_transport"], 1, 1, 1, :)./10^6) +scatter!(ax9, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["drake_transport"], 1, 1, 1, $n)[]/10^6])) +ylims!(ax9, 120, 180) + +lines!(ax10, times./(365days), interior(fds_diagnostics["mean_T"], 1, 1, 1, :)) +scatter!(ax10, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_T"], 1, 1, 1, $n)[]])) +lines!(ax11, times./(365days), interior(fds_diagnostics["mean_S"], 1, 1, 1, :)) +scatter!(ax11, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_S"], 1, 1, 1, $n)[]])) + +title = @lift string(floor(times[$n]./365days))*" years, "*string(mod(times[$n]/(365days/12), 12))*" months" +Label(fig[-1, :], title, fontsize=18) -# ![](one_degree_global_ocean_surface.mp4) -=# \ No newline at end of file +record(fig, "orca2.mp4", 1:Nt, framerate = 10) do i + n[] = i +end \ No newline at end of file From 23ea14ebe2473b510eb111f8d15077f98dce4096 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 14 Apr 2026 15:54:34 +0100 Subject: [PATCH 23/37] add example to literate --- docs/make.jl | 1 + examples/{coarse_example.jl => orca2.jl} | 0 2 files changed, 1 insertion(+) rename examples/{coarse_example.jl => orca2.jl} (100%) diff --git a/docs/make.jl b/docs/make.jl index f3f332f8d..261177c0c 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", true), 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/coarse_example.jl b/examples/orca2.jl similarity index 100% rename from examples/coarse_example.jl rename to examples/orca2.jl From 838edf209fdc73a7d6f38a4c84862905d852b1db Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 14 Apr 2026 15:55:22 +0100 Subject: [PATCH 24/37] add deps for plotting --- examples/orca2.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/orca2.jl b/examples/orca2.jl index 769d8ac53..51b89574b 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -173,6 +173,8 @@ sea_ice.output_writers[:surface] = JLD2Writer(sea_ice.model, sea_ice_outputs; # We are ready to press the big red button and run the simulation. run!(simulation) +using CairoMakie, GeoMakie + 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") From 820ab5a466ffcdad424de8b6ef0c7d1485688f39 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 14 Apr 2026 15:56:04 +0100 Subject: [PATCH 25/37] tuen off always build --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 261177c0c..1fef163cb 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -32,7 +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", true), + 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), From 2cb0840ac4709031eee110850b99dd889a39a24d Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Wed, 15 Apr 2026 19:21:20 +0100 Subject: [PATCH 26/37] fix diagnostics --- examples/orca2.jl | 106 +++++++++++++++++++++++++++------------------- 1 file changed, 62 insertions(+), 44 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index 51b89574b..f15f7284a 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -51,7 +51,7 @@ vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() free_surface = SplitExplicitFreeSurface(grid; substeps=20*3) momentum_advection = VectorInvariant() -tracer_advection = UpwindBiased(order = 3) +tracer_advection = WENO(order=5) dates = DateTime(1993, 1, 1) : Month(1) : DateTime(1993, 11, 1) mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(-100, 0)) @@ -63,9 +63,7 @@ ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surfac closure=(eddy_closure, vertical_mixing, horizontal_viscosity, vertical_diffusivity), forcing = (; S = FS)) -sea_ice = sea_ice_simulation(grid, ocean; - advection=WENO(order=3, minimum_buffer_upwind_order=1), - dynamics = nothing) +sea_ice = sea_ice_simulation(grid, ocean; dynamics = nothing) # ### Initial condition @@ -146,8 +144,7 @@ diagnostics = (drake_transport = Field(Integral(view(u, 108, 20:33, :), dims = ( vertical_slices = (atlantic = view(ocean.model.tracers.T, 125, :, :), - pacific = view(ocean.model.tracers.T, 66, :, :), - amoc = view(v, 100:135, 80, :)) + pacific = view(ocean.model.tracers.T, 66, :, :)) fname_suffix = "orca2" @@ -233,66 +230,87 @@ end times = fds_diagnostics["drake_transport"].times +title = @lift string(floor(times[$n]./365days))*" years, "*string(mod(times[$n]/(365days/12), 12))*" months" + λ, φ = nodes(so) fig = Figure(size=(1200, 1200)) -ax = GeoAxis(fig[1:2, 1:2]; dest="+proj=denoy", limits = (-180, 180, -80, 80))#Axis(fig[1, 1]) -ax2 = GeoAxis(fig[1:2, 3:4]; dest="+proj=denoy", limits = (-180, 180, -80, 80))#Axis(fig[1, 3]) +ax = GeoAxis(fig[1:2, 1:4]; dest="+proj=moll", limits = (-180, 180, -80, 80), title) +ax2 = GeoAxis(fig[3:4, 1:4]; dest="+proj=moll", limits = (-180, 180, -80, 80)) -ax3 = GeoAxis(fig[3, 1]; dest="+proj=ortho +lat_0=90", limits = (-180, 180, 40, 90)) -ax4 = GeoAxis(fig[3, 2]; dest="+proj=ortho +lat_0=-90", limits = (-180, 180, -90, -40)) +ax3 = GeoAxis(fig[1, 5:6]; dest="+proj=stere +lat_0=90", limits = (-180, 180, 40, 90)) +ax4 = GeoAxis(fig[2, 5:6]; dest="+proj=stere +lat_0=-90", limits = (-180, 180, -90, -40)) -ax5 = GeoAxis(fig[3, 3]; dest="+proj=ortho +lat_0=90", limits = (-180, 180, 40, 90)) -ax6 = GeoAxis(fig[3, 4]; dest="+proj=ortho +lat_0=-90", limits = (-180, 180, -90, -40)) +ax5 = GeoAxis(fig[3, 5:6]; dest="+proj=stere +lat_0=90", limits = (-180, 180, 40, 90)) +ax6 = GeoAxis(fig[4, 5:6]; dest="+proj=stere +lat_0=-90", limits = (-180, 180, -90, -40)) -ax7 = Axis(fig[4, 1:2], title = "Atlantic", ylabel = "Depth (m)", xlabel = "Latitude (°)") -ax8 = Axis(fig[4, 3:4], title = "Pacific", ylabel = "Depth (m)", xlabel = "Latitude (°)") +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 (°)") -ax9 = Axis(fig[6, 1:2], title = "Drake transport (Sverdrops)") -ax10 = Axis(fig[6, 3:4], ylabel = "Temperature (°C)", title = "Global Mean", xlabel = "Year") -ax11 = Axis(fig[6, 3:4], yaxisposition = :right, ylabel = "Salinity (psu)", xlabel = "Year") -hidexdecorations!(ax11) +markersize = [[(8*Δx[i, j]/maximum(Δx), 5*Δy[i, j]/maximum(Δy)) for i in 1:Nx, j in 1:Ny]...]; -sco = scatter!(ax, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) -sci = scatter!(ax, [λ...], [φ...], color=ice_plt, colormap=:greys, markersize = 3, colorrange = (0, 4), marker = :rect) +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, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) -scatter!(ax3, [λ...], [φ...], color=ice_plt, colormap=:greys, markersize = 3, colorrange = (0, 4), marker = :rect) +scatter!(ax3, [λ...], [φ...], color=ice_plt, colormap=:greys; markersize, colorrange = (0, 4), marker = :rect) scatter!(ax4, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) -scatter!(ax4, [λ...], [φ...], color=ice_plt, colormap=:greys, markersize = 3, colorrange = (0, 4), marker = :rect) +scatter!(ax4, [λ...], [φ...], color=ice_plt, colormap=:greys; markersize, colorrange = (0, 4), marker = :rect) -Colorbar(fig[0, 1:2], sco, label = "Ocean Surface Speed (m/s)", vertical = false) -Colorbar(fig[5, 1:2], sci, label = "Sea Ice Effective Thickness (m)", vertical = false) +Colorbar(fig[1:2, 0], sco, label = "Ocean Surface Speed (m/s)") +Colorbar(fig[1:2, 7], sci, label = "Sea Ice Effective Thickness (m)") -scT = scatter!(ax2, [λ...], [φ...], color=T_plt, colormap=:magma, markersize = 3, colorrange = (-1, 32), marker = :rect) -scatter!(ax5, [λ...], [φ...], color=T_plt, colormap=:magma, markersize = 3, colorrange = (-1, 32), marker = :rect) -scatter!(ax6, [λ...], [φ...], color=T_plt, colormap=:magma, markersize = 3, colorrange = (-1, 32), marker = :rect) -scTi = scatter!(ax2, [λ...], [φ...], color=ice_T_plt, colormap=:vik, markersize = 3, colorrange = (-2, 0), marker = :rect) -scatter!(ax5, [λ...], [φ...], color=ice_T_plt, colormap=:vik, markersize = 3, colorrange = (-2, 0), marker = :rect) -scatter!(ax6, [λ...], [φ...], color=ice_T_plt, colormap=:vik, markersize = 3, colorrange = (-2, 0), marker = :rect) +scT = scatter!(ax2, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect) +scatter!(ax5, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect) +scatter!(ax6, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect) +scTi = scatter!(ax2, [λ...], [φ...], color=ice_T_plt, colormap=:vik; markersize, colorrange = (-2, 0), marker = :rect) +scatter!(ax5, [λ...], [φ...], color=ice_T_plt, colormap=:vik; markersize, colorrange = (-2, 0), marker = :rect) +scatter!(ax6, [λ...], [φ...], color=ice_T_plt, colormap=:vik; markersize, colorrange = (-2, 0), marker = :rect) -Colorbar(fig[0, 3:4], scT, label = "Sea Surface Temperature (°C)", vertical = false) -Colorbar(fig[5, 3:4], scTi, label = "Ice Surface Temperature (°C)", vertical = false) +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) -lines!(ax9, times./(365days), interior(fds_diagnostics["drake_transport"], 1, 1, 1, :)./10^6) -scatter!(ax9, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["drake_transport"], 1, 1, 1, $n)[]/10^6])) -ylims!(ax9, 120, 180) +#=record(fig, "orca2.mp4", 1:Nt, framerate = 10) do i + n[] = i +end=# +save("orca2_snapshot.png", fig) +nothing #hide -lines!(ax10, times./(365days), interior(fds_diagnostics["mean_T"], 1, 1, 1, :)) -scatter!(ax10, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_T"], 1, 1, 1, $n)[]])) -lines!(ax11, times./(365days), interior(fds_diagnostics["mean_S"], 1, 1, 1, :)) -scatter!(ax11, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_S"], 1, 1, 1, $n)[]])) +# ![](orca2_snapshot.png) -title = @lift string(floor(times[$n]./365days))*" years, "*string(mod(times[$n]/(365days/12), 12))*" months" -Label(fig[-1, :], title, fontsize=18) +fig = Figure() -record(fig, "orca2.mp4", 1:Nt, framerate = 10) do i - n[] = i -end \ No newline at end of file +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) +scatter!(ax, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["drake_transport"], 1, 1, 1, $n)[]/10^6])) +ylims!(ax, 120, 180) + +lines!(ax2, times./(365days), interior(fds_diagnostics["mean_T"], 1, 1, 1, :)) +scatter!(ax2, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_T"], 1, 1, 1, $n)[]])) +lines!(ax3, times./(365days), interior(fds_diagnostics["mean_S"], 1, 1, 1, :), color = Makie.wong_colors()[2]) +scatter!(ax3, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_S"], 1, 1, 1, $n)[]]), color = Makie.wong_colors()[2]) + +lines!(ax4, times./(365days), interior(fds_diagnostics["total_TKE"], 1, 1, 1, :)) +scatter!(ax4, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["total_TKE"], 1, 1, 1, $n)[]])) +lines!(ax5, times./(365days), interior(fds_diagnostics["mean_SSH"], 1, 1, 1, :), color = Makie.wong_colors()[2]) +scatter!(ax5, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_SSH"], 1, 1, 1, $n)[]]), color = Makie.wong_colors()[2]) + +save("orca2_diagnostics.png", fig) +nothing #hide + +# ![](orca2_diagnostics.png) \ No newline at end of file From f946955043e6672ff9283f37d25e6e6848049b3b Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 17 Apr 2026 17:47:24 +0100 Subject: [PATCH 27/37] remove GeoMakie dependency --- examples/orca2.jl | 101 +++++++++++++++++++++++++++------------------- 1 file changed, 60 insertions(+), 41 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index f15f7284a..09e876486 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -1,10 +1,21 @@ -# # [One-degree global ocean--sea ice simulation](@id one-degree-ocean-seaice) +# # [Coarse global ocean--sea ice simulation](@id coarse-degree-ocean-seaice) # -# This example configures a global ocean--sea ice simulation at 1ᵒ horizontal resolution with -# realistic bathymetry and a few closures including the "Gent-McWilliams" `IsopycnalSkewSymmetricDiffusivity`. -# The simulation is forced by repeat-year JRA55 atmospheric reanalysis -# and initialized by temperature, salinity, sea ice concentration, and sea ice thickness -# from the ECCO state estimate. +# 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. @@ -18,15 +29,12 @@ using Statistics using CUDA # ### Grid and Bathymetry - -# We start by constructing an underlying TripolarGrid at ~1 degree resolution, - const data_path = has_cuda_gpu() ? "/cephfs/home/js2430/store/Global/data" : "data" #"/home/js2430/rds/hpc-work/GlobalOceanBioME/data"# arch = GPU() Nz = 30 z = ExponentialDiscretization(Nz, -5500, 0; scale = 1240) -grid = ORCATripolarGrid(arch; dataset = ORCA2(), z, Nz, remove_closed_basins = true)# ORCATripolarGrid +grid = ORCATripolarGrid(arch; dataset = ORCA2(), z, Nz, remove_closed_basins = true, halo = (5, 5, 4))# ORCATripolarGrid # ### Closures # @@ -47,10 +55,16 @@ vertical_diffusivity = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretiz using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=500, κ_symmetric=500) -vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() +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=20*3) -momentum_advection = VectorInvariant() +momentum_advection = WENOVectorInvariant(order=5) tracer_advection = WENO(order=5) dates = DateTime(1993, 1, 1) : Month(1) : DateTime(1993, 11, 1) @@ -95,7 +109,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(2920), # With Runge-Kutta 3rd order time-stepping we can safely use a timestep of 20 minutes. coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=60minutes, stop_time=10*365days) +simulation = Simulation(coupled_model; Δt=90minutes, stop_time=10*365days) # ### A progress messenger # @@ -170,7 +184,7 @@ sea_ice.output_writers[:surface] = JLD2Writer(sea_ice.model, sea_ice_outputs; # We are ready to press the big red button and run the simulation. run!(simulation) -using CairoMakie, GeoMakie +using CairoMakie fds_surface = FieldDataset("surface_"*fname_suffix*".jld2", backend = OnDisk()) fds_ice = FieldDataset("sea_ice_"*fname_suffix*".jld2", backend = OnDisk()) @@ -186,11 +200,7 @@ Nt = length(fds_surface["u"]) uoₙ = Field{Face, Center, Nothing}(grid) voₙ = Field{Center, Face, Nothing}(grid) -uiₙ = Field{Face, Center, Nothing}(grid) -viₙ = Field{Center, Face, Nothing}(grid) - so = Field(sqrt(uoₙ^2 + voₙ^2)) -si = Field(sqrt(uiₙ^2 + viₙ^2)) n = Observable(Nt) @@ -230,44 +240,48 @@ end times = fds_diagnostics["drake_transport"].times -title = @lift string(floor(times[$n]./365days))*" years, "*string(mod(times[$n]/(365days/12), 12))*" months" +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 = GeoAxis(fig[1:2, 1:4]; dest="+proj=moll", limits = (-180, 180, -80, 80), title) -ax2 = GeoAxis(fig[3:4, 1:4]; dest="+proj=moll", limits = (-180, 180, -80, 80)) +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 = GeoAxis(fig[1, 5:6]; dest="+proj=stere +lat_0=90", limits = (-180, 180, 40, 90)) -ax4 = GeoAxis(fig[2, 5:6]; dest="+proj=stere +lat_0=-90", limits = (-180, 180, -90, -40)) +ax3 = PolarAxis(fig[1, 5:6], rgridvisible = false, thetagridvisible = false) +ax4 = PolarAxis(fig[2, 5:6], rgridvisible = false, thetagridvisible = false) -ax5 = GeoAxis(fig[3, 5:6]; dest="+proj=stere +lat_0=90", limits = (-180, 180, 40, 90)) -ax6 = GeoAxis(fig[4, 5:6]; dest="+proj=stere +lat_0=-90", limits = (-180, 180, -90, -40)) +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 (°)") -markersize = [[(8*Δx[i, j]/maximum(Δx), 5*Δy[i, j]/maximum(Δy)) for i in 1:Nx, j in 1:Ny]...]; - 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, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) -scatter!(ax3, [λ...], [φ...], 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, [λ...], [φ...], color=spd_plt, colormap=:deep, nan_color=:lightgray, markersize = 3, colorrange = (0, 0.5), marker = :rect) -scatter!(ax4, [λ...], [φ...], color=ice_plt, colormap=:greys; markersize, colorrange = (0, 4), marker = :rect) +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)") -scT = scatter!(ax2, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect) -scatter!(ax5, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect) -scatter!(ax6, [λ...], [φ...], color=T_plt, colormap=:magma; markersize, colorrange = (-1, 32), marker = :rect) -scTi = scatter!(ax2, [λ...], [φ...], color=ice_T_plt, colormap=:vik; markersize, colorrange = (-2, 0), marker = :rect) -scatter!(ax5, [λ...], [φ...], color=ice_T_plt, colormap=:vik; markersize, colorrange = (-2, 0), marker = :rect) -scatter!(ax6, [λ...], [φ...], color=ice_T_plt, colormap=:vik; markersize, colorrange = (-2, 0), marker = :rect) +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=:roma; markersize, colorrange = (-2, 0), marker = :rect) +scatter!(ax5, π/180 .* [λ...], 90 .-[φ...], color=ice_T_plt, colormap=:roma; markersize = polar_markersize, colorrange = (-2, 0), marker = :rect) +scatter!(ax6, π/180 .* [λ...], [φ...] .+ 90, color=ice_T_plt, colormap=:roma; markersize = polar_markersize, colorrange = (-2, 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)") @@ -278,13 +292,18 @@ Colorbar(fig[3:4, 7], scTi, label = "Ice Surface Temperature (°C)") 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) -#=record(fig, "orca2.mp4", 1:Nt, framerate = 10) do i +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=# -save("orca2_snapshot.png", fig) +end + nothing #hide -# ![](orca2_snapshot.png) +# ![](orca2.mp4) fig = Figure() From fdb5b5d9487f0d0102bd2f1ef2a745e2988e4e5c Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 17 Apr 2026 18:00:10 +0100 Subject: [PATCH 28/37] remove excess methods --- src/DataWrangling/ORCA/ORCA.jl | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/src/DataWrangling/ORCA/ORCA.jl b/src/DataWrangling/ORCA/ORCA.jl index eaec494ac..060200498 100644 --- a/src/DataWrangling/ORCA/ORCA.jl +++ b/src/DataWrangling/ORCA/ORCA.jl @@ -31,15 +31,6 @@ end struct ORCA1 end -default_download_directory(::ORCA1) = download_ORCA_cache -reversed_vertical_axis(::ORCA1) = false -longitude_interfaces(::ORCA1) = (-180, 180) -latitude_interfaces(::ORCA1) = (-80, 90) - -all_dates(::ORCA1, args...) = nothing -first_date(::ORCA1, args...) = nothing -last_date(::ORCA1, args...) = nothing - const ORCA1Metadatum = Metadatum{<:ORCA1} ORCA1_variable_names = Dict( @@ -92,15 +83,6 @@ default_south_rows_to_remove(::ORCA1) = 35 # ORCA2 struct ORCA2 end -default_download_directory(::ORCA2) = download_ORCA_cache -reversed_vertical_axis(::ORCA2) = false -longitude_interfaces(::ORCA2) = (-180, 180) -latitude_interfaces(::ORCA2) = (-80, 90) - -all_dates(::ORCA2, args...) = nothing -first_date(::ORCA2, args...) = nothing -last_date(::ORCA2, args...) = nothing - const ORCA2Metadatum = Metadatum{<:ORCA2} ORCA2_variable_names = Dict( @@ -145,4 +127,15 @@ end default_south_rows_to_remove(::ORCA2) = 0 +const ORCAs = Union{ORCA1, ORCA2} + +default_download_directory(::ORCAs) = download_ORCA_cache +reversed_vertical_axis(::ORCAs) = false +longitude_interfaces(::ORCAs) = (-180, 180) +latitude_interfaces(::ORCAs) = (-80, 90) + +all_dates(::ORCAs, args...) = nothing +first_date(::ORCAs, args...) = nothing +last_date(::ORCAs, args...) = nothing + end # module From 643e8f5f13e91c94a304c25ba83953ff4f8e2d6c Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 17 Apr 2026 18:42:39 +0100 Subject: [PATCH 29/37] accdentally committed the wrong file --- examples/orca2.jl | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index 09e876486..3d00d07e4 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -29,12 +29,10 @@ using Statistics using CUDA # ### Grid and Bathymetry -const data_path = has_cuda_gpu() ? "/cephfs/home/js2430/store/Global/data" : "data" #"/home/js2430/rds/hpc-work/GlobalOceanBioME/data"# - 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))# ORCATripolarGrid +grid = ORCATripolarGrid(arch; dataset = ORCA2(), z, Nz, remove_closed_basins = true, halo = (5, 5, 4)) # ### Closures # @@ -69,7 +67,7 @@ tracer_advection = WENO(order=5) dates = DateTime(1993, 1, 1) : Month(1) : DateTime(1993, 11, 1) mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(-100, 0)) -salinity = Metadata(:salinity; dates, dataset=ECCO4DarwinMonthly(), dir = data_path) +salinity = Metadata(:salinity; dates, dataset=ECCO4DarwinMonthly()) rate = 1/6days FS = DatasetRestoring(salinity, grid; mask, rate) @@ -85,10 +83,10 @@ sea_ice = sea_ice_simulation(grid, ocean; dynamics = nothing) date = DateTime(1993, 1, 1) dataset = ECCO4Monthly() -ecco_temperature = Metadatum(:temperature; date, dataset, dir = data_path) -ecco_salinity = Metadatum(:salinity; date, dataset, dir = data_path) -ecco_sea_ice_thickness = Metadatum(:sea_ice_thickness; date, dataset, dir = data_path) -ecco_sea_ice_concentration = Metadatum(:sea_ice_concentration; date, dataset, dir = data_path) +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) @@ -97,9 +95,8 @@ set!(sea_ice.model, h=ecco_sea_ice_thickness, ℵ=ecco_sea_ice_concentration) # We force the simulation with a JRA55-do atmospheric reanalysis. radiation = Radiation(arch) -atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(2920), - include_rivers_and_icebergs = true, - dir = data_path) +atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(), + include_rivers_and_icebergs = true) # ### Coupled simulation @@ -109,7 +106,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(2920), # With Runge-Kutta 3rd order time-stepping we can safely use a timestep of 20 minutes. coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=90minutes, stop_time=10*365days) +simulation = Simulation(coupled_model; Δt=90minutes, stop_time=5*365days) # ### A progress messenger # From 67f5180ce418c8c81d4b700f17bd0110277bc281 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sat, 18 Apr 2026 12:40:35 +0100 Subject: [PATCH 30/37] typo --- examples/orca2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index 3d00d07e4..823ff9c0f 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -249,7 +249,7 @@ polar_markersize = [[(12*Δx[i, j]/maximum(Δx), 12*Δy[i, j]/maximum(Δy) .* si fig = Figure(size=(1200, 1200)) -ax = Axis(fig[1:2, 1:4], limits = (-180, 180, -80, 80), backgroundcolor = :lightgray, xgridvisible = false, ygridvisible = false, title) +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) From 4c760a1f28e2e24e25885279e80e2db9b680c707 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Sun, 19 Apr 2026 14:57:12 +0100 Subject: [PATCH 31/37] fix diagnostics --- examples/orca2.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index 823ff9c0f..a63b43790 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -313,20 +313,15 @@ hidexdecorations!(ax3) hidexdecorations!(ax5) lines!(ax, times./(365days), interior(fds_diagnostics["drake_transport"], 1, 1, 1, :)./10^6) -scatter!(ax, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["drake_transport"], 1, 1, 1, $n)[]/10^6])) ylims!(ax, 120, 180) lines!(ax2, times./(365days), interior(fds_diagnostics["mean_T"], 1, 1, 1, :)) -scatter!(ax2, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_T"], 1, 1, 1, $n)[]])) lines!(ax3, times./(365days), interior(fds_diagnostics["mean_S"], 1, 1, 1, :), color = Makie.wong_colors()[2]) -scatter!(ax3, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_S"], 1, 1, 1, $n)[]]), color = Makie.wong_colors()[2]) - lines!(ax4, times./(365days), interior(fds_diagnostics["total_TKE"], 1, 1, 1, :)) -scatter!(ax4, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["total_TKE"], 1, 1, 1, $n)[]])) lines!(ax5, times./(365days), interior(fds_diagnostics["mean_SSH"], 1, 1, 1, :), color = Makie.wong_colors()[2]) -scatter!(ax5, (@lift [times[$n]/(365days)]), (@lift [interior(fds_diagnostics["mean_SSH"], 1, 1, 1, $n)[]]), color = Makie.wong_colors()[2]) save("orca2_diagnostics.png", fig) + nothing #hide # ![](orca2_diagnostics.png) \ No newline at end of file From c88585c9e208adfb27ab5861735ed25ef9483807 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 20 Apr 2026 13:46:48 +0100 Subject: [PATCH 32/37] change colormap --- examples/orca2.jl | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index a63b43790..e44e653f6 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -273,12 +273,23 @@ scatter!(ax4, π/180 .* [λ...], [φ...] .+ 90, color=ice_plt, colormap=:greys; 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=:roma; markersize, colorrange = (-2, 0), marker = :rect) -scatter!(ax5, π/180 .* [λ...], 90 .-[φ...], color=ice_T_plt, colormap=:roma; markersize = polar_markersize, colorrange = (-2, 0), marker = :rect) -scatter!(ax6, π/180 .* [λ...], [φ...] .+ 90, color=ice_T_plt, colormap=:roma; markersize = polar_markersize, colorrange = (-2, 0), marker = :rect) +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)") From ffc81025f7a193fdb6de056ce08a38afbc92b6ce Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 24 Apr 2026 15:53:55 +0100 Subject: [PATCH 33/37] remove ORCA12 from tests --- test/test_orca_tripolar_grid.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_orca_tripolar_grid.jl b/test/test_orca_tripolar_grid.jl index b36f6e6dc..e1ff15c48 100644 --- a/test/test_orca_tripolar_grid.jl +++ b/test/test_orca_tripolar_grid.jl @@ -10,9 +10,9 @@ 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") +orca_datasets = (ORCA1, ORCA2)#, ORCA12) +orca_sizes = ((362, 332), (180, 148))#, (4322, 3606)) +bathymetry_names = ("Bathymetry", "bathy_metry")#, "bottom_height") for (n, DATASET) in enumerate(orca_datasets) dataset = DATASET() From 2a94e873c585c34a6411fea1adf7e347e5b11fd4 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Fri, 24 Apr 2026 20:15:44 +0100 Subject: [PATCH 34/37] Update orca2.jl --- examples/orca2.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index e44e653f6..4ac83795f 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -103,7 +103,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(), # 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 20 minutes. +# 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=5*365days) @@ -181,6 +181,8 @@ sea_ice.output_writers[:surface] = JLD2Writer(sea_ice.model, sea_ice_outputs; # 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()) @@ -273,7 +275,6 @@ scatter!(ax4, π/180 .* [λ...], [φ...] .+ 90, color=ice_plt, colormap=:greys; 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) @@ -335,4 +336,4 @@ save("orca2_diagnostics.png", fig) nothing #hide -# ![](orca2_diagnostics.png) \ No newline at end of file +# ![](orca2_diagnostics.png) From a3f0c779a8f410057d6e2ba5e7b238d65644f330 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 27 Apr 2026 15:17:46 +0100 Subject: [PATCH 35/37] Apply suggestions from code review Co-authored-by: Simone Silvestri --- examples/orca2.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index 4ac83795f..125e87c84 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -61,7 +61,7 @@ vertical_mixing = CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization( maximum_tke_diffusivity = 1, maximum_viscosity = 1) -free_surface = SplitExplicitFreeSurface(grid; substeps=20*3) +free_surface = SplitExplicitFreeSurface(grid; substeps=60) momentum_advection = WENOVectorInvariant(order=5) tracer_advection = WENO(order=5) @@ -106,7 +106,7 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(), # 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=5*365days) +simulation = Simulation(coupled_model; Δt=90minutes, stop_time=2.5*365days) # ### A progress messenger # From d4f5abdb826ec36ad858ff7964dde820b539e552 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Mon, 27 Apr 2026 18:09:56 +0100 Subject: [PATCH 36/37] missed test line --- test/test_orca_tripolar_grid.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_orca_tripolar_grid.jl b/test/test_orca_tripolar_grid.jl index 847de3992..1867b906a 100644 --- a/test/test_orca_tripolar_grid.jl +++ b/test/test_orca_tripolar_grid.jl @@ -1,4 +1,5 @@ include("runtests_setup.jl") +include("download_utils.jl") using NumericalEarth using NumericalEarth.DataWrangling: download_dataset, metadata_path From b588697d0dc07b8866261e912ac510598c536ac0 Mon Sep 17 00:00:00 2001 From: Jago Strong-Wright Date: Tue, 28 Apr 2026 17:06:26 +0100 Subject: [PATCH 37/37] restoring under ice --- examples/orca2.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/orca2.jl b/examples/orca2.jl index 125e87c84..9e79519fb 100644 --- a/examples/orca2.jl +++ b/examples/orca2.jl @@ -65,11 +65,11 @@ 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) -mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(-100, 0)) -salinity = Metadata(:salinity; dates, dataset=ECCO4DarwinMonthly()) -rate = 1/6days -FS = DatasetRestoring(salinity, grid; mask, rate) +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),