Merge branch 'main' into ss/one-degree-take3

Project.toml

@@ -2,7 +2,7 @@ name = "ClimaOcean"
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
 license = "MIT"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

examples/ecco_inspect_temperature_salinity.jl

@@ -0,0 +1,87 @@
+using Oceananigans
+using Oceananigans.ImmersedBoundaries: mask_immersed_field!
+
+using GLMakie
+using Printf
+using ClimaOcean
+using ClimaOcean.DataWrangling.ECCO: ECCO_field, ECCOFieldTimeSeries
+using CFTime
+using Dates
+
+arch = CPU()
+Nx = 360 ÷ 4
+Ny = 160 ÷ 4
+
+z = ClimaOcean.DataWrangling.ECCO.ECCO_z
+z = z[20:end]
+Nz = length(z) - 1
+
+grid = LatitudeLongitudeGrid(arch; z,
+                             size = (Nx, Ny, Nz),
+                             latitude  = (-80, 80),
+                             longitude = (0, 360))
+
+bottom_height = regrid_bathymetry(grid; 
+                                  minimum_depth = 10,
+                                  interpolation_passes = 5,
+                                  major_basins = 1)
+
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
+
+T = CenterField(grid)
+S = CenterField(grid)
+
+using SeawaterPolynomials: TEOS10EquationOfState
+using Oceananigans.BuoyancyModels: buoyancy
+
+equation_of_state = TEOS10EquationOfState()
+sb = SeawaterBuoyancy(; equation_of_state)
+tracers = (T=T, S=S)
+b = Field(buoyancy(sb, grid, tracers))
+
+start = DateTimeProlepticGregorian(1993, 1, 1)
+stop  = DateTimeProlepticGregorian(1999, 1, 1)
+dates = range(start; stop, step=Month(1))
+
+Tmeta = ECCOMetadata(:temperature; dates)
+Smeta = ECCOMetadata(:salinity; dates)
+
+Tt = ECCOFieldTimeSeries(Tmeta, grid; time_indices_in_memory=length(dates))
+St = ECCOFieldTimeSeries(Smeta, grid; time_indices_in_memory=length(dates))
+
+fig = Figure(size=(900, 1050))
+
+axT = Axis(fig[1, 1])
+axS = Axis(fig[2, 1])
+axb = Axis(fig[3, 1])
+
+Nt = length(dates)
+grid = T.grid
+Nz = size(grid, 3)
+kslider = Slider(fig[1:3, 0], range=1:Nz, startvalue=Nz, horizontal=false)
+nslider = Slider(fig[4, 1:2], range=1:Nt, startvalue=1)
+k = kslider.value
+n = nslider.value
+
+Tk = @lift view(Tt[$n], :, :, $k)
+Sk = @lift view(St[$n], :, :, $k)
+
+Δb = @lift begin
+    parent(T) .= parent(Tt[$n])
+    parent(S) .= parent(St[$n])
+    compute!(b)
+    mask_immersed_field!(b, NaN)
+    Δb = interior(b, :, :, Nz) .- interior(b, :, :, $k)
+    Δb
+end
+
+hmT = heatmap!(axT, Tk, nan_color=:lightgray, colorrange=(-2, 30), colormap=:thermal)
+hmS = heatmap!(axS, Sk, nan_color=:lightgray, colorrange=(31, 37), colormap=:haline)
+hmb = heatmap!(axb, Δb, nan_color=:lightgray, colorrange=(0, 1e-3), colormap=:magma)
+
+Colorbar(fig[1, 2], hmT, label="Temperature (ᵒC)")
+Colorbar(fig[2, 2], hmS, label="Salinity (g kg⁻¹)")
+Colorbar(fig[3, 2], hmb, label="Buoyancy difference (m s⁻²)")
+
+display(fig)
+

examples/ecco_mixed_layer_depth.jl

