Typo in generate_surface_fluxes.jl + add mediterranean example

Project.toml

@@ -19,6 +19,7 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OrthogonalSphericalShellGrids = "c2be9673-fb75-4747-82dc-aa2bb9f4aed0"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"

docs/make.jl

@@ -17,7 +17,7 @@ to_be_literated = [
     "generate_bathymetry.jl",
     "generate_surface_fluxes.jl",
     "single_column_os_papa_simulation.jl",
-    # "mediterranean_simulation_with_ecco_restoring.jl",
+    "mediterranean_simulation_with_ecco_restoring.jl",
     "near_global_ocean_simulation.jl"
 ]
 
@@ -45,7 +45,7 @@ pages = [
         "Generate bathymetry" => "literated/generate_bathymetry.md",
         "Surface fluxes" => "literated/generate_surface_fluxes.md",
         "Single column simulation" => "literated/single_column_os_papa_simulation.md",
-        # "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md",
+        "Mediterranean simulation with ECCO restoring" => "literated/mediterranean_simulation_with_ecco_restoring.md",
         "Near-global Ocean simulation" => "literated/near_global_ocean_simulation.md",
         ],
 

examples/generate_surface_fluxes.jl

@@ -18,7 +18,7 @@ using CairoMakie
 
 # # Computing fluxes on the ECCO2 grid
 #
-# We start by building the ECCO2 grid, using `ECCO_bottom_height` to identify the bottom height.
+# We start by building the ECCO2 grid, using `ECCO_immersed_grid` function.
 
 grid = ECCO_immersed_grid()
 

examples/mediterranean_simulation_with_ecco_restoring.jl

@@ -42,8 +42,8 @@ z_faces = stretched_vertical_faces(depth = 5000,
                                    stretching = PowerLawStretching(1.070), 
                                    surface_layer_height = 50)
 
-Nx = 15 * Int(λ₂ - λ₁) # 1/15 th of a degree resolution
-Ny = 15 * Int(φ₂ - φ₁) # 1/15 th of a degree resolution
+Nx = 15 * Int(λ₂ - λ₁) # 1/15th of a degree resolution
+Ny = 15 * Int(φ₂ - φ₁) # 1/15th of a degree resolution
 Nz = length(z_faces) - 1
 
 grid = LatitudeLongitudeGrid(GPU();
@@ -60,26 +60,24 @@ grid = LatitudeLongitudeGrid(GPU();
 # are adjusted to refine the bathymetry representation.
 
 bottom_height = regrid_bathymetry(grid, 
-                                  height_above_water = 1,
                                   minimum_depth = 10,
-                                  interpolation_passes = 25,
-                                  connected_regions_allowed = 1)
+                                  interpolation_passes = 10,
+                                  major_basins = 1)
 
 grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
 
 # ## Downloading ECCO data
 #
 # The model is initialized with temperature and salinity fields from the ECCO dataset, 
-# using the function `ECCO_restoring_forcing` to apply restoring forcings for these tracers. 
+# using the function `ECCORestoring` to apply restoring forcings for these tracers. 
 # This allows us to nudge the model towards realistic temperature and salinity profiles.
 
-dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 12, 1)
+@inline surface_mask(x, y, z, t) = z > - 50 
 
-temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
-salinity    = ECCOMetadata(:salinity,    dates, ECCO4Monthly())
+dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 12, 1)
 
-FT = ECCO_restoring_forcing(temperature; architecture = GPU(), timescale = 2days)
-FS = ECCO_restoring_forcing(salinity;    architecture = GPU(), timescale = 2days)
+FT = ECCORestoring(GPU(), :temperature; dates, mask=surface_mask, version=ECCO4Monthly(), rate=1/5days)
+FS = ECCORestoring(GPU(), :salinity;    dates, mask=surface_mask, version=ECCO4Monthly(), rate=1/5days)
 
 # Constructing the Simulation
 #
@@ -94,13 +92,17 @@ ocean = ocean_simulation(grid; forcing = (T = FT, S = FS))
 # In this case, our ECCO dataset has access to a temperature and a salinity
 # field, so we initialize temperature T and salinity S from ECCO.
 
-set!(ocean.model, T = temperature[1], S = salinity[1])
+set!(ocean.model, T=ECCOMetadata(:temperature; dates=dates[1]), 
+                  S=ECCOMetadata(:salinity;    dates=dates[1]))
 
 fig = Figure()
 ax  = Axis(fig[1, 1])
