Merge #513

513: move env_thermo precip arrays to aux r=charleskawczynski a=trontrytel

Peel off from #506 - moves environment precipitation related data structures to aux. 

(does not affect RiCO)

Co-authored-by: Anna Jaruga <[email protected]>

integration_tests/utils/main.jl

@@ -112,6 +112,8 @@ cent_aux_vars_edmf(FT, n_up) = (;
             Hvar = FT(0),
             QTvar = FT(0),
             HQTcov = FT(0),
+            qt_tendency_rain_formation = FT(0),
+            θ_liq_ice_tendency_rain_formation = FT(0),
         ),
         θ_liq_ice_tendency_rain_evap = FT(0),
         qt_tendency_rain_evap = FT(0),

src/EDMF_Environment.jl

@@ -1,11 +1,13 @@
-function update_env_precip_tendencies(en_thermo::EnvironmentThermodynamics, state, k, qt_tendency, θ_liq_ice_tendency)
+function update_env_precip_tendencies(state, k, qt_tendency, θ_liq_ice_tendency, qr_tendency, qs_tendency)
 
     aux_en = center_aux_environment(state)
     tendencies_pr = center_tendencies_precipitation(state)
-    en_thermo.qt_tendency_rain_formation[k] = qt_tendency * aux_en.area[k]
-    tendencies_pr.qr[k] += -en_thermo.qt_tendency_rain_formation[k]
-    en_thermo.θ_liq_ice_tendency_rain_formation[k] = θ_liq_ice_tendency * aux_en.area[k]
 
+    aux_en.qt_tendency_rain_formation[k] = qt_tendency * aux_en.area[k]
+    aux_en.θ_liq_ice_tendency_rain_formation[k] = θ_liq_ice_tendency * aux_en.area[k]
+
+    tendencies_pr.qr[k] += qr_tendency
+    tendencies_pr.qs[k] += qs_tendency
     return
 end
 
@@ -52,7 +54,14 @@ function sgs_mean(en_thermo::EnvironmentThermodynamics, grid, state, en, precip,
             ts,
         )
         update_sat_unsat(en_thermo, state, k, ts)
-        update_env_precip_tendencies(en_thermo, state, k, mph.qt_tendency, mph.θ_liq_ice_tendency)
+        update_env_precip_tendencies(
+            state,
+            k,
+            mph.qt_tendency,
+            mph.θ_liq_ice_tendency,
+            mph.qr_tendency,
+            mph.qs_tendency,
+        )
     end
     return
 end
@@ -76,7 +85,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
     # arrays for storing quadarature points and ints for labeling items in the arrays
     # a python dict would be nicer, but its 30% slower than this (for python 2.7. It might not be the case for python 3)
     env_len = 8
-    src_len = 6
+    src_len = 8
 
     sqpi_inv = 1.0 / sqrt(π)
     sqrt2 = sqrt(2.0)
@@ -89,7 +98,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
     inner_src = zeros(src_len)
     outer_src = zeros(src_len)
     i_ql, i_qi, i_T, i_cf, i_qt_sat, i_qt_unsat, i_T_sat, i_T_unsat = 1:env_len
-    i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH = 1:src_len
+    i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH, i_Sqr, i_Sqs = 1:src_len
 
     @inbounds for k in real_center_indices(grid)
         if (
@@ -194,6 +203,8 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
                     end
                     # products for variance and covariance source terms
                     inner_src[i_Sqt] += mph.qt_tendency * weights[m_h] * sqpi_inv
+                    inner_src[i_Sqr] += mph.qr_tendency * weights[m_h] * sqpi_inv
+                    inner_src[i_Sqs] += mph.qs_tendency * weights[m_h] * sqpi_inv
                     inner_src[i_SH] += mph.θ_liq_ice_tendency * weights[m_h] * sqpi_inv
                     inner_src[i_Sqt_H] += mph.qt_tendency * h_hat * weights[m_h] * sqpi_inv
                     inner_src[i_Sqt_qt] += mph.qt_tendency * qt_hat * weights[m_h] * sqpi_inv
@@ -210,7 +221,14 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
             end
 
             # update environmental variables
-            update_env_precip_tendencies(en_thermo, state, k, outer_src[i_Sqt], outer_src[i_SH])
+            update_env_precip_tendencies(
+                state,
+                k,
+                outer_src[i_Sqt],
+                outer_src[i_SH],
+                outer_src[i_Sqr],
+                outer_src[i_Sqs],
+            )
 
             # update cloudy/dry variables for buoyancy in TKE
             aux_en.cloud_fraction[k] = outer_env[i_cf]
@@ -250,7 +268,14 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
                 dt,
                 ts,
             )
