simplify initial exchange

NEWS.md

@@ -6,6 +6,14 @@ ClimaCoupler.jl Release Notes
 
 ### ClimaCoupler features
 
+#### Simplify initial component model exchange. PR[#1305](https://github.com/CliMA/ClimaCoupler.jl/pull/1305)
+
+Surface humidity is now computed from the atmosphere state and surface temperature,
+rather than sometimes computing it and sometimes retrieving it from a model cache.
+This allows us to simplify the initial component exchange, since we don't need
+to `step!` component models to get them to compute humidity. Without `step!`,
+we can also remove `reinit!`.
+
 #### Turbulent energy flux is split into LHF, SHF. PR[#1309](https://github.com/CliMA/ClimaCoupler.jl/pull/1309)
 Previously, we have exchanged the combined turbulent energy flux `F_turb_energy`;
 This PR splits up the exchanged quantity into `F_lh` and `F_sh`.

config/ci_configs/amip_target_topo_diagedmf_shortrun.yml

@@ -1,4 +1,5 @@
 FLOAT_TYPE: "Float32"
+albedo_model: "CouplerAlbedo"
 apply_limiter: false
 atmos_config_file: "config/longrun_configs/longrun_aquaplanet_allsky_diagedmf_0M.yml"
 co2: "maunaloa"

docs/src/interfacer.md

@@ -269,7 +269,6 @@ get_field(sim::AbstractSurfaceStub, ::Val{:beta}) = sim.cache.beta
 get_field(sim::AbstractSurfaceStub, ::Val{:roughness_buoyancy}) = sim.cache.z0b
 get_field(sim::AbstractSurfaceStub, ::Val{:roughness_momentum}) = sim.cache.z0m
 get_field(sim::AbstractSurfaceStub, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}}) = sim.cache.α