-heatmap!(ax, interior(model.tracers.T, :, :, Nz), colorrange = (10, 20), colormap = :thermal)
+heatmap!(ax, view(model.tracers.T, :, :, Nz), colorrange = (10, 20), colormap = :thermal)
 ax  = Axis(fig[1, 2])
-heatmap!(ax, interior(model.tracers.S, :, :, Nz), colorrange = (35, 40), colormap = :haline)
+heatmap!(ax, view(model.tracers.S, :, :, Nz), colorrange = (35, 40), colormap = :haline)
+
+save("initial_conditions.png", fig)
+# ![](initial_conditions.png)
 
 function progress(sim) 
     u, v, w = sim.model.velocities
@@ -129,17 +131,12 @@ ocean.stop_iteration = 1000
 run!(ocean)
 
 # ## Run the real simulation
-#
-# Now that the solution has adjusted to the bathymetry we can ramp up the time
-# step size. We use a `TimeStepWizard` to automatically adapt to a CFL of 0.2.
-
-wizard = TimeStepWizard(; cfl = 0.2, max_Δt = 10minutes, max_change = 1.1)
 
-ocean.callbacks[:wizard] = Callback(wizard, IterationInterval(10))
+# Let's reset the maximum number of iterations and we can not increase the time step size
 
-# Let's reset the maximum number of iterations
+ocean.Δt = 3minutes
 ocean.stop_iteration = Inf
-ocean.stop_time = 200days
+ocean.stop_time = 100days
 
 ocean.output_writers[:surface_fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers);
                                                          indices = (:, :, Nz),
@@ -164,45 +161,23 @@ S_series = FieldTimeSeries("med_surface_field.jld2", "S")
 c_series = FieldTimeSeries("med_surface_field.jld2", "c")
 iter = Observable(1)
 
-u = @lift begin
-    f = interior(u_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-v = @lift begin
-    f = interior(v_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-T = @lift begin
-    f = interior(T_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-S = @lift begin
-    f = interior(S_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
-c = @lift begin
-    f = interior(c_series[$iter], :, :, 1)
-    f[f .== 0] .= NaN
-    f
-end
+u = @lift(u_series[$iter])
+v = @lift(v_series[$iter])
+T = @lift(T_series[$iter])
+S = @lift(S_series[$iter])
 
 fig = Figure()
 ax  = Axis(fig[1, 1], title = "surface zonal velocity ms⁻¹")
-heatmap!(u)
+heatmap!(ax, u)
 ax  = Axis(fig[1, 2], title = "surface meridional velocity ms⁻¹")
-heatmap!(v)
+heatmap!(ax, v)
 ax  = Axis(fig[2, 1], title = "surface temperature ᵒC")
-heatmap!(T)
+heatmap!(ax, T)
 ax  = Axis(fig[2, 2], title = "surface salinity psu")
-heatmap!(S)
-ax  = Axis(fig[2, 3], title = "passive tracer -")
-heatmap!(c)
+heatmap!(ax, S)
 
 CairoMakie.record(fig, "mediterranean_video.mp4", 1:length(u_series.times); framerate = 5) do i
     @info "recording iteration $i"
     iter[] = i    
-end
+end
+# ![](mediterranean_video.mp4)

src/Bathymetry.jl

@@ -17,8 +17,6 @@ using Oceananigans.BoundaryConditions
 using KernelAbstractions: @kernel, @index
 using JLD2
 
-using ClimaOcean
-
 using NCDatasets
 using Downloads
 using Printf
@@ -92,19 +90,26 @@ function regrid_bathymetry(target_grid;
                            major_basins = Inf) # Allow an `Inf` number of ``lakes''
 
     filepath = joinpath(dir, filename)
-    fileurl  = joinpath(url, filename)
-
-    @root begin # perform all this only on rank 0, aka the "root" rank
-        if !isfile(filepath)
-            try 
-                Downloads.download(fileurl, filepath; progress=download_progress, verbose=true)
-            catch 
-                cmd = `wget --no-check-certificate -O $filepath $fileurl`
-                @root run(cmd)
-            end
+
+    if isfile(filepath)
+        @info "Regridding bathymetry from existing file $filepath."
+    else
+        @info "Downloading bathymetry..."
+        if !ispath(dir)
+            @info "Making bathymetry directory $dir..."
+            mkdir(dir)
+        end
+
+        fileurl = joinpath(url, filename)
+
+        try 
+            Downloads.download(fileurl, filepath; progress=download_progress, verbose=true)
+        catch 
+            cmd = `wget --no-check-certificate -O $filepath $fileurl`
+            run(cmd)
         end
     end
-    
+
     dataset = Dataset(filepath)
 
     FT = eltype(target_grid)