Merge #2053

2053: trying to get rid of MSLP overwrite r=trontrytel a=trontrytel Co-authored-by: Anna Jaruga <ajaruga@caltech.edu>
CliMA · Sep 6, 2023 · 01e1f4b · 01e1f4b
2 parents 7b777bc + 634842f
commit 01e1f4b
Show file tree

Hide file tree

Showing 6 changed files with 18 additions and 14 deletions.
diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl
@@ -170,7 +170,7 @@ end
 
 if config.parsed_args["check_conservation"]
     FT = Spaces.undertype(axes(sol.u[end].c.ρ))
-    @test sum(sol.u[1].c.ρ) ≈ sum(sol.u[end].c.ρ) rtol = 25 * eps(FT)
+    @test sum(sol.u[1].c.ρ) ≈ sum(sol.u[end].c.ρ) rtol = 50 * eps(FT)
     @test sum(sol.u[1].c.ρe_tot) +
           (p.net_energy_flux_sfc[][] - p.net_energy_flux_toa[][]) ≈
           sum(sol.u[end].c.ρe_tot) rtol = 30 * eps(FT)

diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl
@@ -1 +1 @@
-122
+123
diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl
@@ -105,7 +105,7 @@ function (initial_condition::AgnesiHProfile)(params)
         (; x, z) = local_geometry.coordinates
         cp_d = CAP.cp_d(params)
         cv_d = CAP.cv_d(params)
-        p_0 = CAP.MSLP(params)
+        p_0 = CAP.p_ref_theta(params)
         R_d = CAP.R_d(params)
         T_0 = CAP.T_0(params)
         # auxiliary quantities
@@ -141,7 +141,7 @@ function (initial_condition::ScharProfile)(params)
         R_d = CAP.R_d(params)
         cp_d = CAP.cp_d(params)
         cv_d = CAP.cv_d(params)
-        p₀ = CAP.MSLP(params)
+        p₀ = CAP.p_ref_theta(params)
         (; x, z) = local_geometry.coordinates
         θ₀ = FT(280.0)
         buoy_freq = FT(0.01)
@@ -192,7 +192,7 @@ function (initial_condition::DryDensityCurrentProfile)(params)
         θ_c = FT(-15)
         cp_d = CAP.cp_d(params)
         cv_d = CAP.cv_d(params)
-        p_0 = CAP.MSLP(params)
+        p_0 = CAP.p_ref_theta(params)
         R_d = CAP.R_d(params)
         T_0 = CAP.T_0(params)
 
@@ -250,7 +250,7 @@ function (initial_condition::RisingThermalBubbleProfile)(params)
         θ_c = FT(0.5)
         cp_d = CAP.cp_d(params)
         cv_d = CAP.cv_d(params)
-        p_0 = CAP.MSLP(params)
+        p_0 = CAP.p_ref_theta(params)
         R_d = CAP.R_d(params)
         T_0 = CAP.T_0(params)
 

diff --git a/src/parameterized_tendencies/radiation/held_suarez.jl b/src/parameterized_tendencies/radiation/held_suarez.jl
@@ -37,6 +37,7 @@ function forcing_tendency!(Yₜ, Y, p, t, colidx, ::HeldSuarezForcing)
     cv_d = FT(CAP.cv_d(params))
     day = FT(CAP.day(params))
     MSLP = FT(CAP.MSLP(params))
+    p_ref_theta = FT(CAP.p_ref_theta(params))
     grav = FT(CAP.grav(params))
     ΔT_y_dry = FT(CAP.ΔT_y_dry(params))
     ΔT_y_wet = FT(CAP.ΔT_y_wet(params))
@@ -79,14 +80,17 @@ function forcing_tendency!(Yₜ, Y, p, t, colidx, ::HeldSuarezForcing)
                 T_min,
                 (
                     T_equator - ΔT_y * abs2(sin(ᶜφ[colidx])) -
-                    Δθ_z * log(ᶜp[colidx] / MSLP) * abs2(cos(ᶜφ[colidx]))
-                ) * fast_pow(ᶜp[colidx] / MSLP, κ_d),
+                    Δθ_z *
+                    log(ᶜp[colidx] / p_ref_theta) *
+                    abs2(cos(ᶜφ[colidx]))
+                ) * fast_pow(ᶜp[colidx] / p_ref_theta, κ_d),
             )
         )
 
     @. Yₜ.c.uₕ[colidx] -= (k_f * ᶜheight_factor[colidx]) * Y.c.uₕ[colidx]
     if :ρθ in propertynames(Y.c)
-        @. Yₜ.c.ρθ[colidx] -= ᶜΔρT[colidx] * fast_pow((MSLP / ᶜp[colidx]), κ_d)
+        @. Yₜ.c.ρθ[colidx] -=
+            ᶜΔρT[colidx] * fast_pow((p_ref_theta / ᶜp[colidx]), κ_d)
     elseif :ρe_tot in propertynames(Y.c)
         @. Yₜ.c.ρe_tot[colidx] -= ᶜΔρT[colidx] * cv_d
     end

diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
@@ -161,15 +161,14 @@ function create_parameter_set(config::AtmosConfig)
     dt = FT(CA.time_to_seconds(parsed_args["dt"]))
 
     return if CA.is_column_edmf(parsed_args)
-        overrides = (; MSLP = 100000.0, τ_precip = dt)
+        overrides = (; τ_precip = dt)
         create_climaatmos_parameter_set(toml_dict, parsed_args, overrides)
     elseif CA.is_column_without_edmf(parsed_args)
         overrides = (; τ_precip = dt)
         create_climaatmos_parameter_set(toml_dict, parsed_args, overrides)
     else
         overrides = (;
             R_d = 287.0,
-            MSLP = 1.0e5,
             grav = 9.80616,
             Omega = 7.29212e-5,
             planet_radius = 6.371229e6,

diff --git a/src/prognostic_equations/implicit/wfact.jl b/src/prognostic_equations/implicit/wfact.jl
@@ -89,7 +89,7 @@ function Wfact!(W, Y, p, dtγ, t, colidx)
     κ_d = FT(CAP.kappa_d(params))
     cv_d = FT(CAP.cv_d(params))
     T_tri = FT(CAP.T_triple(params))
-    MSLP = FT(CAP.MSLP(params))
+    p_ref_theta = FT(CAP.p_ref_theta(params))
 
     dtγ_ref[] = dtγ
 
@@ -207,12 +207,13 @@ function Wfact!(W, Y, p, dtγ, t, colidx)
         # If we ignore the dependence of pressure on moisture,
         # ∂(ᶠgradᵥ(ᶜp - ᶜp_ref))/∂(ᶜρθ) =
         #     ᶠgradᵥ_stencil(
-        #         R_d / (1 - κ_d) * (ᶜρθ * R_d / MSLP)^(κ_d / (1 - κ_d))
+        #         R_d / (1 - κ_d) * (ᶜρθ * R_d / p_ref_theta)^(κ_d / (1 - κ_d))
         #     )
         ᶜρθ = Y.c.ρθ
         @. ∂ᶠ𝕄ₜ∂ᶜ𝔼[colidx] = map_get_data(
             -1 / ᶠinterp(ᶜρ[colidx]) * ᶠgradᵥ_stencil(
-                R_d / (1 - κ_d) * (ᶜρθ[colidx] * R_d / MSLP)^(κ_d / (1 - κ_d)),
+                R_d / (1 - κ_d) *
+                (ᶜρθ[colidx] * R_d / p_ref_theta)^(κ_d / (1 - κ_d)),
             ),
         )