Add a tidal potential to the barotropic forcing

examples/near_global_ocean_simulation.jl

@@ -104,7 +104,8 @@ radiation = Radiation(arch)
 # The number of snapshots that are loaded into memory is determined by
 # the `backend`. Here, we load 41 snapshots at a time into memory.
 
-atmosphere = JRA55PrescribedAtmosphere(arch; backend=JRA55NetCDFBackend(41))
+tidal_potential = FieldTimeSeries("tidal_potential.jld2", "Φ")
+atmosphere = JRA55PrescribedAtmosphere(arch; tidal_potential, backend=JRA55NetCDFBackend(41))
 
 # ## The coupled simulation
 

experiments/ocean_forced_by_tides.jl

@@ -0,0 +1,73 @@
+
+using ClimaOcean
+using Oceananigans
+using OrthogonalSphericalShellGrids
+using Oceananigans.Units
+using Printf
+
+import Oceananigans.BuoyancyFormulations: ρ′ 
+import ClimaOcean.OceanSeaIceModels: reference_density, heat_capacity
+
+struct ZeroBuoyancy{R}
+    reference_density :: R
+end
+
+@inline ρ′(i, j, k, grid, ::ZeroBuoyancy, T, S) = zero(grid)
+@inline reference_density(r::ZeroBuoyancy) = r.reference_density
+@inline heat_capacity(r::ZeroBuoyancy) = 3990.0
+
+# all passive tracers
+equation_of_state = ZeroBuoyancy(1026.0)
+
+# A barotropic ocean
+grid = TripolarGrid(size = (700, 400, 1), halo = (7, 7, 7), z = (-10, 0))
+bottom_height = regrid_bathymetry(grid)
+grid  = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)
+ocean = ocean_simulation(grid; closure=nothing, bottom_drag_coefficient=0, tracer_advection=nothing, equation_of_state)
+
+# A prescribed tidal forcing
+Φ = FieldTimeSeries("tidal_potential_jra55.jld2", "Φ")
+atmos = PrescribedAtmosphere(Φ.grid, Φ.times; tidal_potential=Φ)
+set!(ocean.model, T=first(atmos.tracers.T), S=35)
+
+# neutral radiation
+radiation=Radiation(ocean_albedo=1, ocean_emissivity=0)
+
+# No fluxes
+struct ZeroFluxes{N, B}
+    solver_stop_criteria :: N
+    bulk_velocity :: B
+end
+
+import ClimaOcean.OceanSeaIceModels.InterfaceComputations: compute_interface_state, ComponentInterfaces
+using ClimaOcean.OceanSeaIceModels.InterfaceComputations: zero_interface_state
+
+@inline compute_interface_state(::ZeroFluxes, args...) = zero_interface_state(Float64)
+
+interfaces = ComponentInterfaces(atmos, ocean, nothing;
+                                 atmosphere_ocean_flux_formulation=ZeroFluxes(nothing, ClimaOcean.OceanSeaIceModels.InterfaceComputations.WindVelocity()),
+                                 radiation)
+
+barotropic_earth_model = OceanSeaIceModel(ocean, nothing; atmosphere=atmos, radiation, interfaces)
+
+barotropic_earth = Simulation(barotropic_earth_model, Δt=20minutes, stop_time=2days)
+
+wall_time = Ref(0.0)
+
+function progress(barotropic_earth)
+    clock   = barotropic_earth.model.clock
+    u, v, w = barotropic_earth.model.ocean.model.velocities
+
+    maxu = maximum(abs, u)
+    maxv = maximum(abs, v)
+    maxw = maximum(abs, w)
+
+    @info @sprintf("Iteration: %d, time: %s, wall_time: %s, max(|u|, |v|, |w|): %.2e %.2e %.2e\n",
+                   clock.iteration, prettytime(clock.time), prettytime(1e-9 * (time_ns() - wall_time[])), maxu, maxv, maxw)
+
+    wall_time[] = time_ns()
+end
+
+add_callback!(barotropic_earth, progress, IterationInterval(100))
+
+run!(barotropic_earth)

src/DataWrangling/JRA55.jl

@@ -652,7 +652,8 @@ Return a `PrescribedAtmosphere` representing JRA55 reanalysis data.
 function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon();
                                    backend = nothing,
                                    time_indexing = Cyclical(),
