Allow attaching a prognostic atmosphere take 2

src/OceanSeaIceModels/InterfaceComputations/atmosphere_ocean_fluxes.jl

@@ -1,7 +1,8 @@
 using Oceananigans.Operators: intrinsic_vector
 using Oceananigans.Grids: inactive_node
 using ClimaOcean.OceanSeaIceModels.PrescribedAtmospheres: thermodynamics_parameters, 
-                                                          reference_height
+                                                          reference_height,
+                                                          boundary_layer_height
 
 function compute_atmosphere_ocean_fluxes!(coupled_model)
     ocean = coupled_model.ocean
@@ -27,7 +28,7 @@ function compute_atmosphere_ocean_fluxes!(coupled_model)
                        Qs = atmosphere_fields.Qs.data,
                        Qℓ = atmosphere_fields.Qℓ.data,
                        Mp = atmosphere_fields.Mp.data,
-                       h_bℓ = atmosphere.boundary_layer_height)
+                       h_bℓ = boundary_layer_height(atmosphere))
 
     flux_formulation = coupled_model.interfaces.atmosphere_ocean_interface.flux_formulation
     interface_fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes

src/OceanSeaIceModels/InterfaceComputations/interpolate_atmospheric_state.jl

@@ -6,7 +6,6 @@ function interpolate_atmospheric_state!(coupled_model)
     ocean = coupled_model.ocean
     atmosphere = coupled_model.atmosphere
     atmosphere_grid = atmosphere.grid
-    sea_ice = coupled_model.sea_ice
 
     # Basic model properties
     grid = ocean.model.grid

src/OceanSeaIceModels/OceanSeaIceModels.jl

@@ -71,6 +71,10 @@ const NoAtmosphereModel = OceanSeaIceModel{<:Any, Nothing}
 compute_atmosphere_ocean_fluxes!(::NoAtmosphereModel) = nothing
 compute_atmosphere_sea_ice_fluxes!(::NoAtmosphereModel) = nothing
 
+const PrescribedAtmosphereModel = OceanSeaIceModel{<:Any, <:PrescribedAtmosphere}
+
+compute_net_atmosphere_fluxes!(::PrescribedAtmosphereModel) = nothing
+
 # "No sea ice" implementation
 const NoSeaIceModel = Union{OceanSeaIceModel{Nothing}, FreezingLimitedCoupledModel}
 

src/OceanSeaIceModels/PrescribedAtmospheres.jl

@@ -4,11 +4,12 @@ using Oceananigans.Grids: grid_name
 using Oceananigans.Utils: prettysummary
 using Oceananigans.Fields: Center
 using Oceananigans.OutputReaders: FieldTimeSeries, update_field_time_series!, extract_field_time_series
+using Oceananigans.TimeSteppers: tick!
 
 using Adapt
 using Thermodynamics.Parameters: AbstractThermodynamicsParameters
 
-import Oceananigans.Models: update_model_field_time_series!
+import Oceananigans.TimeSteppers: time_step!
 
 import Thermodynamics.Parameters:
     gas_constant,   #
@@ -287,8 +288,9 @@ const PATP = PrescribedAtmosphereThermodynamicsParameters
 ##### Prescribed atmosphere (as opposed to dynamically evolving / prognostic)
 #####
 
-struct PrescribedAtmosphere{FT, M, G, U, P, C, F, I, R, TP, TI}
+struct PrescribedAtmosphere{FT, M, G, C, U, P, C, F, I, R, TP, TI}
     grid :: G
+    clock :: Clock{C}
     metadata :: M
     velocities :: U
     pressure :: P
@@ -347,11 +349,26 @@ function default_atmosphere_pressure(grid, times)
     return pa
 end
 
