Updates diagnostic for computing mixed layer depth

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

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)

src/ClimaOcean.jl

@@ -78,7 +78,7 @@ include("VerticalGrids.jl")
 include("InitialConditions/InitialConditions.jl")
 include("DataWrangling/DataWrangling.jl")
 include("Bathymetry.jl")
-include("Diagnostics.jl")
+include("Diagnostics/Diagnostics.jl")
 include("OceanSimulations/OceanSimulations.jl")
 
 using .VerticalGrids