-                                   reference_height = 10,  # meters
+                                   surface_layer_height = 10,  # meters
+                                   tidal_potential = nothing,
                                    include_rivers_and_icebergs = false,
                                    other_kw...)
 
@@ -701,7 +702,7 @@ function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon();
 
     FT = eltype(ua)
     boundary_layer_height = convert(FT, 600)
-    reference_height = convert(FT, reference_height)
+    surface_layer_height = convert(FT, surface_layer_height)
     thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT)
     grid = ua.grid
     metadata = JRA55Data()
@@ -714,10 +715,11 @@ function JRA55PrescribedAtmosphere(architecture::AA, time_indices=Colon();
                                 tracers,
                                 freshwater_flux,
                                 auxiliary_freshwater_flux,
+                                tidal_potential,
                                 downwelling_radiation,
                                 thermodynamics_parameters,
                                 times,
-                                reference_height,
+                                surface_layer_height,
                                 boundary_layer_height)
 end
 

src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl

@@ -1,18 +1,18 @@
 using Oceananigans.Operators: intrinsic_vector
-using Oceananigans.Grids: _node
+using Oceananigans.Grids: _node, AbstractGrid   
 using Oceananigans.Fields: FractionalIndices
 using Oceananigans.OutputReaders: TimeInterpolator
 
 using ...OceanSimulations: forcing_barotropic_potential
 
+using ClimaOcean.OceanSimulations: BarotropicPotentialForcing
 using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: PrescribedAtmosphere
 import ClimaOcean.OceanSeaIceModels: interpolate_atmosphere_state!
 
 # TODO: move to PrescribedAtmospheres
 """Interpolate the atmospheric state onto the ocean / sea-ice grid."""
 function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosphere, coupled_model)
     ocean = coupled_model.ocean
-    atmosphere_grid = atmosphere.grid
 
     # Basic model properties
     grid = ocean.model.grid
@@ -40,13 +40,12 @@ function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosph
 
     # Extract info for time-interpolation
     u = atmosphere.velocities.u # for example
-    atmosphere_times = u.times
     atmosphere_backend = u.backend
     atmosphere_time_indexing = u.time_indexing
 
-    atmosphere_fields = coupled_model.interfaces.exchanger.exchange_atmosphere_state
-    space_fractional_indices = coupled_model.interfaces.exchanger.atmosphere_exchanger
-    exchange_grid = coupled_model.interfaces.exchanger.exchange_grid
+    atmosphere_fields = interfaces.exchanger.exchange_atmosphere_state
+    space_fractional_indices = interfaces.exchanger.atmosphere_exchanger
+    exchange_grid = interfaces.exchanger.exchange_grid
 
     # Simplify NamedTuple to reduce parameter space consumption.
     # See https://github.com/CliMA/ClimaOcean.jl/issues/116.
@@ -106,12 +105,37 @@ function interpolate_atmosphere_state!(interfaces, atmosphere::PrescribedAtmosph
                 auxiliary_time_indexing)
     end
 
-    # Set ocean barotropic pressure forcing
-    barotropic_potential = forcing_barotropic_potential(ocean)
-    ρₒ = coupled_model.interfaces.ocean_properties.reference_density
-    if !isnothing(barotropic_potential)
-        parent(barotropic_potential) .= parent(atmosphere_data.p) ./ ρₒ
+    # barotropic potential in PrescribedAtmosphere?
+    atmosphere_pressure = atmosphere.pressure.data
+    tidal_potential = atmosphere.tidal_potential
+
+    if !isnothing(tidal_potential)
+        tidal_potential_data = tidal_potential.data
+    else
+        tidal_potential_data = nothing
+    end
+
+    # Which forcing is this going to be?
+    u_forcing = ocean.model.forcing.u
+    barotropic_potential = if u_forcing isa BarotropicPotentialForcing
+        u_forcing.potential
+    else
+        n = findfirst(x -> x isa BarotropicPotentialForcing, u_forcing)
+        u_forcing[n].potential
     end
+
+    launch!(arch, grid, kernel_parameters,
+            _compute_barotropic_potential!,
+            barotropic_potential,
+            space_fractional_indices,
+            time_interpolator,
+            atmosphere_pressure,
+            tidal_potential_data,
+            interfaces.ocean_properties.reference_density,
+            atmosphere_backend,
+            atmosphere_time_indexing)
+
+    return nothing
 end
 
 @inline get_fractional_index(i, j, ::Nothing) = nothing
