Add a tidal potential to the barotropic forcing

examples/near_global_ocean_simulation.jl

@@ -104,7 +104,10 @@ 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", "Φ"; architecture=GPU(), backend=InMemory(41))
+tidal_potential = FieldTimeSeries("tidal_potential.jld2", "Φ"; architecture=GPU(), backend=InMemory(41), boundary_conditions=FieldBoundaryConditions(tidal_potential.grid, (Center, Center, Nothing)))
+Oceananigans.BoundaryConditions.fill_halo_regions!(tidal_potential)
+atmosphere = JRA55PrescribedAtmosphere(arch; tidal_potential, backend=JRA55NetCDFBackend(41))
 
 # ## The coupled simulation
 
@@ -116,7 +119,7 @@ coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation)
 # We then create a coupled simulation. We start with a small-ish time step of 90 seconds.
 # We run the simulation for 10 days with this small-ish time step.
 
-simulation = Simulation(coupled_model; Δt=90, stop_time=10days)
+simulation = Simulation(coupled_model; Δt=60, stop_time=30days)
 
 # We define a callback function to monitor the simulation's progress,
 
@@ -175,7 +178,7 @@ run!(simulation)
 
 # After the initial spin up of 10 days, we can increase the time-step and run for longer.
 
-simulation.stop_time = 60days
+simulation.stop_time = 120days
 simulation.Δt = 10minutes
 run!(simulation)
 

examples/one_degree_simulation.jl

@@ -93,7 +93,7 @@ tracer_advection   = Centered()
 ocean = ocean_simulation(grid;
                          momentum_advection,
                          tracer_advection,
-                         closure,
+#                          closure,
                          forcing,
                          free_surface)
 

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, forcing_barotropic_potential
 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,32 @@ 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?
+    barotropic_potential = forcing_barotropic_potential(ocean)
+
+    launch!(arch, grid, kernel_parameters,
+            _compute_barotropic_potential!,
+            barotropic_potential,
+            grid,
+            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 +188,41 @@ end
     end
 end
 
+@inline interpolate_tidal_potential(::Nothing,       grid, args) = zero(grid)
+@inline interpolate_tidal_potential(tidal_potential, grid, args) = interp_atmos_time_series(tidal_potential, args...)
+
+# Fallback
+@kernel _compute_barotropic_potential!(::Nothing, args...) = nothing
+
+@kernel function _compute_barotropic_potential!(barotropic_potential,
+                                                grid,
+                                                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 = interpolate_tidal_potential(tidal_potential, grid, atmos_args)
+
+    @inbounds barotropic_potential[i, j, 1] = pa / ρₒ + Φt / ρₒ
+end
+
 @kernel function _interpolate_auxiliary_freshwater_flux!(freshwater_flux,
                                                          interface_grid,
                                                          clock,
@@ -198,6 +252,9 @@ end
 ##### Utility for interpolating tuples of fields
 #####
 
+# Assumption: a Nothing object interpolates to zero!!
+@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

@@ -3,6 +3,7 @@ module PrescribedAtmospheres
 using Oceananigans.Grids: grid_name
 using Oceananigans.Utils: prettysummary, Time
 using Oceananigans.Fields: Center
+using Oceananigans.BoundaryConditions: FieldBoundaryConditions
 using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series
 using Oceananigans.TimeSteppers: Clock, tick!
 
@@ -289,7 +290,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 +299,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
@@ -320,32 +322,37 @@ function Base.show(io::IO, pa::PrescribedAtmosphere)
 end
 
 function default_atmosphere_velocities(grid, times)
-    ua = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    va = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+    bcs = FieldBoundaryConditions(grid, (Center, Center, Nothing))
+    ua  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
+    va  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
     return (u=ua, v=va)
 end
 
 function default_atmosphere_tracers(grid, times)
-    Ta = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    qa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+    bcs = FieldBoundaryConditions(grid, (Center, Center, Nothing))
+    Ta  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
+    qa  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
     parent(Ta) .= 273.15 + 20
     return (T=Ta, q=qa)
 end
 
 function default_downwelling_radiation(grid, times)
-    Qℓ = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    Qs = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+    bcs = FieldBoundaryConditions(grid, (Center, Center, Nothing))
+    Qℓ  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
+    Qs  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
     return TwoBandDownwellingRadiation(shortwave=Qs, longwave=Qℓ)
 end
 
 function default_freshwater_flux(grid, times)
-    rain = FieldTimeSeries{Center, Center, Nothing}(grid, times)
-    snow = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+    bcs  = FieldBoundaryConditions(grid, (Center, Center, Nothing))
+    rain = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
+    snow = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
     return (; rain, snow)
 end
 
 function default_atmosphere_pressure(grid, times)
-    pa = FieldTimeSeries{Center, Center, Nothing}(grid, times)
+    bcs = FieldBoundaryConditions(grid, (Center, Center, Nothing))
+    pa  = FieldTimeSeries{Center, Center, Nothing}(grid, times; boundary_conditions=bcs)
     parent(pa) .= 101325
     return pa
 end
@@ -375,6 +382,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 +399,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 +419,7 @@ function PrescribedAtmosphere(grid, times;
                                 tracers,
                                 freshwater_flux,
                                 auxiliary_freshwater_flux,
+                                tidal_potential,
                                 downwelling_radiation,
                                 thermodynamics_parameters,
                                 times,

src/OceanSimulations/barotropic_potential_forcing.jl

@@ -24,6 +24,15 @@ const YDirectionBPF = BarotropicPotentialForcing{<:YDirection}
 forcing_barotropic_potential(something) = nothing
 forcing_barotropic_potential(f::BarotropicPotentialForcing) = f.potential.data
 
+function forcing_barotropic_potential(tf::Tuple)
+    n = findfirst(f -> f isa BarotropicPotentialForcing, tf)
+    if isnothing(n)
+        return nothing
+    else
+        return forcing_barotropic_potential(tf[n])
+    end
+end
+
 function forcing_barotropic_potential(mf::MultipleForcings)
     n = findfirst(f -> f isa BarotropicPotentialForcing, mf.forcing)
     if isnothing(n)