Merge pull request #2421 from CliMA/zs/tke_pressure

add pressure work term to TKE equation

.buildkite/pipeline.yml

@@ -580,7 +580,7 @@ steps:
         command: >
           julia --color=yes --project=examples examples/hybrid/driver.jl
           --config_file $CONFIG_PATH/diagnostic_edmfx_dycoms_rf01_explicit_box.yml
-        artifact_paths: "diagnostic_edmfx_dycoms_rf01_explicitbox/*"
+        artifact_paths: "diagnostic_edmfx_dycoms_rf01_explicit_box/*"
         agents:
           slurm_mem: 20GB
 

config/model_configs/diagnostic_edmfx_dycoms_rf01_explicit_box.yml

@@ -23,7 +23,7 @@ y_elem: 2
 z_elem: 30 
 z_max: 1500 
 z_stretch: false
-dt: 40secs 
+dt: 20secs 
 t_end: 4hours 
 dt_save_to_disk: 10mins
 toml: [toml/diagnostic_edmfx_box.toml]

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -245,7 +245,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
         ᶜρʲs,
         ᶜentrʲs,
         ᶜdetrʲs,
-        ᶜnh_pressureʲs,
+        ᶠnh_pressure³ʲs,
         ᶜS_q_totʲs,
         ᶜS_e_totʲs_helper,
     ) = p.precomputed
@@ -312,7 +312,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
             ᶜq_totʲ = ᶜq_totʲs.:($j)
             ᶜentrʲ = ᶜentrʲs.:($j)
             ᶜdetrʲ = ᶜdetrʲs.:($j)
-            ᶜnh_pressureʲ = ᶜnh_pressureʲs.:($j)
+            ᶠnh_pressure³ʲ = ᶠnh_pressure³ʲs.:($j)
             ᶜS_q_totʲ = ᶜS_q_totʲs.:($j)
             ᶜS_e_totʲ_helper = ᶜS_e_totʲs_helper.:($j)
 
@@ -336,8 +336,8 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
             tsʲ_prev_level = Fields.field_values(Fields.level(ᶜtsʲ, i - 1))
             entrʲ_prev_level = Fields.field_values(Fields.level(ᶜentrʲ, i - 1))
             detrʲ_prev_level = Fields.field_values(Fields.level(ᶜdetrʲ, i - 1))
-            nh_pressureʲ_prev_level =
-                Fields.field_values(Fields.level(ᶜnh_pressureʲ, i - 1))
+            nh_pressure³ʲ_prev_halflevel =
+                Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i - 1 - half))
             scale_height =
                 CAP.R_d(params) * CAP.T_surf_ref(params) / CAP.grav(params)
             S_q_totʲ_prev_level =
@@ -369,7 +369,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
             )
 
             # TODO: use updraft top instead of scale height
-            @. nh_pressureʲ_prev_level = ᶠupdraft_nh_pressure(
+            @. nh_pressure³ʲ_prev_halflevel = ᶠupdraft_nh_pressure(
                 params,
                 p.atmos.edmfx_nh_pressure,
                 local_geometry_prev_halflevel,
@@ -380,8 +380,8 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                 scale_height,
             )
 
-            nh_pressureʲ_data_prev_level =
-                nh_pressureʲ_prev_level.components.data.:1
+            nh_pressure³ʲ_data_prev_halflevel =
+                nh_pressure³ʲ_prev_halflevel.components.data.:1
 
             # Updraft q_tot sources from precipitation formation
             # To be applied in updraft continuity, moisture and energy
@@ -446,7 +446,7 @@ function set_diagnostic_edmf_precomputed_quantities_do_integral!(Y, p, t)
                     local_geometry_prev_level.J *
                     local_geometry_prev_level.J *
                     2 *
-                    nh_pressureʲ_data_prev_level
+                    nh_pressure³ʲ_data_prev_halflevel
                 )
 
             # get u³ʲ to calculate divergence term for detrainment, 
@@ -644,24 +644,30 @@ Updates the top boundary condition of precomputed quantities stored in `p` for d
 function set_diagnostic_edmf_precomputed_quantities_top_bc!(Y, p, t)
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
     (; ᶜentrʲs, ᶜdetrʲs, ᶜS_q_totʲs, ᶜS_e_totʲs_helper) = p.precomputed
-    (; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs) = p.precomputed
+    (; ᶠu³⁰, ᶠu³ʲs, ᶜuʲs, ᶠnh_pressure³ʲs) = p.precomputed
 
     # set values for the top level
     i_top = Spaces.nlevels(axes(Y.c))
     u³⁰_halflevel = Fields.field_values(Fields.level(ᶠu³⁰, i_top + half))
     @. u³⁰_halflevel = CT3(0)
 
     for j in 1:n
-        ᶠu³ʲ = ᶠu³ʲs.:($j)
         ᶜuʲ = ᶜuʲs.:($j)
         ᶠu³ʲ = ᶠu³ʲs.:($j)
