Merge branch 'main' of github.com:CliMA/ClimaOcean.jl into main

examples/gm_one_degree_calibration/distributed_gm_one_degree_model_calibration.jl

@@ -93,13 +93,13 @@ save_indices = Dict(
 
 eki_iteration = 0
 
-function initialize_output_writers!(sim, save_indices, iteration, rank)
+function initialize_output_writers!(sim, iteration, rank)
     model = sim.model
     T, S  = model.tracers
     for (name, idx) in save_indices
         delete!(sim.output_writers, name)
 
-        output_prefix = prefix * "_eki_iteration" * string(iteration) * "_rank$(rank)"
+        output_prefix = prefix * string(name) * "_eki_iteration" * string(iteration) * "_particle$(rank)"
         sim.output_writers[name] = JLD2OutputWriter(model, (; T, S); dir,
                                                     schedule = TimeInterval(44days),
                                                     filename = output_prefix,
@@ -135,7 +135,7 @@ function initialize_simulation!(sim, parameters)
     set!(S, S₀)
     sim.model.clock.time = start_time
     global eki_iteration += 1
-    initialize_output_writers!(sim, save_indices, eki_iteration, rank)
+    initialize_output_writers!(sim, eki_iteration, rank)
     return nothing
 end
 

examples/gm_one_degree_calibration/gm_one_degree_model_calibration.jl

@@ -66,13 +66,13 @@ save_indices = Dict(
 
 eki_iteration = 0
 
-function initialize_output_writers!(sim, save_indices, iteration, rank)
+function initialize_output_writers!(sim, iteration, rank)
     model = sim.model
     T, S  = model.tracers
     for (name, idx) in save_indices
         delete!(sim.output_writers, name)
 
-        output_prefix = prefix * "_eki_iteration" * string(iteration) * "_rank$(rank)"
+        output_prefix = prefix * string(name) * "_eki_iteration" * string(iteration) * "_particle$(rank)"
         sim.output_writers[name] = JLD2OutputWriter(model, (; T, S); dir,
                                                     schedule = TimeInterval(44days),
                                                     filename = output_prefix,
@@ -105,7 +105,7 @@ function fake_function(sim, p)
 end
 
 for (idx, sim) in enumerate(simulation_ensemble)
-    initialize_output_writers!(sim, save_indices, eki_iteration, idx)
+    initialize_output_writers!(sim, eki_iteration, idx)
     sim.callbacks[:name] = Callback(fake_function, TimeInterval(1000years); parameters = idx)
 end
 
@@ -119,7 +119,7 @@ function initialize_simulation!(sim, parameters)
     set!(S, S₀)
     sim.model.clock.time = start_time
     global eki_iteration += 1
-    initialize_output_writers!(sim, save_indices, eki_iteration, sim.callbacks[:name].parameters)
+    initialize_output_writers!(sim, eki_iteration, sim.callbacks[:name].parameters)
     return nothing
 end
 

src/NearGlobalSimulations.jl

@@ -102,17 +102,45 @@ end
     return cyclic_interpolate(τ₁, τ₂, time)
 end
 
+using Oceananigans.Operators
 using Oceananigans.TurbulenceClosures
 using Oceananigans.Grids: min_Δx, min_Δy
-using Oceananigans.Operators: Δx, Δy
+using Oceananigans.Operators: Δx, Δy, ℑxyz
 
 @inline Δ²(i, j, k, grid, lx, ly, lz)  = (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2))
 @inline grid_dependent_biharmonic_viscosity(i, j, k, grid, lx, ly, lz, clock, fields, λ) = Δ²(i, j, k, grid, lx, ly, lz)^2 / λ
 
+@inline function leith_dynamic_biharmonic_viscosity(i, j, k, grid, lx, ly, lz, clock, fields, p)
+    location = (lx, ly, lz)
+    from_∂xζ = (Center(), Face(), Center())
+    from_∂yζ = (Face(), Center(), Center())
+    from_∂xδ = (Face(), Center(), Center())
+    from_∂yδ = (Center(), Face(), Center())
+
+    ∂xζ = ℑxyz(i, j, k, grid, from_∂xζ, location, ∂xᶜᶠᶜ, ζ₃ᶠᶠᶜ, fields.u, fields.v)
+    ∂yζ = ℑxyz(i, j, k, grid, from_∂yζ, location, ∂yᶠᶜᶜ, ζ₃ᶠᶠᶜ, fields.u, fields.v)
+    ∂xδ = ℑxyz(i, j, k, grid, from_∂xδ, location, ∂xᶠᶜᶜ, div_xyᶜᶜᶜ, fields.u, fields.v)
+    ∂yδ = ℑxyz(i, j, k, grid, from_∂yδ, location, ∂yᶜᶠᶜ, div_xyᶜᶜᶜ, fields.u, fields.v)
+
+    dynamic_visc = sqrt( p.Cζ * (∂xζ^2 + ∂yζ^2) + p.Cδ * (∂xδ^2 + ∂yδ^2) )
+
+    return Δ²(i, j, k, grid, lx, ly, lz)^2.5 * dynamic_visc
+end
+
 geometric_viscosity(formulation, timescale) = ScalarBiharmonicDiffusivity(formulation, ν = grid_dependent_biharmonic_viscosity, 
                                                                           discrete_form = true, 
                                                                           parameters = timescale)
 
+function leith_viscosity(formulation; C_vort = 3.0, C_div = 3.0)
+
+    Cζ = (C_vort / π)^6 / 8
+    Cδ = (C_div  / π)^6 / 8
+
+    return ScalarBiharmonicDiffusivity(formulation, ν = leith_dynamic_biharmonic_viscosity,
+                                       discrete_form = true,
+                                       parameters = (; Cζ, Cδ))
+end
+
 @inline ϕ²(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2
 @inline spᶠᶜᶜ(i, j, k, grid, Φ) = @inbounds sqrt(Φ.u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, ϕ², Φ.v))
 @inline spᶜᶠᶜ(i, j, k, grid, Φ) = @inbounds sqrt(Φ.v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, ϕ², Φ.u))