Merge pull request #3844 from CliMA/as/hypdiff-ref

Enthalpy perturbation in hyperdiffusion tendency

Author: akshaysridhar

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-246
+247
 
 # **README**
 #
@@ -20,6 +20,9 @@
 
 
 #=
+247
+- Hyperdiffusion of enthalpy perturbation
+
 246
 - Prognostic edmf docstrings + refactor + scm_coriolis change
 

src/cache/cache.jl

@@ -129,6 +129,46 @@ function build_cache(
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
 
+    # Reference enthalpy for hyperdiffusion computation
+    ᶜz = Fields.coordinate_field(Y.c).z
+    T_ref = similar(ᶜz)
+    p_ref = similar(ᶜz)
+    q_tot_ref = similar(ᶜz)
+    cv_m = similar(ᶜz)
+    R_m = similar(ᶜz)
+    thermo_params = CAP.thermodynamics_params(params)
+    decay_scale_height = FT(8000)
+
+    # Reference State
+    temp_profile = TD.TemperatureProfiles.DecayingTemperatureProfile{FT}(
+        thermo_params,
+        FT(300),  # surface temperature
+        FT(220),  # top temperature
+        decay_scale_height,  # decay height scale
+    )
+
+    rel_hum_ref =
+        atmos.moisture_model isa DryModel ? FT(0) :
+        @. FT(0.5) * exp(- ᶜz / decay_scale_height)
+
+    T_0 = CAP.T_0(params)
+    grav = CAP.grav(params)
+    @. T_ref = first(temp_profile(thermo_params, ᶜz))
+    @. p_ref = last(temp_profile(thermo_params, ᶜz))
+    @. q_tot_ref =
+        TD.q_vap_from_RH_liquid(thermo_params, p_ref, T_ref, rel_hum_ref)
+    @. cv_m = TD.cv_m(thermo_params, TD.PhasePartition(q_tot_ref, FT(0), FT(0)))
+    @. R_m = TD.gas_constant_air(
+        thermo_params,
+        TD.PhasePartition(q_tot_ref, FT(0), FT(0)),
+    )
+
+    ᶜh_ref = @. TD.total_specific_enthalpy(
+        cv_m * (T_ref - T_0) + grav * ᶜz,
+        R_m,
+        T_ref,
+    )
+
     core = (
         ᶜΦ,
         ᶠgradᵥ_ᶜΦ = ᶠgradᵥ.(ᶜΦ),
@@ -139,6 +179,7 @@ function build_cache(
         surface_ct3_unit = CT3.(
             unit_basis_vector_data.(CT3, sfc_local_geometry)
         ),
+        ᶜh_ref,
     )
     external_forcing = external_forcing_cache(Y, atmos, params, start_date)
     sfc_setup = surface_setup(params)

src/prognostic_equations/hyperdiffusion.jl

@@ -86,6 +86,7 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     n = n_mass_flux_subdomains(turbconv_model)
     diffuse_tke = use_prognostic_tke(turbconv_model)
     (; ᶜp, ᶜspecific) = p.precomputed
+    (; ᶜh_ref) = p.core
     (; ᶜ∇²u, ᶜ∇²specific_energy) = p.hyperdiff
     if turbconv_model isa PrognosticEDMFX
         (; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²mseʲs) = p.hyperdiff
@@ -97,7 +98,7 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
         C123(wcurlₕ(C123(curlₕ(p.precomputed.ᶜu))))
 
     @. ᶜ∇²specific_energy =
-        wdivₕ(gradₕ(specific(Y.c.ρe_tot, Y.c.ρ) + ᶜp / Y.c.ρ))
+        wdivₕ(gradₕ(specific(Y.c.ρe_tot, Y.c.ρ) + ᶜp / Y.c.ρ - ᶜh_ref))
 
     if diffuse_tke
         (; ᶜtke⁰) = p.precomputed