@@ -169,6 +193,34 @@ end
     end
 end
 
+@kernel function _compute_barotropic_potential!(barotropic_potential,
+                                                space_fractional_indices,
+                                                time_interpolator,
+                                                atmos_pressure,
+                                                tidal_potential,
+                                                ocean_reference_density,
+                                                atmos_backend,
+                                                atmos_time_indexing)
+
+    i, j = @index(Global, NTuple)
+
+    ρₒ = ocean_reference_density
+
+    ii = space_fractional_indices.i
+    jj = space_fractional_indices.j
+    fi = get_fractional_index(i, j, ii)
+    fj = get_fractional_index(i, j, jj)
+
+    x_itp = FractionalIndices(fi, fj, nothing)
+    t_itp = time_interpolator
+    atmos_args = (x_itp, t_itp, atmos_backend, atmos_time_indexing)
+
+    pa = interp_atmos_time_series(atmos_pressure,  atmos_args...) # yes this is a re-interpolation
+    Φt = interp_atmos_time_series(tidal_potential, atmos_args...)
+
+    @inbounds barotropic_potential[i, j, 1] = pa / ρₒ + Φt
+end
+
 @kernel function _interpolate_auxiliary_freshwater_flux!(freshwater_flux,
                                                          interface_grid,
                                                          clock,
@@ -198,6 +250,10 @@ end
 ##### Utility for interpolating tuples of fields
 #####
 
+# Assumption: a Nothing object interpolates to zero!!
+@inline interp_atmos_time_series(::Nothing, x_itp::FractionalIndices,    args...) = zero(x_itp.i)
+@inline interp_atmos_time_series(::Nothing, X, time, grid::AbstractGrid, args...) = zero(grid)
+
 # Note: assumes loc = (c, c, nothing) (and the third location should not matter.)
 @inline interp_atmos_time_series(J::AbstractArray, x_itp::FractionalIndices, t_itp, args...) =
     interpolate(x_itp, t_itp, J, args...)

src/OceanSeaIceModels/PrescribedAtmospheres.jl

@@ -289,7 +289,7 @@ const PATP = PrescribedAtmosphereThermodynamicsParameters
 ##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic)
 #####
 
-mutable struct PrescribedAtmosphere{FT, M, G, T, U, P, C, F, I, R, TP, TI}
+mutable struct PrescribedAtmosphere{FT, M, G, T, U, P, C, F, I, ΦT, R, TP, TI}
     grid :: G
     clock :: Clock{T}
     metadata :: M
@@ -298,6 +298,7 @@ mutable struct PrescribedAtmosphere{FT, M, G, T, U, P, C, F, I, R, TP, TI}
     tracers :: C
     freshwater_flux :: F
     auxiliary_freshwater_flux :: I
+    tidal_potential :: ΦT # this really belongs elsewhere, but we put it here for now
     downwelling_radiation :: R
     thermodynamics_parameters :: TP
     times :: TI
@@ -375,6 +376,7 @@ end
                          boundary_layer_height = 600 # meters,
                          thermodynamics_parameters = PrescribedAtmosphereThermodynamicsParameters(FT),
                          auxiliary_freshwater_flux = nothing,
+                         tidal_potential = nothing,
                          velocities            = default_atmosphere_velocities(grid, times),
                          tracers               = default_atmosphere_tracers(grid, times),
                          pressure              = default_atmosphere_pressure(grid, times),
@@ -391,6 +393,7 @@ function PrescribedAtmosphere(grid, times;
                               boundary_layer_height = convert(eltype(grid), 600),
                               thermodynamics_parameters = nothing,
                               auxiliary_freshwater_flux = nothing,
+                              tidal_potential = nothing,
                               velocities            = default_atmosphere_velocities(grid, times),
                               tracers               = default_atmosphere_tracers(grid, times),
                               pressure              = default_atmosphere_pressure(grid, times),
@@ -410,6 +413,7 @@ function PrescribedAtmosphere(grid, times;
                                 tracers,
                                 freshwater_flux,
                                 auxiliary_freshwater_flux,
+                                tidal_potential,
                                 downwelling_radiation,
                                 thermodynamics_parameters,
                                 times,