Update coupler stepping logic so that fluxes are exchanged as often as possible

config/ci_configs/cmip_oceananigans_climaseaice.yml

@@ -3,7 +3,7 @@ albedo_model: "CouplerAlbedo"
 atmos_config_file: "config/atmos_configs/climaatmos_edonly.yml"
 coupler_toml: ["toml/amip_edonly.toml"]
 dt_atmos: "120secs" # 2 minutes
-dt_cpl: "1800secs" # 30 minutes
+dt_cpl: "360secs" # 6 minutes
 dt_land: "360secs" # 6 minutes
 dt_ocean: "1800secs" # 30 minutes
 dt_seaice: "1800secs" # 30 minutes

config/ci_configs/cmip_oceananigans_climaseaice_bucket.yml

@@ -4,7 +4,7 @@ atmos_config_file: "config/atmos_configs/climaatmos_edonly.yml"
 bucket_albedo_type: "map_temporal"
 coupler_toml: ["toml/amip_edonly.toml"]
 dt_atmos: "120secs" # 2 minutes
-dt_cpl: "1800secs" # 30 minutes
+dt_cpl: "360secs" # 6 minutes
 dt_land: "360secs" # 6 minutes
 dt_ocean: "1800secs" # 30 minutes
 dt_seaice: "1800secs" # 30 minutes

ext/ClimaCouplerCMIPExt/clima_seaice.jl

@@ -23,12 +23,14 @@ It contains the following objects:
 - `ice_properties::IP`: A NamedTuple of sea ice properties, including melting speed, Stefan-Boltzmann constant,
     and the Celsius to Kelvin conversion constant.
 """
-struct ClimaSeaIceSimulation{SIM, A, REMAP, NT, IP} <: Interfacer.AbstractSeaIceSimulation
+struct ClimaSeaIceSimulation{SIM, A, REMAP, NT, IP, MDT} <:
+       Interfacer.AbstractSeaIceSimulation
     ice::SIM
     area_fraction::A
     remapping::REMAP
     ocean_ice_interface::NT
     ice_properties::IP
+    model_Δt::MDT
 end
 
 """
@@ -130,6 +132,9 @@ function ClimaSeaIceSimulation(
     interface_temperature = OC.Field{OC.Center, OC.Center, Nothing}(grid)
     interface_salinity = OC.Field{OC.Center, OC.Center, Nothing}(grid)
 
+    # Initialize model_Δt so that time stepping works properly
+    model_Δt = float(dt)
+
     # Initialize nonzero sea ice if start date provided
     if !isnothing(start_date)
         sic_metadata = CO.DataWrangling.Metadatum(
@@ -207,6 +212,7 @@ function ClimaSeaIceSimulation(
         remapping,
         ocean_ice_interface,
         ice_properties,
+        model_Δt,
     )
 
     # Ensure ocean temperature is above freezing where there is sea ice
@@ -270,8 +276,12 @@ end
 ###############################################################################
 
 # Timestep the simulation forward to time `t`
-Interfacer.step!(sim::ClimaSeaIceSimulation, t) =
-    OC.time_step!(sim.ice, float(t) - sim.ice.model.clock.time)
+function Interfacer.step!(sim::ClimaSeaIceSimulation, t)
+    Δt = float(t) - sim.ice.model.clock.time
+    if isapprox(Δt, sim.model_Δt) || Δt > sim.model_Δt
+        OC.time_step!(sim.ice, Δt)
+    end
+end
 
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:area_fraction}) = sim.area_fraction
 Interfacer.get_field(sim::ClimaSeaIceSimulation, ::Val{:ice_concentration}) =

ext/ClimaCouplerCMIPExt/oceananigans.jl

@@ -18,13 +18,14 @@ It contains the following objects:
 - `remapping::REMAP`: Objects needed to remap from the exchange (spectral) grid to Oceananigans spaces.
 - `ice_concentration::SIC`: An Oceananigans Field representing the sea ice concentration on the ocean/sea ice grid.
 """
-struct OceananigansSimulation{SIM, A, OPROP, REMAP, SIC} <:
+struct OceananigansSimulation{SIM, A, OPROP, REMAP, SIC, MDT} <:
        Interfacer.AbstractOceanSimulation
     ocean::SIM
     area_fraction::A
     ocean_properties::OPROP
     remapping::REMAP
     ice_concentration::SIC
+    model_Δt::MDT
 end
 
 """
