Merge #329

329: Add coupler's (combined) `rho_sfc` and `q_sfc_sat` calculation r=LenkaNovak a=LenkaNovak



Co-authored-by: Valeria Barra <[email protected]>

.buildkite/pipeline.yml

@@ -189,7 +189,7 @@ steps:
         artifact_paths: "experiments/AMIP/modular/output/slabplanet/target_params_in_slab_artifacts/total_energy*.png"
 
       # - label: "Moist earth with slab surface - test: monin allsky sponge idealinsol infreq_dt_cpl"
-        # command: "julia --color=yes --project=experiments/AMIP/modular/ experiments/AMIP/modular/coupler_driver_modular.jl --FLOAT_TYPE Float64 --enable_threading true --coupled true --surface_setup PrescribedSurface --moist equil --vert_diff true --rad allskywithclear --rayleigh_sponge true --alpha_rayleigh_uh 0 --alpha_rayleigh_w 10 --energy_check true --mode_name slabplanet --t_end 10days --dt_save_to_sol 3600secs --dt_cpl 21600 --dt 200secs --dt_rad 6hours --mono_surface true --h_elem 4 --precip_model 0M --run_name target_params_in_slab_test1 --job_id target_params_in_slab_test1" # Unconverged SF (reproduced locally)
+        # command: "julia --color=yes --project=experiments/AMIP/modular/ experiments/AMIP/modular/coupler_driver_modular.jl --FLOAT_TYPE Float64 --enable_threading true --coupled true --surface_setup PrescribedSurface --moist equil --vert_diff true --rad allskywithclear --rayleigh_sponge true --alpha_rayleigh_uh 0 --alpha_rayleigh_w 10 --energy_check true --mode_name slabplanet --t_end 10days --dt_save_to_sol 3600secs --dt_cpl 21600 --dt 200secs --dt_rad 6hours --mono_surface true --h_elem 4 --precip_model 0M --run_name target_params_in_slab_test1 --job_id target_params_in_slab_test1" # Unconverged SF (reproduced locally); works with 200s dt_cpl
         # artifact_paths: "experiments/AMIP/modular/output/slabplanet/target_params_in_slab_test1_artifacts/total_energy*.png"
 
       - label: "Moist earth with slab surface - test: bulk allsky sponge realinsol infreq_dt_cpl"

experiments/AMIP/modular/components/atmosphere/climaatmos_init.jl

@@ -3,6 +3,7 @@ import ClimaAtmos
 using ClimaAtmos: RRTMGPI
 import ClimaCoupler.FluxCalculator: compute_atmos_turbulent_fluxes!
 using ClimaCore: Fields.level, Geometry
+import ClimaCoupler.FieldExchanger: get_thermo_params
 
 # the clima atmos `integrator` is now defined
 struct ClimaAtmosSimulation{P, Y, D, I} <: AtmosModelSimulation
@@ -81,22 +82,23 @@ https://clima.github.io/SurfaceFluxes.jl/dev/SurfaceFluxes/#Monin-Obukhov-Simila
 """
 struct CoupledMoninObukhov end
 """
