Merge #337

337: Enable `q_sfc` calculation from each surface model r=LenkaNovak a=LenkaNovak



Co-authored-by: LenkaNovak <[email protected]>

Project.toml

@@ -22,6 +22,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f"
 TempestRemap_jll = "8573a8c5-1df0-515e-a024-abad257ee284"
 TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
+Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]

experiments/AMIP/modular/Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "fcf142f7db31edf6d14176a318e75e653d1a7ab7"
+project_hash = "6450a469b8b81ae10b8bee0bf984f08dba4cd63e"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "dcfdf328328f2645531c4ddebf841228aef74130"
@@ -235,7 +235,7 @@ uuid = "d934ef94-cdd4-4710-83d6-720549644b70"
 version = "0.3.9"
 
 [[deps.ClimaCoupler]]
-deps = ["ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "UnPack"]
+deps = ["ClimaAtmos", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "ClimaLSM", "Dates", "DocStringExtensions", "Insolation", "JLD2", "NCDatasets", "OrdinaryDiffEq", "Plots", "PrettyTables", "SciMLBase", "Statistics", "SurfaceFluxes", "TempestRemap_jll", "TerminalLoggers", "Thermodynamics", "UnPack"]
 path = "../../.."
 uuid = "4ade58fe-a8da-486c-bd89-46df092ec0c7"
 version = "0.1.0"

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

@@ -1,7 +1,7 @@
 # atmos_init: for ClimaAtmos pre-AMIP interface
 import ClimaAtmos
 using ClimaAtmos: RRTMGPI
-import ClimaCoupler.FluxCalculator: compute_atmos_turbulent_fluxes!
+import ClimaCoupler.FluxCalculator: compute_atmos_turbulent_fluxes!, calculate_surface_air_density
 using ClimaCore: Fields.level, Geometry
 import ClimaCoupler.FieldExchanger: get_thermo_params
 
@@ -29,6 +29,9 @@ function atmos_init(::Type{FT}, parsed_args::Dict) where {FT}
     @. integrator.p.sfc_conditions.ρ_flux_q_tot = Geometry.Covariant3Vector(FT(0.0))
     @. integrator.p.sfc_conditions.ρ_flux_uₕ.components = zeros(axes(integrator.p.sfc_conditions.ρ_flux_uₕ.components))
     parent(integrator.p.ᶠradiation_flux) .= parent(zeros(axes(integrator.p.ᶠradiation_flux)))
+    integrator.p.col_integrated_rain .= FT(0)
+    integrator.p.col_integrated_snow .= FT(0)
+
 
     ClimaAtmosSimulation(integrator.p.params, Y, spaces, integrator)
 end
@@ -45,6 +48,9 @@ get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_energy_flux}) =
 get_field(sim::ClimaAtmosSimulation, ::Val{:turbulent_moisture_flux}) =
     Geometry.WVector.(sim.integrator.p.sfc_conditions.ρ_flux_q_tot)
 
+get_field(sim::ClimaAtmosSimulation, ::Val{:thermo_state_int}) = Spaces.level(sim.integrator.p.ᶜts, 1)
+
+
 function update_field!(sim::ClimaAtmosSimulation, ::Val{:surface_temperature}, field)
     sim.integrator.p.radiation_model.surface_temperature .= RRTMGPI.field2array(field)
 end
@@ -111,48 +117,51 @@ in it. We do not want `new_p` to live in the atmospheric model permanently, beca
 trigger flux calculation during Atmos `step!`. We only want to trigger this once per coupling
 timestep from ClimaCoupler.
 """
+
 function compute_atmos_turbulent_fluxes!(atmos_sim::ClimaAtmosSimulation, csf)
 
-    p = atmos_sim.integrator.p
+    function modified_atmos_p(atmos_sim, csf_sfc)
 
-    thermo_params = get_thermo_params(atmos_sim)
-    ts_int = get_field(atmos_sim, Val(:thermo_state_int))
+        p = atmos_sim.integrator.p
 
-    # 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, q_vap = q_sfc),
-    )
+        coupler_sfc_setup = coupler_surface_setup(CoupledMoninObukhov(), p; csf_sfc = csf_sfc)
+
+        p_names = propertynames(p)
+        p_values = map(x -> x == :sfc_setup ? coupler_sfc_setup : getproperty(p, x), p_names) # TODO: use merge here
 
-    p_names = propertynames(p)
-    p_values = map(x -> x == :sfc_setup ? coupler_sfc_setup : getproperty(p, x), p_names)
+        (; zip(p_names, p_values)...)
+    end
 
-    new_p = (; zip(p_names, p_values)...)
+    p = atmos_sim.integrator.p
 
+    csf_sfc = (; T = csf.T_S, z0m = csf.z0m_S, z0b = csf.z0b_S, beta = csf.beta, q_vap = csf.q_sfc)
+    new_p = modified_atmos_p(atmos_sim, csf_sfc)
     ClimaAtmos.SurfaceConditions.update_surface_conditions!(atmos_sim.integrator.u, new_p, atmos_sim.integrator.t)
-
     p.sfc_conditions .= new_p.sfc_conditions
-
 end
 
+"""
+    get_thermo_params(sim::ClimaAtmosSimulation)
+
+Returns the thermodynamic parameters from the atmospheric model simulation object.
+"""
 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)
+
+"""
+    calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::Fields.Field)
+
+Extension for this  to to calculate surface density.
+"""
+function calculate_surface_air_density(atmos_sim::ClimaAtmosSimulation, T_S::Fields.Field)
+    thermo_params = get_thermo_params(atmos_sim)
+    ts_int = get_field(atmos_sim, Val(:thermo_state_int))
+    extrapolate_ρ_to_sfc.(Ref(thermo_params), ts_int, swap_space!(ones(axes(ts_int.ρ)), T_S))
+end
 
 """
     extrapolate_ρ_to_sfc(thermo_params, ts_int, T_sfc)
 
