Merge #512

512: Make some rain variables ClimaCore fields r=charleskawczynski a=charleskawczynski

This PR makes some rain variables ClimaCore fields, and removes the `UpdraftThermodynamics` struct. Peeled off from #506.

Co-authored-by: Charles Kawczynski <[email protected]>

integration_tests/utils/main.jl

@@ -82,6 +82,8 @@ cent_aux_vars_up(FT) = (;
     area = FT(0),
     q_tot = FT(0),
     θ_liq_ice = FT(0),
+    θ_liq_ice_tendency_precip_formation = FT(0),
+    qt_tendency_precip_formation = FT(0),
 )
 cent_aux_vars_edmf(FT, n_up) = (;
     turbconv = (;
@@ -95,6 +97,8 @@ cent_aux_vars_edmf(FT, n_up) = (;
             q_ice = FT(0),
             T = FT(0),
             cloud_fraction = FT(0),
+            θ_liq_ice_tendency_precip_formation = FT(0),
+            qt_tendency_precip_formation = FT(0),
         ),
         up = ntuple(i -> cent_aux_vars_up(FT), n_up),
         en = (;

src/EDMF_Environment.jl

@@ -1,10 +1,11 @@
 function update_env_precip_tendencies(en_thermo::EnvironmentThermodynamics, state, k, qt_tendency, θ_liq_ice_tendency)
 
     aux_en = center_aux_environment(state)
+    aux_bulk = center_aux_bulk(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_bulk.qt_tendency_precip_formation[k] = qt_tendency * aux_en.area[k]
+    tendencies_pr.qr[k] += -aux_bulk.qt_tendency_precip_formation[k]
+    aux_bulk.θ_liq_ice_tendency_precip_formation[k] = θ_liq_ice_tendency * aux_en.area[k]
 
     return
 end

src/EDMF_Updrafts.jl

@@ -2,7 +2,6 @@
 Computes tendencies to qt and θ_liq_ice due to precipitation formation
 """
 function compute_precipitation_formation_tendencies(
-    up_thermo::UpdraftThermodynamics,
     grid,
     state,
     up::UpdraftVariables,
@@ -13,6 +12,7 @@ function compute_precipitation_formation_tendencies(
     p0_c = center_ref_state(state).p0
     ρ0_c = center_ref_state(state).ρ0
     aux_up = center_aux_updrafts(state)
+    aux_bulk = center_aux_bulk(state)
     prog_pr = center_prog_precipitation(state)
     tendencies_pr = center_tendencies_precipitation(state)
 
@@ -32,13 +32,19 @@ function compute_precipitation_formation_tendencies(
                 dt,
                 ts_up,
             )
-            up_thermo.qt_tendency_rain_formation[i, k] = mph.qt_tendency * aux_up[i].area[k]
-            up_thermo.θ_liq_ice_tendency_rain_formation[i, k] = mph.θ_liq_ice_tendency * aux_up[i].area[k]
+            aux_up[i].qt_tendency_precip_formation[k] = mph.qt_tendency * aux_up[i].area[k]
+            aux_up[i].θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_tendency * aux_up[i].area[k]
         end
     end
     # TODO - to be deleted once we sum all tendencies elsewhere
-    up_thermo.θ_liq_ice_tendency_rain_formation_tot .= up_sum(up_thermo.θ_liq_ice_tendency_rain_formation)
-    up_thermo.qt_tendency_rain_formation_tot .= up_sum(up_thermo.qt_tendency_rain_formation)
-    parent(tendencies_pr.qr) .+= -up_thermo.qt_tendency_rain_formation_tot
+    for k in real_center_indices(grid)
+        aux_bulk.θ_liq_ice_tendency_precip_formation[k] = 0
+        aux_bulk.qt_tendency_precip_formation[k] = 0
+        for i in 1:(up.n_updrafts)
+            aux_bulk.θ_liq_ice_tendency_precip_formation[k] += aux_up[i].θ_liq_ice_tendency_precip_formation[k]
+            aux_bulk.qt_tendency_precip_formation[k] += aux_up[i].qt_tendency_precip_formation[k]
+        end
+        tendencies_pr.qr[k] += -aux_bulk.qt_tendency_precip_formation[k]
+    end
     return
 end

src/Turbulence_PrognosticTKE.jl

@@ -28,6 +28,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
     aux_en = center_aux_environment(state)
     aux_en_f = face_aux_environment(state)
     aux_up = center_aux_updrafts(state)
+    aux_bulk = center_aux_bulk(state)
     ρ0_f = face_ref_state(state).ρ0
     p0_c = center_ref_state(state).p0
     α0_c = center_ref_state(state).α0
@@ -94,12 +95,12 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
             tendencies_gm.v[k] += gm_V_nudge_k
         end
         tendencies_gm.q_tot[k] +=
-            edmf.UpdThermo.qt_tendency_rain_formation_tot[k] +
-            edmf.EnvThermo.qt_tendency_rain_formation[k] +
+            aux_bulk.qt_tendency_precip_formation[k] +
+            edmf.EnvThermo.qt_tendency_precip_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_bulk.θ_liq_ice_tendency_precip_formation[k] +
+            edmf.EnvThermo.θ_liq_ice_tendency_precip_formation[k] +
             aux_tc.θ_liq_ice_tendency_rain_evap[k]
     end
 
@@ -237,7 +238,6 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca
     up = edmf.UpdVar
     en = edmf.EnvVar
     param_set = parameter_set(gm)
-    up_thermo = edmf.UpdThermo
     en_thermo = edmf.EnvThermo
     n_updrafts = up.n_updrafts
     prog_gm = center_prog_grid_mean(state)
@@ -261,7 +261,7 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca
     parent(tendencies_up_f) .= 0
     parent(tendencies_en) .= 0
     parent(tendencies_pr) .= 0
-    compute_precipitation_formation_tendencies(up_thermo, grid, state, edmf.UpdVar, edmf.Precip, TS.dt, param_set) # causes division error in dry bubble first time step
+    compute_precipitation_formation_tendencies(grid, state, edmf.UpdVar, edmf.Precip, TS.dt, param_set) # causes division error in dry bubble first time step
     microphysics(en_thermo, grid, state, en, edmf.Precip, TS.dt, param_set) # saturation adjustment + rain creation
     if edmf.Precip.precipitation_model == "clima_1m"
         compute_rain_evap_tendencies(edmf.PrecipPhys, grid, state, gm, TS)
@@ -485,7 +485,6 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G
     kf_surf = kf_surface(grid)
 
     up = edmf.UpdVar
-    up_thermo = edmf.UpdThermo
     en = edmf.EnvVar
     aux_up = center_aux_updrafts(state)
     aux_en = center_aux_environment(state)
@@ -516,13 +515,13 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G
             adv = upwind_advection_scalar(ρ0_c, aux_up[i].area, aux_up_f[i].w, aux_up[i].θ_liq_ice, grid, k)
             entr = entr_term * aux_en.θ_liq_ice[k]
             detr = detr_term * aux_up[i].θ_liq_ice[k]
-            rain = ρ0_c[k] * up_thermo.θ_liq_ice_tendency_rain_formation[i, k]
+            rain = ρ0_c[k] * aux_up[i].θ_liq_ice_tendency_precip_formation[k]
             tendencies_up[i].ρaθ_liq_ice[k] = -adv + entr - detr + rain
 
             adv = upwind_advection_scalar(ρ0_c, aux_up[i].area, aux_up_f[i].w, aux_up[i].q_tot, grid, k)
             entr = entr_term * aux_en.q_tot[k]
             detr = detr_term * aux_up[i].q_tot[k]
-            rain = ρ0_c[k] * up_thermo.qt_tendency_rain_formation[i, k]
+            rain = ρ0_c[k] * aux_up[i].qt_tendency_precip_formation[k]
             tendencies_up[i].ρaq_tot[k] = -adv + entr - detr + rain
         end
     end

src/diagnostics.jl

@@ -126,6 +126,9 @@ function io_dictionary_aux(state)
 
         "updraft_cloud_fraction" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).bulk.cloud_fraction),
 
+        "updraft_qt_precip" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_bulk(state).qt_tendency_precip_formation),
+        "updraft_thetal_precip" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_bulk(state).θ_liq_ice_tendency_precip_formation),
+
     )
     return io_dict
 end
@@ -166,7 +169,8 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
     aux_tc_f = face_aux_turbconv(state)
     aux_gm_f = face_aux_grid_mean(state)
     prog_pr = center_prog_precipitation(state)
-    a_up_bulk = center_aux_turbconv(state).bulk.area
+    aux_bulk = center_aux_bulk(state)
+    a_up_bulk = aux_bulk.area
     kc_toa = kc_top_of_atmos(grid)
     gm.cloud_base = grid.zc[kc_toa]
     gm.cloud_top = 0.0
@@ -176,7 +180,6 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
     en = edmf.EnvVar
     precip = edmf.Precip
     en_thermo = edmf.EnvThermo
-    up_thermo = edmf.UpdThermo
     n_updrafts = up.n_updrafts
     diag_tc = center_diagnostics_turbconv(state)
     diag_tc_f = face_diagnostics_turbconv(state)
@@ -263,7 +266,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]) *
+                (en_thermo.qt_tendency_precip_formation[k] + aux_bulk.qt_tendency_precip_formation[k]) *
                 ρ0_c[k] *
                 grid.Δz / rho_cloud_liq *
                 3.6 *

src/dycore_api.jl

@@ -69,6 +69,7 @@ face_aux_turbconv(state) = aux_turbconv(state, FaceField())
 center_aux_updrafts(state) = aux_turbconv(state, CentField()).up
 face_aux_updrafts(state) = aux_turbconv(state, FaceField()).up
 center_aux_environment(state) = aux_turbconv(state, CentField()).en
+center_aux_bulk(state) = aux_turbconv(state, CentField()).bulk
 center_aux_environment_2m(state) = aux_turbconv(state, CentField()).en_2m
 face_aux_environment(state) = aux_turbconv(state, FaceField()).en
 

src/io.jl

@@ -129,8 +129,6 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats)
     add_profile(Stats, "diffusive_tendency_qt")
     add_profile(Stats, "total_flux_h")
     add_profile(Stats, "total_flux_qt")
-    add_profile(Stats, "updraft_qt_precip")
-    add_profile(Stats, "updraft_thetal_precip")
     # Diff mixing lengths: Ignacio
     add_profile(Stats, "ed_length_scheme")
     add_profile(Stats, "mixing_length_ratio")
@@ -155,8 +153,6 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti
     write_profile(Stats, "diffusive_tendency_qt", edmf.diffusive_tendency_qt)
     write_profile(Stats, "total_flux_h", edmf.massflux_h .+ edmf.diffusive_flux_h)
     write_profile(Stats, "total_flux_qt", edmf.massflux_qt .+ edmf.diffusive_flux_qt)
-    write_profile(Stats, "updraft_qt_precip", edmf.UpdThermo.qt_tendency_rain_formation_tot)
-    write_profile(Stats, "updraft_thetal_precip", edmf.UpdThermo.θ_liq_ice_tendency_rain_formation_tot)
 
     #Different mixing lengths : Ignacio
     write_profile(Stats, "ed_length_scheme", edmf.mls)

src/types.jl

@@ -250,31 +250,7 @@ function GridMeanVariables(namelist, grid::Grid, param_set::PS) where {PS}
 end
 
 
-struct UpdraftThermodynamics{A1, A2}
-    n_updraft::Int
-    θ_liq_ice_tendency_rain_formation::A2
-    qt_tendency_rain_formation::A2
-    θ_liq_ice_tendency_rain_formation_tot::A1
-    qt_tendency_rain_formation_tot::A1
-    function UpdraftThermodynamics(n_updraft::Int, grid::Grid)
-        # tendencies from each updraft
-        θ_liq_ice_tendency_rain_formation = center_field(grid, n_updraft)
-        qt_tendency_rain_formation = center_field(grid, n_updraft)
-        # tendencies from all updrafts
-        θ_liq_ice_tendency_rain_formation_tot = center_field(grid)
-        qt_tendency_rain_formation_tot = center_field(grid)
-
-        A1 = typeof(θ_liq_ice_tendency_rain_formation_tot)
-        A2 = typeof(θ_liq_ice_tendency_rain_formation)
-        return new{A1, A2}(
-            n_updraft,
-            θ_liq_ice_tendency_rain_formation,
-            qt_tendency_rain_formation,
-            θ_liq_ice_tendency_rain_formation_tot,
-            qt_tendency_rain_formation_tot,
-        )
-    end
-end
+struct UpdraftThermodynamics end
 
 Base.@kwdef mutable struct EnvironmentVariables
     cloud_base::Float64 = 0
@@ -304,8 +280,8 @@ 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
+    qt_tendency_precip_formation::A1
+    θ_liq_ice_tendency_precip_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")
@@ -324,8 +300,8 @@ 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)
+        qt_tendency_precip_formation = center_field(grid)
+        θ_liq_ice_tendency_precip_formation = center_field(grid)
         A1 = typeof(qt_unsat)
         return new{A1}(
             quadrature_order,
@@ -341,8 +317,8 @@ struct EnvironmentThermodynamics{A1}
             Hvar_rain_dt,
             QTvar_rain_dt,
             HQTcov_rain_dt,
-            qt_tendency_rain_formation,
-            θ_liq_ice_tendency_rain_formation,
+            qt_tendency_precip_formation,
+            θ_liq_ice_tendency_precip_formation,
         )
     end
 end
@@ -512,7 +488,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2}
     Precip::PrecipVariables
     PrecipPhys::PrecipPhysics
     UpdVar::UpdraftVariables
-    UpdThermo::UpdraftThermodynamics
     EnvVar::EnvironmentVariables
     EnvThermo::EnvironmentThermodynamics
     entr_sc::A2
@@ -612,8 +587,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2}
 
         # Create the updraft variable class (major diagnostic and prognostic variables)
         UpdVar = UpdraftVariables(n_updrafts, namelist, grid)
-        # Create the class for updraft thermodynamics
-        UpdThermo = UpdraftThermodynamics(n_updrafts, grid)
 
         # Create the environment variable class (major diagnostic and prognostic variables)
         EnvVar = EnvironmentVariables(namelist, grid)
@@ -738,7 +711,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2}
             Precip,
             PrecipPhys,
             UpdVar,
-            UpdThermo,
             EnvVar,
             EnvThermo,
             entr_sc,