-    coupler_surface_setup(::CoupledMoninObukhov, p, csf_sfc = (; T = nothing, z0m = nothing, z0b = nothing, beta = nothing))
+    coupler_surface_setup(::CoupledMoninObukhov, p, csf_sfc = (; T = nothing, z0m = nothing, z0b = nothing, beta = nothing, q_vap = nothing))
 
 Sets up `surface_setup` as a `Fields.Field` of `SurfaceState`s.
 """
 function coupler_surface_setup(
     ::CoupledMoninObukhov,
     p;
-    csf_sfc = (; T = nothing, z0m = nothing, z0b = nothing, beta = nothing),
+    csf_sfc = (; T = nothing, z0m = nothing, z0b = nothing, beta = nothing, q_vap = nothing),
 )
 
-    surface_state(z0m, z0b, T, beta) = ClimaAtmos.SurfaceConditions.SurfaceState(;
+    surface_state(z0m, z0b, T, beta, q_vap) = ClimaAtmos.SurfaceConditions.SurfaceState(;
         parameterization = ClimaAtmos.SurfaceConditions.MoninObukhov(; z0m, z0b),
         T,
         beta,
+        q_vap,
     )
-    surface_state_field = @. surface_state(csf_sfc.z0m, csf_sfc.z0b, csf_sfc.T, csf_sfc.beta)
+    surface_state_field = @. surface_state(csf_sfc.z0m, csf_sfc.z0b, csf_sfc.T, csf_sfc.beta, csf_sfc.q_vap)
     return surface_state_field
 end
 
@@ -113,10 +115,24 @@ function compute_atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf)
 
     p = atmos_sim.integrator.p
 
+    thermo_params = get_thermo_params(atmos_sim)
+    ts_int = get_field(atmos_sim, Val(:thermo_state_int))
+
+    # surface density is needed for q_sat and requires atmos and sfc states, so it is calculated and saved in the coupler
+    parent(csf.ρ_sfc) .=
+        parent(extrapolate_ρ_to_sfc.(Ref(thermo_params), ts_int, swap_space!(ones(axes(ts_int.ρ)), csf.T_S)))
+
+    # Surface-specific vapor specific humidity, q_sfc, is calculated as a bulk quantity here, assuming a saturated liquid surface.
+    # (atmos has been doing this until this PR, so here we're moving that calculation to the coupler)
+    # The current method is ok for the time being because q_sfc is a passive variable in Land and nonexistent in the slab components
+    # TODO (next PR) - use land's q_sfc
+    # To approximately account for undersaturation of the surface, we can use beta to adjust evaporation to appoximate undersaturation.
+    q_sfc = TD.q_vap_saturation_generic.(thermo_params, csf.T_S, csf.ρ_sfc, TD.Liquid())
+
     coupler_sfc_setup = coupler_surface_setup(
         CoupledMoninObukhov(),
         p;
-        csf_sfc = (; T = csf.T_S, z0m = csf.z0m_S, z0b = csf.z0b_S, beta = csf.beta),
+        csf_sfc = (; T = csf.T_S, z0m = csf.z0m_S, z0b = csf.z0b_S, beta = csf.beta, q_vap = q_sfc),
     )
 
     p_names = propertynames(p)
@@ -126,13 +142,21 @@ function compute_atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf)
 
     ClimaAtmos.SurfaceConditions.update_surface_conditions!(atmos_sim.integrator.u, new_p, atmos_sim.integrator.t)
 
-    # p.sfc_conditions.ρ_flux_h_tot .= new_p.sfc_conditions.ρ_flux_h_tot
-    # p.sfc_conditions.ρ_flux_q_tot .= new_p.sfc_conditions.ρ_flux_q_tot
-    # p.sfc_conditions.ρ_flux_uₕ .= new_p.sfc_conditions.ρ_flux_uₕ
-    # p.sfc_conditions.ts .= new_p.sfc_conditions.ts
-    # p.sfc_conditions.buoyancy_flux .= new_p.sfc_conditions.buoyancy_flux
-    # p.sfc_conditions.obukhov_length .= new_p.sfc_conditions.obukhov_length
-
     p.sfc_conditions .= new_p.sfc_conditions
 
 end
+
+get_thermo_params(sim::ClimaAtmosSimulation) = CAP.thermodynamics_params(sim.integrator.p.params)
+get_field(sim::ClimaAtmosSimulation, ::Val{:thermo_state_int}) = Spaces.level(sim.integrator.p.ᶜts, 1)
+
+"""
+    extrapolate_ρ_to_sfc(thermo_params, ts_int, T_sfc)
+
+Uses the ideal gas law and hydrostatic balance to extrapolate for surface density.
+"""
+function extrapolate_ρ_to_sfc(thermo_params, ts_in, T_sfc)
+    T_int = TD.air_temperature(thermo_params, ts_in)
+    Rm_int = TD.gas_constant_air(thermo_params, ts_in)
+    ρ_air = TD.air_density(thermo_params, ts_in)
+    ρ_air * (T_sfc / T_int)^(TD.cv_m(thermo_params, ts_in) / Rm_int)
+end

experiments/AMIP/modular/components/land/bucket_utils.jl

@@ -1,83 +1,11 @@
-"""
-    get_land_temp(slab_sim::BucketSimulation)
-
-Returns the surface temperature of the earth;
-a method for the bucket model
-when used as the land model.
-"""
-function get_land_temp(slab_sim::BucketSimulation)
-    return ClimaLSM.surface_temperature(
-        slab_sim.model,
-        slab_sim.integrator.u,
-        slab_sim.integrator.p,
-        slab_sim.integrator.t,
-    )
-end
-
-
-"""
-    get_land_temp_from_state(land_sim, u)
-
-Returns the surface temperature of the earth, computed
-from the state u.
-"""
-function get_land_temp_from_state(land_sim, u)
-    return ClimaLSM.surface_temperature(land_sim.model, u, land_sim.integrator.p, land_sim.integrator.t)
-end
-
-"""
-    get_land_roughness(slab_sim::BucketSimulation)
-
-Returns the roughness length parameters of the bucket;
-a method for the bucket model
-when used as the land model.
-"""
-function get_land_roughness(slab_sim::BucketSimulation)
-    return slab_sim.model.parameters.z_0m, slab_sim.model.parameters.z_0b
-end
-
-"""
-   land_albedo(slab_sim::BucketSimulation)
-
-Returns the surface albedo of the earth;
-a method for the bucket model
-when used as the land model.
-"""
-function land_albedo(slab_sim::BucketSimulation)
-    return ClimaLSM.surface_albedo(slab_sim.model, slab_sim.integrator.u, slab_sim.integrator.p)
-end
-
-"""
-   land_beta(slab_sim::BucketSimulation)
-
-Returns the surface evaporative scaling factor over land;
-a method for the bucket model when used as the land model.
-Note that this is slightly different from the coupler's β,
-which includes the scaling factor over non-land surfaces.
-"""
-function land_beta(slab_sim::BucketSimulation)
-    return ClimaLSM.surface_evaporative_scaling(slab_sim.model, slab_sim.integrator.u, slab_sim.integrator.p)
-end
-
-
-"""
-    get_land_q(slab_sim::Bucketimulation, _...)
-
-Returns the surface specific humidity of the earth;
-a method for the bucket
-when used as the land model.
-"""
-function get_land_q(slab_sim::BucketSimulation, _...)
-    return ClimaLSM.surface_specific_humidity(slab_sim.model, slab_sim.integrator.u, slab_sim.integrator.p)
-end
 
 """
     get_bucket_energy(bucket_sim)
 
 Returns the volumetric internal energy of the bucket land model.
 """
 function get_land_energy(bucket_sim::BucketSimulation, e_per_area)
-
+    # required by ConservationChecker
     e_per_area .= zeros(axes(bucket_sim.integrator.u.bucket.W))
     soil_depth = FT = eltype(bucket_sim.integrator.u.bucket.W)
     ClimaCore.Fields.bycolumn(axes(bucket_sim.integrator.u.bucket.T)) do colidx
@@ -92,6 +20,15 @@ function get_land_energy(bucket_sim::BucketSimulation, e_per_area)
     return e_per_area
 end
 
+"""
+    get_land_temp_from_state(land_sim, u)
+Returns the surface temperature of the earth, computed
+from the state u.
+"""
+function get_land_temp_from_state(land_sim, u)
+    # required by viz_explorer.jl
+    return ClimaLSM.surface_temperature(land_sim.model, u, land_sim.integrator.p, land_sim.integrator.t)
+end
 
 """
     make_lsm_domain(
@@ -144,7 +81,8 @@ function make_lsm_domain(
 end
 
 # required by Interfacer
-get_field(sim::BucketSimulation, ::Val{:surface_temperature}) = sim.integrator.p.bucket.T_sfc
+get_field(sim::BucketSimulation, ::Val{:surface_temperature}) =
+    ClimaLSM.surface_temperature(sim.model, sim.integrator.u, sim.integrator.p, sim.integrator.t)
 get_field(sim::BucketSimulation, ::Val{:roughness_momentum}) = sim.model.parameters.z_0m
 get_field(sim::BucketSimulation, ::Val{:roughness_buoyancy}) = sim.model.parameters.z_0b
 get_field(sim::BucketSimulation, ::Val{:beta}) =
@@ -157,9 +95,6 @@ get_field(sim::BucketSimulation, ::Val{:area_fraction}) = sim.area_fraction
 # The surface air density 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.
-function update_field!(sim::BucketSimulation, ::Val{:air_density}, field)
-    parent(sim.integrator.p.bucket.ρ_sfc) .= parent(field)
-end
 
 function update_field!(sim::BucketSimulation, ::Val{:turbulent_energy_flux}, field)
     parent(sim.integrator.p.bucket.turbulent_energy_flux) .= parent(field)

experiments/AMIP/modular/coupler_driver_modular.jl

@@ -382,16 +382,24 @@ cs = CoupledSimulation{FT}(
 =#
 
 ## share states and fluxes between models
+# 1) import states into the coupler
 turbulent_fluxes = CombinedAtmosGrid()
 update_surface_fractions!(cs)
-import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes)
-update_sim!(cs.model_sims.atmos_sim, cs.fields, turbulent_fluxes) # would be good to rm dep in cs
-parsed_args["ode_algo"] == "ARS343" ? step!(atmos_sim, Δt_cpl) : nothing
-compute_combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes) # here computed using atmos functions
+import_combined_surface_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta
+# 2) import states into the atmos model
+update_sim!(cs.model_sims.atmos_sim, cs.fields, turbulent_fluxes)
+# 3) initiate atmos states
+parsed_args["ode_algo"] == "ARS343" ? step!(atmos_sim, Δt_cpl) : nothing # TODO: this should be in the integrator init
+# 4) calculate fluxes and update `sfc_conditions` (needed for radiation) and save in atmos cache
+compute_combined_turbulent_fluxes!(cs.model_sims, cs.fields, turbulent_fluxes)
+# 5) reset atmos model's states and step with the correct `sfc_conditions` to caclucate radiative fluxes
+reinit!(atmos_sim)
+parsed_args["ode_algo"] == "ARS343" ? step!(atmos_sim, Δt_cpl) : nothing # calcualte F_rad from surface conditions
+# 6) export fluxes from atmos model
 import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, turbulent_fluxes)
+# 7) update surface models
 update_model_sims!(cs.model_sims, cs.fields, turbulent_fluxes)
-
-## reinitialize (TODO: avoid with interfaces)
+# 8) reinitialize all models
 reinit_model_sims!(cs.model_sims)
 
 #=