-Uses the ideal gas law and hydrostatic balance to extrapolate for surface density.
+Uses the ideal gas law and hydrostatic balance to extrapolate for surface density pointwise.
 """
 function extrapolate_ρ_to_sfc(thermo_params, ts_in, T_sfc)
     T_int = TD.air_temperature(thermo_params, ts_in)

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

@@ -210,10 +210,10 @@ function bucket_init(
     # Initial conditions with no moisture
     Y, p, coords = initialize(model)
     anomaly = false
-    anomaly_tropics = true
+    anomaly_tropics = false
     hs_sfc = false
     Y.bucket.T = map(coords.subsurface) do coord
-        T_sfc_0 = FT(270.0)
+        T_sfc_0 = FT(285.0)
         radlat = coord.lat / FT(180) * pi
         ΔT = FT(0)
         if anomaly == true

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

@@ -83,6 +83,8 @@ end
 # required by Interfacer
 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{:surface_humidity}) =
+    ClimaLSM.surface_specific_humidity(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}) =
@@ -113,6 +115,10 @@ function update_field!(sim::BucketSimulation, ::Val{:snow_precipitation}, field)
     parent(sim.integrator.p.bucket.P_snow) .= parent(field)
 end
 
+function update_field!(sim::BucketSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.bucket.ρ_sfc) .= parent(field)
+end
+
 step!(sim::BucketSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
 
 reinit!(sim::BucketSimulation) = reinit!(sim.integrator)

experiments/AMIP/modular/components/ocean/slab_ocean_init.jl

@@ -36,14 +36,21 @@ function slab_ocean_space_init(::Type{FT}, space, p) where {FT}
 
     # initial condition
     T_sfc = map(coords) do coord
+
         T_sfc_0 = FT(p.T_init) #- FT(275) # close to the average of T_1 in atmos
-        anom_ampl = FT(0)
+        anomaly = false
+        anomaly_tropics = false
+        anom = FT(0)
         radlat = coord.lat / FT(180) * pi
-        lat_0 = FT(60) / FT(180) * pi
-        lon_0 = FT(-90) / FT(180) * pi
-        radlon = coord.long / FT(180) * pi
-        stdev = FT(5) / FT(180) * pi
-        anom = anom_ampl * exp(-((radlat - lat_0)^2 / 2stdev^2 + (radlon - lon_0)^2 / 2stdev^2))
+        if anomaly == true
+            lat_0 = FT(60) / FT(180) * pi
+            lon_0 = FT(-90) / FT(180) * pi
+            radlon = coord.long / FT(180) * pi
+            stdev = FT(5) / FT(180) * pi
+            anom = anom_ampl * exp(-((radlat - lat_0)^2 / 2stdev^2 + (radlon - lon_0)^2 / 2stdev^2))
+        elseif anomaly_tropics == true
+            anom = FT(20 * cos(radlat)^4)
+        end
         T_sfc = T_sfc_0 + anom
     end
 
@@ -59,14 +66,16 @@ function slab_ocean_rhs!(dY, Y, cache, t)
     FT = eltype(Y.T_sfc)
     rhs = @. -(F_turb_energy + F_radiative) / (p.h * p.ρ * p.c)
     parent(dY.T_sfc) .= parent(rhs .* Regridder.binary_mask.(area_fraction, threshold = eps(FT)))
+
+    @. cache.q_sfc = TD.q_vap_saturation_generic.(cache.thermo_params, Y.T_sfc, cache.ρ_sfc, TD.Liquid())
 end
 
 """
     ocean_init(::Type{FT}; tspan, dt, saveat, space, area_fraction, stepper = Euler()) where {FT}
 
 Initializes the `DiffEq` problem, and creates a Simulation-type object containing the necessary information for `step!` in the coupling loop.
 """
-function ocean_init(::Type{FT}; tspan, dt, saveat, space, area_fraction, stepper = Euler()) where {FT}
+function ocean_init(::Type{FT}; tspan, dt, saveat, space, area_fraction, thermo_params, stepper = Euler()) where {FT}
 
     params = OceanSlabParameters(FT(20), FT(1500.0), FT(800.0), FT(280.0), FT(1e-3), FT(1e-5), FT(0.06))
 
@@ -75,7 +84,10 @@ function ocean_init(::Type{FT}; tspan, dt, saveat, space, area_fraction, stepper
         params = params,
         F_turb_energy = ClimaCore.Fields.zeros(space),
         F_radiative = ClimaCore.Fields.zeros(space),
+        q_sfc = ClimaCore.Fields.zeros(space),
+        ρ_sfc = ClimaCore.Fields.zeros(space),
         area_fraction = area_fraction,
+        thermo_params = thermo_params,
     )
     problem = OrdinaryDiffEq.ODEProblem(slab_ocean_rhs!, Y, tspan, cache)
     integrator = OrdinaryDiffEq.init(problem, stepper, dt = dt, saveat = saveat)
@@ -92,6 +104,7 @@ clean_sst(SST, _info) = (swap_space!(zeros(axes(_info.land_fraction)), SST) .+ f
 
 # required by Interfacer
 get_field(sim::SlabOceanSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
+get_field(sim::SlabOceanSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
 get_field(sim::SlabOceanSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
 get_field(sim::SlabOceanSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
 get_field(sim::SlabOceanSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
@@ -108,6 +121,9 @@ end
 function update_field!(sim::SlabOceanSimulation, ::Val{:radiative_energy_flux}, field)
     parent(sim.integrator.p.F_radiative) .= parent(field)
 end
+function update_field!(sim::SlabOceanSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.ρ_sfc) .= parent(field)
+end
 
 step!(sim::SlabOceanSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
 

experiments/AMIP/modular/components/ocean/slab_seaice_init.jl

@@ -73,23 +73,28 @@ function ice_rhs!(du, u, p, _)
     unphysical = @. Regridder.binary_mask.(T_freeze - (Y.T_sfc + FT(rhs) * p.dt), threshold = FT(0)) .*
        Regridder.binary_mask.(area_fraction, threshold = eps(FT))
     parent(dY.T_sfc) .= parent(rhs .* unphysical)
+
+    @. p.q_sfc = TD.q_vap_saturation_generic.(p.thermo_params, Y.T_sfc, p.ρ_sfc, TD.Ice())
 end
 
 """
     ice_init(::Type{FT}; tspan, dt, saveat, space, ice_fraction, stepper = Euler()) where {FT}
 
 Initializes the `DiffEq` problem, and creates a Simulation-type object containing the necessary information for `step!` in the coupling loop.
 """
-function ice_init(::Type{FT}; tspan, saveat, dt, space, area_fraction, stepper = Euler()) where {FT}
+function ice_init(::Type{FT}; tspan, saveat, dt, space, area_fraction, thermo_params, stepper = Euler()) where {FT}
 
     params = IceSlabParameters(FT(2), FT(900.0), FT(2100.0), FT(271.2), FT(1e-3), FT(1e-5), FT(271.2), FT(2.0), FT(0.8))
 
     Y = slab_ice_space_init(FT, space, params)
     additional_cache = (;
         F_turb_energy = ClimaCore.Fields.zeros(space),
         F_radiative = ClimaCore.Fields.zeros(space),
+        q_sfc = ClimaCore.Fields.zeros(space),
+        ρ_sfc = ClimaCore.Fields.zeros(space),
         area_fraction = area_fraction,
         dt = dt,
+        thermo_params = thermo_params,
     )
 
     problem = OrdinaryDiffEq.ODEProblem(ice_rhs!, Y, tspan, (; additional_cache..., params = params))
@@ -112,6 +117,7 @@ get_ice_fraction(h_ice::FT, mono::Bool, threshold = 0.5) where {FT} =
 
 # required by Interfacer
 get_field(sim::PrescribedIceSimulation, ::Val{:surface_temperature}) = sim.integrator.u.T_sfc
+get_field(sim::PrescribedIceSimulation, ::Val{:surface_humidity}) = sim.integrator.p.q_sfc
 get_field(sim::PrescribedIceSimulation, ::Val{:roughness_momentum}) = sim.integrator.p.params.z0m
 get_field(sim::PrescribedIceSimulation, ::Val{:roughness_buoyancy}) = sim.integrator.p.params.z0b
 get_field(sim::PrescribedIceSimulation, ::Val{:beta}) = convert(eltype(sim.integrator.u), 1.0)
@@ -128,6 +134,9 @@ end
 function update_field!(sim::PrescribedIceSimulation, ::Val{:radiative_energy_flux}, field)
     parent(sim.integrator.p.F_radiative) .= parent(field)
 end
+function update_field!(sim::PrescribedIceSimulation, ::Val{:air_density}, field)
+    parent(sim.integrator.p.ρ_sfc) .= parent(field)
+end
 
 step!(sim::PrescribedIceSimulation, t) = step!(sim.integrator, t - sim.integrator.t, true)
 reinit!(sim::PrescribedIceSimulation) = reinit!(sim.integrator)