diff --git a/.buildkite/Manifest-v1.11.toml b/.buildkite/Manifest-v1.11.toml index be7fb8df4c..b75678b418 100644 --- a/.buildkite/Manifest-v1.11.toml +++ b/.buildkite/Manifest-v1.11.toml @@ -2918,9 +2918,9 @@ version = "1.0.2" [[deps.Thermodynamics]] deps = ["ForwardDiff", "Random", "RootSolvers"] -git-tree-sha1 = "ed3bda4ef913084e1e028d4fb05e26f1b69356f6" +git-tree-sha1 = "3958121b3ed912e0858960b682de6e98cb4949dd" uuid = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" -version = "0.15.4" +version = "0.15.7" weakdeps = ["ClimaParams"] [deps.Thermodynamics.extensions] diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index fc4efd60ef..36451f97ca 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -609,10 +609,10 @@ steps: # --config_file $CONFIG_PATH/baroclinic_wave_equil_amd.yml # --job_id baroclinic_wave_equil_amd - artifact_paths: "baroclinic_wave_equil_amd/output_active/*" - soft_fail: true - agents: - slurm_constraint: icelake|cascadelake|skylake|epyc + # artifact_paths: "baroclinic_wave_equil_amd/output_active/*" + # soft_fail: true + # agents: + # slurm_constraint: icelake|cascadelake|skylake|epyc - label: ":computer: held suarez" key: held_suarez diff --git a/NEWS.md b/NEWS.md index bd6d9935d9..9a3bcd7c28 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,9 @@ ClimaAtmos.jl Release Notes main ------- +- [#4231](https://github.com/CliMA/ClimaAtmos.jl/pull/4231) [badge-💥breaking] removes grid-scale +thermo state, including ᶜts in p.precomputed.sfc_conditions. + - [#4211](https://github.com/CliMA/ClimaAtmos.jl/pull/4211) add experimental methods to remove negative microphysical condensate values @@ -10,6 +13,7 @@ v0.34.0 ------- - [#4198](https://github.com/CliMA/ClimaAtmos.jl/pull/4198) [badge-💥breaking] modifies surface conditions to use SurfaceFluxes v0.15. + - [#4220](https://github.com/CliMA/ClimaAtmos.jl/pull/4220) modifies `SphereGrid` to use spacefillingcurve. v0.33.2 diff --git a/Project.toml b/Project.toml index 9e32bde03e..cd1a5ec380 100644 --- a/Project.toml +++ b/Project.toml @@ -73,7 +73,7 @@ SparseMatrixColorings = "0.4.20" StaticArrays = "1.9" Statistics = "1" SurfaceFluxes = "0.15" -Thermodynamics = "0.15.3" +Thermodynamics = "0.15.7" UnrolledUtilities = "0.1.9" YAML = "0.4" julia = "1.9" diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 88dc414dfa..b0a97fea4b 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -1051,6 +1051,7 @@ AquaplanetPlots = Union{ Val{:gpu_aquaplanet_dyamond_summer}, Val{:edonly_edmfx_aquaplanet}, Val{:mpi_sphere_aquaplanet_rhoe_equil_clearsky}, + Val{:aquaplanet_equil_allsky_gw_raw}, Val{:aquaplanet_nonequil_allsky_gw_res}, Val{:aquaplanet_nonequil_allsky_gw_res_2M}, Val{:rcemipii_sphere_diagnostic_edmfx}, diff --git a/reproducibility_tests/ref_counter.jl b/reproducibility_tests/ref_counter.jl index e3a3945edf..28bcc87371 100644 --- a/reproducibility_tests/ref_counter.jl +++ b/reproducibility_tests/ref_counter.jl @@ -1,4 +1,4 @@ -303 +304 # **README** # @@ -21,6 +21,9 @@ #= +304 +- Remove grid-scale thermo state from precomputed quantities and uses new thermodynamics functions. + 303 - Remove thermo state from initial conditions, which changes the behavior of prognostic EDMF case, possibly from the change in `enthalpy` function. diff --git a/src/cache/cloud_fraction.jl b/src/cache/cloud_fraction.jl index d9ee47627e..942fc0d889 100644 --- a/src/cache/cloud_fraction.jl +++ b/src/cache/cloud_fraction.jl @@ -13,7 +13,7 @@ end Diagnose horizontal covariances based on vertical gradients (i.e. taking turbulence production as the only term) """ -function compute_covariance(Y, p, thermo_params, ᶜts) +function compute_covariance(Y, p, thermo_params) coeff = CAP.diagnostic_covariance_coeff(p.params) turbconv_model = p.atmos.turbconv_model (; ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed @@ -23,10 +23,20 @@ function compute_covariance(Y, p, thermo_params, ᶜts) if isnothing(turbconv_model) if p.atmos.call_cloud_diagnostics_per_stage isa CallCloudDiagnosticsPerStage - @. ᶜgradᵥ_q_tot = - ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts))) - @. ᶜgradᵥ_θ_liq_ice = - ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts))) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed + @. ᶜgradᵥ_q_tot = ᶜgradᵥ(ᶠinterp(ᶜq_tot_safe)) + @. ᶜgradᵥ_θ_liq_ice = ᶜgradᵥ( + ᶠinterp( + TD.liquid_ice_pottemp( + thermo_params, + ᶜT, + Y.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ), + ), + ) end end # Reminder that gradients need to be precomputed when using compute_gm_mixing_length @@ -99,14 +109,16 @@ end """ function compute_cloud_fraction_quadrature_diagnostics( - thermo_params, SG_quad, ts, q′q′, θ′θ′, θ′q′ + thermo_params, SG_quad, p_c, q_mean, θ_mean, q′q′, θ′θ′, θ′q′ ) where: - thermo params - thermodynamics parameters - SG_quad is a struct containing information about Gaussian quadrature order, sampling point values and weights - - ts is the thermodynamic state + - p_c is the air pressure + - q_mean is the total specific humidity + - θ_mean is the liquid-ice potential temperature - q′q′, θ′θ′, θ′q′ are the covariances of q_tot and liquid ice potential temperature The function imposes additional limits on the quadrature points @@ -117,15 +129,13 @@ computed as a sum over quadrature points with weights. function compute_cloud_fraction_quadrature_diagnostics( thermo_params, SG_quad::SGSQuadrature, - ts, + p_c, + q_mean, + θ_mean, q′q′, θ′θ′, θ′q′, ) - # Grab mean pressure, liquid ice potential temperature and total specific humidity - p_c = TD.air_pressure(thermo_params, ts) - q_mean = TD.total_specific_humidity(thermo_params, ts) - θ_mean = TD.liquid_ice_pottemp(thermo_params, ts) # Return physical values based on quadrature points and limited covarainces function get_x_hat(χ1, χ2) @@ -158,7 +168,7 @@ function compute_cloud_fraction_quadrature_diagnostics( function f(x1_hat, x2_hat) FT = eltype(thermo_params) - _ts = thermo_state(thermo_params; p = p_c, θ = x1_hat, q_tot = x2_hat) + _ts = TD.PhaseEquil_pθq(thermo_params, p_c, x1_hat, x2_hat) hc = TD.has_condensate(thermo_params, _ts) cf = hc ? FT(1) : FT(0) # cloud fraction @@ -192,23 +202,22 @@ NVTX.@annotate function set_cloud_fraction!( moist_model::Union{EquilMoistModel, NonEquilMoistModel}, ::GridScaleCloud, ) - (; ᶜts) = p.precomputed - thermo_params = CAP.thermodynamics_params(p.params) + (; ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed FT = eltype(p.params) if moist_model isa EquilMoistModel @. p.precomputed.cloud_diagnostics_tuple = make_cloud_fraction_named_tuple( - ifelse(TD.has_condensate(thermo_params, ᶜts), FT(1), FT(0)), - TD.PhasePartition(thermo_params, ᶜts).liq, - TD.PhasePartition(thermo_params, ᶜts).ice, + ifelse(TD.has_condensate(ᶜq_liq_rai + ᶜq_ice_sno), FT(1), FT(0)), + ᶜq_liq_rai, + ᶜq_ice_sno, ) else q_liq = @. lazy(specific(Y.c.ρq_liq, Y.c.ρ)) q_ice = @. lazy(specific(Y.c.ρq_ice, Y.c.ρ)) @. p.precomputed.cloud_diagnostics_tuple = make_cloud_fraction_named_tuple( - ifelse(TD.has_condensate(thermo_params, ᶜts), FT(1), FT(0)), + ifelse(TD.has_condensate(ᶜq_liq_rai + ᶜq_ice_sno), FT(1), FT(0)), q_liq, q_ice, ) @@ -222,18 +231,40 @@ NVTX.@annotate function set_cloud_fraction!( ) thermo_params = CAP.thermodynamics_params(p.params) turbconv_model = p.atmos.turbconv_model + (; ᶜp, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed - ᶜts = turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶜts⁰ : p.precomputed.ᶜts + # For PrognosticEDMFX, use environment state; otherwise use grid-scale + if turbconv_model isa PrognosticEDMFX + ᶜts⁰ = p.precomputed.ᶜts⁰ + ᶜp_env = @. lazy(TD.air_pressure(thermo_params, ᶜts⁰)) + ᶜq_mean = @. lazy(TD.total_specific_humidity(thermo_params, ᶜts⁰)) + ᶜθ_mean = @. lazy(TD.liquid_ice_pottemp(thermo_params, ᶜts⁰)) + else + ᶜp_env = ᶜp + ᶜq_mean = ᶜq_tot_safe + ᶜθ_mean = @. lazy( + TD.liquid_ice_pottemp( + thermo_params, + ᶜT, + Y.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ), + ) + end # Compute covariance based on the gradients of q_tot and theta_liq_ice - ᶜq′q′, ᶜθ′θ′, ᶜθ′q′ = compute_covariance(Y, p, thermo_params, ᶜts) + ᶜq′q′, ᶜθ′θ′, ᶜθ′q′ = compute_covariance(Y, p, thermo_params) # Compute SGS cloud fraction diagnostics based on environment quadrature points ... @. p.precomputed.cloud_diagnostics_tuple = compute_cloud_fraction_quadrature_diagnostics( thermo_params, qc.SG_quad, - ᶜts, + ᶜp_env, + ᶜq_mean, + ᶜθ_mean, ᶜq′q′, ᶜθ′θ′, ᶜθ′q′, @@ -247,18 +278,21 @@ NVTX.@annotate function set_cloud_fraction!( qc, thermo_params, turbconv_model, - ᶜts, + ᶜp_env, + ᶜq_mean, + ᶜθ_mean, ) end # ... weight by environment area fraction if using PrognosticEDMFX (assumed 1 otherwise) ... if turbconv_model isa PrognosticEDMFX ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, p.atmos.turbconv_model)) + ᶜρ⁰ = @. lazy(TD.air_density(thermo_params, p.precomputed.ᶜts⁰)) @. p.precomputed.cloud_diagnostics_tuple *= NamedTuple{(:cf, :q_liq, :q_ice)}( tuple( - draft_area(ᶜρa⁰, TD.air_density(thermo_params, ᶜts)), - draft_area(ᶜρa⁰, TD.air_density(thermo_params, ᶜts)), - draft_area(ᶜρa⁰, TD.air_density(thermo_params, ᶜts)), + draft_area(ᶜρa⁰, ᶜρ⁰), + draft_area(ᶜρa⁰, ᶜρ⁰), + draft_area(ᶜρa⁰, ᶜρ⁰), ), ) end @@ -294,7 +328,9 @@ function set_ml_cloud_fraction!( ml_cloud::MLCloud, thermo_params, turbconv_model, - ᶜts, + ᶜp_env, + ᶜq_mean, + ᶜθ_mean, ) FT = eltype(p.params) ᶜmixing_length_field = p.scratch.ᶜtemp_scalar @@ -325,8 +361,9 @@ function set_ml_cloud_fraction!( ᶜmixing_length_field, ᶜ∇q, ᶜ∇θ, - specific.(Y.c.ρq_tot, Y.c.ρ), - ᶜts, + ᶜq_mean, + ᶜp_env, + ᶜθ_mean, thermo_params, ) end @@ -337,20 +374,18 @@ function compute_ml_cloud_fraction( ᶜ∇q, ᶜ∇θ, q_tot, - ᶜts, + p_c, + θli, thermo_params, ) FT = eltype(thermo_params) - # Saturation state at current thermodynamic state - q_sat = TD.q_vap_saturation(thermo_params, ᶜts) - - # Liquid–ice potential temperature at current thermodynamic state - θli = TD.liquid_ice_pottemp(thermo_params, ᶜts) - + # Saturation state at current thermodynamic state (via θ, p, q_tot) + ts_current = TD.PhaseEquil_pθq(thermo_params, p_c, θli, q_tot) + q_sat = TD.q_vap_saturation(thermo_params, ts_current) # distance to saturation in temperature space Δθli, θli_sat, dqsatdθli = - saturation_distance(q_tot, q_sat, ᶜts, θli, thermo_params, FT(0.1)) + saturation_distance(q_tot, q_sat, p_c, θli, thermo_params, FT(0.1)) # form the pi groups π_1 = (q_sat - q_tot) / q_sat @@ -359,17 +394,16 @@ function compute_ml_cloud_fraction( π_4 = (ᶜ∇θ * ᶜmixing_length_field) / θli_sat return apply_cf_nn(nn_model, π_1, π_2, π_3, π_4) - end -function saturation_distance(q_tot, q_sat, ᶜts, θli, thermo_params, Δθli_fd) +function saturation_distance(q_tot, q_sat, p_c, θli, thermo_params, Δθli_fd) # Perturbed thermodynamic states for finite-difference ts_perturbed = TD.PhaseEquil_pθq( thermo_params, - ᶜts.p, - θli .+ Δθli_fd, - ᶜts.q_tot, + p_c, + θli + Δθli_fd, + q_tot, ) q_sat_perturbed = TD.q_vap_saturation(thermo_params, ts_perturbed) diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index 6f4046344a..c0550de195 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -162,7 +162,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( FT = eltype(Y) n = n_mass_flux_subdomains(turbconv_model) (; ᶜΦ) = p.core - (; ᶜp, ᶠu³, ᶜK) = p.precomputed + (; ᶜp, ᶠu³, ᶜK, ᶜh_tot) = p.precomputed (; ustar, obukhov_length, buoyancy_flux, ρ_flux_h_tot, ρ_flux_q_tot) = p.precomputed.sfc_conditions (; ᶜρaʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs, ᶜtsʲs, ᶜρʲs) = p.precomputed @@ -172,11 +172,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( thermo_params = CAP.thermodynamics_params(params) turbconv_params = CAP.turbconv_params(params) - ᶜts = p.precomputed.ᶜts #TODO replace ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) ρ_int_level = Fields.field_values(Fields.level(Y.c.ρ, 1)) uₕ_int_level = Fields.field_values(Fields.level(Y.c.uₕ, 1)) @@ -450,7 +447,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( (; params) = p (; dt) = p (; ᶜΦ, ᶜgradᵥ_ᶠΦ) = p.core - (; ᶜp, ᶠu³, ᶜts, ᶜK) = p.precomputed + (; ᶜp, ᶠu³, ᶜT, ᶜh_tot, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜK) = p.precomputed (; ᶜρaʲs, ᶠu³ʲs, @@ -496,8 +493,6 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( # integral ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) ᶜtke = @. lazy(specific(Y.c.ρtke, Y.c.ρ)) @@ -533,7 +528,9 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( K_prev_level = Fields.field_values(Fields.level(ᶜK, i - 1)) h_tot_prev_level = Fields.field_values(Fields.level(ᶜh_tot, i - 1)) q_tot_prev_level = Fields.field_values(Fields.level(ᶜq_tot, i - 1)) - ts_prev_level = Fields.field_values(Fields.level(ᶜts, i - 1)) + T_prev_level = Fields.field_values(Fields.level(ᶜT, i - 1)) + q_liq_rai_prev_level = Fields.field_values(Fields.level(ᶜq_liq_rai, i - 1)) + q_ice_sno_prev_level = Fields.field_values(Fields.level(ᶜq_ice_sno, i - 1)) p_prev_level = Fields.field_values(Fields.level(ᶜp, i - 1)) z_prev_level = Fields.field_values(Fields.level(ᶜz, i - 1)) dz_prev_level = Fields.field_values(Fields.level(ᶜdz, i - 1)) @@ -685,7 +682,14 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( u³_prev_halflevel, local_geometry_prev_halflevel, ), - TD.relative_humidity(thermo_params, ts_prev_level), + TD.relative_humidity( + thermo_params, + T_prev_level, + p_prev_level, + q_tot_prev_level, + q_liq_rai_prev_level, + q_ice_sno_prev_level, + ), FT(0), tke_prev_level, p.atmos.edmfx_model.entr_model, @@ -757,11 +761,12 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( # 0-moment microphysics: sink of q_tot from precipitation removal if microphysics_model isa Microphysics0Moment @. S_q_totʲ_prev_level = q_tot_0M_precipitation_sources( - thermo_params, microphys_0m_params, dt, q_totʲ_prev_level, - tsʲ_prev_level, + TD.total_specific_humidity(thermo_params, tsʲ_prev_level), + TD.liquid_specific_humidity(thermo_params, tsʲ_prev_level), + TD.ice_specific_humidity(thermo_params, tsʲ_prev_level), ) # 1-moment microphysics: cloud water (liquid and ice) and # precipitation (rain and snow) sources. q_tot is constant, because @@ -769,6 +774,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( elseif moisture_model isa NonEquilMoistModel && microphysics_model isa Microphysics1Moment # Rain formation from the updrafts + Tʲ_prev_level = @. lazy(TD.air_temperature(thermo_params, tsʲ_prev_level)) compute_precipitation_sources!( Sᵖ_prev_level, Sᵖ_snow_prev_level, @@ -782,7 +788,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( q_iceʲ_prev_level, q_raiʲ_prev_level, q_snoʲ_prev_level, - tsʲ_prev_level, + Tʲ_prev_level, dt, microphys_1m_params, thermo_params, @@ -798,7 +804,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( q_iceʲ_prev_level, q_raiʲ_prev_level, q_snoʲ_prev_level, - tsʲ_prev_level, + Tʲ_prev_level, dt, microphys_1m_params, thermo_params, @@ -813,7 +819,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( q_raiʲ_prev_level, q_snoʲ_prev_level, ρʲ_prev_level, - TD.air_temperature(thermo_params, tsʲ_prev_level), + Tʲ_prev_level, dt, ) @. S_q_iceʲ_prev_level += cloud_sources( @@ -825,7 +831,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( q_raiʲ_prev_level, q_snoʲ_prev_level, ρʲ_prev_level, - TD.air_temperature(thermo_params, tsʲ_prev_level), + Tʲ_prev_level, dt, ) @@ -885,7 +891,10 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( u³ʲ_prev_halflevel, local_geometry_prev_halflevel, u³_prev_halflevel, - ts_prev_level, + T_prev_level, + q_tot_prev_level, + q_liq_rai_prev_level, + q_ice_sno_prev_level, FT(0), entrʲ_prev_level, vert_div_level, @@ -986,6 +995,11 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( mseʲ_prev_level, ) if microphysics_model isa Microphysics0Moment + ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, tsʲ_prev_level)) + ᶜq_liq_raiʲ = + @. lazy(TD.liquid_specific_humidity(thermo_params, tsʲ_prev_level)) + ᶜq_ice_snoʲ = + @. lazy(TD.ice_specific_humidity(thermo_params, tsʲ_prev_level)) @. ρaʲu³ʲ_datamse += microphysics_sources( local_geometry_halflevel.J, local_geometry_prev_level.J, @@ -993,7 +1007,9 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( S_q_totʲ_prev_level * e_tot_0M_precipitation_sources_helper( thermo_params, - tsʲ_prev_level, + ᶜTʲ, + ᶜq_liq_raiʲ, + ᶜq_ice_snoʲ, Φ_prev_level, ), ) @@ -1334,7 +1350,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures! t, ) (; params) = p - (; ᶜu, ᶜts, ᶠu³) = p.precomputed + (; ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶠu³) = p.precomputed (; ustar) = p.precomputed.sfc_conditions (; ᶜlinear_buoygrad, ᶜstrain_rate_norm) = p.precomputed (; ρtke_flux) = p.precomputed @@ -1347,10 +1363,15 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_closures! @. ᶜu⁰ = C123(Y.c.uₕ) + ᶜinterp(C123(ᶠu³⁰)) end # Set here, but used in a different function + ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) @. ᶜlinear_buoygrad = buoyancy_gradients( BuoyGradMean(), thermo_params, - ᶜts, + ᶜT, + Y.c.ρ, + ᶜq_tot, + ᶜq_liq_rai, + ᶜq_ice_sno, p.precomputed.cloud_diagnostics_tuple.cf, C3, p.precomputed.ᶜgradᵥ_q_tot, @@ -1398,19 +1419,19 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita t, microphysics_model::Microphysics0Moment, ) - thermo_params = CAP.thermodynamics_params(p.params) microphys_0m_params = CAP.microphysics_0m_params(p.params) (; dt) = p - (; ᶜts, ᶜSqₜᵖ⁰) = p.precomputed + (; ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜSqₜᵖ⁰) = p.precomputed # Environment precipitation sources (to be applied to grid mean) ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources( - thermo_params, microphys_0m_params, dt, ᶜq_tot, - ᶜts, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, ) return nothing end @@ -1425,7 +1446,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita cloud_params = CAP.microphysics_cloud_params(p.params) (; dt) = p - (; ᶜts, ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed + (; ᶜT, ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed ᶜSᵖ = p.scratch.ᶜtemp_scalar ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2 @@ -1443,7 +1464,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita specific.(Y.c.ρq_ice, Y.c.ρ), specific.(Y.c.ρq_rai, Y.c.ρ), specific.(Y.c.ρq_sno, Y.c.ρ), - ᶜts, + ᶜT, dt, microphys_1m_params, thermo_params, @@ -1459,7 +1480,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita specific.(Y.c.ρq_ice, Y.c.ρ), specific.(Y.c.ρq_rai, Y.c.ρ), specific.(Y.c.ρq_sno, Y.c.ρ), - ᶜts, + ᶜT, dt, microphys_1m_params, thermo_params, @@ -1474,7 +1495,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita specific(Y.c.ρq_rai, Y.c.ρ), specific(Y.c.ρq_sno, Y.c.ρ), Y.c.ρ, - TD.air_temperature(thermo_params, ᶜts), + ᶜT, dt, ) @. ᶜSqᵢᵖ⁰ += cloud_sources( @@ -1486,7 +1507,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita specific(Y.c.ρq_rai, Y.c.ρ), specific(Y.c.ρq_sno, Y.c.ρ), Y.c.ρ, - TD.air_temperature(thermo_params, ᶜts), + ᶜT, dt, ) return nothing diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl index 98bc705649..1574f346d5 100644 --- a/src/cache/precipitation_precomputed_quantities.jl +++ b/src/cache/precipitation_precomputed_quantities.jl @@ -56,7 +56,7 @@ function set_precipitation_velocities!( microphysics_model::Microphysics1Moment, _, ) - (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed + (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜT, ᶜu) = p.precomputed (; ᶜΦ) = p.core cmc = CAP.microphysics_cloud_params(p.params) cmp = CAP.microphysics_1m_params(p.params) @@ -99,10 +99,10 @@ function set_precipitation_velocities!( ) / Y.c.ρ @. ᶜwₕhₜ = Geometry.WVector( - ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) + - ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) + - ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) + - ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))), + ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) + + ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) + + ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) + + ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))), ) / Y.c.ρ return nothing end @@ -284,7 +284,7 @@ function set_precipitation_velocities!( microphysics_model::Microphysics2Moment, _, ) - (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₙₗ, ᶜwₙᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed + (; ᶜwₗ, ᶜwᵢ, ᶜwᵣ, ᶜwₛ, ᶜwₙₗ, ᶜwₙᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜT, ᶜu) = p.precomputed (; ᶜΦ) = p.core cmc = CAP.microphysics_cloud_params(p.params) @@ -359,10 +359,10 @@ function set_precipitation_velocities!( ) / Y.c.ρ @. ᶜwₕhₜ = Geometry.WVector( - ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) + - ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) + - ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) + - ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))), + ᶜwₗ * Y.c.ρq_liq * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) + + ᶜwᵢ * Y.c.ρq_ice * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) + + ᶜwᵣ * Y.c.ρq_rai * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))) + + ᶜwₛ * Y.c.ρq_sno * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₛ, ᶜu))), ) / Y.c.ρ return nothing end @@ -555,7 +555,7 @@ function set_precipitation_velocities!( Y, p, ::NonEquilMoistModel, ::Microphysics2MomentP3, ) ## liquid quantities (2M warm rain) - (; ᶜwₗ, ᶜwᵣ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜts, ᶜu) = p.precomputed + (; ᶜwₗ, ᶜwᵣ, ᶜwnₗ, ᶜwnᵣ, ᶜwₜqₜ, ᶜwₕhₜ, ᶜT, ᶜu) = p.precomputed (; ᶜΦ) = p.core (; ρ, ρq_liq, ρn_liq, ρq_rai, ρn_rai) = Y.c @@ -606,9 +606,9 @@ function set_precipitation_velocities!( @. ᶜwₜqₜ = Geometry.WVector(ᶜwₗ * ρq_liq + ᶜwᵢ * ρq_ice + ᶜwᵣ * ρq_rai) / ρ @. ᶜwₕhₜ = Geometry.WVector( - ᶜwₗ * ρq_liq * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) + - ᶜwᵢ * ρq_ice * (Iᵢ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) + - ᶜwᵣ * ρq_rai * (Iₗ(thp, ᶜts) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))), + ᶜwₗ * ρq_liq * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwₗ, ᶜu))) + + ᶜwᵢ * ρq_ice * (Iᵢ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵢ, ᶜu))) + + ᶜwᵣ * ρq_rai * (Iₗ(thp, ᶜT) + ᶜΦ + $(Kin(ᶜwᵣ, ᶜu))), ) / ρ return nothing end @@ -624,22 +624,23 @@ sources from the sub-domains. set_precipitation_cache!(Y, p, _, _) = nothing function set_precipitation_cache!(Y, p, ::Microphysics0Moment, _) (; params, dt) = p - (; ᶜts) = p.precomputed + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precomputed (; ᶜΦ) = p.core cm_params = CAP.microphysics_0m_params(params) thermo_params = CAP.thermodynamics_params(params) @. ᶜS_ρq_tot = Y.c.ρ * q_tot_0M_precipitation_sources( - thermo_params, cm_params, dt, Y.c.ρq_tot / Y.c.ρ, - ᶜts, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, ) @. ᶜS_ρe_tot = ᶜS_ρq_tot * - e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, ᶜΦ) + e_tot_0M_precipitation_sources_helper(thermo_params, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶜΦ) return nothing end function set_precipitation_cache!( @@ -653,7 +654,7 @@ function set_precipitation_cache!( (; ᶜΦ) = p.core (; ᶜSqₜᵖ⁰, ᶜSqₜᵖʲs, ᶜρaʲs) = p.precomputed (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precomputed - (; ᶜts, ᶜtsʲs) = p.precomputed + (; ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶜtsʲs) = p.precomputed thermo_params = CAP.thermodynamics_params(p.params) n = n_mass_flux_subdomains(p.atmos.turbconv_model) @@ -663,15 +664,20 @@ function set_precipitation_cache!( @. ᶜS_ρe_tot = ᶜSqₜᵖ⁰ * ρ * - e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts, ᶜΦ) + e_tot_0M_precipitation_sources_helper(thermo_params, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶜΦ) for j in 1:n + ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, ᶜtsʲs.:($$j))) + ᶜq_liq_raiʲ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜtsʲs.:($$j))) + ᶜq_ice_snoʲ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜtsʲs.:($$j))) @. ᶜS_ρq_tot += ᶜSqₜᵖʲs.:($$j) * ᶜρaʲs.:($$j) @. ᶜS_ρe_tot += ᶜSqₜᵖʲs.:($$j) * ᶜρaʲs.:($$j) * e_tot_0M_precipitation_sources_helper( thermo_params, - ᶜtsʲs.:($$j), + ᶜTʲ, + ᶜq_liq_raiʲ, + ᶜq_ice_snoʲ, ᶜΦ, ) end @@ -693,18 +699,32 @@ function set_precipitation_cache!( n = n_mass_flux_subdomains(p.atmos.turbconv_model) @. ᶜS_ρq_tot = ᶜSqₜᵖ⁰ * ᶜρa⁰ + ᶜT⁰ = @. lazy(TD.air_temperature(thermo_params, ᶜts⁰)) + ᶜq_liq_rai⁰ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜts⁰)) + ᶜq_ice_sno⁰ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜts⁰)) @. ᶜS_ρe_tot = ᶜSqₜᵖ⁰ * ᶜρa⁰ * - e_tot_0M_precipitation_sources_helper(thermo_params, ᶜts⁰, ᶜΦ) + e_tot_0M_precipitation_sources_helper( + thermo_params, + ᶜT⁰, + ᶜq_liq_rai⁰, + ᶜq_ice_sno⁰, + ᶜΦ, + ) for j in 1:n + ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, ᶜtsʲs.:($$j))) + ᶜq_liq_raiʲ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜtsʲs.:($$j))) + ᶜq_ice_snoʲ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜtsʲs.:($$j))) @. ᶜS_ρq_tot += ᶜSqₜᵖʲs.:($$j) * Y.c.sgsʲs.:($$j).ρa @. ᶜS_ρe_tot += ᶜSqₜᵖʲs.:($$j) * Y.c.sgsʲs.:($$j).ρa * e_tot_0M_precipitation_sources_helper( thermo_params, - ᶜtsʲs.:($$j), + ᶜTʲ, + ᶜq_liq_raiʲ, + ᶜq_ice_snoʲ, ᶜΦ, ) end @@ -712,7 +732,7 @@ function set_precipitation_cache!( end function set_precipitation_cache!(Y, p, ::Microphysics1Moment, _) (; dt) = p - (; ᶜts, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed + (; ᶜT, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) @@ -744,7 +764,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics1Moment, _) ᶜq_ice, ᶜq_rai, ᶜq_sno, - ᶜts, + ᶜT, dt, cmp, thp, @@ -761,7 +781,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics1Moment, _) ᶜq_ice, ᶜq_rai, ᶜq_sno, - ᶜts, + ᶜT, dt, cmp, thp, @@ -790,8 +810,7 @@ function set_precipitation_cache!( end function set_precipitation_cache!(Y, p, ::Microphysics2Moment, _) (; dt) = p - (; ᶜts) = p.precomputed - (; ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed + (; ᶜT, ᶜSqₗᵖ, ᶜSqᵢᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ) = p.precomputed (; ᶜSnₗᵖ, ᶜSnᵣᵖ) = p.precomputed ᶜSᵖ = p.scratch.ᶜtemp_scalar @@ -818,7 +837,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics2Moment, _) lazy.(specific.(Y.c.ρq_ice, Y.c.ρ)), lazy.(specific.(Y.c.ρq_rai, Y.c.ρ)), lazy.(specific.(Y.c.ρq_sno, Y.c.ρ)), - ᶜts, + ᶜT, dt, cmp, thp, @@ -857,7 +876,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics2MomentP3, ::Nothing) # NOTE: the above function sets `ᶜSqᵢᵖ` to `0`. For P3, need to update `ᶜSqᵢᵖ` below!! ### Icy processes (P3) - (; ᶜScoll, ᶜts, ᶜlogλ) = p.precomputed + (; ᶜScoll, ᶜT, ᶜlogλ) = p.precomputed # get thermodynamics and microphysics params (; params) = p @@ -874,7 +893,7 @@ function set_precipitation_cache!(Y, p, ::Microphysics2MomentP3, ::Nothing) # compute warm precipitation sources on the grid mean (based on SB2006 2M scheme) compute_cold_precipitation_sources_P3!( - ᶜScoll, params_2mp3, thermo_params, ᶜY_reduced, ᶜts, ᶜlogλ, + ᶜScoll, params_2mp3, thermo_params, ᶜY_reduced, ᶜT, ᶜlogλ, ) return nothing @@ -894,8 +913,7 @@ function set_precipitation_surface_fluxes!( p, microphysics_model::Microphysics0Moment, ) - ᶜT = p.scratch.ᶜtemp_scalar - (; ᶜts) = p.precomputed # assume ᶜts has been updated + (; ᶜT) = p.precomputed (; ᶜS_ρq_tot, ᶜS_ρe_tot) = p.precomputed (; surface_rain_flux, surface_snow_flux) = p.precomputed (; col_integrated_precip_energy_tendency) = p.conservation_check @@ -909,7 +927,6 @@ function set_precipitation_surface_fluxes!( thermo_params = CAP.thermodynamics_params(p.params) T_freeze = TD.Parameters.T_freeze(thermo_params) FT = eltype(p.params) - @. ᶜT = TD.air_temperature(thermo_params, ᶜts) ᶜ3d_rain = @. lazy(ifelse(ᶜT >= T_freeze, ᶜS_ρq_tot, FT(0))) ᶜ3d_snow = @. lazy(ifelse(ᶜT < T_freeze, ᶜS_ρq_tot, FT(0))) Operators.column_integral_definite!(surface_rain_flux, ᶜ3d_rain) diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index ed521f9f10..f6cd3be19c 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -16,7 +16,10 @@ The following grid-scale quantities are treated implicitly and are precomputed: - `ᶜu`: covariant velocity on cell centers - `ᶠu`: contravariant velocity on cell faces - `ᶜK`: kinetic energy on cell centers - - `ᶜts`: thermodynamic state on cell centers + - `ᶜT`: air temperature on cell centers + - `ᶜq_tot_safe`: total water specific humidity on cell centers + - `ᶜq_liq_rai`: liquid water specific humidity on cell centers + - `ᶜq_ice_sno`: ice specific humidity on cell centers - `ᶜp`: air pressure on cell centers If the `turbconv_model` is `PrognosticEDMFX`, there also two SGS versions of every quantity except for `ᶜp` (which is shared across all subdomains): @@ -40,9 +43,29 @@ function implicit_precomputed_quantities(Y, atmos) ᶠu³ = similar(Y.f, CT3{FT}), ᶠu = similar(Y.f, CT123{FT}), ᶜK = similar(Y.c, FT), - ᶜts = similar(Y.c, TST), + ᶜT = similar(Y.c, FT), + ᶜh_tot = similar(Y.c, FT), ᶜp = similar(Y.c, FT), ) + # Moisture-related quantities depend on moisture model: + # - EquilMoistModel: allocate fields + thermo_state cache for saturation adjustment + # - Ohters: allocate fields only + thermo_state_type_gs = @NamedTuple{T::FT, q_tot_safe::FT, q_liq_rai::FT, q_ice_sno::FT} + moist_quantities = + if moisture_model isa EquilMoistModel + (; + ᶜq_tot_safe = similar(Y.c, FT), + ᶜq_liq_rai = similar(Y.c, FT), + ᶜq_ice_sno = similar(Y.c, FT), + ᶜthermo_state = similar(Y.c, thermo_state_type_gs), + ) + else # DryModel or NonEquilMoistModel + (; + ᶜq_tot_safe = similar(Y.c, FT), + ᶜq_liq_rai = similar(Y.c, FT), + ᶜq_ice_sno = similar(Y.c, FT), + ) + end sgs_quantities = (;) prognostic_sgs_quantities = turbconv_model isa PrognosticEDMFX ? @@ -60,7 +83,12 @@ function implicit_precomputed_quantities(Y, atmos) ᶜρʲs = similar(Y.c, NTuple{n, FT}), ᶠnh_pressure₃_dragʲs = similar(Y.f, NTuple{n, C3{FT}}), ) : (;) - return (; gs_quantities..., sgs_quantities..., prognostic_sgs_quantities...) + return (; + gs_quantities..., + moist_quantities..., + sgs_quantities..., + prognostic_sgs_quantities..., + ) end """ @@ -379,65 +407,22 @@ function add_sgs_ᶜK!(ᶜK, Y, ᶜρa⁰, ᶠu₃⁰, turbconv_model) return nothing end -function thermo_state( - thermo_params; - ρ = nothing, - p = nothing, - θ = nothing, - e_int = nothing, - q_tot = nothing, - q_pt = nothing, -) - get_ts(ρ::Real, ::Nothing, ::Nothing, e_int::Real, ::Nothing, ::Nothing) = - TD.PhaseDry_ρe(thermo_params, ρ, e_int) - get_ts(ρ::Real, ::Nothing, ::Nothing, e_int::Real, q_tot::Real, ::Nothing) = - TD.PhaseEquil_ρeq( - thermo_params, - ρ, - e_int, - q_tot, - 3, - eltype(thermo_params)(0.003), - ) - get_ts(ρ::Real, ::Nothing, ::Nothing, e_int::Real, ::Nothing, q_pt) = - TD.PhaseNonEquil(thermo_params, e_int, ρ, q_pt) - get_ts(::Nothing, p::Real, θ::Real, ::Nothing, q_tot::Real, ::Nothing) = - TD.PhaseEquil_pθq(thermo_params, p, θ, q_tot) - return get_ts(ρ, p, θ, e_int, q_tot, q_pt) -end - -# GPU-compatible thermo_vars using multiple dispatch. Each method returns a -# NamedTuple with the thermodynamic variables needed for the corresponding -# moisture model. +# Combined getter function for all thermodynamic state variables. +# Returns a NamedTuple with T, q_tot_safe, q_liq_rai, q_ice_sno. +# This avoids redundant saturation_adjustment calls for EquilMoistModel. -function thermo_vars(::DryModel, ::Any, ᶜY, K, Φ) - e_int = specific(ᶜY.ρe_tot, ᶜY.ρ) - K - Φ - return (; e_int) -end - -function thermo_vars(::EquilMoistModel, ::Any, ᶜY, K, Φ) - e_int = specific(ᶜY.ρe_tot, ᶜY.ρ) - K - Φ - q_tot = specific(ᶜY.ρq_tot, ᶜY.ρ) - return (; e_int, q_tot) -end - -function thermo_vars(::NonEquilMoistModel, ::Any, ᶜY, K, Φ) - e_int = specific(ᶜY.ρe_tot, ᶜY.ρ) - K - Φ - q_pt_args = (; - q_tot = specific(ᶜY.ρq_tot, ᶜY.ρ), - q_liq = specific(ᶜY.ρq_liq, ᶜY.ρ) + specific(ᶜY.ρq_rai, ᶜY.ρ), - q_ice = specific(ᶜY.ρq_ice, ᶜY.ρ) + specific(ᶜY.ρq_sno, ᶜY.ρ), +function thermo_state_gs(thermo_params, ρ, ρe_tot, ρq_tot, K, Φ) + e_int = specific(ρe_tot, ρ) - K - Φ + q_tot_safe = max(0, specific(ρq_tot, ρ)) + sa_result = TD.saturation_adjustment(thermo_params, TD.ρe(), ρ, e_int, q_tot_safe) + return (; + T = sa_result.T, + q_tot_safe, + q_liq_rai = sa_result.q_liq, + q_ice_sno = sa_result.q_ice, ) - return (; e_int, q_pt = TD.PhasePartition(q_pt_args...)) end -ts_gs(thermo_params, moisture_model, microphysics_model, ᶜY, K, Φ, ρ) = - thermo_state( - thermo_params; - thermo_vars(moisture_model, microphysics_model, ᶜY, K, Φ)..., - ρ, - ) - function eddy_diffusivity_coefficient_H(D₀, H, z_sfc, z) return D₀ * exp(-(z - z_sfc) / H) end @@ -464,13 +449,12 @@ elsewhere, but doing it here ensures that it occurs whenever the precomputed quantities are updated. """ NVTX.@annotate function set_implicit_precomputed_quantities!(Y, p, t) - (; turbconv_model, moisture_model, microphysics_model) = p.atmos + (; turbconv_model, moisture_model) = p.atmos (; ᶜΦ) = p.core - (; ᶜu, ᶠu³, ᶠu, ᶜK, ᶜts, ᶜp) = p.precomputed + (; ᶜu, ᶠu³, ᶠu, ᶜK, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜh_tot, ᶜp) = p.precomputed ᶠuₕ³ = p.scratch.ᶠtemp_CT3 n = n_mass_flux_subdomains(turbconv_model) thermo_params = CAP.thermodynamics_params(p.params) - thermo_args = (thermo_params, moisture_model, microphysics_model) @. ᶠuₕ³ = $compute_ᶠuₕ³(Y.c.uₕ, Y.c.ρ) @@ -488,8 +472,8 @@ NVTX.@annotate function set_implicit_precomputed_quantities!(Y, p, t) # quantities of the form ᶜρaχ⁰ / ᶜρ⁰ and ᶜρaχʲ / ᶜρʲ to ᶜK, rather than # quantities of the form ᶜρaχ⁰ / ᶜρ and ᶜρaχʲ / ᶜρ. However, we cannot # compute ᶜρ⁰ and ᶜρʲ without first computing ᶜts⁰ and ᶜtsʲ, both of - # which depend on the value of ᶜp, which in turn depends on ᶜts. Since - # ᶜts depends on ᶜK, this + # which depend on the value of ᶜp, which in turn depends on ᶜT. Since + # ᶜT depends on ᶜK, this # means that the amount by which ᶜK needs to be incremented is a # function of ᶜK itself. So, unless we run a nonlinear solver here, this # circular dependency will prevent us from computing the exact value of @@ -498,8 +482,38 @@ NVTX.@annotate function set_implicit_precomputed_quantities!(Y, p, t) # @. ᶜK += Y.c.ρtke / Y.c.ρ # TODO: We should think more about these increments before we use them. end - @. ᶜts = ts_gs(thermo_args..., Y.c, ᶜK, ᶜΦ, Y.c.ρ) - @. ᶜp = TD.air_pressure(thermo_params, ᶜts) + if moisture_model isa EquilMoistModel + # Compute thermodynamic state variables using combined getter function. + # This avoids redundant saturation_adjustment calls for EquilMoistModel. + (; ᶜthermo_state) = p.precomputed + @. ᶜthermo_state = + thermo_state_gs(thermo_params, Y.c.ρ, Y.c.ρe_tot, Y.c.ρq_tot, ᶜK, ᶜΦ) + @. ᶜT = ᶜthermo_state.T + @. ᶜq_tot_safe = ᶜthermo_state.q_tot_safe + @. ᶜq_liq_rai = ᶜthermo_state.q_liq_rai + @. ᶜq_ice_sno = ᶜthermo_state.q_ice_sno + else # DryModel or NonEquilMoistModel + # For DryModel: q values are set to zero + # For NonEquilMoistModel: q values are computed from state variables + if moisture_model isa DryModel + @. ᶜq_tot_safe = zero(eltype(ᶜT)) + @. ᶜq_liq_rai = zero(eltype(ᶜT)) + @. ᶜq_ice_sno = zero(eltype(ᶜT)) + else # NonEquilMoistModel + @. ᶜq_tot_safe = max(0, specific(Y.c.ρq_tot, Y.c.ρ)) + @. ᶜq_liq_rai = + max(0, specific(Y.c.ρq_liq, Y.c.ρ) + specific(Y.c.ρq_rai, Y.c.ρ)) + @. ᶜq_ice_sno = + max(0, specific(Y.c.ρq_ice, Y.c.ρ) + specific(Y.c.ρq_sno, Y.c.ρ)) + end + ᶜe_int = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ) - ᶜK - ᶜΦ) + @. ᶜT = + TD.air_temperature(thermo_params, ᶜe_int, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + end + ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) + @. ᶜh_tot = + TD.total_enthalpy(thermo_params, ᶜe_tot, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + @. ᶜp = TD.air_pressure(thermo_params, ᶜT, Y.c.ρ, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) if turbconv_model isa PrognosticEDMFX set_prognostic_edmf_precomputed_quantities_draft!(Y, p, ᶠuₕ³, t) @@ -521,7 +535,6 @@ current state `Y`. This is only called before each evaluation of NVTX.@annotate function set_explicit_precomputed_quantities!(Y, p, t) (; turbconv_model, moisture_model, cloud_model) = p.atmos (; call_cloud_diagnostics_per_stage) = p.atmos - (; ᶜts) = p.precomputed thermo_params = CAP.thermodynamics_params(p.params) FT = eltype(p.params) @@ -530,10 +543,20 @@ NVTX.@annotate function set_explicit_precomputed_quantities!(Y, p, t) end if turbconv_model isa AbstractEDMF - @. p.precomputed.ᶜgradᵥ_q_tot = - ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts))) - @. p.precomputed.ᶜgradᵥ_θ_liq_ice = - ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts))) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed + @. p.precomputed.ᶜgradᵥ_q_tot = ᶜgradᵥ(ᶠinterp(ᶜq_tot_safe)) + @. p.precomputed.ᶜgradᵥ_θ_liq_ice = ᶜgradᵥ( + ᶠinterp( + TD.liquid_ice_pottemp( + thermo_params, + ᶜT, + Y.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ), + ), + ) end # The buoyancy gradient depends on the cloud fraction, and the cloud fraction diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index d39ee59e21..a76864730d 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -175,7 +175,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos FT = eltype(params) n = n_mass_flux_subdomains(turbconv_model) - (; ᶜu, ᶜp, ᶠu³, ᶜts, ᶠu³⁰, ᶜts⁰) = p.precomputed + (; ᶜu, ᶜp, ᶠu³, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno, ᶜts⁰) = p.precomputed (; ᶜlinear_buoygrad, ᶜstrain_rate_norm, ρtke_flux) = p.precomputed (; ᶜuʲs, @@ -192,10 +192,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos ᶜz = Fields.coordinate_field(Y.c).z z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half) - ᶜdz = Fields.Δz_field(axes(Y.c)) ᶜlg = Fields.local_geometry_field(Y.c) - ᶠlg = Fields.local_geometry_field(Y.f) - ᶜρa⁰ = @. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model)) ᶜtke = @. lazy(specific(Y.c.ρtke, Y.c.ρ)) ᶜvert_div = p.scratch.ᶜtemp_scalar @@ -301,10 +298,15 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_explicit_clos (; ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice, cloud_diagnostics_tuple) = p.precomputed + ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) @. ᶜlinear_buoygrad = buoyancy_gradients( # TODO - do we need to modify buoyancy gradients based on NonEq + 1M tracers? BuoyGradMean(), thermo_params, - ᶜts, + ᶜT, + Y.c.ρ, + ᶜq_tot, + ᶜq_liq_rai, + ᶜq_ice_sno, cloud_diagnostics_tuple.cf, C3, ᶜgradᵥ_q_tot, @@ -362,16 +364,22 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation n = n_mass_flux_subdomains(p.atmos.turbconv_model) for j in 1:n @. ᶜSqₜᵖʲs.:($$j) = q_tot_0M_precipitation_sources( - thp, cmp, dt, Y.c.sgsʲs.:($$j).q_tot, - ᶜtsʲs.:($$j), + TD.total_specific_humidity(thp, ᶜtsʲs.:($$j)), + TD.liquid_specific_humidity(thp, ᶜtsʲs.:($$j)), + TD.ice_specific_humidity(thp, ᶜtsʲs.:($$j)), ) end # sources from the environment ᶜq_tot⁰ = ᶜspecific_env_value(@name(q_tot), Y, p) - @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources(thp, cmp, dt, ᶜq_tot⁰, ᶜts⁰) + @. ᶜSqₜᵖ⁰ = q_tot_0M_precipitation_sources( + cmp, dt, ᶜq_tot⁰, + TD.total_specific_humidity(thp, ᶜts⁰), + TD.liquid_specific_humidity(thp, ᶜts⁰), + TD.ice_specific_humidity(thp, ᶜts⁰), + ) return nothing end NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation!( @@ -426,6 +434,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ) # Precipitation sources and sinks from the updrafts + Tʲ = @. lazy(TD.air_temperature(thp, ᶜtsʲs.:($$j))) compute_precipitation_sources!( ᶜSᵖ, ᶜSᵖ_snow, @@ -439,7 +448,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation Y.c.sgsʲs.:($j).q_ice, Y.c.sgsʲs.:($j).q_rai, Y.c.sgsʲs.:($j).q_sno, - ᶜtsʲs.:($j), + Tʲ, dt, cmp, thp, @@ -454,7 +463,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation Y.c.sgsʲs.:($j).q_ice, Y.c.sgsʲs.:($j).q_rai, Y.c.sgsʲs.:($j).q_sno, - ᶜtsʲs.:($j), + Tʲ, dt, cmp, thp, @@ -469,7 +478,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation Y.c.sgsʲs.:($$j).q_rai, Y.c.sgsʲs.:($$j).q_sno, ᶜρʲs.:($$j), - TD.air_temperature(thp, ᶜtsʲs.:($$j)), + Tʲ, dt, ) @. ᶜSqᵢᵖʲs.:($$j) += cloud_sources( @@ -481,7 +490,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation Y.c.sgsʲs.:($$j).q_rai, Y.c.sgsʲs.:($$j).q_sno, ᶜρʲs.:($$j), - TD.air_temperature(thp, ᶜtsʲs.:($$j)), + Tʲ, dt, ) end @@ -493,6 +502,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_rai⁰ = ᶜspecific_env_value(@name(q_rai), Y, p) ᶜq_sno⁰ = ᶜspecific_env_value(@name(q_sno), Y, p) ᶜρ⁰ = @. lazy(TD.air_density(thp, ᶜts⁰)) + T⁰ = @. lazy(TD.air_temperature(thp, ᶜts⁰)) compute_precipitation_sources!( ᶜSᵖ, ᶜSᵖ_snow, @@ -506,7 +516,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_ice⁰, ᶜq_rai⁰, ᶜq_sno⁰, - ᶜts⁰, + T⁰, dt, cmp, thp, @@ -521,7 +531,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_ice⁰, ᶜq_rai⁰, ᶜq_sno⁰, - ᶜts⁰, + T⁰, dt, cmp, thp, @@ -536,7 +546,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_rai⁰, ᶜq_sno⁰, ᶜρ⁰, - TD.air_temperature(thp, ᶜts⁰), + T⁰, dt, ) @. ᶜSqᵢᵖ⁰ += cloud_sources( @@ -548,7 +558,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_rai⁰, ᶜq_sno⁰, ᶜρ⁰, - TD.air_temperature(thp, ᶜts⁰), + T⁰, dt, ) return nothing @@ -665,6 +675,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ) # Precipitation sources and sinks from the updrafts + Tʲ = @. lazy(TD.air_temperature(thp, ᶜtsʲs.:($$j))) compute_warm_precipitation_sources_2M!( ᶜSᵖ, ᶜS₂ᵖ, @@ -680,7 +691,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation Y.c.sgsʲs.:($j).q_ice, Y.c.sgsʲs.:($j).q_rai, Y.c.sgsʲs.:($j).q_sno, - ᶜtsʲs.:($j), + Tʲ, dt, cm2p, thp, @@ -724,7 +735,8 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation max(0, w_component.(Geometry.WVector.(ᶜuʲs.:($$j)))), (cm2p,), thp, - ᶜtsʲs.:($$j), + TD.air_temperature(thp, ᶜtsʲs.:($$j)), + TD.air_pressure(thp, ᶜtsʲs.:($$j)), dt, ) end @@ -738,6 +750,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_rai⁰ = ᶜspecific_env_value(@name(q_rai), Y, p) ᶜq_sno⁰ = ᶜspecific_env_value(@name(q_sno), Y, p) ᶜρ⁰ = @. lazy(TD.air_density(thp, ᶜts⁰)) + T⁰ = @. lazy(TD.air_temperature(thp, ᶜts⁰)) compute_warm_precipitation_sources_2M!( ᶜSᵖ, ᶜS₂ᵖ, @@ -753,7 +766,7 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation ᶜq_ice⁰, ᶜq_rai⁰, ᶜq_sno⁰, - ᶜts⁰, + T⁰, dt, cm2p, thp, @@ -797,7 +810,8 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation w_component.(Geometry.WVector.(ᶜu⁰)), (cm2p,), thp, - ᶜts⁰, + TD.air_temperature(thp, ᶜts⁰), + TD.air_pressure(thp, ᶜts⁰), dt, ) return nothing diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index b4e92f2293..063ba0f831 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -109,13 +109,23 @@ end NVTX.@annotate function cloud_fraction_model_callback!(integrator) Y = integrator.u p = integrator.p - (; ᶜts, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = p.precomputed + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜgradᵥ_q_tot, ᶜgradᵥ_θ_liq_ice) = + p.precomputed thermo_params = CAP.thermodynamics_params(p.params) if isnothing(p.atmos.turbconv_model) - @. ᶜgradᵥ_q_tot = - ᶜgradᵥ(ᶠinterp(TD.total_specific_humidity(thermo_params, ᶜts))) - @. ᶜgradᵥ_θ_liq_ice = - ᶜgradᵥ(ᶠinterp(TD.liquid_ice_pottemp(thermo_params, ᶜts))) + @. ᶜgradᵥ_q_tot = ᶜgradᵥ(ᶠinterp(ᶜq_tot_safe)) + @. ᶜgradᵥ_θ_liq_ice = ᶜgradᵥ( + ᶠinterp( + TD.liquid_ice_pottemp( + thermo_params, + ᶜT, + Y.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ), + ), + ) end set_cloud_fraction!(Y, p, p.atmos.moisture_model, p.atmos.cloud_model) end diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 03dce1738e..cd6c59e473 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -21,7 +21,7 @@ # ::T, # ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} # thermo_params = CAP.thermodynamics_params(cache.params) -# out .= TD.relative_humidity.(thermo_params, cache.ᶜts) +# out .= TD.relative_humidity.(thermo_params, cache.precomputed.ᶜT, state.c.ρ, cache.precomputed.ᶜq_tot_safe, cache.precomputed.ᶜq_liq_rai, cache.precomputed.ᶜq_ice_sno)) # end # # 2. Define a function that has the correct signature and calls this function @@ -129,11 +129,10 @@ add_diagnostic_variable!( standard_name = "air_temperature", units = "K", compute! = (out, state, cache, time) -> begin - thermo_params = CAP.thermodynamics_params(cache.params) if isnothing(out) - return TD.air_temperature.(thermo_params, cache.precomputed.ᶜts) + return cache.precomputed.ᶜT else - out .= TD.air_temperature.(thermo_params, cache.precomputed.ᶜts) + out .= cache.precomputed.ᶜT end end, ) @@ -148,10 +147,26 @@ add_diagnostic_variable!( units = "K", compute! = (out, state, cache, time) -> begin thermo_params = CAP.thermodynamics_params(cache.params) + (; ᶜT, ᶜp, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = cache.precomputed if isnothing(out) - return TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts) + return TD.potential_temperature.( + thermo_params, + ᶜT, + ᶜp, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ) else - out .= TD.dry_pottemp.(thermo_params, cache.precomputed.ᶜts) + out .= + TD.potential_temperature.( + thermo_params, + ᶜT, + ᶜp, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ) end end, ) @@ -165,10 +180,11 @@ add_diagnostic_variable!( units = "m^2 s^-2", compute! = (out, state, cache, time) -> begin thermo_params = CAP.thermodynamics_params(cache.params) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = cache.precomputed if isnothing(out) - return TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜts) + return TD.enthalpy.(thermo_params, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) else - out .= TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜts) + out .= TD.enthalpy.(thermo_params, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) end end, ) @@ -399,10 +415,26 @@ function compute_hur!( moisture_model::T, ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} thermo_params = CAP.thermodynamics_params(cache.params) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = cache.precomputed if isnothing(out) - return TD.relative_humidity.(thermo_params, cache.precomputed.ᶜts) + return TD.relative_humidity.( + thermo_params, + ᶜT, + state.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ) else - out .= TD.relative_humidity.(thermo_params, cache.precomputed.ᶜts) + out .= + TD.relative_humidity.( + thermo_params, + ᶜT, + state.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ) end end @@ -559,18 +591,11 @@ function compute_hussfc!( time, moisture_model::T, ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} - thermo_params = CAP.thermodynamics_params(cache.params) + # q_vap_sfc is the total specific humidity at the surface (no liquid/ice) if isnothing(out) - return TD.total_specific_humidity.( - thermo_params, - cache.precomputed.sfc_conditions.ts, - ) + return copy(cache.precomputed.sfc_conditions.q_vap_sfc) else - out .= - TD.total_specific_humidity.( - thermo_params, - cache.precomputed.sfc_conditions.ts, - ) + out .= cache.precomputed.sfc_conditions.q_vap_sfc end end @@ -593,18 +618,10 @@ add_diagnostic_variable!( units = "K", comments = "Temperature of the lower boundary of the atmosphere", compute! = (out, state, cache, time) -> begin - thermo_params = CAP.thermodynamics_params(cache.params) if isnothing(out) - return TD.air_temperature.( - thermo_params, - cache.precomputed.sfc_conditions.ts, - ) + return copy(cache.precomputed.sfc_conditions.T_sfc) else - out .= - TD.air_temperature.( - thermo_params, - cache.precomputed.sfc_conditions.ts, - ) + out .= cache.precomputed.sfc_conditions.T_sfc end end, ) @@ -619,18 +636,11 @@ add_diagnostic_variable!( units = "K", comments = "Temperature at the bottom cell center of the atmosphere", compute! = (out, state, cache, time) -> begin - thermo_params = CAP.thermodynamics_params(cache.params) + T_level = Fields.level(cache.precomputed.ᶜT, 1) if isnothing(out) - return TD.air_temperature.( - thermo_params, - Fields.level(cache.precomputed.ᶜts, 1), - ) + return copy(T_level) else - out .= - TD.air_temperature.( - thermo_params, - Fields.level(cache.precomputed.ᶜts, 1), - ) + out .= T_level end end, ) @@ -1323,25 +1333,18 @@ add_diagnostic_variable!( ### function compute_dsevi!(out, state, cache, time) thermo_params = CAP.thermodynamics_params(cache.params) + ᶜT = cache.precomputed.ᶜT if isnothing(out) out = zeros(axes(Fields.level(state.f, half))) cp = CAP.cp_d(cache.params) dse = cache.scratch.ᶜtemp_scalar - @. dse = - state.c.ρ * ( - cp * TD.air_temperature(thermo_params, cache.precomputed.ᶜts) + - cache.core.ᶜΦ - ) + @. dse = state.c.ρ * (TD.dry_static_energy(thermo_params, ᶜT, cache.core.ᶜΦ)) Operators.column_integral_definite!(out, dse) return out else cp = CAP.cp_d(cache.params) dse = cache.scratch.ᶜtemp_scalar - @. dse = - state.c.ρ * ( - cp * TD.air_temperature(thermo_params, cache.precomputed.ᶜts) + - cache.core.ᶜΦ - ) + @. dse = state.c.ρ * (TD.dry_static_energy(thermo_params, ᶜT, cache.core.ᶜΦ)) Operators.column_integral_definite!(out, dse) end end @@ -1370,7 +1373,6 @@ function compute_clvi!( time, moisture_model::T, ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} - thermo_params = CAP.thermodynamics_params(cache.params) if isnothing(out) out = zeros(axes(Fields.level(state.f, half))) cloud_cover = cache.scratch.ᶜtemp_scalar @@ -1454,6 +1456,9 @@ function compute_hurvi!( moisture_model::T, ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} thermo_params = CAP.thermodynamics_params(cache.params) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = cache.precomputed + # Vapor specific humidity = q_tot - q_liq - q_ice + ᶜq_vap = @. lazy(ᶜq_tot_safe - ᶜq_liq_rai - ᶜq_ice_sno) if isnothing(out) out = zeros(axes(Fields.level(state.f, half))) # compute vertical integral of saturation specific humidity @@ -1462,16 +1467,11 @@ function compute_hurvi!( sat = cache.scratch.ᶜtemp_scalar @. sat = state.c.ρ * - TD.q_vap_saturation(thermo_params, cache.precomputed.ᶜts) + TD.q_vap_saturation(thermo_params, ᶜT, state.c.ρ, ᶜq_liq_rai, ᶜq_ice_sno) Operators.column_integral_definite!(sat_vi, sat) # compute saturation-weighted vertical integral of specific humidity - hur = cache.scratch.ᶜtemp_scalar - @. hur = TD.vapor_specific_humidity( - CAP.thermodynamics_params(cache.params), - cache.precomputed.ᶜts, - ) hur_weighted = cache.scratch.ᶜtemp_scalar_2 - @. hur_weighted = state.c.ρ * hur / sat_vi + @. hur_weighted = state.c.ρ * ᶜq_vap / sat_vi Operators.column_integral_definite!(out, hur_weighted) return out else @@ -1481,16 +1481,11 @@ function compute_hurvi!( sat = cache.scratch.ᶜtemp_scalar @. sat = state.c.ρ * - TD.q_vap_saturation(thermo_params, cache.precomputed.ᶜts) + TD.q_vap_saturation(thermo_params, ᶜT, state.c.ρ, ᶜq_liq_rai, ᶜq_ice_sno) Operators.column_integral_definite!(sat_vi, sat) # compute saturation-weighted vertical integral of specific humidity - hur = cache.scratch.ᶜtemp_scalar - @. hur = TD.vapor_specific_humidity( - CAP.thermodynamics_params(cache.params), - cache.precomputed.ᶜts, - ) hur_weighted = cache.scratch.ᶜtemp_scalar_2 - @. hur_weighted = state.c.ρ * hur / sat_vi + @. hur_weighted = state.c.ρ * ᶜq_vap / sat_vi Operators.column_integral_definite!(out, hur_weighted) end end @@ -1498,7 +1493,6 @@ end add_diagnostic_variable!( short_name = "hurvi", long_name = "Relative Humidity Saturation-Weighted Vertical Integral", - standard_name = "relative_humidity_vi", units = "kg m^-2", comments = "Integrated relative humidity over the vertical column", compute! = compute_hurvi!, @@ -1520,18 +1514,12 @@ function compute_husv!( time, moisture_model::T, ) where {T <: Union{EquilMoistModel, NonEquilMoistModel}} - thermo_params = CAP.thermodynamics_params(cache.params) + # Vapor specific humidity = q_tot - q_liq - q_ice + (; ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = cache.precomputed if isnothing(out) - return TD.vapor_specific_humidity.( - CAP.thermodynamics_params(cache.params), - cache.precomputed.ᶜts, - ) + return @. ᶜq_tot_safe - ᶜq_liq_rai - ᶜq_ice_sno else - out .= - TD.vapor_specific_humidity.( - CAP.thermodynamics_params(cache.params), - cache.precomputed.ᶜts, - ) + out .= @. ᶜq_tot_safe - ᶜq_liq_rai - ᶜq_ice_sno end end @@ -1552,7 +1540,6 @@ add_diagnostic_variable!( add_diagnostic_variable!( short_name = "uapredicted", long_name = "Predicted Eastward Wind", - standard_name = "predicted_eastward_wind", units = "m s^-1", comments = "Predicted steady-state eastward (zonal) wind component", compute! = (out, state, cache, time) -> begin @@ -1567,7 +1554,6 @@ add_diagnostic_variable!( add_diagnostic_variable!( short_name = "vapredicted", long_name = "Predicted Northward Wind", - standard_name = "predicted_northward_wind", units = "m s^-1", comments = "Predicted steady-state northward (meridional) wind component", compute! = (out, state, cache, time) -> begin @@ -1582,7 +1568,6 @@ add_diagnostic_variable!( add_diagnostic_variable!( short_name = "wapredicted", long_name = "Predicted Upward Air Velocity", - standard_name = "predicted_upward_air_velocity", units = "m s^-1", comments = "Predicted steady-state vertical wind component", compute! = (out, state, cache, time) -> begin @@ -1597,7 +1582,6 @@ add_diagnostic_variable!( add_diagnostic_variable!( short_name = "uaerror", long_name = "Error of Eastward Wind", - standard_name = "error_eastward_wind", units = "m s^-1", comments = "Error of steady-state eastward (zonal) wind component", compute! = (out, state, cache, time) -> begin @@ -1615,7 +1599,6 @@ add_diagnostic_variable!( add_diagnostic_variable!( short_name = "vaerror", long_name = "Error of Northward Wind", - standard_name = "error_northward_wind", units = "m s^-1", comments = "Error of steady-state northward (meridional) wind component", compute! = (out, state, cache, time) -> begin @@ -1633,7 +1616,6 @@ add_diagnostic_variable!( add_diagnostic_variable!( short_name = "waerror", long_name = "Error of Upward Air Velocity", - standard_name = "error_upward_air_velocity", units = "m s^-1", comments = "Error of steady-state vertical wind component", compute! = (out, state, cache, time) -> begin @@ -1655,12 +1637,22 @@ function compute_cape!(out, state, cache, time) thermo_params = lazy.(CAP.thermodynamics_params(cache.params)) g = lazy.(TD.Parameters.grav(thermo_params)) - # Get surface parcel properties - surface_thermal_state = lazy.(cache.precomputed.sfc_conditions.ts) - surface_q = - lazy.(TD.total_specific_humidity.(thermo_params, surface_thermal_state)) + # Get surface parcel properties from sfc_conditions + # At the surface, q_tot ≈ q_vap (no condensate) + surface_q = lazy.(cache.precomputed.sfc_conditions.q_vap_sfc) + surface_T = lazy.(cache.precomputed.sfc_conditions.T_sfc) + # Use lowest level pressure as approximate surface pressure + surface_p = lazy.(Fields.level(ᶠinterp.(cache.precomputed.ᶜp), half)) + # Compute liquid-ice potential temperature at surface (no condensate, so q_liq=q_ice=0) surface_θ_liq_ice = - lazy.(TD.liquid_ice_pottemp.(thermo_params, surface_thermal_state)) + lazy.( + TD.liquid_ice_pottemp.( + thermo_params, + surface_T, + surface_p, + TD.PhasePartition.(surface_q), + ), + ) # Create parcel thermodynamic states at each level based on energy & moisture at surface parcel_ts_moist = @@ -1675,8 +1667,11 @@ function compute_cape!(out, state, cache, time) # Calculate virtual temperatures for parcel & environment parcel_Tv = lazy.(TD.virtual_temperature.(thermo_params, parcel_ts_moist)) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = cache.precomputed env_Tv = - lazy.(TD.virtual_temperature.(thermo_params, cache.precomputed.ᶜts)) + lazy.( + TD.virtual_temperature.(thermo_params, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + ) # Calculate buoyancy from the difference in virtual temperatures ᶜbuoyancy = cache.scratch.ᶜtemp_scalar @@ -1713,11 +1708,20 @@ add_diagnostic_variable!( function compute_mslp!(out, state, cache, time) thermo_params = CAP.thermodynamics_params(cache.params) g = TD.Parameters.grav(thermo_params) - ts_level = Fields.level(cache.precomputed.ᶜts, 1) - R_m_surf = @. lazy(TD.gas_constant_air(thermo_params, ts_level)) + q_tot_safe_level = Fields.level(cache.precomputed.ᶜq_tot_safe, 1) + q_liq_rai_level = Fields.level(cache.precomputed.ᶜq_liq_rai, 1) + q_ice_sno_level = Fields.level(cache.precomputed.ᶜq_ice_sno, 1) + R_m_surf = @. lazy( + TD.gas_constant_air( + thermo_params, + q_tot_safe_level, + q_liq_rai_level, + q_ice_sno_level, + ), + ) p_level = Fields.level(cache.precomputed.ᶜp, 1) - t_level = @. lazy(TD.air_temperature(thermo_params, ts_level)) + t_level = Fields.level(cache.precomputed.ᶜT, 1) z_level = Fields.level(Fields.coordinate_field(state.c.ρ).z, 1) # Reduce to mean sea level using hypsometric formulation with lapse rate adjustment @@ -1795,17 +1799,11 @@ add_diagnostic_variable!( # Covariances (3d) ### function compute_covariance_diagnostics!(out, state, cache, time, type) - turbconv_model = cache.atmos.turbconv_model thermo_params = CAP.thermodynamics_params(cache.params) - if isa(turbconv_model, PrognosticEDMFX) - ᶜts = cache.precomputed.ᶜts⁰ - else - ᶜts = cache.precomputed.ᶜts - end # Reuse central compute_covariance function (ᶜq′q′, ᶜθ′θ′, ᶜθ′q′) = compute_covariance( - state, cache, thermo_params, ᶜts, + state, cache, thermo_params, ) result = if type == :qt_qt diff --git a/src/diagnostics/edmfx_diagnostics.jl b/src/diagnostics/edmfx_diagnostics.jl index bd6dfe51ef..d9843ab909 100644 --- a/src/diagnostics/edmfx_diagnostics.jl +++ b/src/diagnostics/edmfx_diagnostics.jl @@ -190,9 +190,9 @@ function compute_haup!( ) thermo_params = CAP.thermodynamics_params(cache.params) if isnothing(out) - return TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜtsʲs.:1) + return TD.enthalpy.(thermo_params, cache.precomputed.ᶜtsʲs.:1) else - out .= TD.specific_enthalpy.(thermo_params, cache.precomputed.ᶜtsʲs.:1) + out .= TD.enthalpy.(thermo_params, cache.precomputed.ᶜtsʲs.:1) end end diff --git a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl index 1fd87b8657..6ef1cb34f8 100644 --- a/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl +++ b/src/parameterized_tendencies/gravity_wave_drag/non_orographic_gravity_wave.jl @@ -139,8 +139,7 @@ function non_orographic_gravity_wave_compute_tendency!( ::NonOrographicGravityWave, ) #unpack - ᶜT = p.scratch.ᶜtemp_scalar - (; ᶜts) = p.precomputed + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed (; params) = p (; ᶜdTdz, @@ -163,12 +162,11 @@ function non_orographic_gravity_wave_compute_tendency!( grav = CAP.grav(params) # compute buoyancy frequency - @. ᶜT = TD.air_temperature(thermo_params, ᶜts) - ᶜdTdz .= Geometry.WVector.(ᶜgradᵥ.(ᶠinterp.(ᶜT))).components.data.:1 + ᶜcp_m = @. lazy(TD.cp_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno)) @. ᶜbuoyancy_frequency = - (grav / ᶜT) * (ᶜdTdz + grav / TD.cp_m(thermo_params, ᶜts)) + (grav / ᶜT) * (ᶜdTdz + grav / ᶜcp_m) ᶜbuoyancy_frequency = @. ifelse( ᶜbuoyancy_frequency < FT(2.5e-5), FT(sqrt(2.5e-5)), diff --git a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl index adf5025924..ab7186fc03 100644 --- a/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl +++ b/src/parameterized_tendencies/gravity_wave_drag/orographic_gravity_wave.jl @@ -78,9 +78,8 @@ function orographic_gravity_wave_cache(Y, ogw::OrographicGravityWave) end function orographic_gravity_wave_tendency!(Yₜ, Y, p, t, ::OrographicGravityWave) - ᶜT = p.scratch.ᶜtemp_scalar (; params) = p - (; ᶜts, ᶜp) = p.precomputed + (; ᶜp, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed (; ᶜdTdz) = p.orographic_gravity_wave (; topo_k_pbl, @@ -107,7 +106,6 @@ function orographic_gravity_wave_tendency!(Yₜ, Y, p, t, ::OrographicGravityWav ᶠz = Fields.coordinate_field(Y.f).z # get PBL info - @. ᶜT = TD.air_temperature(thermo_params, ᶜts) Fields.bycolumn(axes(Y.c.ρ)) do colidx parent(topo_k_pbl[colidx]) .= get_pbl(ᶜp[colidx], ᶜT[colidx], ᶜz[colidx], grav, cp_d) @@ -115,7 +113,9 @@ function orographic_gravity_wave_tendency!(Yₜ, Y, p, t, ::OrographicGravityWav # buoyancy frequency at cell centers parent(ᶜdTdz) .= parent(Geometry.WVector.(ᶜgradᵥ.(ᶠinterp.(ᶜT)))) - ᶜN = @. (grav / ᶜT) * (ᶜdTdz + grav / TD.cp_m(thermo_params, ᶜts)) # this is actually ᶜN^2 + ᶜcp_m = @. lazy(TD.cp_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno)) + #TODO: make a function to compute buoyancy frequency + ᶜN = @. (grav / ᶜT) * (ᶜdTdz + grav / ᶜcp_m) # this is actually ᶜN^2 @. ᶜN = ifelse(ᶜN < eps(FT), sqrt(eps(FT)), sqrt(abs(ᶜN))) # to avoid small numbers # prepare physical uv input variables for gravity_wave_forcing() diff --git a/src/parameterized_tendencies/les_sgs_models/anisotropic_minimum_dissipation.jl b/src/parameterized_tendencies/les_sgs_models/anisotropic_minimum_dissipation.jl index 88fa8cd058..6d054eb636 100644 --- a/src/parameterized_tendencies/les_sgs_models/anisotropic_minimum_dissipation.jl +++ b/src/parameterized_tendencies/les_sgs_models/anisotropic_minimum_dissipation.jl @@ -59,9 +59,7 @@ function horizontal_amd_tendency!(Yₜ, Y, p, t, les::AnisotropicMinimumDissipat (; atmos, precomputed, scratch, params) = p FT = eltype(Y) c_amd = les.c_amd - grav = CAP.grav(params) - thermo_params = CAP.thermodynamics_params(params) - (; ᶜu, ᶠu³, ᶜts) = precomputed + (; ᶜu, ᶠu³) = precomputed (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_strain, ᶠtemp_strain) = scratch (; ᶜtemp_scalar, ᶠtemp_scalar, ᶠtemp_scalar_2, ᶜtemp_UVW, ᶠtemp_UVW) = scratch @@ -128,8 +126,7 @@ function horizontal_amd_tendency!(Yₜ, Y, p, t, les::AnisotropicMinimumDissipat @. Yₜ.f.u₃ -= C3(wdivₕ(ᶠρ * ᶠτ_amd) / ᶠρ) ## Total energy tendency - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = precomputed ∇h_tot = @. lazy(Geometry.project(axis_uvw, gradₕ(ᶜh_tot))) ∂̂h_tot = @. lazy(Δ_h * ∇h_tot) ᶜD_amd = @. ᶜtemp_scalar = max( @@ -200,10 +197,6 @@ in `remaining_tendencies.jl` """ function vertical_amd_tendency!(Yₜ, Y, p, t, les::AnisotropicMinimumDissipation) FT = eltype(Y) - (; sfc_temp_C3) = p.scratch - (; ᶜts, sfc_conditions) = p.precomputed - (; ρ_flux_uₕ, ρ_flux_h_tot) = sfc_conditions - thermo_params = CAP.thermodynamics_params(p.params) c_amd = les.c_amd @@ -224,7 +217,7 @@ function vertical_amd_tendency!(Yₜ, Y, p, t, les::AnisotropicMinimumDissipatio ### AMD ### - (; ᶜu, ᶠu³, ᶜts) = p.precomputed + (; ᶜu, ᶠu³) = p.precomputed (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_strain, ᶠtemp_strain) = p.scratch (; ᶜtemp_scalar, ᶠtemp_scalar, ᶜtemp_UVW, ᶠtemp_UVW) = p.scratch @@ -301,8 +294,7 @@ function vertical_amd_tendency!(Yₜ, Y, p, t, les::AnisotropicMinimumDissipatio @. Yₜ.f.u₃ -= C3(ᶠdivᵥ(Y.c.ρ * ᶜτ_amd) / ᶠρ) ## Total energy tendency - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed # TODO: Fix @lazy broadcast (components access) ∇h_tot = @. lazy(Geometry.project(axis_uvw, ᶠgradᵥ_scalar(ᶜh_tot))) ∂̂h_tot = @. lazy(ᶠΔ_z * ∇h_tot) diff --git a/src/parameterized_tendencies/les_sgs_models/constant_horizontal_diffusion.jl b/src/parameterized_tendencies/les_sgs_models/constant_horizontal_diffusion.jl index f6e3766a94..6adc8ff67c 100644 --- a/src/parameterized_tendencies/les_sgs_models/constant_horizontal_diffusion.jl +++ b/src/parameterized_tendencies/les_sgs_models/constant_horizontal_diffusion.jl @@ -18,18 +18,11 @@ function horizontal_constant_diffusion_tendency!( FT = eltype(Y) thermo_params = CAP.thermodynamics_params(p.params) (; ᶜtemp_scalar) = p.scratch - (; ᶜts) = p.precomputed ᶜD = @. ᶜtemp_scalar = FT(chd.D) # Total energy diffusion - ᶜh_tot = @. lazy( - TD.total_specific_enthalpy( - thermo_params, - ᶜts, - specific(Y.c.ρe_tot, Y.c.ρ), - ), - ) + (; ᶜh_tot) = p.precomputed @. Yₜ.c.ρe_tot += wdivₕ(Y.c.ρ * ᶜD * gradₕ(ᶜh_tot)) # Tracer diffusion diff --git a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl index 7a63dfe195..dae7ba946b 100644 --- a/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl +++ b/src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl @@ -7,24 +7,34 @@ import ClimaCore.Operators as Operators import ClimaCore: Geometry """ - lilly_stratification_correction(p, ᶜS) + lilly_stratification_correction(Y, p, ᶜS) Return a lazy representation of the Lilly stratification correction factor based on the local Richardson number. # Arguments +- `Y`: The model state. - `p`: The model parameters, e.g. `AtmosCache`. - `ᶜS`: The cell-centered strain rate tensor. """ -function lilly_stratification_correction(p, ᶜS) - (; ᶜts) = p.precomputed +function lilly_stratification_correction(Y, p, ᶜS) + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed (; ᶜtemp_scalar) = p.scratch grav = CAP.grav(p.params) Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params)) thermo_params = CAP.thermodynamics_params(p.params) FT = eltype(Pr_t) # Stratification correction - ᶜθ_v = @. lazy(TD.virtual_pottemp(thermo_params, ᶜts)) + ᶜθ_v = @. lazy( + TD.virtual_pottemp( + thermo_params, + ᶜT, + Y.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ), + ) ᶜ∇ᵥθ = @. ᶜtemp_scalar = Geometry.WVector(ᶜgradᵥ(ᶠinterp(ᶜθ_v))).components.data.:1 ᶜN² = @. lazy(grav / ᶜθ_v * ᶜ∇ᵥθ) ᶜS_norm = strain_rate_norm(ᶜS, Geometry.WAxis()) @@ -72,7 +82,7 @@ function set_smagorinsky_lilly_precomputed_quantities!(Y, p, model) ax_xy = is_smagorinsky_UVW_coupled(model) ? Geometry.UVWAxis() : Geometry.UVAxis() ax_z = is_smagorinsky_UVW_coupled(model) ? Geometry.UVWAxis() : Geometry.WAxis() - ᶜfb = lilly_stratification_correction(p, ᶜS) + ᶜfb = lilly_stratification_correction(Y, p, ᶜS) if is_smagorinsky_UVW_coupled(model) ᶜL_h = ᶜL_v = @. lazy(c_smag * cbrt(Δx * Δy * ᶜΔz) * ᶜfb) else @@ -102,7 +112,7 @@ vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, ::Nothing) = nothing function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, model::SmagorinskyLilly) is_smagorinsky_horizontal(model) || return nothing - (; ᶜts, ᶜS, ᶠS, ᶜνₜ_h, ᶜD_h) = p.precomputed + (; ᶜS, ᶠS, ᶜνₜ_h, ᶜD_h) = p.precomputed (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶜtemp_scalar, ᶠtemp_scalar) = p.scratch thermo_params = CAP.thermodynamics_params(p.params) ᶜρ = Y.c.ρ @@ -120,8 +130,7 @@ function horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, model::Smagorinsk @. Yₜ.f.u₃ -= C3(wdivₕ(ᶠρ * ᶠτ_smag) / ᶠρ) ## Total energy tendency - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, ᶜρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed @. Yₜ.c.ρe_tot += wdivₕ(ᶜρ * ᶜD_h * gradₕ(ᶜh_tot)) ## Tracer diffusion and associated mass changes @@ -139,10 +148,9 @@ end function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, model::SmagorinskyLilly) is_smagorinsky_vertical(model) || return nothing FT = eltype(Y) - (; ᶜts, ᶜS, ᶠS, ᶜνₜ_v) = p.precomputed + (; ᶜS, ᶠS, ᶜνₜ_v) = p.precomputed (; ᶜtemp_UVWxUVW, ᶠtemp_UVWxUVW, ᶠtemp_scalar, ᶠtemp_scalar_2) = p.scratch Pr_t = CAP.Prandtl_number_0(CAP.turbconv_params(p.params)) - thermo_params = CAP.thermodynamics_params(p.params) ᶜρ = Y.c.ρ ᶠρ = @. ᶠtemp_scalar = ᶠinterp(ᶜρ) @@ -172,8 +180,7 @@ function vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, model::SmagorinskyL @. Yₜ.f.u₃ -= C3(ᶠdivᵥ(ᶜρ * ᶜτ_smag) / ᶠρ) ## Total energy tendency - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, ᶜρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρχ(-(ᶠρ * ᶠD_smag * ᶠgradᵥ(ᶜh_tot))) ## Tracer diffusion and associated mass changes diff --git a/src/parameterized_tendencies/microphysics/cloud_condensate.jl b/src/parameterized_tendencies/microphysics/cloud_condensate.jl index fd0680a52a..7c73b52160 100644 --- a/src/parameterized_tendencies/microphysics/cloud_condensate.jl +++ b/src/parameterized_tendencies/microphysics/cloud_condensate.jl @@ -29,14 +29,12 @@ function cloud_condensate_tendency!( ::Microphysics1Moment, _, ) - (; ᶜts) = p.precomputed + (; ᶜT) = p.precomputed (; params, dt) = p FT = eltype(params) thp = CAP.thermodynamics_params(params) cmc = CAP.microphysics_cloud_params(params) - Tₐ = @. lazy(TD.air_temperature(thp, ᶜts)) - @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources( cmc.liquid, @@ -47,7 +45,7 @@ function cloud_condensate_tendency!( specific(Y.c.ρq_rai, Y.c.ρ), specific(Y.c.ρq_sno, Y.c.ρ), Y.c.ρ, - Tₐ, + ᶜT, dt, ) @. Yₜ.c.ρq_ice += @@ -60,7 +58,7 @@ function cloud_condensate_tendency!( specific(Y.c.ρq_rai, Y.c.ρ), specific(Y.c.ρq_sno, Y.c.ρ), Y.c.ρ, - Tₐ, + ᶜT, dt, ) end @@ -73,13 +71,11 @@ function cloud_condensate_tendency!( ::Microphysics2Moment, _, ) - (; ᶜts, ᶜu) = p.precomputed + (; ᶜT, ᶜp, ᶜu) = p.precomputed (; params, dt) = p thp = CAP.thermodynamics_params(params) cmp = CAP.microphysics_2m_params(params) - Tₐ = @. lazy(TD.air_temperature(thp, ᶜts)) - @. Yₜ.c.ρq_liq += Y.c.ρ * cloud_sources( cmp.liquid, @@ -90,7 +86,7 @@ function cloud_condensate_tendency!( specific(Y.c.ρq_rai, Y.c.ρ), specific(Y.c.ρq_sno, Y.c.ρ), Y.c.ρ, - Tₐ, + ᶜT, dt, ) @. Yₜ.c.ρq_ice += @@ -103,7 +99,7 @@ function cloud_condensate_tendency!( specific(Y.c.ρq_rai, Y.c.ρ), specific(Y.c.ρq_sno, Y.c.ρ), Y.c.ρ, - Tₐ, + ᶜT, dt, ) @@ -137,7 +133,8 @@ function cloud_condensate_tendency!( w_component.(Geometry.WVector.(ᶜu)), (cmp,), thp, - ᶜts, + ᶜT, + ᶜp, dt, ) end diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index 5f6366ce3a..7cc790adb8 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -12,9 +12,6 @@ import CloudMicrophysics.Parameters as CMP # Define some aliases and functions to make the code more readable const Tₐ = TD.air_temperature -const Pₐ = TD.air_pressure -const PP = TD.PhasePartition -const qᵥ = TD.vapor_specific_humidity # Clip any specific humidity function clip(q) @@ -29,7 +26,7 @@ function limit(q, dt, n::Int) end """ - cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ, dt) + cloud_sources(cm_params, thp, qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T, dt) - cm_params - CloudMicrophysics parameters struct for cloud water or ice condensate - thp - Thermodynamics parameters struct @@ -39,7 +36,7 @@ end - qᵣ - rain specific humidity - qₛ - snow specific humidity - ρ - air density - - Tₐ - air temperature + - T - air temperature - dt - model time step Returns the condensation/evaporation or deposition/sublimation rate for @@ -140,50 +137,51 @@ function cloud_sources( end """ - q_tot_0M_precipitation_sources(thp, cmp, dt, qₜ, ts) + q_tot_0M_precipitation_sources(cmp, dt, qₜ, q_tot, q_liq, q_ice) - - thp, cmp - structs with thermodynamic and microphysics parameters + - cmp - struct with microphysics parameters - dt - model time step - qₜ - total water specific humidity - - ts - thermodynamic state (see Thermodynamics.jl package for details) + - q_tot, q_liq, q_ice - phase partition quantities Returns the qₜ source term due to precipitation formation defined as Δm_tot / (m_dry + m_tot) for the 0-moment scheme """ -function q_tot_0M_precipitation_sources(thp, cmp::CMP.Parameters0M, dt, qₜ, ts) +function q_tot_0M_precipitation_sources(cmp::CMP.Parameters0M, dt, qₜ, q_tot, q_liq, q_ice) return -triangle_inequality_limiter( - -CM0.remove_precipitation(cmp, PP(thp, ts)), + -CM0.remove_precipitation(cmp, TD.PhasePartition(q_tot, q_liq, q_ice)), qₜ / dt, ) end """ - e_tot_0M_precipitation_sources_helper(thp, ts, Φ) + e_tot_0M_precipitation_sources_helper(thp, T, q_liq, q_ice, Φ) - thp - set with thermodynamics parameters - - ts - thermodynamic state (see td package for details) + - T - temperature + - q_liq, q_ice - liquid and ice specific humidities - Φ - geopotential Returns the total energy source term multiplier from precipitation formation for the 0-moment scheme """ -function e_tot_0M_precipitation_sources_helper(thp, ts, Φ) +function e_tot_0M_precipitation_sources_helper(thp, T, q_liq, q_ice, Φ) - λ = TD.liquid_fraction(thp, ts) - Iₗ = TD.internal_energy_liquid(thp, ts) - Iᵢ = TD.internal_energy_ice(thp, ts) + λ = TD.liquid_fraction(thp, T, q_liq, q_ice) + Iₗ = TD.internal_energy_liquid(thp, T) + Iᵢ = TD.internal_energy_ice(thp, T) return λ * Iₗ + (1 - λ) * Iᵢ + Φ end """ - compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₗᵖ, Sqᵢᵖ, Sqᵣᵖ, Sqₛᵖ, ρ, qₜ, qₗ, qᵢ, qᵣ, qₛ, ts, dt, mp, thp) + compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₗᵖ, Sqᵢᵖ, Sqᵣᵖ, Sqₛᵖ, ρ, qₜ, qₗ, qᵢ, qᵣ, qₛ, T, dt, mp, thp) - Sᵖ, Sᵖ_snow - temporary containters to help compute precipitation source terms - Sqₗᵖ, Sqᵢᵖ, Sqᵣᵖ, Sqₛᵖ - cached storage for precipitation source terms - ρ - air density - qₜ, qₗ, qᵢ, qᵣ, qₛ - total water, liquid and ice, rain and snow specific humidity - - ts - thermodynamic state (see td package for details) + - T - air temperature - dt - model time step - thp, cmp - structs with thermodynamic and microphysics parameters @@ -205,7 +203,7 @@ function compute_precipitation_sources!( qᵢ, qᵣ, qₛ, - ts, + T, dt, mp, thp, @@ -262,15 +260,15 @@ function compute_precipitation_sources!( ) # if T < T_freeze cloud droplets freeze to become snow # else the snow melts and both cloud water and snow become rain - α(thp, ts) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze) + α(thp, T) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, T) * (T - mp.ps.T_freeze) @. Sᵖ_snow = ifelse( - Tₐ(thp, ts) < mp.ps.T_freeze, + T < mp.ps.T_freeze, Sᵖ, - FT(-1) * triangle_inequality_limiter(Sᵖ * α(thp, ts), limit(qₛ, dt, 5)), + FT(-1) * triangle_inequality_limiter(Sᵖ * α(thp, T), limit(qₛ, dt, 5)), ) @. Sqₛᵖ += Sᵖ_snow @. Sqₗᵖ -= Sᵖ - @. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow) + @. Sqᵣᵖ += ifelse(T < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow) # accretion: q_ice + q_rai -> q_sno @. Sᵖ = triangle_inequality_limiter( @@ -291,7 +289,7 @@ function compute_precipitation_sources!( # accretion: q_rai + q_sno -> q_rai or q_sno @. Sᵖ = ifelse( - Tₐ(thp, ts) < mp.ps.T_freeze, + T < mp.ps.T_freeze, triangle_inequality_limiter( CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₛ, qᵣ, ρ), limit(qᵣ, dt, 5), @@ -309,13 +307,13 @@ function compute_precipitation_sources!( end """ - compute_precipitation_sinks!(Sᵖ, Sqᵣᵖ, Sqₛᵖ, ρ, qₜ, qₗ, qᵢ, qᵣ, qₛ, ts, dt, mp, thp) + compute_precipitation_sinks!(Sᵖ, Sqᵣᵖ, Sqₛᵖ, ρ, qₜ, qₗ, qᵢ, qᵣ, qₛ, T, dt, mp, thp) - Sᵖ - a temporary containter to help compute precipitation source terms - Sqᵣᵖ, Sqₛᵖ - cached storage for precipitation source terms - ρ - air density - qₜ, qₗ, qᵢ, qᵣ, qₛ - total water, cloud liquid and ice, rain and snow specific humidities - - ts - thermodynamic state (see td package for details) + - T - air temperature - dt - model time step - thp, cmp - structs with thermodynamic and microphysics parameters @@ -334,11 +332,12 @@ function compute_precipitation_sinks!( qᵢ, qᵣ, qₛ, - ts, + T, dt, mp, thp, ) + # TODO: maybe use q_tot_safe, q_liq_rai, q_ice_sno instead of qₜ, qₗ, qᵢ, qᵣ, qₛ FT = eltype(thp) sps = (mp.ps, mp.tv.snow, mp.aps, thp) rps = (mp.pr, mp.tv.rain, mp.aps, thp) @@ -346,15 +345,15 @@ function compute_precipitation_sinks!( #! format: off # evaporation: q_rai -> q_vap @. Sᵖ = -triangle_inequality_limiter( - -CM1.evaporation_sublimation(rps..., qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ(thp, ts)), + -CM1.evaporation_sublimation(rps..., qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T), limit(qᵣ, dt, 5), - limit(qᵥ(thp, ts), dt, 5), + limit(TD.vapor_specific_humidity(qₜ, qₗ + qᵢ, qᵣ + qₛ), dt, 5), ) @. Sqᵣᵖ += Sᵖ # melting: q_sno -> q_rai @. Sᵖ = triangle_inequality_limiter( - CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)), + CM1.snow_melt(sps..., qₛ, ρ, T), limit(qₛ, dt, 5), limit(qᵣ, dt, 5), ) @@ -362,11 +361,11 @@ function compute_precipitation_sinks!( @. Sqₛᵖ -= Sᵖ # deposition/sublimation: q_vap <-> q_sno - @. Sᵖ = CM1.evaporation_sublimation(sps..., qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, Tₐ(thp, ts)) + @. Sᵖ = CM1.evaporation_sublimation(sps..., qₜ, qₗ, qᵢ, qᵣ, qₛ, ρ, T) @. Sᵖ = ifelse( Sᵖ > FT(0), - triangle_inequality_limiter(Sᵖ, limit(qᵥ(thp, ts), dt, 5), limit(qₛ, dt, 5)), - -triangle_inequality_limiter(FT(-1) * Sᵖ, limit(qₛ, dt, 5), limit(qᵥ(thp, ts), dt, 5)), + triangle_inequality_limiter(Sᵖ, limit(TD.vapor_specific_humidity(qₜ, qₗ + qᵣ, qᵢ + qₛ), dt, 5), limit(qₛ, dt, 5)), + -triangle_inequality_limiter(FT(-1) * Sᵖ, limit(qₛ, dt, 5), limit(TD.vapor_specific_humidity(qₜ, qₗ + qᵣ, qᵢ + qₛ), dt, 5)), ) @. Sqₛᵖ += Sᵖ #! format: on @@ -463,7 +462,8 @@ end w, cmp, thermo_params, - ts, + T, + p, dt, ) @@ -484,7 +484,8 @@ velocity. The result is returned as a tendency (per second) of liquid droplet nu - `w`: Vertical velocity [m/s]. - `cmp`: Microphysics parameters - `thermo_params`: Thermodynamic parameters for computing saturation, pressure, temperature, etc. -- `ts`: Thermodynamic state (e.g., prognostic variables) used for evaluating the phase partition. +- `T`: Air temperature [K]. +- `p`: Air pressure [Pa]. - `dt`: Time step (s) over which the activation tendency is applied. # Returns @@ -502,7 +503,8 @@ function aerosol_activation_sources( w, cmp, thermo_params, - ts, + T, + p, dt, ) @@ -510,8 +512,6 @@ function aerosol_activation_sources( air_params = cmp.aps arg_params = cmp.arg aerosol_params = cmp.aerosol - T = Tₐ(thermo_params, ts) - p = Pₐ(thermo_params, ts) S = CMTDI.supersaturation_over_liquid(thermo_params, qₜ, qₗ, qᵢ, ρ, T) n_aer = seasalt_num + sulfate_num if S < 0 || n_aer < eps(FT) @@ -564,14 +564,14 @@ function aerosol_activation_sources( end """ - compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, nₗ, nᵣ, qₜ, qₗ, qᵢ, qᵣ, qₛ, ts, dt, sb, thp) + compute_warm_precipitation_sources_2M!(Sᵖ, S₂ᵖ, Snₗᵖ, Snᵣᵖ, Sqₗᵖ, Sqᵣᵖ, ρ, nₗ, nᵣ, qₜ, qₗ, qᵢ, qᵣ, qₛ, T, dt, sb, thp) - `Sᵖ`, `S₂ᵖ` - temporary containters to help compute precipitation source terms - `Snₗᵖ`, `Snᵣᵖ`, `Sqₗᵖ`, `Sqᵣᵖ` - cached storage for precipitation source terms - `ρ` - air density - `nₗ`, `nᵣ` - cloud liquid and rain number concentration per mass [1 / kg of moist air] - `qₜ`, `qₗ`, `qᵢ`, `qᵣ`, `qₛ` - total water, cloud liquid, cloud ice, rain and snow specific humidity - - `ts` - thermodynamic state (see td package for details) + - `T` - air temperature - `dt` - model time step - `thp`, `mp` - structs with thermodynamic and microphysics parameters @@ -593,7 +593,7 @@ function compute_warm_precipitation_sources_2M!( qᵢ, qᵣ, qₛ, - ts, + T, dt, mp, thp, @@ -690,10 +690,10 @@ function compute_warm_precipitation_sources_2M!( qₛ, ρ, ρ * nᵣ, - Tₐ(thp, ts), + T, ).evap_rate_1, limit(qᵣ, dt, 5), - limit(qᵥ(thp, ts), dt, 5), + limit(TD.vapor_specific_humidity(qₜ, qₗ + qᵣ, qᵢ + qₛ), dt, 5), ) @. Sqᵣᵖ += Sᵖ @@ -711,7 +711,7 @@ function compute_warm_precipitation_sources_2M!( qₛ, ρ, ρ * nᵣ, - Tₐ(thp, ts), + T, ).evap_rate_0 / ρ, limit(nᵣ, dt, 5), ) @@ -768,7 +768,7 @@ function compute_cold_precipitation_sources_P3!( params_2mp3, # Parameters for 2M and P3 schemes, see `get_microphysics_2m_p3_parameters` thermo_params, # An instance of `Thermodynamics.Parameters.ThermodynamicsParameters` ᶜY_reduced, # A reduced set of prognostic variables needed for P3 sources - ᶜts, # Thermodynamic state + ᶜT, # Air temperature from cache ᶜlogλ, # Logarithm of the P3 distribution slope parameter ) @@ -784,7 +784,7 @@ function compute_cold_precipitation_sources_P3!( warm.sb.pdf_c, warm.sb.pdf_r, ρq_liq, ρn_liq, ρq_rai, ρn_rai, warm.aps, thermo_params, (cold.velocity_params,), - ρ, Tₐ(thermo_params, ᶜts), + ρ, ᶜT, ) return nothing diff --git a/src/parameterized_tendencies/radiation/held_suarez.jl b/src/parameterized_tendencies/radiation/held_suarez.jl index 411cd263ba..5c9ea33352 100644 --- a/src/parameterized_tendencies/radiation/held_suarez.jl +++ b/src/parameterized_tendencies/radiation/held_suarez.jl @@ -44,15 +44,14 @@ end Base.Broadcast.broadcastable(x::HeldSuarezForcingParams) = tuple(x) function compute_ΔρT( - thermo_params::TDP.ThermodynamicsParameters, - ts_surf::TD.ThermodynamicState, + T_sfc::FT, ρ::FT, p::FT, lat::FT, z_surface::FT, s::HeldSuarezForcingParams, ) where {FT} - σ = compute_σ(thermo_params, z_surface, p, ts_surf, s) + σ = compute_σ(z_surface, p, T_sfc, s) k_a = 1 / (40 * s.day) k_s = 1 / (4 * s.day) @@ -71,36 +70,24 @@ function compute_ΔρT( end function compute_σ( - thermo_params::TDP.ThermodynamicsParameters, z_surface::FT, p::FT, - ts_surf::TD.ThermodynamicState, + T_sfc::FT, s::HeldSuarezForcingParams, ) where {FT} - p / ( - s.MSLP * exp( - -s.grav * z_surface / s.R_d / - TD.air_temperature(thermo_params, ts_surf), - ) - ) + p / (s.MSLP * exp(-s.grav * z_surface / s.R_d / T_sfc)) end height_factor(σ::FT, σ_b::FT) where {FT} = max(0, (σ - σ_b) / (1 - σ_b)) -height_factor( - thermo_params::TDP.ThermodynamicsParameters, - z_surface::FT, - p::FT, - ts_surf::TD.ThermodynamicState, - s::HeldSuarezForcingParams, -) where {FT} = - height_factor(compute_σ(thermo_params, z_surface, p, ts_surf, s), s.σ_b) +height_factor(z_surface::FT, p::FT, T_sfc::FT, s::HeldSuarezForcingParams) where {FT} = + height_factor(compute_σ(z_surface, p, T_sfc, s), s.σ_b) function held_suarez_forcing_tendency_ρe_tot( ᶜρ, ᶜuₕ, ᶜp, params, - ts_surf, + T_sfc, moisture_model, forcing, ) @@ -120,7 +107,6 @@ function held_suarez_forcing_tendency_ρe_tot( grav = FT(CAP.grav(params)) Δθ_z = FT(CAP.Δθ_z(params)) T_min = FT(CAP.T_min_hs(params)) - thermo_params = CAP.thermodynamics_params(params) σ_b = CAP.σ_b(params) k_f = 1 / day @@ -143,15 +129,7 @@ function held_suarez_forcing_tendency_ρe_tot( ) return @. lazy( - -compute_ΔρT( - thermo_params, - ts_surf, - ᶜρ, - ᶜp, - lat, - z_surface, - hs_params, - ) * cv_d, + -compute_ΔρT(T_sfc, ᶜρ, ᶜp, lat, z_surface, hs_params) * cv_d, ) end @@ -159,7 +137,7 @@ function held_suarez_forcing_tendency_uₕ( ᶜuₕ, ᶜp, params, - ts_surf, + T_sfc, moisture_model, forcing, ) @@ -177,7 +155,6 @@ function held_suarez_forcing_tendency_uₕ( grav = FT(CAP.grav(params)) Δθ_z = FT(CAP.Δθ_z(params)) T_min = FT(CAP.T_min_hs(params)) - thermo_params = CAP.thermodynamics_params(params) σ_b = CAP.σ_b(params) k_f = 1 / day @@ -199,10 +176,5 @@ function held_suarez_forcing_tendency_uₕ( MSLP, ) - return @. lazy( - -( - k_f * - height_factor(thermo_params, z_surface, ᶜp, ts_surf, hs_params) - ) * ᶜuₕ, - ) + return @. lazy(-(k_f * height_factor(z_surface, ᶜp, T_sfc, hs_params)) * ᶜuₕ) end diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl index c3ab8ebc60..ced7c02c63 100644 --- a/src/parameterized_tendencies/radiation/radiation.jl +++ b/src/parameterized_tendencies/radiation/radiation.jl @@ -462,21 +462,17 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationDYCOMS) (; params) = p (; ᶜκρq, ∫_0_∞_κρq, ᶠ∫_0_z_κρq, isoline_z_ρ_ρq, ᶠradiation_flux) = p.radiation - (; ᶜts) = p.precomputed - thermo_params = CAP.thermodynamics_params(params) + (; ᶜq_liq_rai) = p.precomputed cp_d = CAP.cp_d(params) FT = Spaces.undertype(axes(Y.c)) NT = NamedTuple{(:z, :ρ, :ρq_tot), NTuple{3, FT}} ᶜz = Fields.coordinate_field(Y.c).z ᶠz = Fields.coordinate_field(Y.f).z - # TODO: According to the paper, we should replace liquid_specific_humidity - # with TD.mixing_ratios(thermo_params, ᶜts).liq, but this wouldn't - # match the original code from TurbulenceConvection. + # TODO: According to the paper, we should replace ᶜq_liq_rai + # with mixing ratio. @. ᶜκρq = - radiation_mode.kappa * - Y.c.ρ * - TD.liquid_specific_humidity(thermo_params, ᶜts) + radiation_mode.kappa * Y.c.ρ * ᶜq_liq_rai Operators.column_integral_definite!(∫_0_∞_κρq, ᶜκρq) @@ -540,10 +536,11 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationTRMM_LBA) thermo_params = CAP.thermodynamics_params(params) ᶜdTdt_rad = p.radiation.ᶜdTdt_rad ᶜρ = Y.c.ρ - ᶜts_gm = p.precomputed.ᶜts + (; ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed zc = Fields.coordinate_field(axes(ᶜρ)).z @. ᶜdTdt_rad = rad(FT(t), zc) - @. Yₜ.c.ρe_tot += ᶜρ * TD.cv_m(thermo_params, ᶜts_gm) * ᶜdTdt_rad + @. Yₜ.c.ρe_tot += + ᶜρ * TD.cv_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) * ᶜdTdt_rad return nothing end @@ -555,11 +552,10 @@ radiation_model_cache(Y, radiation_mode::RadiationISDAC; args...) = (;) # Don't function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationISDAC) (; F₀, F₁, κ) = radiation_mode (; params, precomputed) = p - (; ᶜts) = precomputed - thermo_params = CAP.thermodynamics_params(params) + (; ᶜq_liq_rai) = precomputed ᶜρq = p.scratch.ᶜtemp_scalar - @. ᶜρq = Y.c.ρ * TD.liquid_specific_humidity(thermo_params, ᶜts) + @. ᶜρq = Y.c.ρ * ᶜq_liq_rai LWP_zₜ = p.scratch.temp_field_level # column integral of LWP (zₜ = top-of-domain) Operators.column_integral_definite!(LWP_zₜ, ᶜρq) diff --git a/src/parameterized_tendencies/radiation/update_inputs.jl b/src/parameterized_tendencies/radiation/update_inputs.jl index 57567f1ba1..6cc06969ea 100644 --- a/src/parameterized_tendencies/radiation/update_inputs.jl +++ b/src/parameterized_tendencies/radiation/update_inputs.jl @@ -47,21 +47,19 @@ function update_temperature_pressure!((; u, p, t)::I) where {I} T_min = CAP.optics_lookup_temperature_min(params) T_max = CAP.optics_lookup_temperature_max(params) - (; ᶜts, ᶜp, sfc_conditions) = p.precomputed + (; ᶜT, sfc_conditions) = p.precomputed model = p.radiation.rrtmgp_model # update surface temperature - sfc_ts = sfc_conditions.ts - sfc_T = Fields.array2field(model.surface_temperature, axes(sfc_ts)) - @. sfc_T = TD.air_temperature(thermo_params, sfc_ts) + sfc_T = Fields.array2field(model.surface_temperature, axes(sfc_conditions.T_sfc)) + @. sfc_T = sfc_conditions.T_sfc # update layer pressure model.center_pressure .= Fields.field2array(p.precomputed.ᶜp) - # compute layer temperature - ᶜT = Fields.array2field(model.center_temperature, axes(u.c)) + # compute layer temperature (clamped to RRTMGP bounds) + rrtmgp_ᶜT = Fields.array2field(model.center_temperature, axes(u.c)) # TODO: move this to RRTMGP - @. ᶜT = - min(max(TD.air_temperature(thermo_params, ᶜts), FT(T_min)), FT(T_max)) + @. rrtmgp_ᶜT = min(max(ᶜT, FT(T_min)), FT(T_max)) return nothing end @@ -75,7 +73,7 @@ function update_relative_humidity!((; u, p, t)::I) where {I} (; rrtmgp_model) = p.radiation thermo_params = CAP.thermodynamics_params(p.params) FT = eltype(thermo_params) - (; ᶜts) = p.precomputed + (; ᶜT, ᶜp, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed ᶜrh = Fields.array2field(rrtmgp_model.center_relative_humidity, axes(u.c)) ᶜvmr_h2o = Fields.array2field( rrtmgp_model.center_volume_mixing_ratio_h2o, @@ -94,7 +92,8 @@ function update_relative_humidity!((; u, p, t)::I) where {I} # temporarily store ᶜq_tot in ᶜvmr_h2o ᶜq_tot = ᶜvmr_h2o @. ᶜq_tot = - max_relative_humidity * TD.q_vap_saturation(thermo_params, ᶜts) + max_relative_humidity * + TD.q_vap_saturation(thermo_params, ᶜT, u.c.ρ, ᶜq_liq_rai, ᶜq_ice_sno) # filter ᶜq_tot so that it is monotonically decreasing with z for i in 2:Spaces.nlevels(axes(ᶜq_tot)) @@ -109,9 +108,22 @@ function update_relative_humidity!((; u, p, t)::I) where {I} else @. ᶜvmr_h2o = TD.vol_vapor_mixing_ratio( thermo_params, - TD.PhasePartition(thermo_params, ᶜts), + TD.PhasePartition(ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno), + ) + @. ᶜrh = min( + max( + TD.relative_humidity( + thermo_params, + ᶜT, + ᶜp, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ), + 0, + ), + 1, ) - @. ᶜrh = min(max(TD.relative_humidity(thermo_params, ᶜts), 0), 1) end return nothing diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 3ff710a3db..e67ad8cd47 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -36,7 +36,7 @@ Modifies `Yₜ.c.ρ`, `Yₜ.c.ρe_tot`, `Yₜ.c.uₕ`, and EDMFX-related fields NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t) n = n_mass_flux_subdomains(p.atmos.turbconv_model) (; ᶜΦ) = p.core - (; ᶜu, ᶜK, ᶜp, ᶜts) = p.precomputed + (; ᶜu, ᶜK, ᶜp, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed (; params) = p thermo_params = CAP.thermodynamics_params(params) cp_d = thermo_params.cp_d @@ -55,8 +55,7 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t) end end - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed @. Yₜ.c.ρe_tot -= split_divₕ(Y.c.ρ * ᶜu, ᶜh_tot) if p.atmos.turbconv_model isa PrognosticEDMFX @@ -72,10 +71,12 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t) @. Yₜ.c.ρtke -= split_divₕ(Y.c.ρ * ᶜu, ᶜtke) end - ᶜΦ_r = @. lazy(phi_r(thermo_params, ᶜts)) - ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts)) - ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜts)) - ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts)) + (; ᶜq_tot_safe) = p.precomputed + ᶜΦ_r = @. lazy(phi_r(thermo_params, ᶜp)) + ᶜθ_v = p.scratch.ᶜtemp_scalar + @. ᶜθ_v = theta_v(thermo_params, ᶜT, ᶜp, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜp)) + ᶜΠ = @. lazy(TD.exner_given_pressure(thermo_params, ᶜp)) ᶜθ_v_diff = @. lazy(ᶜθ_v - ᶜθ_vr) # split form pressure gradient: 0.5 * cp_d * [θv ∇Π + ∇(θv Π) - Π∇θv] @. Yₜ.c.uₕ -= C12( @@ -220,8 +221,8 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) point_type = eltype(Fields.coordinate_field(Y.c)) (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J - (; ᶜf³, ᶠf¹², ᶜΦ) = p.core - (; ᶜu, ᶠu³, ᶜK, ᶜts) = p.precomputed + (; ᶜf³, ᶠf¹²) = p.core + (; ᶜu, ᶠu³, ᶜK) = p.precomputed (; edmfx_mse_q_tot_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing (; ᶜuʲs, ᶜKʲs, ᶠKᵥʲs) = n > 0 ? p.precomputed : all_nothing (; energy_q_tot_upwinding, tracer_upwinding) = p.atmos.numerics @@ -270,8 +271,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) # ... and upwinding correction of energy and total water. # (The central advection of energy and total water is done implicitly.) if energy_q_tot_upwinding != Val(:none) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed vtt = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, dt, energy_q_tot_upwinding) vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, dt, Val(:none)) @. Yₜ.c.ρe_tot += vtt - vtt_central diff --git a/src/prognostic_equations/eddy_diffusion_closures.jl b/src/prognostic_equations/eddy_diffusion_closures.jl index f9a6eb9d82..1ca959e0f1 100644 --- a/src/prognostic_equations/eddy_diffusion_closures.jl +++ b/src/prognostic_equations/eddy_diffusion_closures.jl @@ -14,7 +14,12 @@ import SurfaceFluxes.UniversalFunctions as UF thermo_params, # Arguments for the first method (most commonly called): - ts::TD.ThermodynamicState, + T, # Air temperature [K] + ρ, # Air density [kg/m³] + q_tot, # Total specific humidity [kg/kg] + q_liq, # Liquid specific humidity [kg/kg] + q_ice, # Ice specific humidity [kg/kg] + cf, # Cloud fraction ::Type{C3}, # Covariant3 vector type, for projecting gradients ∂qt∂z::AbstractField, # Vertical gradient of total specific humidity ∂θli∂z::AbstractField, # Vertical gradient of liquid-ice potential temperature @@ -36,23 +41,22 @@ fraction. The calculation involves: thermodynamic variables (`∂θᵥ/∂z`, `∂θₗᵢ/∂z`, `∂qₜ/∂z`), obtained from the input fields after projection. 3. Blending the resulting unsaturated and saturated buoyancy gradients based on - the environmental cloud fraction derived from the thermodynamic state `ts`. - -The dispatch on `EnvBuoyGradVars` (used internally by the first method after -constructing it from `ts` and projected gradients) occurs during its construction. -The analytical solutions employed are consistent for both mean fields and -conditional fields derived from assumed distributions over conserved thermodynamic -variables. + the environmental cloud fraction. Arguments: - `closure`: The environmental buoyancy gradient closure type (e.g., `BuoyGradMean`). - `thermo_params`: Thermodynamic parameters from `CLIMAParameters`. -- `ts`: Center-level thermodynamic state of the environment. +- `T`: Air temperature [K] +- `ρ`: Air density [kg/m³] +- `q_tot`: Total specific humidity [kg/kg] +- `q_liq`: Liquid specific humidity [kg/kg] +- `q_ice`: Ice specific humidity [kg/kg] +- `cf`: Cloud fraction - `C3`: The `ClimaCore.Geometry.Covariant3Vector` type, used for projecting input vertical gradients. - `∂qt∂z`: Field of vertical gradients of total specific humidity. - `∂θli∂z`: Field of vertical gradients of liquid-ice potential temperature. - `local_geometry`: Field of local geometry at cell centers, used for gradient projection. -The second method takes a precomputed `EnvBuoyGradVars` object instead of `ts` and gradient fields. +The second method takes a precomputed `EnvBuoyGradVars` object instead of T, ρ, q_tot, q_liq, q_ice and gradient fields. Returns: - `∂b∂z`: The mean vertical buoyancy gradient [s⁻²], as a field of scalars. @@ -62,7 +66,11 @@ function buoyancy_gradients end function buoyancy_gradients( ebgc::AbstractEnvBuoyGradClosure, thermo_params, - ts, + T, + ρ, + q_tot, + q_liq, + q_ice, cf, ::Type{C3}, ∂qt∂z, @@ -73,7 +81,11 @@ function buoyancy_gradients( ebgc, thermo_params, EnvBuoyGradVars( - ts, + T, + ρ, + max(q_tot, 0), + max(q_liq, 0), + max(q_ice, 0), cf, projected_vector_buoy_grad_vars( C3, @@ -96,18 +108,13 @@ function buoyancy_gradients( Rv_over_Rd = TDP.Rv_over_Rd(thermo_params) R_v = TDP.R_v(thermo_params) - ts = bg_model.ts - ∂b∂θv = g / TD.virtual_pottemp(thermo_params, ts) - - T = TD.air_temperature(thermo_params, ts) - λ = TD.liquid_fraction(thermo_params, ts) - lh = - λ * TD.latent_heat_vapor(thermo_params, T) + - (1 - λ) * TD.latent_heat_sublim(thermo_params, T) - cp_m = TD.cp_m(thermo_params, ts) - q_sat = TD.q_vap_saturation(thermo_params, ts) - q_tot = TD.total_specific_humidity(thermo_params, ts) - θ = TD.dry_pottemp(thermo_params, ts) + (; T, ρ, q_tot, q_liq, q_ice) = bg_model + ∂b∂θv = g / TD.virtual_pottemp(thermo_params, T, ρ, q_tot, q_liq, q_ice) + + lh = TD.latent_heat(thermo_params, T, q_liq, q_ice) + cp_m = TD.cp_m(thermo_params, q_tot, q_liq, q_ice) + q_sat = TD.q_vap_saturation(thermo_params, T, ρ, q_liq, q_ice) + θ = TD.potential_temperature(thermo_params, T, ρ, q_tot, q_liq, q_ice) ∂b∂θli_unsat = ∂b∂θv * (1 + (Rv_over_Rd - 1) * q_tot) ∂b∂qt_unsat = ∂b∂θv * (Rv_over_Rd - 1) * θ ∂b∂θli_sat = ( diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl index 3e82cd1b7e..8f99a4ab76 100644 --- a/src/prognostic_equations/edmfx_entr_detr.jl +++ b/src/prognostic_equations/edmfx_entr_detr.jl @@ -255,7 +255,10 @@ function detrainment_from_thermo_state( u³ʲ_prev_halflevel, local_geometry_prev_halflevel, u³_prev_halflevel, - ts_prev_level, + T_prev_level, + q_tot_prev_level, + q_liq_prev_level, + q_ice_prev_level, ᶜbuoy⁰, entrʲ_prev_level, vert_div_level, @@ -284,7 +287,14 @@ function detrainment_from_thermo_state( local_geometry_prev_halflevel, ), get_physical_w(u³_prev_halflevel, local_geometry_prev_halflevel), - TD.relative_humidity(thermo_params, ts_prev_level), + TD.relative_humidity( + thermo_params, + T_prev_level, + p_prev_level, + q_tot_prev_level, + q_liq_prev_level, + q_ice_prev_level, + ), FT(0), entrʲ_prev_level, vert_div_level, @@ -601,7 +611,7 @@ function edmfx_first_interior_entr_tendency!( ) (; params, dt) = p - (; ᶜK, ᶜρʲs, ᶜentrʲs, ᶜts) = p.precomputed + (; ᶜK, ᶜρʲs, ᶜentrʲs) = p.precomputed FT = eltype(params) n = n_mass_flux_subdomains(p.atmos.turbconv_model) @@ -637,8 +647,7 @@ function edmfx_first_interior_entr_tendency!( Fields.local_geometry_field(Fields.level(Y.f, Fields.half)), ) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed ᶜh_tot_int_val = Fields.field_values(Fields.level(ᶜh_tot, 1)) ᶜK_int_val = Fields.field_values(Fields.level(ᶜK, 1)) ᶜmse⁰ = ᶜspecific_env_mse(Y, p) diff --git a/src/prognostic_equations/edmfx_precipitation.jl b/src/prognostic_equations/edmfx_precipitation.jl index 03097de319..62d842c98c 100644 --- a/src/prognostic_equations/edmfx_precipitation.jl +++ b/src/prognostic_equations/edmfx_precipitation.jl @@ -48,11 +48,16 @@ function edmfx_precipitation_tendency!( @. Yₜ.c.sgsʲs.:($$j).ρa += Y.c.sgsʲs.:($$j).ρa * ᶜSqₜᵖʲs.:($$j) + ᶜTʲ = @. lazy(TD.air_temperature(thermo_params, ᶜtsʲs.:($$j))) + ᶜq_liq_raiʲ = @. lazy(TD.liquid_specific_humidity(thermo_params, ᶜtsʲs.:($$j))) + ᶜq_ice_snoʲ = @. lazy(TD.ice_specific_humidity(thermo_params, ᶜtsʲs.:($$j))) @. Yₜ.c.sgsʲs.:($$j).mse += ᶜSqₜᵖʲs.:($$j) * ( e_tot_0M_precipitation_sources_helper( thermo_params, - ᶜtsʲs.:($$j), + ᶜTʲ, + ᶜq_liq_raiʲ, + ᶜq_ice_snoʲ, ᶜΦ, ) - TD.internal_energy(thermo_params, ᶜtsʲs.:($$j)) ) diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl index 4692a8efc6..ffc2f841aa 100644 --- a/src/prognostic_equations/edmfx_sgs_flux.jl +++ b/src/prognostic_equations/edmfx_sgs_flux.jl @@ -42,7 +42,7 @@ function edmfx_sgs_mass_flux_tendency!( (; edmfx_sgsflux_upwinding, edmfx_tracer_upwinding) = p.atmos.numerics (; ᶠu³) = p.precomputed (; ᶠu³ʲs, ᶜKʲs, ᶜρʲs) = p.precomputed - (; ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜts) = p.precomputed + (; ᶠu³⁰, ᶜK⁰, ᶜts⁰) = p.precomputed (; dt) = p thermo_params = CAP.thermodynamics_params(p.params) @@ -57,8 +57,7 @@ function edmfx_sgs_mass_flux_tendency!( # [best after removal of precomputed quantities] ᶠu³_diff = p.scratch.ᶠtemp_CT3 ᶜa_scalar = p.scratch.ᶜtemp_scalar - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed for j in 1:n @. ᶠu³_diff = ᶠu³ʲs.:($$j) - ᶠu³ @. ᶜa_scalar = @@ -189,7 +188,7 @@ function edmfx_sgs_mass_flux_tendency!( n = n_mass_flux_subdomains(turbconv_model) (; edmfx_sgsflux_upwinding, edmfx_tracer_upwinding) = p.atmos.numerics (; ᶠu³) = p.precomputed - (; ᶜρaʲs, ᶜρʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs, ᶜts) = p.precomputed + (; ᶜρaʲs, ᶜρʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs) = p.precomputed (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J FT = eltype(Y) @@ -197,8 +196,7 @@ function edmfx_sgs_mass_flux_tendency!( if p.atmos.edmfx_model.sgs_mass_flux isa Val{true} thermo_params = CAP.thermodynamics_params(p.params) # energy - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed ᶠu³_diff = p.scratch.ᶠtemp_CT3 ᶜa_scalar = p.scratch.ᶜtemp_scalar for j in 1:n @@ -378,8 +376,7 @@ function edmfx_sgs_diffusive_flux_tendency!( FT = Spaces.undertype(axes(Y.c)) (; dt, params) = p turbconv_params = CAP.turbconv_params(params) - thermo_params = CAP.thermodynamics_params(params) - (; ᶜu, ᶜts) = p.precomputed + (; ᶜu) = p.precomputed (; ρtke_flux) = p.precomputed ᶠgradᵥ = Operators.GradientC2F() ᶜtke = @. lazy(specific(Y.c.ρtke, Y.c.ρ)) @@ -412,8 +409,7 @@ function edmfx_sgs_diffusive_flux_tendency!( top = Operators.SetValue(C3(FT(0))), bottom = Operators.SetValue(C3(FT(0))), ) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠρaK_h * ᶠgradᵥ(ᶜh_tot))) if use_prognostic_tke(turbconv_model) diff --git a/src/prognostic_equations/forcing/external_forcing.jl b/src/prognostic_equations/forcing/external_forcing.jl index 94973278d7..b9b65471d5 100644 --- a/src/prognostic_equations/forcing/external_forcing.jl +++ b/src/prognostic_equations/forcing/external_forcing.jl @@ -318,7 +318,7 @@ function external_forcing_tendency!( # horizontal advection, vertical fluctuation, nudging, subsidence (need to add), (; params) = p thermo_params = CAP.thermodynamics_params(params) - (; ᶜts) = p.precomputed + (; ᶜT) = p.precomputed (; ᶜdTdt_fluc, ᶜdqtdt_fluc, @@ -338,13 +338,12 @@ function external_forcing_tendency!( @. ᶜuₕ_nudge = C12(Geometry.UVVector(ᶜu_nudge, ᶜv_nudge), ᶜlg) @. Yₜ.c.uₕ -= (Y.c.uₕ - ᶜuₕ_nudge) * ᶜinv_τ_wind - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed # nudging tendency ᶜdTdt_nudging = p.scratch.ᶜtemp_scalar ᶜdqtdt_nudging = p.scratch.ᶜtemp_scalar_2 @. ᶜdTdt_nudging = - -(TD.air_temperature(thermo_params, ᶜts) - ᶜT_nudge) * ᶜinv_τ_scalar + -(ᶜT - ᶜT_nudge) * ᶜinv_τ_scalar @. ᶜdqtdt_nudging = -(specific(Y.c.ρq_tot, Y.c.ρ) - ᶜqt_nudge) * ᶜinv_τ_scalar @@ -360,9 +359,9 @@ function external_forcing_tendency!( # total energy @. Yₜ.c.ρe_tot += Y.c.ρ * ( - TD.cv_m(thermo_params, ᶜts) * ᶜdTdt_sum + + TD.cv_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) * ᶜdTdt_sum + ( - cv_v * (TD.air_temperature(thermo_params, ᶜts) - T_0) + Lv_0 - + cv_v * (ᶜT - T_0) + Lv_0 - R_v * T_0 ) * ᶜdqtdt_sum ) @@ -561,7 +560,7 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::ISDACForcing) FT = Spaces.undertype(axes(Y.c)) (; params) = p thermo_params = CAP.thermodynamics_params(params) - (; ᶜts, ᶜp) = p.precomputed + (; ᶜp, ᶜT) = p.precomputed ᶜinv_τ_scalar = APL.ISDAC_inv_τ_scalar(FT) # s⁻¹ ᶜinv_τ_wind = APL.ISDAC_inv_τ_wind(FT) # s⁻¹ @@ -588,7 +587,7 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::ISDACForcing) ᶜdTdt_nudging = p.scratch.ᶜtemp_scalar ᶜdqtdt_nudging = p.scratch.ᶜtemp_scalar_2 @. ᶜdTdt_nudging = - -(TD.air_temperature(thermo_params, ᶜts) - ta_ISDAC(ᶜp, ᶜz)) * + -(ᶜT - ta_ISDAC(ᶜp, ᶜz)) * ᶜinv_τ_scalar(ᶜz) @. ᶜdqtdt_nudging = -(specific(Y.c.ρq_tot, Y.c.ρ) - q_tot(ᶜz)) * ᶜinv_τ_scalar(ᶜz) @@ -598,11 +597,12 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::ISDACForcing) cv_v = TD.Parameters.cv_v(thermo_params) R_v = TD.Parameters.R_v(thermo_params) # total energy + (; ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed @. Yₜ.c.ρe_tot += Y.c.ρ * ( - TD.cv_m(thermo_params, ᶜts) * ᶜdTdt_nudging + + TD.cv_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) * ᶜdTdt_nudging + ( - cv_v * (TD.air_temperature(thermo_params, ᶜts) - T_0) + Lv_0 - + cv_v * (ᶜT - T_0) + Lv_0 - R_v * T_0 ) * ᶜdqtdt_nudging ) diff --git a/src/prognostic_equations/forcing/large_scale_advection.jl b/src/prognostic_equations/forcing/large_scale_advection.jl index 7f86040ece..7778c10a85 100644 --- a/src/prognostic_equations/forcing/large_scale_advection.jl +++ b/src/prognostic_equations/forcing/large_scale_advection.jl @@ -8,7 +8,7 @@ import ClimaCore.Spaces as Spaces import ClimaCore.Fields as Fields """ - large_scale_advection_tendency_ρq_tot(ᶜρ, thermo_params, ᶜts, t, ls_adv) + large_scale_advection_tendency_ρq_tot(ᶜρ, thermo_params, ᶜp, t, ls_adv) Computes the tendency of the total water content (`ρq_tot`) due to prescribed large-scale advection of total specific humidity (`q_tot`). @@ -18,14 +18,14 @@ is not active), it returns `NullBroadcasted()`, resulting in no tendency. Otherwise, it retrieves a profile function `prof_dqtdt` from `ls_adv`. This function provides the prescribed advective tendency of `q_tot` (i.e., -``(\\partial q_{tot}/\\partial t)_{LS_adv}``) as a function of the thermodynamic -state (`ᶜts`), time (`t`), and height (`ᶜz`). The final tendency for `ρq_tot` +``(\\partial q_{tot}/\\partial t)_{LS_adv}``) as a function of pressure (`ᶜp`), +time (`t`), and height (`ᶜz`). The final tendency for `ρq_tot` is then computed as ``ᶜρ * (\\partial q_{tot}/\\partial t)_{LS_adv}``. Arguments: - `ᶜρ`: Cell-center air density field. - `thermo_params`: Thermodynamic parameters. -- `ᶜts`: Cell-center thermodynamic state field. +- `ᶜp`: Cell-center pressure field. - `t`: Current simulation time. - `ls_adv`: `LargeScaleAdvection` object containing profile functions for tendencies, or another type if large-scale advection is inactive. @@ -38,19 +38,23 @@ Returns: function large_scale_advection_tendency_ρq_tot( ᶜρ, thermo_params, - ᶜts, + ᶜT, + ᶜp, + q_tot, + q_liq, + q_ice, t, ls_adv, ) ls_adv isa LargeScaleAdvection || return NullBroadcasted() (; prof_dTdt, prof_dqtdt) = ls_adv ᶜz = Fields.coordinate_field(axes(ᶜρ)).z - ᶜdqtdt_hadv = @. lazy(prof_dqtdt(thermo_params, ᶜts, t, ᶜz)) + ᶜdqtdt_hadv = @. lazy(prof_dqtdt(thermo_params, ᶜp, t, ᶜz)) return @. lazy(ᶜρ * ᶜdqtdt_hadv) end """ - large_scale_advection_tendency_ρe_tot(ᶜρ, thermo_params, ᶜts, t, ls_adv) + large_scale_advection_tendency_ρe_tot(ᶜρ, thermo_params, ᶜT, ᶜp, q_tot, q_liq, q_ice, t, ls_adv) Computes the tendency of total energy (`ρe_tot`) due to prescribed large-scale advection of temperature (`T`) and total specific humidity (`q_tot`). @@ -71,7 +75,9 @@ vapor for this energy calculation). Arguments: - `ᶜρ`: Cell-center air density field. - `thermo_params`: Thermodynamic parameters. -- `ᶜts`: Cell-center thermodynamic state field. +- `ᶜT`: Cell-center temperature field. +- `ᶜp`: Cell-center pressure field. +- `q_tot`, `q_liq`, `q_ice`: Specific humidity fields. - `t`: Current simulation time. - `ls_adv`: `LargeScaleAdvection` object containing profile functions for tendencies, or another type if large-scale advection is inactive. @@ -84,22 +90,26 @@ Returns: function large_scale_advection_tendency_ρe_tot( ᶜρ, thermo_params, - ᶜts, + ᶜT, + ᶜp, + q_tot, + q_liq, + q_ice, t, ls_adv, ) ls_adv isa LargeScaleAdvection || return NullBroadcasted() (; prof_dTdt, prof_dqtdt) = ls_adv z = Fields.coordinate_field(axes(ᶜρ)).z - ᶜdTdt_hadv = @. lazy(prof_dTdt(thermo_params, ᶜts, t, z)) - ᶜdqtdt_hadv = @. lazy(prof_dqtdt(thermo_params, ᶜts, t, z)) + ᶜdTdt_hadv = @. lazy(prof_dTdt(thermo_params, ᶜp, t, z)) + ᶜdqtdt_hadv = @. lazy(prof_dqtdt(thermo_params, ᶜp, t, z)) # Moisture advection term does not contain potential energy because # it's just horizontal advection of specific humidity return @. lazy( ᶜρ * ( - TD.cv_m(thermo_params, ᶜts) * ᶜdTdt_hadv + - TD.internal_energy_vapor(thermo_params, ᶜts) * ᶜdqtdt_hadv + TD.cv_m(thermo_params, q_tot, q_liq, q_ice) * ᶜdTdt_hadv + + TD.internal_energy_vapor(thermo_params, ᶜT) * ᶜdqtdt_hadv ), ) end diff --git a/src/prognostic_equations/forcing/subsidence.jl b/src/prognostic_equations/forcing/subsidence.jl index 9bab691600..e4ed517471 100644 --- a/src/prognostic_equations/forcing/subsidence.jl +++ b/src/prognostic_equations/forcing/subsidence.jl @@ -90,7 +90,7 @@ function subsidence_tendency!(Yₜ, Y, p, t, ::Subsidence) (; moisture_model) = p.atmos subsidence_profile = p.atmos.subsidence.prof thermo_params = CAP.thermodynamics_params(p.params) - (; ᶜts) = p.precomputed + (; ᶜh_tot) = p.precomputed ᶠz = Fields.coordinate_field(axes(Y.f)).z ᶠlg = Fields.local_geometry_field(Y.f) @@ -99,28 +99,29 @@ function subsidence_tendency!(Yₜ, Y, p, t, ::Subsidence) subsidence_profile(ᶠz) * CT3(unit_basis_vector_data(CT3, ᶠlg)) # LS Subsidence - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) - ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) subsidence!(Yₜ.c.ρe_tot, Y.c.ρ, ᶠsubsidence³, ᶜh_tot, Val{:first_order}()) - subsidence!(Yₜ.c.ρq_tot, Y.c.ρ, ᶠsubsidence³, ᶜq_tot, Val{:first_order}()) - if moisture_model isa NonEquilMoistModel - ᶜq_liq = @. lazy(specific(Y.c.ρq_liq, Y.c.ρ)) - subsidence!( - Yₜ.c.ρq_liq, - Y.c.ρ, - ᶠsubsidence³, - ᶜq_liq, - Val{:first_order}(), - ) - ᶜq_ice = @. lazy(specific(Y.c.ρq_ice, Y.c.ρ)) - subsidence!( - Yₜ.c.ρq_ice, - Y.c.ρ, - ᶠsubsidence³, - ᶜq_ice, - Val{:first_order}(), - ) + + if !(moisture_model isa DryModel) + ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) + subsidence!(Yₜ.c.ρq_tot, Y.c.ρ, ᶠsubsidence³, ᶜq_tot, Val{:first_order}()) + if moisture_model isa NonEquilMoistModel + ᶜq_liq = @. lazy(specific(Y.c.ρq_liq, Y.c.ρ)) + subsidence!( + Yₜ.c.ρq_liq, + Y.c.ρ, + ᶠsubsidence³, + ᶜq_liq, + Val{:first_order}(), + ) + ᶜq_ice = @. lazy(specific(Y.c.ρq_ice, Y.c.ρ)) + subsidence!( + Yₜ.c.ρq_ice, + Y.c.ρ, + ᶠsubsidence³, + ᶜq_ice, + Val{:first_order}(), + ) + end end return nothing diff --git a/src/prognostic_equations/gm_sgs_closures.jl b/src/prognostic_equations/gm_sgs_closures.jl index 6bee98a9c3..3527160cb1 100644 --- a/src/prognostic_equations/gm_sgs_closures.jl +++ b/src/prognostic_equations/gm_sgs_closures.jl @@ -54,7 +54,7 @@ This function performs several steps: Arguments: - `ᶜmixing_length`: Output `ClimaCore.Field` where the computed mixing length will be stored. - `Y`: The current state vector (containing `Y.c.uₕ`). -- `p`: Cache containing parameters (`p.params`), precomputed fields (e.g., `ᶜts`, +- `p`: Cache containing parameters (`p.params`), precomputed fields (e.g., `ᶜT`, `ᶠu³`, vertical gradients of thermodynamic variables), and scratch space. Modifies `ᶜmixing_length` in place. Also modifies fields in `p.precomputed` @@ -67,13 +67,26 @@ NVTX.@annotate function compute_gm_mixing_length(Y, p) ᶜdz = Fields.Δz_field(axes(Y.c)) ᶜlg = Fields.local_geometry_field(Y.c) - (; ᶜts, ᶠu³, ᶜlinear_buoygrad, ᶜstrain_rate_norm, cloud_diagnostics_tuple) = + (; + ᶜT, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ᶠu³, + ᶜlinear_buoygrad, + ᶜstrain_rate_norm, + cloud_diagnostics_tuple, + ) = p.precomputed @. ᶜlinear_buoygrad = buoyancy_gradients( BuoyGradMean(), thermo_params, - ᶜts, + ᶜT, + Y.c.ρ, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, cloud_diagnostics_tuple.cf, C3, p.precomputed.ᶜgradᵥ_q_tot, diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl index 1fcf1d6646..d74b3505e1 100644 --- a/src/prognostic_equations/hyperdiffusion.jl +++ b/src/prognostic_equations/hyperdiffusion.jl @@ -91,7 +91,6 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t) (; hyperdiff, turbconv_model) = p.atmos (; params) = p (; ᶜΦ) = p.core - (; ᶜts) = p.precomputed thermo_params = CAP.thermodynamics_params(params) isnothing(hyperdiff) && return nothing @@ -108,7 +107,7 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t) # Grid scale hyperdiffusion @. ᶜ∇²u = C123(wgradₕ(divₕ(ᶜu))) - C123(wcurlₕ(C123(curlₕ(ᶜu)))) - ᶜh_ref = @. lazy(h_dr(thermo_params, ᶜts, ᶜΦ)) + ᶜh_ref = @. lazy(h_dr(thermo_params, ᶜp, ᶜΦ)) @. ᶜ∇²specific_energy = wdivₕ(gradₕ(specific(Y.c.ρe_tot, Y.c.ρ) + ᶜp / Y.c.ρ - ᶜh_ref)) diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 07d7dc55a9..940f859148 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -112,11 +112,9 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) ᶜJ = Fields.local_geometry_field(axes(Y.c)).J ᶠJ = Fields.local_geometry_field(axes(Y.f)).J (; ᶠgradᵥ_ᶜΦ) = p.core - (; ᶠu³, ᶜp, ᶜts) = p.precomputed + (; ᶠu³, ᶜp, ᶜh_tot, ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed thermo_params = CAP.thermodynamics_params(params) cp_d = CAP.cp_d(params) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠu³) @@ -197,10 +195,11 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) vertical_advection_of_water_tendency!(Yₜ, Y, p, t) # This is equivalent to grad_v(Φ) + grad_v(p) / ρ - ᶜΦ_r = @. lazy(phi_r(thermo_params, ᶜts)) - ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts)) - ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜts)) - ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts)) + ᶜΦ_r = @. lazy(phi_r(thermo_params, ᶜp)) + ᶜθ_v = p.scratch.ᶜtemp_scalar + @. ᶜθ_v = theta_v(thermo_params, ᶜT, ᶜp, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜp)) + ᶜΠ = @. lazy(TD.exner_given_pressure(thermo_params, ᶜp)) @. Yₜ.f.u₃ -= ᶠgradᵥ_ᶜΦ - ᶠgradᵥ(ᶜΦ_r) + cp_d * (ᶠinterp(ᶜθ_v - ᶜθ_vr)) * ᶠgradᵥ(ᶜΠ) diff --git a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl index bc8f775d86..d171a6965b 100644 --- a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl +++ b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl @@ -428,7 +428,8 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) (; matrix) = cache (; params) = p (; ᶜΦ, ᶠgradᵥ_ᶜΦ) = p.core - (; ᶠu³, ᶜK, ᶜts, ᶜp) = p.precomputed + (; ᶜu, ᶠu³, ᶜK, ᶜp, ᶜT, ᶜh_tot) = p.precomputed + (; ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed (; ∂ᶜK_∂ᶜuₕ, ∂ᶜK_∂ᶠu₃, @@ -466,12 +467,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) Δcv_i = FT(CAP.cp_i(params) - CAP.cv_v(params)) e_int_v0 = FT(CAP.e_int_v0(params)) e_int_s0 = FT(CAP.e_int_i0(params)) + e_int_v0 - # This term appears a few times in the Jacobian, and is technically - # minus ∂e_int_∂q_tot - ∂e_int_∂q_tot = T_0 * (Δcv_v - R_d) - e_int_v0 thermo_params = CAP.thermodynamics_params(params) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) ᶜρ = Y.c.ρ ᶜuₕ = Y.c.uₕ @@ -485,9 +481,9 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) ᶜkappa_m = p.scratch.ᶜtemp_scalar @. ᶜkappa_m = - TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts) + TD.gas_constant_air(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) / + TD.cv_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) - ᶜT = @. lazy(TD.air_temperature(thermo_params, ᶜts)) ᶜ∂p∂ρq_tot = p.scratch.ᶜtemp_scalar_2 @. ᶜ∂p∂ρq_tot = ᶜkappa_m * (-e_int_v0 - R_d * T_0 - Δcv_v * (ᶜT - T_0)) + ΔR_v * ᶜT @@ -537,8 +533,9 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)] ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)] - ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts)) - ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts)) + ᶜθ_v = p.scratch.ᶜtemp_scalar_3 + @. ᶜθ_v = theta_v(thermo_params, ᶜT, ᶜp, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + ᶜΠ = @. lazy(TD.exner_given_pressure(thermo_params, ᶜp)) # In implicit tendency, we use the new pressure-gradient formulation (PGF) and gravitational acceleration: # grad(p) / ρ + grad(Φ) = cp_d * θ_v * grad(Π) + grad(Φ). # Here below, we use the old formulation of (grad(Φ) + grad(p) / ρ). @@ -660,8 +657,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) p.scratch.ᶜbidiagonal_adjoint_matrix_c3 ⋅ p.scratch.ᶠband_matrix_wvec ⋅ DiagonalMatrixRow( - e_int_func(thermo_params, p.precomputed.ᶜts) + ᶜΦ + - $(Kin(ᶜwₚ, p.precomputed.ᶜu)), + e_int_func(thermo_params, ᶜT) + ᶜΦ + $(Kin(ᶜwₚ, ᶜu)), ) end end @@ -1257,10 +1253,15 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) ## grid-mean ρe_tot ᶜkappa_m = p.scratch.ᶜtemp_scalar @. ᶜkappa_m = - TD.gas_constant_air(thermo_params, ᶜts) / - TD.cv_m(thermo_params, ᶜts) + TD.gas_constant_air( + thermo_params, + ᶜq_tot_safe, + ᶜq_liq_rai, + ᶜq_ice_sno, + ) / + TD.cv_m(thermo_params, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) + - ᶜT = @. lazy(TD.air_temperature(thermo_params, ᶜts)) ᶜ∂p∂ρq_tot = p.scratch.ᶜtemp_scalar_2 @. ᶜ∂p∂ρq_tot = ᶜkappa_m * (-e_int_v0 - R_d * T_0 - Δcv_v * (ᶜT - T_0)) + ΔR_v * ᶜT diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index 6e01d385f6..c32dc1235c 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -144,10 +144,9 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) (; ls_adv, scm_coriolis) = p.atmos (; params) = p thermo_params = CAP.thermodynamics_params(params) - (; ᶜp, sfc_conditions, ᶜts, ᶜK) = p.precomputed + (; ᶜp, ᶜK, ᶜT, ᶜh_tot, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed + (; sfc_conditions) = p.precomputed - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) vst_uₕ = viscous_sponge_tendency_uₕ(ᶜuₕ, viscous_sponge) vst_u₃ = viscous_sponge_tendency_u₃(ᶠu₃, viscous_sponge) vst_ρe_tot = viscous_sponge_tendency_ρe_tot(ᶜρ, ᶜh_tot, viscous_sponge) @@ -203,11 +202,11 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) end # For HeldSuarezForcing, the radiation_mode is used as the forcing parameter forcing = radiation_mode isa HeldSuarezForcing ? radiation_mode : nothing - hs_args = (ᶜuₕ, ᶜp, params, sfc_conditions.ts, moisture_model, forcing) + hs_args = (ᶜuₕ, ᶜp, params, sfc_conditions.T_sfc, moisture_model, forcing) hs_tendency_uₕ = held_suarez_forcing_tendency_uₕ(hs_args...) hs_tendency_ρe_tot = held_suarez_forcing_tendency_ρe_tot(ᶜρ, hs_args...) edmf_cor_tend_uₕ = scm_coriolis_tendency_uₕ(ᶜuₕ, scm_coriolis) - lsa_args = (ᶜρ, thermo_params, ᶜts, t, ls_adv) + lsa_args = (ᶜρ, thermo_params, ᶜT, ᶜp, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, t, ls_adv) bc_lsa_tend_ρe_tot = large_scale_advection_tendency_ρe_tot(lsa_args...) # TODO: fuse, once we fix diff --git a/src/prognostic_equations/surface_flux.jl b/src/prognostic_equations/surface_flux.jl index ea52352b88..4e12e37de4 100644 --- a/src/prognostic_equations/surface_flux.jl +++ b/src/prognostic_equations/surface_flux.jl @@ -106,7 +106,7 @@ function surface_flux_tendency!(Yₜ, Y, p, t) FT = eltype(Y) (; params) = p (; turbconv_model) = p.atmos - (; sfc_conditions, ᶜts) = p.precomputed + (; sfc_conditions, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed thermo_params = CAP.thermodynamics_params(params) if !disable_momentum_vertical_diffusion(p.atmos.vertical_diffusion) @@ -114,8 +114,7 @@ function surface_flux_tendency!(Yₜ, Y, p, t) @. Yₜ.c.uₕ -= btt end - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed btt = boundary_tendency_scalar(ᶜh_tot, sfc_conditions.ρ_flux_h_tot) @. Yₜ.c.ρe_tot -= btt diff --git a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl index 7f236d81aa..f1d7faa739 100644 --- a/src/prognostic_equations/vertical_diffusion_boundary_layer.jl +++ b/src/prognostic_equations/vertical_diffusion_boundary_layer.jl @@ -90,7 +90,7 @@ function vertical_diffusion_boundary_layer_tendency!( (; vertical_diffusion) = p.atmos α_vert_diff_tracer = CAP.α_vert_diff_tracer(p.params) thermo_params = CAP.thermodynamics_params(p.params) - (; ᶜu, ᶜp, ᶜts) = p.precomputed + (; ᶜu, ᶜp, ᶜT, ᶜq_liq_rai, ᶜq_ice_sno) = p.precomputed ᶠgradᵥ = Operators.GradientC2F() # apply BCs to ᶜdivᵥ, which wraps ᶠgradᵥ ᶜK_h = p.scratch.ᶜtemp_scalar if vertical_diffusion isa DecayWithHeightDiffusion @@ -114,8 +114,7 @@ function vertical_diffusion_boundary_layer_tendency!( top = Operators.SetValue(C3(0)), bottom = Operators.SetValue(C3(0)), ) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜh_tot) = p.precomputed @. Yₜ.c.ρe_tot -= ᶜdivᵥ_ρe_tot(-(ᶠinterp(Y.c.ρ) * ᶠinterp(ᶜK_h) * ᶠgradᵥ(ᶜh_tot))) ᶜρχₜ_diffusion = p.scratch.ᶜtemp_scalar_2 diff --git a/src/prognostic_equations/water_advection.jl b/src/prognostic_equations/water_advection.jl index 1f366038ac..55624f9b5f 100644 --- a/src/prognostic_equations/water_advection.jl +++ b/src/prognostic_equations/water_advection.jl @@ -36,7 +36,7 @@ function vertical_advection_of_water_tendency!(Yₜ, Y, p, t) (; params) = p (; ᶜΦ) = p.core - (; ᶜu, ᶜts) = p.precomputed + (; ᶜu, ᶜT) = p.precomputed thp = CAP.thermodynamics_params(params) ᶜJ = Fields.local_geometry_field(Y.c).J @@ -71,7 +71,7 @@ function vertical_advection_of_water_tendency!(Yₜ, Y, p, t) @. Yₜ.c.ρe_tot -= ᶜprecipdivᵥ( ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( Geometry.WVector(-(ᶜw)) * specific(ᶜρq, Y.c.ρ) * - (e_int_func(thp, ᶜts) + ᶜΦ + $(Kin(ᶜw, ᶜu))), + (e_int_func(thp, ᶜT) + ᶜΦ + $(Kin(ᶜw, ᶜu))), ), ) end @@ -108,7 +108,7 @@ function vertical_advection_of_water_tendency!(Yₜ, Y, p, t) Geometry.WVector(-(ᶜwʲ)) * draft_area(Y.c.sgsʲs.:(1).ρa, ᶜρʲs.:(1)) * ᶜqʲ * ( - e_int_func(thp, ᶜtsʲs.:(1)) - e_int_func(thp, ᶜts) - + e_int_func(thp, ᶜtsʲs.:(1)) - e_int_func(thp, ᶜT) - $(Kin(ᶜw, ᶜu)) ), ), @@ -119,7 +119,7 @@ function vertical_advection_of_water_tendency!(Yₜ, Y, p, t) -1 * ᶜprecipdivᵥ( ᶠinterp(ᶜρ⁰ * ᶜJ) / ᶠJ * ᶠright_bias( Geometry.WVector(-(ᶜwaq⁰)) * - (e_int_func(thp, ᶜts⁰) - e_int_func(thp, ᶜts) - $(Kin(ᶜw, ᶜu))), + (e_int_func(thp, ᶜts⁰) - e_int_func(thp, ᶜT) - $(Kin(ᶜw, ᶜu))), ), ) end diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 9d13e3b64e..90a4ab5bde 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -464,9 +464,11 @@ function get_large_scale_advection_model(parsed_args, ::Type{FT}) where {FT} end # See https://clima.github.io/AtmosphericProfilesLibrary.jl/dev/ # for which functions accept which arguments. - prof_dqtdt = (thermo_params, ᶜts, t, z) -> prof_dqtdt₀(z) - prof_dTdt = (thermo_params, ᶜts, t, z) -> - prof_dTdt₀(TD.exner(thermo_params, ᶜts), z) + # TODO: do not assume dry air? + prof_dqtdt = (thermo_params, ᶜp, t, z) -> prof_dqtdt₀(z) + prof_dTdt = + (thermo_params, ᶜp, t, z) -> + prof_dTdt₀(TD.exner_given_pressure(thermo_params, ᶜp), z) return LargeScaleAdvection(prof_dTdt, prof_dqtdt) end diff --git a/src/solver/types.jl b/src/solver/types.jl index a112a4fda3..50fd16c7b1 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -340,20 +340,28 @@ Base.broadcastable(x::BuoyGradMean) = tuple(x) Variables used in the environmental buoyancy gradient computation. """ -Base.@kwdef struct EnvBuoyGradVars{FT, TS} - ts::TS +Base.@kwdef struct EnvBuoyGradVars{FT} + T::FT + ρ::FT + q_tot::FT + q_liq::FT + q_ice::FT cf::FT ∂qt∂z::FT ∂θli∂z::FT end function EnvBuoyGradVars( - ts::TD.ThermodynamicState, + T, + ρ, + q_tot, + q_liq, + q_ice, cf, ∂qt∂z_∂θli∂z, ) (; ∂qt∂z, ∂θli∂z) = ∂qt∂z_∂θli∂z - return EnvBuoyGradVars(ts, cf, ∂qt∂z, ∂θli∂z) + return EnvBuoyGradVars(T, ρ, q_tot, q_liq, q_ice, cf, ∂qt∂z, ∂θli∂z) end Base.eltype(::EnvBuoyGradVars{FT}) where {FT} = FT @@ -371,7 +379,6 @@ function MixingLength(master, wall, tke, buoy, l_grid) return MixingLength(promote(master, wall, tke, buoy, l_grid)...) end - abstract type AbstractEDMF end """ diff --git a/src/surface_conditions/surface_conditions.jl b/src/surface_conditions/surface_conditions.jl index 5fd48feb09..b57e95c451 100644 --- a/src/surface_conditions/surface_conditions.jl +++ b/src/surface_conditions/surface_conditions.jl @@ -14,12 +14,16 @@ function update_surface_conditions!(Y, p, t) ) int_local_geometry_values = Fields.field_values(Fields.level(Fields.local_geometry_field(Y.c), 1)) - (; ᶜts, ᶜu, sfc_conditions) = p.precomputed + (; ᶜT, ᶜq_tot_safe, ᶜq_liq_rai, ᶜq_ice_sno, ᶜu, sfc_conditions) = p.precomputed (; params, sfc_setup, atmos) = p thermo_params = CAP.thermodynamics_params(params) surface_fluxes_params = CAP.surface_fluxes_params(params) surface_temp_params = CAP.surface_temp_params(params) - int_ts_values = Fields.field_values(Fields.level(ᶜts, 1)) + int_T_values = Fields.field_values(Fields.level(ᶜT, 1)) + int_ρ_values = Fields.field_values(Fields.level(Y.c.ρ, 1)) + int_q_tot_values = Fields.field_values(Fields.level(ᶜq_tot_safe, 1)) + int_q_liq_values = Fields.field_values(Fields.level(ᶜq_liq_rai, 1)) + int_q_ice_values = Fields.field_values(Fields.level(ᶜq_ice_sno, 1)) int_u_values = Fields.field_values(Fields.level(ᶜu, 1)) int_z_values = Fields.field_values(Fields.level(Fields.coordinate_field(Y.c).z, 1)) sfc_conditions_values = Fields.field_values(sfc_conditions) @@ -37,7 +41,11 @@ function update_surface_conditions!(Y, p, t) @. sfc_conditions_values = surface_state_to_conditions( wrapped_sfc_setup, sfc_local_geometry_values, - int_ts_values, + int_T_values, + int_ρ_values, + int_q_tot_values, + int_q_liq_values, + int_q_ice_values, projected_vector_data(CT1, int_u_values, int_local_geometry_values), projected_vector_data(CT2, int_u_values, int_local_geometry_values), int_z_values, @@ -81,16 +89,12 @@ function set_dummy_surface_conditions!(p) (; params, atmos) = p (; sfc_conditions) = p.precomputed FT = eltype(params) - thermo_params = CAP.thermodynamics_params(params) + @. sfc_conditions.T_sfc = FT(300) + @. sfc_conditions.q_vap_sfc = FT(0) @. sfc_conditions.ustar = FT(0.2) @. sfc_conditions.obukhov_length = FT(1e-4) @. sfc_conditions.buoyancy_flux = FT(0) - if atmos.moisture_model isa DryModel - @. sfc_conditions.ts = TD.PhaseDry_ρT(thermo_params, FT(1), FT(300)) - else - @. sfc_conditions.ts = TD.PhaseNonEquil_ρTq( - thermo_params, FT(1), FT(300), TD.PhasePartition(FT(0)), - ) + if !(atmos.moisture_model isa DryModel) @. sfc_conditions.ρ_flux_q_tot = C3(FT(0)) end @. sfc_conditions.ρ_flux_h_tot = C3(FT(0)) @@ -109,10 +113,14 @@ ifelsenothing(x::Nothing, default) = default surface_state_to_conditions( wrapped_sfc_setup, surface_local_geometry, - interior_ts, - interior_u, - interior_v, - interior_z, + T_int, + ρ_int, + q_tot_int, + q_liq_int, + q_ice_int, + u_int, + v_int, + z_int, thermo_params, surface_fluxes_params, surface_temp_params, @@ -127,10 +135,14 @@ first interior point. Implements the assumptions listed for `SurfaceState`. function surface_state_to_conditions( wrapped_sfc_setup::WSS, surface_local_geometry, - interior_ts, - interior_u, - interior_v, - interior_z, + T_int, + ρ_int, + q_tot_int, + q_liq_int, + q_ice_int, + u_int, + v_int, + z_int, thermo_params, surface_fluxes_params, surface_temp_params, @@ -138,11 +150,11 @@ function surface_state_to_conditions( sfc_temp_var, t, ) where {WSS} - surf_state = surface_state(wrapped_sfc_setup, surface_local_geometry, interior_z, t) + surf_state = surface_state(wrapped_sfc_setup, surface_local_geometry, z_int, t) parameterization = surf_state.parameterization (; coordinates) = surface_local_geometry Φ_sfc = geopotential(SFP.grav(surface_fluxes_params), coordinates.z) - Δz = interior_z - coordinates.z + Δz = z_int - coordinates.z FT = eltype(thermo_params) (!isnothing(surf_state.q_vap) && atmos.moisture_model isa DryModel) && @@ -160,14 +172,9 @@ function surface_state_to_conditions( u = ifelsenothing(surf_state.u, FT(0)) v = ifelsenothing(surf_state.v, FT(0)) - u_int = SA.SVector(interior_u, interior_v) - u_sfc = SA.SVector(u, v) + uv_int = SA.SVector(u_int, v_int) + uv_sfc = SA.SVector(u, v) - ρ_int = TD.air_density(thermo_params, interior_ts) - T_int = TD.air_temperature(thermo_params, interior_ts) - q_tot_int = TD.total_specific_humidity(thermo_params, interior_ts) - q_liq_int = TD.liquid_specific_humidity(thermo_params, interior_ts) - q_ice_int = TD.ice_specific_humidity(thermo_params, interior_ts) ρ_sfc = SF.surface_density( surface_fluxes_params, T_int, @@ -180,14 +187,11 @@ function surface_state_to_conditions( ) if atmos.moisture_model isa DryModel q_vap = 0 - ts = TD.PhaseDry_ρT(thermo_params, ρ_sfc, T_sfc) else # Assume that the surface is water with saturated air directly # above it. - q_vap_sat = TD.q_vap_saturation_generic(thermo_params, T_sfc, ρ_sfc, TD.Liquid()) + q_vap_sat = TD.q_vap_saturation(thermo_params, T_sfc, ρ_sfc, TD.Liquid()) q_vap = ifelsenothing(surf_state.q_vap, q_vap_sat) - q = TD.PhasePartition(q_vap) - ts = TD.PhaseNonEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q) end gustiness = ifelsenothing(surf_state.gustiness, FT(1)) @@ -220,9 +224,10 @@ function surface_state_to_conditions( "q_flux cannot be specified when using a DryModel", ) end - ρ = TD.air_density(thermo_params, ts) - shf = θ_flux * ρ * TD.cp_m(thermo_params, ts) - lhf = q_flux * ρ * TD.latent_heat_vapor(thermo_params, ts) + # Use ρ_sfc directly and compute cp_m from primitives + q_sfc = TD.PhasePartition(q_vap) + shf = θ_flux * ρ_sfc * TD.cp_m(thermo_params, q_sfc) + lhf = q_flux * ρ_sfc * TD.latent_heat_vapor(thermo_params, T_sfc) end flux_specs = SF.FluxSpecs(ustar = parameterization.ustar, shf = shf, lhf = lhf) config = SF.default_surface_flux_config(FT) @@ -232,9 +237,9 @@ function surface_state_to_conditions( return atmos_surface_conditions( surface_fluxes_params, SF.surface_fluxes(surface_fluxes_params, T_int, q_tot_int, q_liq_int, q_ice_int, - ρ_int, T_sfc, q_vap, Φ_sfc, Δz, 0, u_int, u_sfc, nothing, config, + ρ_int, T_sfc, q_vap, Φ_sfc, Δz, 0, uv_int, uv_sfc, nothing, config, UF.PointValueScheme(), nothing, flux_specs), - ts, + ρ_sfc, surface_local_geometry, ) end @@ -286,16 +291,15 @@ function surface_temperature( end """ - atmos_surface_conditions(surface_conditions, ts, surface_local_geometry) + atmos_surface_conditions(surface_conditions, ρ_sfc, surface_local_geometry) -Adds local geometry information to the `SurfaceFluxes.SurfaceFluxConditions` struct -along with information about the thermodynamic state. The resulting values are the -ones actually used by ClimaAtmos operator boundary conditions. +Adds local geometry information to the `SurfaceFluxes.SurfaceFluxConditions` struct. +The resulting values are the ones actually used by ClimaAtmos operator boundary conditions. """ function atmos_surface_conditions( surface_fluxes_params, surface_conditions, - ts, + ρ_sfc, surface_local_geometry, ) (; ustar, L_MO, ρτxz, ρτyz, shf, lhf, evaporation, T_sfc, q_vap_sfc) = @@ -303,8 +307,6 @@ function atmos_surface_conditions( # surface normal z = surface_normal(surface_local_geometry) - thermo_params = SFP.thermodynamics_params(surface_fluxes_params) - ρ_sfc = TD.air_density(thermo_params, ts) buoy_flux = SF.buoyancy_flux(surface_fluxes_params, shf, lhf, T_sfc, ρ_sfc, q_vap_sfc) @@ -314,7 +316,8 @@ function atmos_surface_conditions( moisture_flux = (; ρ_flux_q_tot = vector_from_component(evaporation, z)) return (; - ts, + T_sfc, + q_vap_sfc, ustar, obukhov_length = L_MO, buoyancy_flux = buoy_flux, @@ -348,12 +351,11 @@ function surface_conditions_type(atmos, ::Type{FT}) where {FT} # NOTE: Technically ρ_flux_q_tot is not really needed for a dry model, but # SF always has evaporation moisture_flux_names = (:ρ_flux_q_tot,) - names = (:ts, :ustar, :obukhov_length, :buoyancy_flux, :ρ_flux_uₕ, + names = (:T_sfc, :q_vap_sfc, :ustar, :obukhov_length, :buoyancy_flux, :ρ_flux_uₕ, energy_flux_names..., moisture_flux_names..., ) type_tuple = Tuple{ - atmos.moisture_model isa DryModel ? TD.PhaseDry{FT} : TD.PhaseNonEquil{FT}, - FT, FT, FT, + FT, FT, FT, FT, FT, typeof(C3(FT(0)) ⊗ C12(FT(0), FT(0))), ntuple(_ -> C3{FT}, Val(length(energy_flux_names)))..., ntuple(_ -> C3{FT}, Val(length(moisture_flux_names)))..., diff --git a/src/utils/refstate_thermodynamics.jl b/src/utils/refstate_thermodynamics.jl index 8b5b543d8d..30a048449f 100644 --- a/src/utils/refstate_thermodynamics.jl +++ b/src/utils/refstate_thermodynamics.jl @@ -5,52 +5,42 @@ # https://github.com/CliMA/ClimaParams.jl/pull/253 const s_ref = 7 -function dry_exner_function(thermo_params, ts) - p = TD.air_pressure(thermo_params, ts) - return TD.exner_given_pressure( - thermo_params, - p, - ) -end - - -function theta_v(thermo_params, ts) +function theta_v(thermo_params, T, p, q_tot, q_liq, q_ice) R_d = TD.TP.R_d(thermo_params) - T = TD.air_temperature(thermo_params, ts) - R_m = TD.gas_constant_air(thermo_params, ts) - Π = dry_exner_function(thermo_params, ts) + R_m = TD.gas_constant_air(thermo_params, TD.PhasePartition(q_tot, q_liq, q_ice)) + Π = TD.exner_given_pressure(thermo_params, p) return T * R_m / (Π * R_d) end -function air_temperature_reference(thermo_params, ts) +function air_temperature_reference(thermo_params, p) T_min = TD.TP.T_min_ref(thermo_params) T_sfc = TD.TP.T_surf_ref(thermo_params) - Π = dry_exner_function(thermo_params, ts) + Π = TD.exner_given_pressure(thermo_params, p) return T_min + (T_sfc - T_min) * Π^s_ref end -function theta_vr(thermo_params, ts) - T_r = air_temperature_reference(thermo_params, ts) - Π = dry_exner_function(thermo_params, ts) +function theta_vr(thermo_params, p) + T_r = air_temperature_reference(thermo_params, p) + Π = TD.exner_given_pressure(thermo_params, p) return T_r / Π end -function phi_r(thermo_params, ts) +function phi_r(thermo_params, p) cp_d = TD.TP.cp_d(thermo_params) T_min = TD.TP.T_min_ref(thermo_params) T_sfc = TD.TP.T_surf_ref(thermo_params) - Π = dry_exner_function(thermo_params, ts) + Π = TD.exner_given_pressure(thermo_params, p) return -cp_d * (T_min * log(Π) + (T_sfc - T_min) / s_ref * (Π^s_ref - 1)) end -function h_dr(thermo_params, ts, Φ) +function h_dr(thermo_params, p, Φ) T_0 = TD.TP.T_0(thermo_params) cp_d = TD.TP.cp_d(thermo_params) - T_r = air_temperature_reference(thermo_params, ts) + T_r = air_temperature_reference(thermo_params, p) return cp_d * (T_r - T_0) + Φ end diff --git a/src/utils/variable_manipulations.jl b/src/utils/variable_manipulations.jl index a742c2b7fb..e3155cb0dd 100644 --- a/src/utils/variable_manipulations.jl +++ b/src/utils/variable_manipulations.jl @@ -435,10 +435,7 @@ Returns: """ function ᶜspecific_env_mse(Y, p) turbconv_model = p.atmos.turbconv_model - (; ᶜK, ᶜts) = p.precomputed - thermo_params = CAP.thermodynamics_params(p.params) - ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ)) - ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot)) + (; ᶜK, ᶜh_tot) = p.precomputed # grid-scale moist static energy density `ρ * mse`. ᶜρmse = @. lazy(Y.c.ρ * (ᶜh_tot - ᶜK)) diff --git a/test/coupler_compatibility.jl b/test/coupler_compatibility.jl index 879f3f604e..7700008765 100644 --- a/test/coupler_compatibility.jl +++ b/test/coupler_compatibility.jl @@ -86,13 +86,11 @@ const T2 = 290 # temperature to T1 and then to T2. @. sfc_setup.T = FT(T1) CA.set_precomputed_quantities!(Y, p_overwritten, t) - sfc_T = - @. TD.air_temperature(thermo_params, p.precomputed.sfc_conditions.ts) + sfc_T = p.precomputed.sfc_conditions.T_sfc @test all(isequal(T1), parent(sfc_T)) @. sfc_setup.T = FT(T2) CA.set_precomputed_quantities!(Y, p_overwritten, t) - sfc_T = - @. TD.air_temperature(thermo_params, p.precomputed.sfc_conditions.ts) + sfc_T = p.precomputed.sfc_conditions.T_sfc @test all(isequal(T2), parent(sfc_T)) end