+        ᶠnh_pressure³ʲ = ᶠnh_pressure³ʲs.:($j)
         ᶜentrʲ = ᶜentrʲs.:($j)
         ᶜdetrʲ = ᶜdetrʲs.:($j)
         ᶜS_q_totʲ = ᶜS_q_totʲs.:($j)
         ᶜS_e_totʲ_helper = ᶜS_e_totʲs_helper.:($j)
 
         u³ʲ_halflevel = Fields.field_values(Fields.level(ᶠu³ʲ, i_top + half))
         @. u³ʲ_halflevel = CT3(0)
+        nh_pressure³ʲ_halflevel =
+            Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i_top - half))
+        @. nh_pressure³ʲ_halflevel = CT3(0)
+        nh_pressure³ʲ_halflevel =
+            Fields.field_values(Fields.level(ᶠnh_pressure³ʲ, i_top + half))
+        @. nh_pressure³ʲ_halflevel = CT3(0)
 
         entrʲ_level = Fields.field_values(Fields.level(ᶜentrʲ, i_top))
         detrʲ_level = Fields.field_values(Fields.level(ᶜdetrʲ, i_top))

src/cache/precomputed_quantities.jl

@@ -78,6 +78,7 @@ function precomputed_quantities(Y, atmos)
             ᶜρʲs = similar(Y.c, NTuple{n, FT}),
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
+            ᶠnh_pressure₃ʲs = similar(Y.f, NTuple{n, C3{FT}}),
         ) : (;)
     diagnostic_sgs_quantities =
         atmos.turbconv_model isa DiagnosticEDMFX ?
@@ -92,7 +93,7 @@ function precomputed_quantities(Y, atmos)
             ᶜq_totʲs = similar(Y.c, NTuple{n, FT}),
             ᶜentrʲs = similar(Y.c, NTuple{n, FT}),
             ᶜdetrʲs = similar(Y.c, NTuple{n, FT}),
-            ᶜnh_pressureʲs = similar(Y.c, NTuple{n, CT3{FT}}),
+            ᶠnh_pressure³ʲs = similar(Y.f, NTuple{n, CT3{FT}}),
             ᶜS_q_totʲs = similar(Y.c, NTuple{n, FT}),
             ᶜS_q_tot⁰ = similar(Y.c, FT),
             ᶜS_e_totʲs_helper = similar(Y.c, NTuple{n, FT}),

src/prognostic_equations/edmfx_closures.jl

@@ -108,7 +108,7 @@ function edmfx_nh_pressure_tendency!(
     n = n_mass_flux_subdomains(turbconv_model)
     (; params) = p
     (; ᶠgradᵥ_ᶜΦ) = p.core
-    (; ᶜρʲs, ᶠu₃⁰) = p.precomputed
+    (; ᶜρʲs, ᶠnh_pressure₃ʲs, ᶠu₃⁰) = p.precomputed
     FT = eltype(Y)
     ᶜz = Fields.coordinate_field(Y.c).z
     z_sfc = Fields.level(Fields.coordinate_field(Y.f).z, Fields.half)
@@ -131,7 +131,7 @@ function edmfx_nh_pressure_tendency!(
         end
         updraft_top = updraft_top - z_sfc[colidx][]
 
-        @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -= ᶠupdraft_nh_pressure(
+        @. ᶠnh_pressure₃ʲs.:($$j)[colidx] = ᶠupdraft_nh_pressure(
             params,
             p.atmos.edmfx_nh_pressure,
             ᶠlg[colidx],
@@ -144,6 +144,8 @@ function edmfx_nh_pressure_tendency!(
             ᶠu₃⁰[colidx],
             updraft_top,
         )
+
+        @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -= ᶠnh_pressure₃ʲs.:($$j)[colidx]
     end
 end
 

src/prognostic_equations/edmfx_tke.jl

@@ -24,6 +24,21 @@ function edmfx_tke_tendency!(
     (; ᶜK_u, ᶜK_h) = p.precomputed
     (; dt) = p
     ᶜρa⁰ = turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶜρa⁰ : Y.c.ρ
+    nh_pressure3ʲs =
+        turbconv_model isa PrognosticEDMFX ? p.precomputed.ᶠnh_pressure₃ʲs :
+        p.precomputed.ᶠnh_pressure³ʲs
+    ᶜtke_press = p.scratch.ᶜtemp_scalar[colidx]
+    @. ᶜtke_press = 0
+    for j in 1:n
+        ᶜρaʲ_colidx =
+            turbconv_model isa PrognosticEDMFX ? Y.c.sgsʲs.:($j).ρa[colidx] :
+            p.precomputed.ᶜρaʲs.:($j)[colidx]
+        @. ᶜtke_press +=
+            ᶜρaʲ_colidx *
+            adjoint(ᶜinterp.(ᶠu³ʲs.:($$j)[colidx] - ᶠu³⁰[colidx])) *
+            ᶜinterp(C3(nh_pressure3ʲs.:($$j)[colidx]))
+    end
+
 
     if use_prognostic_tke(turbconv_model)
         # shear production
@@ -43,6 +58,7 @@ function edmfx_tke_tendency!(
                 )
         end
         # pressure work
+        @. Yₜ.c.sgs⁰.ρatke[colidx] += ᶜtke_press[colidx]
         # dissipation
         @. Yₜ.c.sgs⁰.ρatke[colidx] -=
             ᶜρa⁰[colidx] * c_d * max(ᶜtke⁰[colidx], 0)^(FT(3) / 2) / max(