Merge branch 'main' into ss/omip-prototype

Author: simone-silvestri

.github/workflows/ci.yml

@@ -39,7 +39,7 @@ jobs:
             arch: aarch64
             version: '1.10'
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@latest
         with:
           version: ${{ matrix.version }}
@@ -76,7 +76,7 @@ jobs:
             arch: aarch64
             version: '1.10'
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@latest
         with:
           version: ${{ matrix.version }}
@@ -112,7 +112,7 @@ jobs:
             arch: aarch64
             version: '1.10'
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@latest
         with:
           version: ${{ matrix.version }}

Project.toml

@@ -2,7 +2,7 @@ name = "ClimaOcean"
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
 license = "MIT"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.8.2"
+version = "0.8.4"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -51,11 +51,11 @@ DataDeps = "0.7"
 DocStringExtensions = "0.9"
 Downloads = "1.6"
 ImageMorphology = "0.4"
-JLD2 = "0.4, 0.5"
+JLD2 = "0.4, 0.5, 0.6"
 KernelAbstractions = "0.9"
 MPI = "0.20"
 NCDatasets = "0.12, 0.13, 0.14"
-Oceananigans = "0.97.6"
+Oceananigans = "0.97.6, 0.98"
 OffsetArrays = "1.14"
 PrecompileTools = "1"
 PythonCall = "0.9"
@@ -65,7 +65,7 @@ SeawaterPolynomials = "0.3.5"
 StaticArrays = "1"
 Statistics = "1.9"
 SurfaceFluxes = "0.11, 0.12"
-Thermodynamics = "0.12, 0.13"
+Thermodynamics = "0.14"
 ZipFile = "0.10"
 julia = "1.10"
 

docs/make.jl

@@ -1,5 +1,6 @@
 using
   ClimaOcean,
+  CUDA,
   Documenter,
   DocumenterCitations,
   Literate
@@ -43,6 +44,8 @@ for file in to_be_literated
     withenv("JULIA_DEBUG" => "Literate") do
         Literate.markdown(filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor(), execute = true)
     end
+    GC.gc()
+    CUDA.reclaim()
 end
 
 #####

examples/one_degree_simulation.jl

@@ -24,11 +24,11 @@ using CUDA
 arch = GPU()
 Nx = 360
 Ny = 180
-Nz = 40
+Nz = 50
+
+depth = 5000meters
+z = ExponentialCoordinate(Nz, -depth, 0; scale = depth/4)
 
-depth = 4000meters
-z = ExponentialCoordinate(Nz, -depth, 0; scale = 0.85*depth)
-z = Oceananigans.Grids.MutableVerticalDiscretization(z)
 underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z)
 
 # Next, we build bathymetry on this grid, using interpolation passes to smooth the bathymetry.
@@ -49,29 +49,29 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height);
 #
 # 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. We also include some explicit horizontal diffusivity.
+# parameterization.
 
-eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=2e3, κ_symmetric=2e3)
-horizontal_viscosity = HorizontalScalarDiffusivity(ν=4000)
+eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=10, κ_symmetric=10)
 vertical_mixing = ClimaOcean.OceanSimulations.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.
+# 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, horizontal_viscosity, vertical_mixing))
+                         timestepper = :SplitRungeKutta3,
+                         closure=(eddy_closure, vertical_mixing))
 
 @info "We've built an ocean simulation with model:"
 @show ocean.model
 
 # ### Sea Ice simulation
-# We also need to build a sea ice simulation. We use the default configuration:
+#
+# 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.
 
@@ -84,8 +84,8 @@ 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_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)
@@ -103,12 +103,10 @@ atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(80),
 # Now we are ready to build the coupled ocean--sea ice model and bring everything
 # together into a `simulation`.
 
-# We use a relatively short time step initially and only run for a few days to
-# avoid numerical instabilities from the initial "shock" of the adjustment of the
-# flow fields.
+# 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=8minutes, stop_time=20days)
+simulation = Simulation(coupled_model; Δt=20minutes, stop_time=365days)
 
 # ### A progress messenger
 #
@@ -141,7 +139,7 @@ function progress(sim)
 end
 
 # And add it as a callback to the simulation.
-add_callback!(simulation, progress, IterationInterval(1000))
+add_callback!(simulation, progress, TimeInterval(5days))
 
 # ### Output
 #