@@ -67,7 +68,7 @@ function OceananigansSimulation(
     tspan,
     output_dir,
     simple_ocean = false,
-    dt = nothing,
+    dt = 1800.0, # 30 minutes
     comms_ctx = ClimaComms.context(),
     coupled_param_dict = CP.create_toml_dict(FT),
     extra_kwargs...,
@@ -156,10 +157,10 @@ function OceananigansSimulation(
         closure = (horizontal_viscosity, vertical_mixing)
     end
 
-    Δt = isnothing(dt) ? CO.OceanSimulations.estimate_maximum_Δt(grid) : float(dt)
+    model_Δt = float(dt)
     ocean = CO.ocean_simulation(
         grid;
-        Δt,
+        Δt = model_Δt,
         timestepper = :SplitRungeKutta3,
         momentum_advection,
         tracer_advection,
@@ -287,6 +288,7 @@ function OceananigansSimulation(
         ocean_properties,
         remapping,
         ice_concentration,
+        model_Δt,
     )
 end
 
@@ -335,9 +337,13 @@ end
 ### Functions required by ClimaCoupler.jl for a AbstractSurfaceSimulation
 ###############################################################################
 
-# Timestep the simulation forward to time `t`
-Interfacer.step!(sim::OceananigansSimulation, t) =
-    OC.time_step!(sim.ocean, float(t) - sim.ocean.model.clock.time)
+# Timestep the simulation forward to time `t`. This may not actually do anything.
+function Interfacer.step!(sim::OceananigansSimulation, t)
+    Δt = float(t) - sim.ocean.model.clock.time
+    if isapprox(Δt, sim.model_Δt) || Δt > sim.model_Δt
+        OC.time_step!(sim.ocean, Δt)
+    end
+end
 
 Interfacer.get_field(sim::OceananigansSimulation, ::Val{:area_fraction}) = sim.area_fraction
 

ext/ClimaCouplerClimaAtmosExt.jl

@@ -523,11 +523,23 @@ function FluxCalculator.update_turbulent_fluxes!(sim::ClimaAtmosSimulation, fiel
 end
 
 # extensions required by FieldExchanger
-Interfacer.step!(sim::ClimaAtmosSimulation, t::Real) =
-    Interfacer.step!(sim.integrator, t - sim.integrator.t, true)
+function Interfacer.step!(sim::ClimaAtmosSimulation, t::Real)
+    Δt = t - sim.integrator.t
+    if Δt < sim.integrator.dt   # don't step if we haven't reached a step boundary
+        # (can happen if the coupler dt is less than this model's)
+        return nothing
+    end
+    Interfacer.step!(sim.integrator, Δt, true)
+end
+
 function Interfacer.step!(sim::ClimaAtmosSimulation, t::ITime)
-    while sim.integrator.t < t
-        Interfacer.step!(sim.integrator)
+    # Don't step if we haven't reached a step boundary
+    # (This can happen if the coupler dt is less than this model's)
+    Δt = float(t) - float(sim.integrator.t)
+    if isapprox(Δt, sim.integrator.dt) || Δt > sim.integrator.dt
+        while sim.integrator.t < t
+            Interfacer.step!(sim.integrator)
+        end
     end
     return nothing
 end

ext/ClimaCouplerClimaLandExt/climaland_bucket.jl

@@ -316,10 +316,14 @@ function Interfacer.update_field!(
 end
 
 function Interfacer.step!(sim::BucketSimulation, t)
-    while float(sim.integrator.t) < float(t)
-        Interfacer.step!(sim.integrator)
+    # Don't step if we haven't reached a step boundary
+    # (This can happen if the coupler dt is less than this model's)
+    Δt = float(t) - float(sim.integrator.t)
+    if isapprox(Δt, sim.integrator.dt) || Δt > sim.integrator.dt
+        while float(sim.integrator.t) < float(t)
+            Interfacer.step!(sim.integrator)
+        end
     end
-    return nothing
 end
 Interfacer.close_output_writers(sim::BucketSimulation) =
     isnothing(sim.output_writer) || close(sim.output_writer)

ext/ClimaCouplerClimaLandExt/climaland_integrated.jl

@@ -381,8 +381,13 @@ function Interfacer.update_field!(
 end
 
 function Interfacer.step!(sim::ClimaLandSimulation, t)
-    while float(sim.integrator.t) < float(t)
-        Interfacer.step!(sim.integrator)
+    # Don't step if we haven't reached a step boundary
+    # (This can happen if the coupler dt is less than this model's)
+    Δt = float(t) - float(sim.integrator.t)
+    if isapprox(Δt, sim.integrator.dt) || Δt > sim.integrator.dt
+        while sim.integrator.t < t
+            Interfacer.step!(sim.integrator)
+        end
     end
     return nothing
 end

src/Input.jl

@@ -636,7 +636,13 @@ function parse_component_dts!(config_dict)
         end
         for key in component_dt_names
             component_dt = Float64(Utilities.time_to_seconds(config_dict[key]))
-            @assert isapprox(Δt_cpl % component_dt, 0.0) "Coupler dt must be divisible by all component dt's\n dt_cpl = $Δt_cpl\n $key = $component_dt"
+            if key == "dt_atmos"
+                # ensure that the coupler dt is an integer multiple of the atmos dt
+                @assert isapprox(Δt_cpl % component_dt, 0.0) "Coupler time step must be an integer multiple of the atmos dt\n dt_cpl = $Δt_cpl\n $key = $component_dt"
+            else
+                # all other (surface) model dts must be divisible by the coupler dt
+                @assert isapprox(component_dt % Δt_cpl, 0.0) "All surface component dts must be divisible by the coupler dt\n $key = $component_dt\n dt_cpl = $Δt_cpl"
+            end
             component_dt_dict[key] = component_dt
         end
     else