Add a tidal potential to the barotropic forcing

examples/near_global_ocean_simulation.jl

@@ -105,7 +105,10 @@ radiation = Radiation(arch)
 # The number of snapshots that are loaded into memory is determined by
 # the `backend`. Here, we load 41 snapshots at a time into memory.
 
-atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
+tidal_potential = FieldTimeSeries("tidal_potential.jld2", "Φ"; architecture=GPU(), backend=InMemory(41))
+tidal_potential = FieldTimeSeries("tidal_potential.jld2", "Φ"; architecture=GPU(), backend=InMemory(41), boundary_conditions=FieldBoundaryConditions(tidal_potential.grid, (Center, Center, Nothing)))
+Oceananigans.BoundaryConditions.fill_halo_regions!(tidal_potential)
+atmosphere = JRA55PrescribedAtmosphere(arch; tidal_potential, backend=JRA55NetCDFBackend(41))
 
 # ## The coupled simulation
 
@@ -117,7 +120,7 @@ coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
 # We then create a coupled simulation. We start with a small-ish time step of 90 seconds.
 # We run the simulation for 10 days with this small-ish time step.
 
-simulation = Simulation(coupled_model; Δt=90, stop_time=10days)
+simulation = Simulation(coupled_model; Δt=60, stop_time=30days)
 
 # We define a callback function to monitor the simulation's progress,
 
@@ -176,7 +179,7 @@ run!(simulation)
 
 # After the initial spin up of 10 days, we can increase the time-step and run for longer.
 
-simulation.stop_time = 60days
+simulation.stop_time = 120days
 simulation.Δt = 10minutes
 run!(simulation)
 

examples/one_degree_simulation.jl

@@ -92,7 +92,7 @@ tracer_advection   = Centered()
 ocean = ocean_simulation(grid;
                          momentum_advection,
                          tracer_advection,
-                         closure,
+#                          closure,
                          forcing,
                          free_surface)
 

experiments/one_degree_simulation/Manifest.toml

@@ -271,13 +271,13 @@ weakdeps = ["SparseArrays"]
 deps = ["Adapt", "CFTime", "CUDA", "ClimaSeaIce", "CubicSplines", "DataDeps", "Dates", "Downloads", "ImageMorphology", "JLD2", "KernelAbstractions", "MPI", "NCDatasets", "Oceananigans", "OffsetArrays", "OrthogonalSphericalShellGrids", "Printf", "Scratch", "SeawaterPolynomials", "StaticArrays", "Statistics", "SurfaceFluxes", "Thermodynamics"]
 path = "../.."
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
-version = "0.4.2"
+version = "0.4.7"
 
 [[deps.ClimaSeaIce]]
 deps = ["Adapt", "KernelAbstractions", "Oceananigans", "RootSolvers", "Roots", "SeawaterPolynomials"]
-git-tree-sha1 = "e53033740cf85c15ddcbe55eadf15aeb6b816d2e"
+git-tree-sha1 = "f480a37fa4933fa1c698b24bde75cbfc27f5ce72"
 uuid = "6ba0ff68-24e6-4315-936c-2e99227c95a4"
-version = "0.1.5"
+version = "0.2.4"
 
 [[deps.CloseOpenIntervals]]
 deps = ["Static", "StaticArrayInterface"]
@@ -333,9 +333,9 @@ version = "0.13.0"
 
 [[deps.CommonDataModel]]
 deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf", "Statistics"]
-git-tree-sha1 = "98d64d5b9e5263884276656a43c45424b3a645c2"
+git-tree-sha1 = "358bf5a7d5c1387b995a43577673290c5d344758"
 uuid = "1fbeeb36-5f17-413c-809b-666fb144f157"
-version = "0.3.7"
+version = "0.3.8"
 
 [[deps.CommonSolve]]
 git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
@@ -485,10 +485,10 @@ uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
 version = "1.15.1"
 
 [[deps.DiskArrays]]
-deps = ["LRUCache", "Mmap", "OffsetArrays"]
-git-tree-sha1 = "5109a9314f8904f96a14b80d40dce5c9972c584f"
+deps = ["LRUCache", "OffsetArrays"]
+git-tree-sha1 = "ef25c513cad08d7ebbed158c91768ae32f308336"
 uuid = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