-            update_env_precip_tendencies(en_thermo, state, k, mph.qt_tendency, mph.θ_liq_ice_tendency)
+            update_env_precip_tendencies(
+                state,
+                k,
+                mph.qt_tendency,
+                mph.θ_liq_ice_tendency,
+                mph.qr_tendency,
+                mph.qs_tendency,
+            )
             update_sat_unsat(en_thermo, state, k, ts)
 
             en_thermo.Hvar_rain_dt[k] = 0.0

src/Turbulence_PrognosticTKE.jl

@@ -95,11 +95,11 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
         end
         tendencies_gm.q_tot[k] +=
             edmf.UpdThermo.qt_tendency_rain_formation_tot[k] +
-            edmf.EnvThermo.qt_tendency_rain_formation[k] +
+            aux_en.qt_tendency_rain_formation[k] +
             aux_tc.qt_tendency_rain_evap[k]
         tendencies_gm.θ_liq_ice[k] +=
             edmf.UpdThermo.θ_liq_ice_tendency_rain_formation_tot[k] +
-            edmf.EnvThermo.θ_liq_ice_tendency_rain_formation[k] +
+            aux_en.θ_liq_ice_tendency_rain_formation[k] +
             aux_tc.θ_liq_ice_tendency_rain_evap[k]
     end
 

src/diagnostics.jl

@@ -263,7 +263,7 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
         # per timestep per EDMF surface area [mm/h]
         if (precip.precipitation_model == "cutoff")
             precip.cutoff_precipitation_rate -=
-                (en_thermo.qt_tendency_rain_formation[k] + up_thermo.qt_tendency_rain_formation_tot[k]) *
+                (aux_en.qt_tendency_rain_formation[k] + up_thermo.qt_tendency_rain_formation_tot[k]) *
                 ρ0_c[k] *
                 grid.Δz / rho_cloud_liq *
                 3.6 *

src/types.jl

@@ -305,8 +305,6 @@ struct EnvironmentThermodynamics{A1}
     Hvar_rain_dt::A1
     QTvar_rain_dt::A1
     HQTcov_rain_dt::A1
-    qt_tendency_rain_formation::A1
-    θ_liq_ice_tendency_rain_formation::A1
     function EnvironmentThermodynamics(namelist, grid::Grid)
         quadrature_order = parse_namelist(namelist, "thermodynamics", "quadrature_order"; default = 3)
         quadrature_type = parse_namelist(namelist, "thermodynamics", "quadrature_type"; default = "gaussian")
@@ -325,8 +323,6 @@ struct EnvironmentThermodynamics{A1}
         QTvar_rain_dt = center_field(grid)
         HQTcov_rain_dt = center_field(grid)
 
-        qt_tendency_rain_formation = center_field(grid)
-        θ_liq_ice_tendency_rain_formation = center_field(grid)
         A1 = typeof(qt_unsat)
         return new{A1}(
             quadrature_order,
@@ -342,8 +338,6 @@ struct EnvironmentThermodynamics{A1}
             Hvar_rain_dt,
             QTvar_rain_dt,
             HQTcov_rain_dt,
-            qt_tendency_rain_formation,
-            θ_liq_ice_tendency_rain_formation,
         )
     end
 end