remove grid-scale thermo state

.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]

.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

NEWS.md

@@ -3,13 +3,17 @@ 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
 
 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

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"

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},

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.

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)