-version = "0.4.10"
+version = "0.3.23"
 
 [[deps.Distances]]
 deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
@@ -624,9 +624,9 @@ version = "0.8.3"
 
 [[deps.FilePathsBase]]
 deps = ["Compat", "Dates"]
-git-tree-sha1 = "2ec417fc319faa2d768621085cc1feebbdee686b"
+git-tree-sha1 = "3bab2c5aa25e7840a4b065805c0cdfc01f3068d2"
 uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
-version = "0.9.23"
+version = "0.9.24"
 weakdeps = ["Mmap", "Test"]
 
     [deps.FilePathsBase.extensions]
@@ -797,10 +797,10 @@ uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
 version = "1.0.2"
 
 [[deps.HDF5_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
-git-tree-sha1 = "38c8874692d48d5440d5752d6c74b0c6b0b60739"
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "OpenSSL_jll", "TOML", "Zlib_jll", "libaec_jll"]
+git-tree-sha1 = "e94f84da9af7ce9c6be049e9067e511e17ff89ec"
 uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
-version = "1.14.2+1"
+version = "1.14.6+0"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
@@ -1273,9 +1273,9 @@ version = "0.20.22"
 
 [[deps.MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "e7159031670cee777cc2840aef7a521c3603e36c"
+git-tree-sha1 = "3aa3210044138a1749dbd350a9ba8680869eb503"
 uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
-version = "4.3.0+0"
+version = "4.3.0+1"
 
 [[deps.MPIPreferences]]
 deps = ["Libdl", "Preferences"]
@@ -1285,9 +1285,9 @@ version = "0.1.11"
 
 [[deps.MPItrampoline_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "97aac4a518b6f01851f8821272780e1ba56fe90d"
+git-tree-sha1 = "ff91ca13c7c472cef700f301c8d752bc2aaff1a8"
 uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
-version = "5.5.2+0"
+version = "5.5.3+0"
 
 [[deps.MacroTools]]
 git-tree-sha1 = "72aebe0b5051e5143a079a4685a46da330a40472"
@@ -1382,9 +1382,9 @@ version = "2023.1.10"
 
 [[deps.NCDatasets]]
 deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "DiskArrays", "NetCDF_jll", "NetworkOptions", "Printf"]
-git-tree-sha1 = "2c9dc92001ac06d432f363f37ff5552954d9947c"
+git-tree-sha1 = "4704758e7dd9c843f0a99b18a26aa2d88dd8b8e6"
 uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
-version = "0.14.6"
+version = "0.14.0"
 
 [[deps.NVTX]]
 deps = ["Colors", "JuliaNVTXCallbacks_jll", "Libdl", "NVTX_jll"]
@@ -1405,10 +1405,10 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 version = "1.1.2"
 
 [[deps.NetCDF_jll]]
-deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenMPI_jll", "XML2_jll", "Zlib_jll", "Zstd_jll", "libzip_jll"]
-git-tree-sha1 = "a8af1798e4eb9ff768ce7fdefc0e957097793f15"
+deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LazyArtifacts", "LibCURL_jll", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "TOML", "XML2_jll", "Zlib_jll", "Zstd_jll", "libaec_jll", "libzip_jll"]
+git-tree-sha1 = "d574803b6055116af212434460adf654ce98e345"
 uuid = "7243133f-43d8-5620-bbf4-c2c921802cf3"
-version = "400.902.209+0"
+version = "401.900.300+0"
 
 [[deps.Netpbm]]
 deps = ["FileIO", "ImageCore", "ImageMetadata"]
@@ -1427,17 +1427,18 @@ version = "0.5.5"
 
 [[deps.Oceananigans]]
 deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "GPUArrays", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "StaticArrays", "Statistics", "StructArrays", "TimesDates"]
-path = "/Users/gregorywagner/Projects/alt/Oceananigans.jl"
+git-tree-sha1 = "9ef8914f4d0b0f3c50cbca229b82ae4b5e854481"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.95.15"
+version = "0.95.23"
 
     [deps.Oceananigans.extensions]
     OceananigansEnzymeExt = "Enzyme"
     OceananigansMakieExt = ["MakieCore", "Makie"]
     OceananigansMetalExt = "Metal"
-    OceananigansReactantExt = ["Reactant", "KernelAbstractions"]
+    OceananigansReactantExt = ["Reactant", "KernelAbstractions", "ConstructionBase"]
 
     [deps.Oceananigans.weakdeps]
+    ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
     Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
     Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
     MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
@@ -1483,9 +1484,9 @@ version = "0.8.1+2"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"]
-git-tree-sha1 = "6c1cf6181ffe0aa33eb33250ca2a60e54a15ea66"
+git-tree-sha1 = "da913f03f17b449951e0461da960229d4a3d1a8c"
 uuid = "fe0851c0-eecd-5654-98d4-656369965a5c"
-version = "5.0.7+0"
+version = "5.0.7+1"
 
 [[deps.OpenSSL]]
 deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"]
@@ -1518,9 +1519,9 @@ version = "1.8.0"
 
 [[deps.OrthogonalSphericalShellGrids]]
 deps = ["Adapt", "Distances", "KernelAbstractions", "MPI", "Oceananigans", "OffsetArrays", "Statistics"]
-git-tree-sha1 = "d11fd9caa8246493ca6da68bf00fef6a86bfa835"
+git-tree-sha1 = "eadf4cc0db734a18b611afdaaf91d607dbdee2d8"
 uuid = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0"
-version = "0.2.1"
+version = "0.2.2"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -2100,9 +2101,9 @@ version = "0.5.0"
 
 [[deps.TimeZones]]
 deps = ["Artifacts", "Dates", "Downloads", "InlineStrings", "Mocking", "Printf", "Scratch", "TZJData", "Unicode", "p7zip_jll"]
-git-tree-sha1 = "38bb1023fb94bfbaf2a29e1e0de4bbba6fe0bf6d"
+git-tree-sha1 = "2c705e96825b66c4a3f25031a683c06518256dd3"
 uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
-version = "1.21.2"
+version = "1.21.3"
 weakdeps = ["RecipesBase"]
 
     [deps.TimeZones.extensions]
@@ -2383,9 +2384,9 @@ version = "2.0.3+0"
 
 [[deps.libpng_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "055a96774f383318750a1a5e10fd4151f04c29c5"
+git-tree-sha1 = "068dfe202b0a05b8332f1e8e6b4080684b9c7700"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.46+0"
+version = "1.6.47+0"
 
 [[deps.libsixel_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "libpng_jll"]

experiments/one_degree_simulation/one_degree_simulation.jl

@@ -5,28 +5,28 @@ using Oceananigans
 using Oceananigans.Units
 using Dates
 using Printf
-using GLMakie
+
+filename = "tidal_potential_jra55.jld2"
+url = "https://www.dropbox.com/scl/fi/5qu30pcuy1z5q0s5xpwaj/" *
+      "tidal_potential_jra55.jld2?rlkey=2597z4vkiuq3stjsry45rpfyd&st=n42stc80&dl=0"
+
+if !isfile(filename)
+    Base.download(url, filename)
+end
 
 Oceananigans.defaults.FloatType = Float64
-arch = CPU()
-Nx = 256
-Ny = 128
-Nz = 32
+arch = GPU()
+Nx = 360 * 4
+Ny = 170 * 4
+Nz = 60
 z_faces = exponential_z_faces(; Nz, depth=6000, h=34)
-underlying_grid = TripolarGrid(arch; size=(Nx, Ny, Nz), z=z_faces)
+underlying_grid = TripolarGrid(arch; size=(Nx, Ny, Nz), halo=(7, 7, 7), z=z_faces)
 
-bottom_height = regrid_bathymetry(underlying_grid; minimum_depth=30, interpolation_passes=20, major_basins=1)
-view(bottom_height, 73:78, 88:89, 1) .= -1000 # open Gibraltar strait 
+bottom_height = regrid_bathymetry(underlying_grid; minimum_depth=20, interpolation_passes=40, major_basins=3)
+#view(bottom_height, 73:78, 88:89, 1) .= -1000 # open Gibraltar strait 
 
 grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); active_cells_map=true)
 
-catke = ClimaOcean.OceanSimulations.default_ocean_closure()
-viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000)
-closure = (catke, viscous_closure)
-
-#gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000)
-#closure = (gm, catke, viscous_closure)
-
 dates = DateTime(1993, 1, 1) : Month(1) : DateTime(1993, 11, 1)
 mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(-100, 0))
 temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly())
@@ -36,19 +36,34 @@ FT = ECCORestoring(temperature, grid; mask, rate)
 FS = ECCORestoring(salinity, grid; mask, rate)
 forcing = (T=FT, S=FS)
 
+#=
+catke = ClimaOcean.OceanSimulations.default_ocean_closure()
+viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000)
+gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000)
+closure = (gm, catke, viscous_closure)
 momentum_advection = VectorInvariant()
 tracer_advection = Centered(order=2)
 free_surface = SplitExplicitFreeSurface(grid; substeps=70)
 ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, forcing)
