Add parameter to epsilon equation so buoyancy flux term may always be positive

src/Grids/rectilinear_grid.jl

@@ -284,6 +284,8 @@ function RectilinearGrid(architecture::AbstractArchitecture = CPU(),
                                        Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
 end
 
+
+
 """ Validate user input arguments to the `RectilinearGrid` constructor. """
 function validate_rectilinear_grid_args(topology, size, halo, FT, extent, x, y, z)
     TX, TY, TZ = topology = validate_topology(topology)
@@ -343,6 +345,20 @@ function Base.show(io::IO, grid::RectilinearGrid, withsummary=true)
                      "└── ", z_summary)
 end
 
+#####
+##### For "column ensemble models"
+#####
+
+struct ColumnEnsembleSize{C<:Tuple{Int, Int}}
+    ensemble :: C
+    Nz :: Int
+    Hz :: Int
+end
+
+ColumnEnsembleSize(; Nz, ensemble=(0, 0), Hz=1) = ColumnEnsembleSize(ensemble, Nz, Hz)
+validate_size(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(e.ensemble[1], e.ensemble[2], e.Nz)
+validate_halo(TX, TY, TZ, size, e::ColumnEnsembleSize) = tuple(0, 0, e.Hz)
+
 #####
 ##### Utilities
 #####
@@ -378,10 +394,16 @@ function constructor_arguments(grid::RectilinearGrid)
 
     # Kwargs
     topo = topology(grid)
-    size = (grid.Nx, grid.Ny, grid.Nz)
-    halo = (grid.Hx, grid.Hy, grid.Hz)
-    size = pop_flat_elements(size, topo)
-    halo = pop_flat_elements(halo, topo)
+
+    if (topo[1] == Flat && grid.Nx > 1) ||
+       (topo[2] == Flat && grid.Ny > 1)
+        size = halo = ColumnEnsembleSize(Nz=grid.Nz, Hz=grid.Hz, ensemble=(grid.Nx, grid.Ny))
+    else
+        size = (grid.Nx, grid.Ny, grid.Nz)
+        halo = (grid.Hx, grid.Hy, grid.Hz)
+        size = pop_flat_elements(size, topo)
+        halo = pop_flat_elements(halo, topo)
+    end
 
     kwargs = Dict(:size => size,
                   :halo => halo,

src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl

@@ -2,13 +2,12 @@ using CUDA: @allowscalar
 
 using Oceananigans: UpdateStateCallsite
 using Oceananigans.Advection: AbstractAdvectionScheme
-using Oceananigans.Grids: Flat, Bounded
+using Oceananigans.Grids: Flat, Bounded, ColumnEnsembleSize
 using Oceananigans.Fields: ZeroField
 using Oceananigans.Coriolis: AbstractRotation
 using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure
 using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVDArray
 
-import Oceananigans.Grids: validate_size, validate_halo
 import Oceananigans.Models: validate_tracer_advection
 import Oceananigans.BoundaryConditions: fill_halo_regions!
 import Oceananigans.TurbulenceClosures: time_discretization, compute_diffusivities!
@@ -102,17 +101,6 @@ end
     return ∇_dot_qᶜ(i, j, k, grid, closure, c, tracer_index, args...)
 end
 
-struct ColumnEnsembleSize{C<:Tuple{Int, Int}}
-    ensemble :: C
-    Nz :: Int
-    Hz :: Int
-end
-
-ColumnEnsembleSize(; Nz, ensemble=(0, 0), Hz=1) = ColumnEnsembleSize(ensemble, Nz, Hz)
-
-validate_size(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(e.ensemble[1], e.ensemble[2], e.Nz)
-validate_halo(TX, TY, TZ, size, e::ColumnEnsembleSize) = tuple(0, 0, e.Hz)
-
 @inline function time_discretization(closure_array::AbstractArray)
     first_closure = @allowscalar first(closure_array) # assumes all closures have same time-discretization
     return time_discretization(first_closure)

src/OutputReaders/field_time_series.jl

@@ -214,16 +214,16 @@ Base.length(backend::PartlyInMemory) = backend.length
 #####
 
 mutable struct FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, P, N, KW} <: AbstractField{LX, LY, LZ, G, ET, 4}
-                   data :: D
-                   grid :: G
-                backend :: K
+    data :: D
+    grid :: G
+    backend :: K
     boundary_conditions :: B
-                indices :: I
-                  times :: χ
-                   path :: P
-                   name :: N
-          time_indexing :: TI
-             reader_kw :: KW
+    indices :: I
+    times :: χ
+    path :: P
+    name :: N
+    time_indexing :: TI
+    reader_kw :: KW
 
     function FieldTimeSeries{LX, LY, LZ}(data::D,
                                          grid::G,
@@ -350,7 +350,8 @@ new_data(FT, grid, loc, indices, ::Nothing) = nothing
 
 # Apparently, not explicitly specifying Int64 in here makes this function
 # fail on x86 processors where `Int` is implied to be `Int32`
-# see ClimaOcean commit 3c47d887659d81e0caed6c9df41b7438e1f1cd52 at https://github.com/CliMA/ClimaOcean.jl/actions/runs/8804916198/job/24166354095)
+# see ClimaOcean commit 3c47d887659d81e0caed6c9df41b7438e1f1cd52 at
+# https://github.com/CliMA/ClimaOcean.jl/actions/runs/8804916198/job/24166354095)
 function new_data(FT, grid, loc, indices, Nt::Union{Int, Int64})
     space_size = total_size(grid, loc, indices)
     underlying_data = zeros(FT, architecture(grid), space_size..., Nt)

src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl

@@ -10,7 +10,8 @@ using CUDA
 Base.@kwdef struct TKEDissipationEquations{FT}
     Cᵋϵ :: FT = 1.92
     Cᴾϵ :: FT = 1.44
-    Cᵇϵ :: FT = -0.65
+    Cᵇϵ⁺ :: FT = -0.65
+    Cᵇϵ⁻ :: FT = -0.65
     Cᵂu★ :: FT = 0.0
     CᵂwΔ :: FT = 0.0
     Cᵂα  :: FT = 0.11 # Charnock parameter
@@ -133,7 +134,11 @@ end
 
     # Patankar trick for ϵ-equation
     Cᵋϵ = closure_ij.tke_dissipation_equations.Cᵋϵ
-    Cᵇϵ = closure_ij.tke_dissipation_equations.Cᵇϵ
+    Cᵇϵ⁺ = closure_ij.tke_dissipation_equations.Cᵇϵ⁺
+    Cᵇϵ⁻ = closure_ij.tke_dissipation_equations.Cᵇϵ⁻
+
+    N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)
+    Cᵇϵ = ifelse(N² ≥ 0, Cᵇϵ⁺, Cᵇϵ⁻) 
 
     Cᵇϵ_wb⁻ = min(Cᵇϵ * wb, zero(grid))
     Cᵇϵ_wb⁺ = max(Cᵇϵ * wb, zero(grid))