@@ -157,27 +155,19 @@ sea_ice_outputs = merge((h = sea_ice.model.ice_thickness,
                          sea_ice.model.velocities)
 
 ocean.output_writers[:surface] = JLD2Writer(ocean.model, ocean_outputs;
-                                            schedule = TimeInterval(5days),
+                                            schedule = TimeInterval(1days),
                                             filename = "ocean_one_degree_surface_fields",
                                             indices = (:, :, grid.Nz),
                                             overwrite_existing = true)
 
 sea_ice.output_writers[:surface] = JLD2Writer(ocean.model, sea_ice_outputs;
-                                              schedule = TimeInterval(5days),
+                                              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.
-
-# After we run for a short time (here we set up the simulation with `stop_time = 20days`),
-# we increase the timestep and run for longer.
-
-run!(simulation)
-
-simulation.Δt = 30minutes
-simulation.stop_time = 365days
 run!(simulation)
 
 # ### A movie

ext/ClimaOceanPythonCallExt.jl

@@ -7,7 +7,7 @@ using Oceananigans
 using Oceananigans.DistributedComputations: @root
 
 using Dates: DateTime
-using ClimaOcean.DataWrangling.Copernicus: CopernicusMetadata
+using ClimaOcean.DataWrangling.Copernicus: CopernicusMetadata, CopernicusMetadatum
 
 import ClimaOcean.DataWrangling: download_dataset
 
@@ -25,7 +25,17 @@ function install_copernicusmarine()
     return cli
 end
 
-function download_dataset(meta::CopernicusMetadata, grid=nothing; skip_existing = true, additional_kw...)
+# Download each date individually, instead of downloading the entire dataset at once.
+# This is useful for a possible extension of the temporal horizon of the dataset.
+function download_dataset(metadata::CopernicusMetadata; kwargs...)
+    paths = Array{String}(undef, length(metadata))
+    for (m, metadatum) in enumerate(metadata)
+        paths[m] = download_dataset(metadatum; kwargs...)
+    end
+    return paths
+end
+
+function download_dataset(meta::CopernicusMetadatum; skip_existing = true, additional_kw...)
     output_directory = meta.dir
     output_filename = ClimaOcean.DataWrangling.metadata_filename(meta)
     output_path = joinpath(output_directory, output_filename)

src/DataWrangling/Copernicus/Copernicus.jl

@@ -6,7 +6,7 @@ using NCDatasets
 using Printf
 
 using Oceananigans.Fields: Center
-using ClimaOcean.DataWrangling: Metadata, Metadatum, metadata_path
+using ClimaOcean.DataWrangling: Metadata, Metadatum
 using Dates: DateTime, Day, Month
 
 import Oceananigans.Fields:
@@ -45,6 +45,8 @@ dataset_name(::GLORYSStatic) = "GLORYSStatic"
 dataset_name(::GLORYSDaily) = "GLORYSDaily"
 dataset_name(::GLORYSMonthly) = "GLORYSMonthly"
 
+Base.size(::CopernicusDataset, variable) = (4320, 2040, 50)
+
 all_dates(::GLORYSStatic, var) = [nothing]
 all_dates(::GLORYSDaily, var) = range(DateTime("1993-01-01"), stop=DateTime("2021-06-30"), step=Day(1))
 all_dates(::GLORYSMonthly, var) = range(DateTime("1993-01-01"), stop=DateTime("2024-12-01"), step=Month(1))
@@ -57,10 +59,11 @@ copernicusmarine_dataset_id(::GLORYSMonthly) = "cmems_mod_glo_phy_my_0.083deg_P1
 struct CMEMSHourlyAnalysis <: CopernicusDataset end
 copernicusmarine_dataset_id(::CMEMSHourlyAnalysis) = "cmems_mod_glo_phy_anfc_0.083deg_PT1H-m"
 
-CopernicusMetadata{D} = Metadata{<:CopernicusDataset, D}
-CopernicusMetadatum = Metadatum{<:CopernicusDataset}
+const CopernicusMetadata{D} = Metadata{<:CopernicusDataset, D}
+const CopernicusMetadatum = Metadatum{<:CopernicusDataset}
 
 Base.size(::CopernicusMetadatum) = (4320, 2040, 50, 1)
+
 reversed_vertical_axis(::CopernicusDataset) = true
 
 available_variables(::CopernicusDataset) = copernicus_dataset_variable_names
@@ -93,6 +96,8 @@ function bbox_strs(c)
     return first, second
 end
 
+colon2dash(s::String) = replace(s, ":" => "-")
+
 function metadata_prefix(metadata::CopernicusMetadata)
     var = copernicus_dataset_variable_names[metadata.name]
     dataset = dataset_name(metadata.dataset)
@@ -109,7 +114,7 @@ function metadata_prefix(metadata::CopernicusMetadata)
     return string(var, "_",
                   dataset, "_",
                   start_date, "_",
-                  end_date, suffix)
+                  end_date, suffix) |> colon2dash
 end
 
 function metadata_filename(metadata::CopernicusMetadata)

src/DataWrangling/JRA55/JRA55_metadata.jl

@@ -118,7 +118,7 @@ JRA55_dataset_variable_names = Dict(
     :northward_velocity              => "vas",      # Northward near-surface wind
 )
 
-JRA55_multiple_year_url = "https://esgf-data2.llnl.gov/thredds/fileServer/user_pub_work/input4MIPs/CMIP6/OMIP/MRI/MRI-JRA55-do-1-5-0/"
+JRA55_multiple_year_url = "http://esgf-node.ornl.gov/thredds/fileServer/user_pub_work/input4MIPs/CMIP6/OMIP/MRI/MRI-JRA55-do-1-5-0/"
 
 JRA55_multiple_year_prefix = Dict(
     :river_freshwater_flux           => "land/day",

src/OceanSeaIceModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl

@@ -14,7 +14,7 @@ using KernelAbstractions.Extras.LoopInfo: @unroll
 using Statistics: norm
 
 import Thermodynamics as AtmosphericThermodynamics
-import Thermodynamics.Parameters: molmass_ratio
+import Thermodynamics.Parameters: Rv_over_Rd
 
 #####
 ##### Bulk turbulent fluxes based on similarity theory
@@ -264,7 +264,7 @@ L_★ = u_★² / ϰ b_★ .
 @inline function buoyancy_scale(θ★, q★, ℂ, 𝒬, g)
     𝒯ₐ = AtmosphericThermodynamics.virtual_temperature(ℂ, 𝒬)
     qₐ = AtmosphericThermodynamics.vapor_specific_humidity(ℂ, 𝒬)
-    ε  = AtmosphericThermodynamics.Parameters.molmass_ratio(ℂ)
+    ε  = AtmosphericThermodynamics.Parameters.Rv_over_Rd(ℂ)
     δ  = ε - 1 # typically equal to 0.608
 
     b★ = g / 𝒯ₐ * (θ★ * (1 + δ * qₐ) + δ * 𝒯ₐ * q★)