+=#
+
+free_surface = SplitExplicitFreeSurface(grid; substeps=70)
+ocean = ocean_simulation(grid; forcing, free_surface, bottom_drag_coefficient=0.01)
 
 set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
                   S=ECCOMetadata(:salinity;    dates=first(dates)))
 
 radiation  = Radiation(arch)
-atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
-coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) 
+tidal_potential = FieldTimeSeries("tidal_potential_jra55.jld2", "Φ"; architecture=GPU(), backend=InMemory(41))
+boundary_conditions = FieldBoundaryConditions(tidal_potential.grid, (Center, Center, Nothing))
 
-simulation = Simulation(coupled_model; Δt=10minutes, stop_iteration=100)
+tidal_potential = FieldTimeSeries("tidal_potential_jra55.jld2", "Φ"; architecture=GPU(), backend=InMemory(41), boundary_conditions) 
+Oceananigans.BoundaryConditions.fill_halo_regions!(tidal_potential)
+atmosphere = JRA55PrescribedAtmosphere(arch; tidal_potential, backend=JRA55NetCDFBackend(41))
+
+# atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
+coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) 
+simulation = Simulation(coupled_model; Δt=1minutes, stop_time=5days)
 
 wall_time = Ref(time_ns())
 
@@ -75,5 +90,23 @@ end
 
 add_callback!(simulation, progress, IterationInterval(10))
 