-get_field(sim::AbstractSurfaceStub, ::Val{:surface_humidity}) = TD.q_vap_saturation_generic.(sim.cache.thermo_params, sim.cache.T_sfc, sim.cache.ρ_sfc, sim.cache.phase)
 get_field(sim::AbstractSurfaceStub, ::Val{:surface_temperature}) = sim.cache.T_sfc
 ```
 and with the corresponding `update_field!` functions

experiments/ClimaEarth/components/atmosphere/climaatmos.jl

@@ -215,7 +215,7 @@ moisture_flux(::CA.DryModel, integrator) = StaticArrays.SVector(eltype(integrato
 moisture_flux(::Union{CA.EquilMoistModel, CA.NonEquilMoistModel}, integrator) =
     CC.Geometry.WVector.(integrator.p.precomputed.sfc_conditions.ρ_flux_q_tot)
 
-ρq_tot(::CA.DryModel, integrator) = StaticArrays.SVector(eltype(integrator.u)(0))
+ρq_tot(::CA.DryModel, integrator) = zeros(axes(integrator.u.c.uₕ))
 ρq_tot(::Union{CA.EquilMoistModel, CA.NonEquilMoistModel}, integrator) = integrator.u.c.ρq_tot
 
 # extensions required by the Interfacer
@@ -270,7 +270,8 @@ function Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:LW_d})
     )
 end
 Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:specific_humidity}) =
-    CC.Fields.level(sim.integrator.u.c.ρq_tot ./ sim.integrator.u.c.ρ, 1)
+    CC.Fields.level(ρq_tot(sim.integrator.p.atmos.moisture_model, sim.integrator) ./ sim.integrator.u.c.ρ, 1)
+Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:air_density}) = CC.Fields.level(sim.integrator.u.c.ρ, 1)
 Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:radiative_energy_flux_sfc}) =
     surface_radiation_flux(sim.integrator.p.atmos.radiation_mode, sim.integrator)
 Interfacer.get_field(sim::ClimaAtmosSimulation, ::Val{:snow_precipitation}) =

experiments/ClimaEarth/components/land/climaland_bucket.jl

@@ -250,8 +250,6 @@ Interfacer.get_field(sim::BucketSimulation, ::Val{:surface_direct_albedo}) =
     CL.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
 Interfacer.get_field(sim::BucketSimulation, ::Val{:surface_diffuse_albedo}) =
     CL.surface_albedo(sim.model, sim.integrator.u, sim.integrator.p)
-Interfacer.get_field(sim::BucketSimulation, ::Val{:surface_humidity}) =
-    CL.surface_specific_humidity(nothing, sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
 Interfacer.get_field(sim::BucketSimulation, ::Val{:surface_temperature}) =
     CL.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
 
@@ -336,16 +334,27 @@ end
 function FluxCalculator.surface_thermo_state(
     sim::BucketSimulation,
     thermo_params::TD.Parameters.ThermodynamicsParameters,
-    thermo_state_int,
+    atmos_sim::Interfacer.AtmosModelSimulation,
 )
-
     T_sfc = Interfacer.get_field(sim, Val(:surface_temperature))
+    FieldExchanger.compute_surface_humidity!(
+        sim.integrator.p.bucket.q_sfc,
+        Interfacer.get_field(atmos_sim, Val(:air_temperature)),
+        Interfacer.get_field(atmos_sim, Val(:specific_humidity)),
+        Interfacer.get_field(atmos_sim, Val(:air_density)),
+        T_sfc,
+        thermo_params,
+    )
     # Note that the surface air density, ρ_sfc, is computed using the atmospheric state at the first level and making ideal gas
     # and hydrostatic balance assumptions. The land model does not compute the surface air density so this is
     # a reasonable stand-in.
-    ρ_sfc = FluxCalculator.extrapolate_ρ_to_sfc.(thermo_params, thermo_state_int, T_sfc) # ideally the # calculate elsewhere, here just getter...
-    q_sfc = Interfacer.get_field(sim, Val(:surface_humidity)) # already calculated in rhs! (cache)
-    @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, q_sfc)
+    ρ_sfc =
+        FluxCalculator.extrapolate_ρ_to_sfc.(
+            thermo_params,
+            Interfacer.get_field(atmos_sim, Val(:thermo_state_int)),
+            T_sfc,
+        )
+    return @. TD.PhaseEquil_ρTq.(thermo_params, ρ_sfc, T_sfc, sim.integrator.p.bucket.q_sfc)
 end
 
 """

experiments/ClimaEarth/components/ocean/prescr_seaice.jl

@@ -155,7 +155,6 @@ Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) =
 Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
 Interfacer.get_field(sim::PrescribedIceSimulation, ::Union{Val{:surface_direct_albedo}, Val{:surface_diffuse_albedo}}) =
     sim.integrator.p.params.α
-Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
 Interfacer.get_field(sim::PrescribedIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
 
 """

experiments/ClimaEarth/components/ocean/slab_ocean.jl

@@ -111,7 +111,6 @@ Interfacer.get_field(sim::SlabOceanSimulation, ::Val{:roughness_buoyancy}) = sim
 Interfacer.get_field(sim::SlabOceanSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
 Interfacer.get_field(sim::SlabOceanSimulation, ::Val{:surface_direct_albedo}) = sim.integrator.p.α_direct
 Interfacer.get_field(sim::SlabOceanSimulation, ::Val{:surface_diffuse_albedo}) = sim.integrator.p.α_diffuse
-Interfacer.get_field(sim::SlabOceanSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
 Interfacer.get_field(sim::SlabOceanSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
 
 """

experiments/ClimaEarth/setup_run.jl

@@ -496,46 +496,34 @@ function CoupledSimulation(config_dict::AbstractDict)
         The concrete steps for proper initialization are:
         =#
 
-        # 1.coupler updates surface model area fractions
+        # 1. Coupler updates surface model area fractions
         FieldExchanger.update_surface_fractions!(cs)
 
-        # 2.surface density (`ρ_sfc`): calculated by the coupler by adiabatically extrapolating atmospheric thermal state to the surface.
-        # For this, we need to import surface and atmospheric fields. The model sims are then updated with the new surface density.
+        # 2. Import atmospheric and surface fields into the coupler fields, then broadcast them back out to all components.
         FieldExchanger.import_combined_surface_fields!(cs.fields, cs.model_sims)
         FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims)
         FieldExchanger.update_model_sims!(cs.model_sims, cs.fields)
 
-        # Set all initial cache values for the land model, now that we have updated drivers
+        # 3. Set all initial cache values for the land model, now that we have updated drivers
         land_set_initial_cache! = CL.make_set_initial_cache(cs.model_sims.land_sim.model)
         land_set_initial_cache!(
             cs.model_sims.land_sim.integrator.p,
             cs.model_sims.land_sim.integrator.u,
             cs.model_sims.land_sim.integrator.t,
         )
 
-        # 3.surface vapor specific humidity (`q_sfc`): step surface models with the new surface density to calculate their respective `q_sfc` internally
-        ## TODO: the q_sfc calculation follows the design of the bucket q_sfc, but it would be neater to abstract this from step! (#331)
-        Interfacer.step!(land_sim, tspan[1] + Δt_cpl)
-        Interfacer.step!(ocean_sim, tspan[1] + Δt_cpl)
-        Interfacer.step!(ice_sim, tspan[1] + Δt_cpl)
-
-        # 4.turbulent fluxes: now we have all information needed for calculating the initial
-        # turbulent surface fluxes
+        # 4. Compute radiative fluxes and update the coupler fields and model simulations with the new fluxes
+        # Any other callbacks that modify a model's cache should be called here as well.
+        if hasradiation(cs.model_sims.atmos_sim.integrator)
+            CA.rrtmgp_model_callback!(cs.model_sims.atmos_sim.integrator)
+            FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims)
+            FieldExchanger.update_model_sims!(cs.model_sims, cs.fields)
+        end
 
-        ## calculate turbulent fluxes in surface models and save the weighted average in coupler fields
+        # 5.turbulent fluxes: Now we have all information needed for calculating the initial
+        # turbulent surface fluxes. Calculate and update turbulent fluxes for each surface model,
+        # and save the weighted average in coupler fields
         FluxCalculator.turbulent_fluxes!(cs.model_sims, cs.fields, cs.boundary_space, cs.thermo_params)
-
-        # Updating only surface temperature because it is required by the RRTGMP callback (
-        # called below at reinit). Turbulent fluxes in atmos are updated in `update_model_sims`.
-        Interfacer.update_field!(atmos_sim, Val(:surface_temperature), cs.fields)
-
-
-        # 5.reinitialize models + radiative flux: prognostic states and time are set to their initial conditions. For atmos, this also triggers the callbacks and sets a nonzero radiation flux (given the new sfc_conditions)
-        FieldExchanger.reinit_model_sims!(cs.model_sims)
-
-        # 6.update all fluxes.
-        FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims)
-        FieldExchanger.update_model_sims!(cs.model_sims, cs.fields)
     end
     return cs
 end
@@ -681,29 +669,25 @@ function step!(cs::CoupledSimulation)
     !isnothing(cs.conservation_checks) && ConservationChecker.check_conservation!(cs)
     ClimaComms.barrier(comms_ctx)
 
-    ## update water albedo from wind at dt_water_albedo
-    ## (this will be extended to a radiation callback from the coupler)
-    TimeManager.maybe_trigger_callback(cs.callbacks.water_albedo, cs)
-
-    ## update the surface fractions for surface models,
-    ## and update all component model simulations with the current fluxes stored in the coupler
-    FieldExchanger.update_surface_fractions!(cs)
-    FieldExchanger.update_model_sims!(cs.model_sims, cs.fields)
-
     ## step component model simulations sequentially for one coupling timestep (Δt_cpl)
     FieldExchanger.step_model_sims!(cs.model_sims, cs.t[])
 
-    ## update the coupler with the new surface properties and calculate the turbulent fluxes
-    FieldExchanger.import_combined_surface_fields!(cs.fields, cs.model_sims) # i.e. T_sfc, surface_albedo, z0, beta
-    ## calculate turbulent fluxes in surfaces and save the weighted average in coupler fields
-    FluxCalculator.turbulent_fluxes!(cs.model_sims, cs.fields, cs.boundary_space, cs.thermo_params)
-
-    Interfacer.update_field!(atmos_sim, Val(:surface_temperature), cs.fields)
-
+    ## update the surface fractions for surface models
+    FieldExchanger.update_surface_fractions!(cs)
 
+    ## update the coupler with the new surface properties
+    FieldExchanger.import_combined_surface_fields!(cs.fields, cs.model_sims)
     ## update the coupler with the new atmospheric properties
-    FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims) # radiative and/or turbulent
+    FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims)
+    ## update the model simulations with the new coupler fields
+    FieldExchanger.update_model_sims!(cs.model_sims, cs.fields)
+
+    ## calculate turbulent fluxes in the coupler and update the model simulations with them
+    FluxCalculator.turbulent_fluxes!(cs.model_sims, cs.fields, cs.boundary_space, cs.thermo_params)
 
+    ## update water albedo from wind at dt_water_albedo
+    ## (this will be extended to a radiation callback from the coupler)
+    TimeManager.maybe_trigger_callback(cs.callbacks.water_albedo, cs)
     ## callback to checkpoint model state
     TimeManager.maybe_trigger_callback(cs.callbacks.checkpoint, cs)
 

experiments/ClimaEarth/user_io/debug_plots.jl

@@ -249,7 +249,7 @@ Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:snow_water_equiv}) = sim.i
 Interfacer.get_field(sim::ClimaLandSimulation, ::Val{:snow_liquid_water}) = sim.integrator.u.snow.S_l
 
 # currently selected plot fields
-plot_field_names(sim::Interfacer.SurfaceModelSimulation) = (:area_fraction, :surface_temperature, :surface_humidity)
+plot_field_names(sim::Interfacer.SurfaceModelSimulation) = (:area_fraction, :surface_temperature)
 plot_field_names(sim::ClimaLandSimulation) = (
     :area_fraction,
     :surface_direct_albedo,
@@ -264,5 +264,5 @@ plot_field_names(sim::ClimaLandSimulation) = (
     :snow_water_equiv,
     :snow_liquid_water,
 )
-plot_field_names(sim::BucketSimulation) = (:area_fraction, :surface_temperature, :surface_humidity, :σS, :Ws, :W)
+plot_field_names(sim::BucketSimulation) = (:area_fraction, :surface_temperature, :σS, :Ws, :W)
 plot_field_names(sim::ClimaAtmosSimulation) = (:w, :ρq_tot, :ρe_tot, :liquid_precipitation, :snow_precipitation)

src/FieldExchanger.jl

@@ -7,6 +7,8 @@ atmospheric and surface component models.
 module FieldExchanger
 
 import ..Interfacer, ..FluxCalculator, ..Utilities
+import Thermodynamics as TD
+import Thermodynamics.Parameters as TDP
 
 export import_atmos_fields!,
     update_surface_fractions!,
@@ -248,4 +250,33 @@ function combine_surfaces!(combined_field, sims, field_name)
     end
 end
 
+"""
+    compute_surface_humidity!(q_sfc, T_atmos, q_atmos, ρ_atmos, T_sfc, thermo_params)
+
+Compute the surface specific humidity based on the atmospheric state and surface temperature.
+The phase of the surface is determined by the surface temperature, and the saturation
+specific humidity is computed accordingly.
+All fields should be on the exchange grid.
+# Arguments
+- `q_sfc`: [CC.Fields.Field] output field for surface specific humidity.
+- `T_atmos`: [CC.Fields.Field] atmospheric temperature.
+- `q_atmos`: [CC.Fields.Field] atmospheric specific humidity.
+- `ρ_atmos`: [CC.Fields.Field] atmospheric air density.
+- `T_sfc`: [CC.Fields.Field] surface temperature.
+- `thermo_params`: [TD.Parameters.ThermodynamicsParameters] the thermodynamic parameters.
+"""
+function compute_surface_humidity!(q_sfc, T_atmos, q_atmos, ρ_atmos, T_sfc, thermo_params)
+    thermo_state_atmos = TD.PhaseEquil_ρTq.(thermo_params, ρ_atmos, T_atmos, q_atmos)
+    ρ_sfc = FluxCalculator.extrapolate_ρ_to_sfc.(thermo_params, thermo_state_atmos, T_sfc)
+
+    T_freeze = TDP.T_freeze(thermo_params)
+    @. q_sfc = ifelse(
+        T_sfc .> T_freeze,
+        TD.q_vap_saturation_generic(thermo_params, T_sfc, ρ_sfc, TD.Liquid()),
+        TD.q_vap_saturation_generic(thermo_params, T_sfc, ρ_sfc, TD.Ice()),
+    )
+    return nothing
+end
+
+
 end # module

src/FluxCalculator.jl

@@ -82,6 +82,11 @@ function turbulent_fluxes!(
         compute_surface_fluxes!(csf, sim, atmos_sim, boundary_space, thermo_params)
     end
 
+    # Update the atmosphere with the fluxes across all surface models
+    # The surface models have already been updated with the fluxes in `compute_surface_fluxes!`
+    # TODO this should be `update_turbulent_fluxes` to match the surface models
+    Interfacer.update_field!(atmos_sim, Val(:turbulent_fluxes), csf)
+
     # TODO: add allowable bounds here, check explicitly that all fluxes are equal
     return nothing
 end
@@ -128,7 +133,7 @@ end
 """
     surface_thermo_state(sim::Interfacer.SurfaceModelSimulation,
                          thermo_params::TD.Parameters.ThermodynamicsParameters,
-                         thermo_state_int)
+                         atmos_sim::Interfacer.AtmosModelSimulation)
 
 Return the surface thermo state the surface model simulation `sim`.
 
@@ -139,17 +144,22 @@ default, is computed assuming a liquid phase).
 function surface_thermo_state(
     sim::Interfacer.SurfaceModelSimulation,
     thermo_params::TD.Parameters.ThermodynamicsParameters,
-    thermo_state_int,
+    atmos_sim::Interfacer.AtmosModelSimulation,
 )
-    FT = eltype(parent(thermo_state_int))
+    FT = eltype(atmos_sim.integrator.u)
 
     T_sfc = Interfacer.get_field(sim, Val(:surface_temperature))
     # Note that the surface air density, ρ_sfc, is computed using the atmospheric state at the first level and making ideal gas
     # and hydrostatic balance assumptions. The land model does not compute the surface air density so this is
     # a reasonable stand-in.
     #
     # NOTE: This allocates! Fix me!
-    ρ_sfc = FluxCalculator.extrapolate_ρ_to_sfc.(thermo_params, thermo_state_int, T_sfc) # ideally the # calculate elsewhere, here just getter...
+    ρ_sfc =
+        FluxCalculator.extrapolate_ρ_to_sfc.(
+            thermo_params,
+            Interfacer.get_field(atmos_sim, Val(:thermo_state_int)),
+            T_sfc,
+        ) # ideally the # calculate elsewhere, here just getter...
 
     # For SurfaceStabs, this is just liquid phase
     q_sfc = TD.q_vap_saturation_generic.(thermo_params, T_sfc, ρ_sfc, TD.Liquid()) # default: saturated liquid surface
@@ -305,7 +315,7 @@ function compute_surface_fluxes!(
     # get area fraction (min = 0, max = 1)
     area_fraction = Interfacer.get_field(sim, Val(:area_fraction))
 
-    thermo_state_sfc = FluxCalculator.surface_thermo_state(sim, thermo_params, thermo_state_int)
+    thermo_state_sfc = FluxCalculator.surface_thermo_state(sim, thermo_params, atmos_sim)
 
     surface_params = FluxCalculator.get_surface_params(atmos_sim)
 

src/Interfacer.jl

@@ -215,7 +215,6 @@ get_field(
         Val{:roughness_momentum},
         Val{:surface_direct_albedo},
         Val{:surface_diffuse_albedo},
-        Val{:surface_humidity},
         Val{:surface_temperature},
     },
 ) = get_field_error(sim, val)

src/surface_stub.jl

@@ -40,8 +40,6 @@ get_field(sim::AbstractSurfaceStub, ::Val{:roughness_buoyancy}) = sim.cache.z0b
 get_field(sim::AbstractSurfaceStub, ::Val{:roughness_momentum}) = sim.cache.z0m
 get_field(sim::AbstractSurfaceStub, ::Val{:surface_direct_albedo}) = sim.cache.α_direct
 get_field(sim::AbstractSurfaceStub, ::Val{:surface_diffuse_albedo}) = sim.cache.α_diffuse
-get_field(sim::AbstractSurfaceStub, ::Val{:surface_humidity}) =
-    TD.q_vap_saturation_generic.(sim.cache.thermo_params, sim.cache.T_sfc, sim.cache.ρ_sfc, sim.cache.phase)
 get_field(sim::AbstractSurfaceStub, ::Val{:surface_temperature}) = sim.cache.T_sfc
 
 """

test/field_exchanger_tests.jl

@@ -43,11 +43,6 @@ Interfacer.get_field(sim::TestSurfaceSimulation1, ::Val{:area_fraction}) = sim.c
 Interfacer.get_field(sim::TestSurfaceSimulation2, ::Val{:area_fraction}) =
     sim.cache_field .* CC.Spaces.undertype(axes(sim.cache_field))(0.5)
 
-Interfacer.get_field(sim::Union{TestSurfaceSimulation1, TestSurfaceSimulation2}, ::Val{:surface_humidity}) =
-    sim.cache_field .* 0
-Interfacer.get_field(sim::Union{TestSurfaceSimulation2, TestSurfaceSimulation2}, ::Val{:surface_humidity}) =
-    sim.cache_field .* 0
-
 Interfacer.reinit!(::TestSurfaceSimulation1) = nothing
 Interfacer.step!(::TestSurfaceSimulation1, _) = nothing