+@inline function time_step!(atmos::PrescribedAtmosphere, Δt) 
+    tick!(atmos.clock, Δt)
+
+    time = Time(atmos.clock.time)
+    ftses = extract_field_time_series(atmos)
+
+    for fts in ftses
+        update_field_time_series!(fts, time)
+    end    
+
+    return nothing
+end
+
 @inline thermodynamics_parameters(atmos::PrescribedAtmosphere) = atmos.thermodynamics_parameters
 @inline reference_height(atmos::PrescribedAtmosphere) = atmos.reference_height
+@inline boundary_layer_height(atmos::PrescribedAtmosphere) = atmos.boundary_layer_height    
 
 """
     PrescribedAtmosphere(grid, times;
+                         clock = Clock{Float64}(time = 0),
                          metadata = nothing,
                          reference_height = 10, # meters
                          boundary_layer_height = 600 # meters,
@@ -367,6 +384,7 @@ Return a representation of a prescribed time-evolving atmospheric
 state with data given at `times`.
 """
 function PrescribedAtmosphere(grid, times;
+                              clock = Clock{Float64}(time = 0),
                               metadata = nothing,  
                               reference_height = convert(eltype(grid), 10),
                               boundary_layer_height = convert(eltype(grid), 600),
@@ -384,6 +402,7 @@ function PrescribedAtmosphere(grid, times;
     end
 
     return PrescribedAtmosphere(grid,
+                                clock,
                                 metadata,
                                 velocities,
                                 pressure,
@@ -397,17 +416,6 @@ function PrescribedAtmosphere(grid, times;
                                 convert(FT, boundary_layer_height))
 end
 
-update_model_field_time_series!(::Nothing, time) = nothing
-
-function update_model_field_time_series!(atmos::PrescribedAtmosphere, time)
-    ftses = extract_field_time_series(atmos)
-    for fts in ftses
-        update_field_time_series!(fts, time)
-    end
-
-    return nothing
-end
-
 struct TwoBandDownwellingRadiation{SW, LW}
     shortwave :: SW
     longwave :: LW

src/OceanSeaIceModels/ocean_sea_ice_model.jl

@@ -1,5 +1,4 @@
 using Oceananigans
-using Oceananigans.Models: update_model_field_time_series!
 using Oceananigans.TimeSteppers: Clock
 using Oceananigans: SeawaterBuoyancy
 

src/OceanSeaIceModels/time_step_ocean_sea_ice_model.jl

@@ -13,6 +13,7 @@ using Printf
 function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_tendencies=true)
     ocean = coupled_model.ocean
     sea_ice = coupled_model.sea_ice
+    atmosphere = coupled_model.atmosphere
 
     # Be paranoid and update state at iteration 0
     coupled_model.clock.iteration == 0 && update_state!(coupled_model, callbacks)
@@ -39,6 +40,9 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_
     ocean.Δt = Δt
     time_step!(ocean)
 
+    # Time step the atmosphere
+    time_step!(atmosphere, Δt)
+
     # TODO:
     # - Store fractional ice-free / ice-covered _time_ for more
     #   accurate flux computation?
@@ -49,16 +53,17 @@ function time_step!(coupled_model::OceanSeaIceModel, Δt; callbacks=[], compute_
 end
 
 function update_state!(coupled_model::OceanSeaIceModel, callbacks=[]; compute_tendencies=true)
-    time = Time(coupled_model.clock.time)
-    update_model_field_time_series!(coupled_model.atmosphere, time)
-    interpolate_atmospheric_state!(coupled_model)
+
+    # This function needs to be specialized to allow different atmospheric models
+    interpolate_atmospheric_state!(coupled_model) 
 
     # Compute interface states
     compute_atmosphere_ocean_fluxes!(coupled_model)
     compute_atmosphere_sea_ice_fluxes!(coupled_model)
     compute_sea_ice_ocean_fluxes!(coupled_model)
 
-    # compute_net_atmosphere_fluxes!(coupled_model)
+    # This function needs to be specialized to allow different atmospheric models
+    compute_net_atmosphere_fluxes!(coupled_model)
     compute_net_ocean_fluxes!(coupled_model)
     #compute_net_sea_ice_fluxes!(coupled_model)