@@ -0,0 +1,81 @@
+using ClimaOcean
+using ClimaOcean.Diagnostics: MixedLayerDepthField
+using ClimaOcean.DataWrangling.ECCO: ECCO_field, ECCOFieldTimeSeries
+using Oceananigans
+using GLMakie
+using Printf
+using CFTime
+using Dates
+
+using SeawaterPolynomials: TEOS10EquationOfState
+using Oceananigans.BuoyancyFormulations: buoyancy
+
+arch = CPU()
+Nx = 360 ÷ 1
+Ny = 160 ÷ 1
+
+z = ClimaOcean.DataWrangling.ECCO.ECCO_z
+z = z[20:end]
+Nz = length(z) - 1
+
+grid = LatitudeLongitudeGrid(arch; z,
+                             size = (Nx, Ny, Nz),
+                             latitude  = (-80, 80),
+                             longitude = (0, 360))
+
+bottom_height = regrid_bathymetry(grid; 
+                                  minimum_depth = 10,
+                                  interpolation_passes = 5,
+                                  major_basins = 1)
+
+grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height))
+
+start = DateTimeProlepticGregorian(1993, 1, 1)
+stop  = DateTimeProlepticGregorian(2003, 1, 1)
+dates = range(start; stop, step=Month(1))
+
+Tmeta = ECCOMetadata(:temperature; dates)
+Smeta = ECCOMetadata(:salinity; dates)
+
+Tt = ECCOFieldTimeSeries(Tmeta, grid; time_indices_in_memory=2)
+St = ECCOFieldTimeSeries(Smeta, grid; time_indices_in_memory=2)
+ht = FieldTimeSeries{Center, Center, Nothing}(grid, Tt.times)
+
+equation_of_state = TEOS10EquationOfState()
+sb = SeawaterBuoyancy(; equation_of_state)
+tracers = (T=Tt[1], S=St[1])
+h = MixedLayerDepthField(sb, grid, tracers)
+
+Nt = length(ht)
+for n = 1:Nt-1
+    local tracers
+    tracers = (T=Tt[n], S=St[n])
+    h.operand.buoyancy_perturbation = buoyancy(sb, grid, tracers)
+    @show n
+    @time compute!(h)
+    parent(ht[n]) .= parent(h)
+end
+
+function titlestr(n)
+    d = dates[n]
+    yr = year(d)
+    mn = monthname(d)
+    return string("ECCO mixed layer depth on ", mn, " ", yr)
+end
+
+fig = Figure(size=(1500, 800))
+axh = Axis(fig[2, 1], xlabel="Longitude", ylabel="Latitude")
+n = Observable(1)
+
+str = @lift titlestr($n)
+Label(fig[1, 1], str, tellwidth=false)
+
+hn = @lift ht[$n]
+hm = heatmap!(axh, hn, colorrange=(0, 500), colormap=:magma, nan_color=:lightgray)
+Colorbar(fig[2, 2], hm, label="Mixed layer depth (m)")
+display(fig)
+
+record(fig, "ecco_mld.mp4", 1:Nt-1, framerate=4) do nn
+    @info "Drawing frame $nn of $Nt..."
+    n[] = nn
+end

experiments/three_degree_simulation/three_degree_simulation.jl

@@ -0,0 +1,72 @@
+using ClimaOcean
+using OrthogonalSphericalShellGrids
+using Oceananigans
+using Oceananigans.Units
+using CFTime
+using Dates
+using Printf
+
+arch = CPU()
+Nx = 120
+Ny = 60
+Nz = 50
+
+z_faces = exponential_z_faces(; Nz, depth=6000, h=34)
+
+underlying_grid = TripolarGrid(arch; size=(Nx, Ny, Nz), z=z_faces)
+bottom_height = regrid_bathymetry(underlying_grid)
+grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height))
+
+gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=4000, κ_symmetric=4000)
+catke = ClimaOcean.OceanSimulations.default_ocean_closure()
+viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=4000)
+
+dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 12, 1)
+temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
+salinity = ECCOMetadata(:salinity, dates, ECCO4Monthly())
+
+restoring_rate  = 1/2days
+mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90))
+FT = ECCORestoring(temperature, grid; mask, rate=restoring_rate)
+FS = ECCORestoring(salinity, grid; mask, rate=restoring_rate)
+
+ocean = ocean_simulation(grid;
+                         momentum_advection =  VectorInvariant(),
+                         tracer_advection = Centered(order=2),
+                         closure = (gm, catke, viscous_closure),
+                         forcing = (T=FT, S=FT),
+                         tracers = (:T, :S, :e))
+
+set!(ocean.model, T=ECCOMetadata(:temperature; dates=first(dates)),
+                  S=ECCOMetadata(:salinity; dates=first(dates)))
+
+radiation  = Radiation(arch)
+atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(20))
+coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) 
+simulation = Simulation(coupled_model; Δt=20minutes, stop_time=30days)
+
+wall_time = Ref(time_ns())
+
+function progress(sim)
+    ocean = sim.model.ocean
+    u, v, w = ocean.model.velocities
+    T = ocean.model.tracers.T
+    Tmax = maximum(interior(T))
+    Tmin = minimum(interior(T))
+    umax = (maximum(abs, interior(u)), maximum(abs, interior(v)), maximum(abs, interior(w)))
+    step_time = 1e-9 * (time_ns() - wall_time[])
+
+    @info @sprintf("Time: %s, n: %d, Δt: %s, max|u|: (%.2e, %.2e, %.2e) m s⁻¹, \
+                   extrema(T): (%.2f, %.2f) ᵒC, wall time: %s \n",
+                   prettytime(sim), iteration(sim), prettytime(sim.Δt),
+                   umax..., Tmax, Tmin, prettytime(step_time))
+
+    wall_time[] = time_ns()
+
+    return nothing
+end
+
+add_callback!(simulation, progress, IterationInterval(10))
+
+run!(simulation)
+

mwe.jl

@@ -0,0 +1,11 @@
+using ClimaOcean
+using NCDatasets
+
+cachepath = ClimaOcean.DataWrangling.JRA55.download_jra55_cache
+filename = "RYF.tas.1990_1991.nc"
+filepath = joinpath(cachepath, filename)
+
+ds = Dataset(filepath)
+Nx, Ny, Nt = size(ds["tas"])
+ds["tas"][1, 1, [Nt, 1]]
+close(ds)