+u, v, w = ocean.model.velocities
+ζ = ∂x(v) - ∂y(u)
+s = @at (Center, Center, Center) sqrt(u^2 + v^2)
+outputs = merge(ocean.model.tracers, ocean.model.velocities, (; ζ, s))
+simulation.output_writers[:surface] = JLD2OutputWriter(ocean.model, outputs;
+                                                       schedule = TimeInterval(3hours),
+                                                       filename = "tidally_forced_surface_fields",
+                                                       indices = (:, :, grid.Nz),
+                                                       overwrite_existing = true,
+                                                       array_type = Array{Float32})
+
+simulation.output_writers[:middepth] = JLD2OutputWriter(ocean.model, outputs;
+                                                        schedule = TimeInterval(3hours),
+                                                        filename = "tidally_forced_deep_fields",
+                                                        indices = (:, :, Int(grid.Nz/2)),
+                                                        overwrite_existing = true,
+                                                        array_type = Array{Float32})
+
 run!(simulation)
 

src/Bathymetry.jl

@@ -205,6 +205,7 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
     Nλt, Nφt = Nt = size(target_grid)
     Nλn, Nφn = Nn = size(native_z)
 
+    #=
     resxt = minimum_xspacing(target_grid)
     resyt = minimum_yspacing(target_grid)
 
@@ -222,10 +223,11 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
                      "is finer than the target grid in both horizontal directions.")
         return target_z
     end
-
+    =#
+
     # Interpolate in passes
-    latitude  = y_domain(native_z.grid)
     longitude = x_domain(native_z.grid)
+    latitude  = y_domain(native_z.grid)
 
     ΔNλ = floor((Nλn - Nλt) / passes)
     ΔNφ = floor((Nφn - Nφt) / passes)
@@ -241,7 +243,6 @@ function interpolate_bathymetry_in_passes(native_z, target_grid;
 
     for pass = 1:passes - 1
         new_size = (Nλ[pass], Nφ[pass], 1)
-
         @debug "Bathymetry interpolation pass $pass with size $new_size"
 
         new_grid = LatitudeLongitudeGrid(architecture(target_grid), Float32,