Periodic Coriolis force for planet with two equators

[rbenamran]

Hi,

I am trying to recreate a planet with two equators model (specifically in the ShallowWaterModel), as found in the supplementary materials from Delplace, Marston, and Venaille, “Topological Origin of Equatorial Waves.”, or also in Warneford and Dellar, “Super- and Sub-Rotating Equatorial Jets in Shallow Water Models of Jovian Atmospheres.”, referred to as the square planet.

I want to set up a model with a periodic Coriolis force in Oceananigans and am having trouble doing so. Is it possible to somehow set up a changing Coriolis force in y such that the sine wave vanishes twice?

Thanks!

[glwagner]

Oh my periodic Coriolis! There's two ways to achieve this. One is to implement the desired Coriolis force as a forcing function. I suspect you'll need to use DiscreteForcing, because you'll have to write the Coriolis stencil carefully to get either enstrophy or energy conserving results. For this option you can check out the docs on forcing functions:

https://clima.github.io/OceananigansDocumentation/stable/model_setup/forcing_functions/

and also our implementation of Coriolis stencils for spherical domains:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/Coriolis/hydrostatic_spherical_coriolis.jl

For example, you might need something like

using Oceananigans.Operators: Δx_vᶜᶠᵃ, Δy_uᶠᶜᵃ, Δxᶠᶜᵃ, Δyᶜᶠᵃ

const f0 = 1e-4
const Ly = 2*pi

@inline f(y) = f0 * sin(4*pi*y/Ly)

@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, v) = f(ynode(Face, j, grid)) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_vᶜᶠᵃ, v)
@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, u) = f(ynode(Face, j, grid)) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_uᶠᶜᵃ, u)

@inline crazy_x_f_cross_U(i, j, k, grid, fields) =
    @inbounds + ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, fields.v) / Δxᶠᶜᵃ(i, j, k, grid)

@inline crazy_y_f_cross_U(i, j, k, grid, fields) =
    @inbounds - ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, fields.u) / Δyᶜᶠᵃ(i, j, k, grid)

u_coriolis_forcing = Forcing(crazy_x_f_cross_U, discrete_form=true)
v_coriolis_forcing = Forcing(crazy_y_f_cross_U, discrete_form=true)

A second possibility is either to add a new coriolis implementation (something like VariableRectilinearCoriolis that takes in a user-defined coriolis parameter f(y)), or to generalize an existing one (like HydrostaticSphericalCoriolis).

[rbenamran]

Thanks for this great answer and the forcing code example. For now I will try to set up the Coriolis through discrete forcing and will look further into the spherical implementation.
A VariableRectilinearCoriolis would be an interesting feature suggestion for the future in my opinion!

[glwagner]

Make sure to check the signs - I think I changed them to be correct, but my code is misleading because the functions crazy_y_f_cross_U actually return the negative of y_f_cross_U (because forcing terms appear on the RHS of the equation, see equation 1 here: https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model).

Let's ponder a bit more what a native implementation might look like. The Coriolis module needs a bit of clean up anyways, since we should support either enstrophy or energy conserving stencils for FPlane and BetaPlane on rectilinear grids in addition to spherical grids. We haven't had a need to pursue this yet because we don't support horizontally-stretched rectilinear grids yet (where the stencil choice becomes relevant). Perhaps when we get around to that refactor, we can generalize the API so the "square planet" can be simulated. Thanks for the idea!