Merge #515

515: Additional data structs for snow r=trontrytel a=trontrytel

The last of the peel offs from #506 Hopefully no changes as well. 🤞 

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

integration_tests/utils/main.jl

@@ -119,9 +119,11 @@ cent_aux_vars_edmf(FT, n_up) = (;
             qt_tendency_precip_formation = FT(0),
             θ_liq_ice_tendency_precip_formation = FT(0),
         ),
-        θ_liq_ice_tendency_rain_evap = FT(0),
-        qt_tendency_rain_evap = FT(0),
-        qr_tendency_rain_evap = FT(0),
+        θ_liq_ice_tendency_precip_sinks = FT(0),
+        qt_tendency_precip_sinks = FT(0),
+        qr_tendency_evap = FT(0),
+        qs_tendency_melt = FT(0),
+        qs_tendency_dep_sub = FT(0),
         qr_tendency_advection = FT(0),
         qs_tendency_advection = FT(0),
         en_2m = (;
@@ -194,7 +196,7 @@ cent_prognostic_vars_edmf(FT, n_up) = (;
     turbconv = (;
         en = cent_prognostic_vars_en(FT),
         up = ntuple(i -> cent_prognostic_vars_up(FT), n_up),
-        pr = (; qr = FT(0), qs = FT(0)),
+        pr = (; q_rai = FT(0), q_sno = FT(0)),
     ),
 )
 # cent_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model

src/EDMF_Environment.jl

@@ -6,8 +6,8 @@ function update_env_precip_tendencies(state, k, qt_tendency, θ_liq_ice_tendency
     aux_en.qt_tendency_precip_formation[k] = qt_tendency * aux_en.area[k]
     aux_en.θ_liq_ice_tendency_precip_formation[k] = θ_liq_ice_tendency * aux_en.area[k]
 
-    tendencies_pr.qr[k] += qr_tendency * aux_en.area[k]
-    tendencies_pr.qs[k] += qs_tendency * aux_en.area[k]
+    tendencies_pr.q_rai[k] += qr_tendency * aux_en.area[k]
+    tendencies_pr.q_sno[k] += qs_tendency * aux_en.area[k]
     return
 end
 
@@ -46,8 +46,8 @@ function sgs_mean(en_thermo::EnvironmentThermodynamics, grid, state, en, precip,
         mph = precipitation_formation(
             param_set,
             precip.precipitation_model,
-            prog_pr.qr[k],
-            prog_pr.qs[k],
+            prog_pr.q_rai[k],
+            prog_pr.q_sno[k],
             aux_en.area[k],
             ρ0_c[k],
             dt,
@@ -180,8 +180,8 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
                     mph = precipitation_formation(
                         param_set,
                         precip.precipitation_model,
-                        prog_pr.qr[k],
-                        prog_pr.qs[k],
+                        prog_pr.q_rai[k],
+                        prog_pr.q_sno[k],
                         aux_en.area[k],
                         ρ0_c[k],
                         dt,
@@ -261,8 +261,8 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p
             mph = precipitation_formation(
                 param_set,
                 precip.precipitation_model,
-                prog_pr.qr[k],
-                prog_pr.qs[k],
+                prog_pr.q_rai[k],
+                prog_pr.q_sno[k],
                 aux_en.area[k],
                 ρ0_c[k],
                 dt,

src/EDMF_Precipitation.jl

@@ -1,7 +1,7 @@
 """
 Computes the qr advection (down) tendency
 """
-function compute_rain_advection_tendencies(::PrecipPhysics, grid, state, gm, TS::TimeStepping)
+function compute_precipitation_advection_tendencies(::PrecipPhysics, grid, state, gm, TS::TimeStepping)
     param_set = parameter_set(gm)
     Δz = grid.Δz
     Δt = TS.dt
@@ -19,8 +19,8 @@ function compute_rain_advection_tendencies(::PrecipPhysics, grid, state, gm, TS:
     term_vel_snow = aux_tc.term_vel_snow
 
     @inbounds for k in real_center_indices(grid)
-        term_vel_rain[k] = CM1.terminal_velocity(param_set, rain_type, ρ_0_c[k], prog_pr.qr[k])
-        term_vel_snow[k] = CM1.terminal_velocity(param_set, snow_type, ρ_0_c[k], prog_pr.qs[k])
+        term_vel_rain[k] = CM1.terminal_velocity(param_set, rain_type, ρ_0_c[k], prog_pr.q_rai[k])
+        term_vel_snow[k] = CM1.terminal_velocity(param_set, snow_type, ρ_0_c[k], prog_pr.q_sno[k])
     end
 
     @inbounds for k in reverse(real_center_indices(grid))
@@ -40,8 +40,8 @@ function compute_rain_advection_tendencies(::PrecipPhysics, grid, state, gm, TS:
 
         ρ_0_cut = ccut_downwind(ρ_0_c, grid, k)
 
-        QR_cut = ccut_downwind(prog_pr.qr, grid, k)
-        QS_cut = ccut_downwind(prog_pr.qs, grid, k)
+        QR_cut = ccut_downwind(prog_pr.q_rai, grid, k)
+        QS_cut = ccut_downwind(prog_pr.q_sno, grid, k)
 
         w_cut_rain = ccut_downwind(term_vel_rain, grid, k)
         w_cut_snow = ccut_downwind(term_vel_snow, grid, k)
@@ -58,16 +58,16 @@ function compute_rain_advection_tendencies(::PrecipPhysics, grid, state, gm, TS:
         aux_tc.qr_tendency_advection[k] = ∇ρQRw / ρ_0_c_k
         aux_tc.qs_tendency_advection[k] = ∇ρQSw / ρ_0_c_k
 
-        tendencies_pr.qr[k] += aux_tc.qr_tendency_advection[k]
-        tendencies_pr.qs[k] += aux_tc.qs_tendency_advection[k]
+        tendencies_pr.q_rai[k] += aux_tc.qr_tendency_advection[k]
+        tendencies_pr.q_sno[k] += aux_tc.qs_tendency_advection[k]
     end
     return
 end
 
 """
 Computes the tendencies to θ_liq_ice, qt and qr due to rain evaporation
 """
-function compute_rain_evap_tendencies(::PrecipPhysics, grid, state, gm, TS::TimeStepping)
+function compute_precipitation_sink_tendencies(::PrecipPhysics, grid, state, gm, TS::TimeStepping)
     param_set = parameter_set(gm)
     Δt = TS.dt
     p0_c = center_ref_state(state).p0
@@ -86,13 +86,20 @@ function compute_rain_evap_tendencies(::PrecipPhysics, grid, state, gm, TS::Time
         q = TD.PhasePartition(ts)
 
         # TODO - move limiters elsewhere
-        qt_tendency =
-            min(prog_pr.qr[k] / Δt, -CM1.evaporation_sublimation(param_set, rain_type, q, prog_pr.qr[k], ρ0_c[k], T_gm))
+        rain_evap =
+            -min(
+                prog_pr.q_rai[k] / Δt,
+                -CM1.evaporation_sublimation(param_set, rain_type, q, prog_pr.q_rai[k], ρ0_c[k], T_gm),
+            )
 
-        aux_tc.qt_tendency_rain_evap[k] = qt_tendency
-        aux_tc.qr_tendency_rain_evap[k] = -qt_tendency
-        tendencies_pr.qr[k] += aux_tc.qr_tendency_rain_evap[k]
-        aux_tc.θ_liq_ice_tendency_rain_evap[k] = θ_liq_ice_helper(ts, qt_tendency)
+        aux_tc.qr_tendency_evap[k] = rain_evap
+        aux_tc.qs_tendency_melt[k] = 0.0
+        aux_tc.qs_tendency_dep_sub[k] = 0.0
+
+        aux_tc.qt_tendency_precip_sinks[k] = -rain_evap
+        aux_tc.θ_liq_ice_tendency_precip_sinks[k] = θ_liq_ice_helper(ts, -rain_evap)
+
+        tendencies_pr.q_rai[k] += rain_evap
     end
     return
 end

src/EDMF_Updrafts.jl

@@ -26,15 +26,18 @@ function compute_precipitation_formation_tendencies(
             mph = precipitation_formation(
                 param_set,
                 precip.precipitation_model,
-                prog_pr.qr[k],
-                prog_pr.qs[k],
+                prog_pr.q_rai[k],
+                prog_pr.q_sno[k],
                 aux_up[i].area[k],
                 ρ0_c[k],
                 dt,
                 ts_up,
             )
             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]
+
+            tendencies_pr.q_rai[k] += mph.qr_tendency * aux_up[i].area[k]
+            tendencies_pr.q_sno[k] += mph.qs_tendency * aux_up[i].area[k]
         end
     end
     # TODO - to be deleted once we sum all tendencies elsewhere
@@ -45,7 +48,6 @@ function compute_precipitation_formation_tendencies(
             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

@@ -97,11 +97,11 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
         tendencies_gm.q_tot[k] +=
             aux_bulk.qt_tendency_precip_formation[k] +
             aux_en.qt_tendency_precip_formation[k] +
-            aux_tc.qt_tendency_rain_evap[k]
+            aux_tc.qt_tendency_precip_sinks[k]
         tendencies_gm.θ_liq_ice[k] +=
             aux_bulk.θ_liq_ice_tendency_precip_formation[k] +
             aux_en.θ_liq_ice_tendency_precip_formation[k] +
-            aux_tc.θ_liq_ice_tendency_rain_evap[k]
+            aux_tc.θ_liq_ice_tendency_precip_sinks[k]
     end
 
     aux_up_f = face_aux_updrafts(state)
@@ -242,7 +242,6 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca
     n_updrafts = up.n_updrafts
     prog_gm = center_prog_grid_mean(state)
     prog_en = center_prog_environment(state)
-    has_precip = edmf.Precip.precipitation_model == "clima_1m"
 
     # Update aux / pre-tendencies filters. TODO: combine these into a function that minimizes traversals
     # Some of these methods should probably live in `compute_tendencies`, when written, but we'll
@@ -264,8 +263,8 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca
     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)
-        compute_rain_advection_tendencies(edmf.PrecipPhys, grid, state, gm, TS)
+        compute_precipitation_sink_tendencies(edmf.PrecipPhys, grid, state, gm, TS)
+        compute_precipitation_advection_tendencies(edmf.PrecipPhys, grid, state, gm, TS)
     end
 
     # compute tendencies
@@ -301,8 +300,9 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca
             prog_up[i].ρaθ_liq_ice[k] += Δt * tendencies_up[i].ρaθ_liq_ice[k]
             prog_up[i].ρaq_tot[k] += Δt * tendencies_up[i].ρaq_tot[k]
         end
-        if has_precip
-            prog_pr.qr[k] += tendencies_pr.qr[k] * TS.dt
+        if edmf.Precip.precipitation_model == "clima_1m"
+            prog_pr.q_rai[k] += tendencies_pr.q_rai[k] * TS.dt
+            prog_pr.q_sno[k] += tendencies_pr.q_sno[k] * TS.dt
         end
     end
 

src/diagnostics.jl

@@ -110,7 +110,9 @@ function io_dictionary_aux(state)
         "diffusive_flux_s" => (; dims = ("zf", "t"), group = "profiles", field = face_aux_grid_mean(state).diffusive_flux_s),
         "total_flux_s" => (; dims = ("zf", "t"), group = "profiles", field = face_aux_grid_mean(state).massflux_s .+ face_aux_grid_mean(state).diffusive_flux_s),
 
-        "qr_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_precipitation(state).qr),
+        "qr_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_precipitation(state).q_rai),
+        "qs_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_precipitation(state).q_sno),
+
         "mixing_length" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).mixing_length),
 
         "nh_pressure" => (; dims = ("zf", "t"), group = "profiles", field = face_diagnostics_turbconv(state).nh_pressure),
@@ -257,10 +259,12 @@ function compute_diagnostics!(edmf, gm, grid, state, Case, TS)
     end
 
     precip.mean_rwp = 0.0
+    precip.mean_swp = 0.0
     precip.cutoff_precipitation_rate = 0.0
 
     @inbounds for k in real_center_indices(grid)
-        precip.mean_rwp += ρ0_c[k] * prog_pr.qr[k] * grid.Δz
+        precip.mean_rwp += ρ0_c[k] * prog_pr.q_rai[k] * grid.Δz
+        precip.mean_swp += ρ0_c[k] * prog_pr.q_sno[k] * grid.Δz
 
         # precipitation rate from cutoff microphysics scheme defined as a total amount of removed water
         # per timestep per EDMF surface area [mm/h]

src/io.jl

@@ -50,12 +50,14 @@ end
 
 function initialize_io(precip::PrecipVariables, Stats::NetCDFIO_Stats)
     add_ts(Stats, "rwp_mean")
+    add_ts(Stats, "swp_mean")
     add_ts(Stats, "cutoff_precipitation_rate")
     return
 end
 
 function io(precip::PrecipVariables, grid, state, Stats::NetCDFIO_Stats)
     write_ts(Stats, "rwp_mean", precip.mean_rwp)
+    write_ts(Stats, "swp_mean", precip.mean_rwp)
 
     #TODO - change to rain rate that depends on rain model choice
     write_ts(Stats, "cutoff_precipitation_rate", precip.cutoff_precipitation_rate)

src/types.jl

@@ -182,6 +182,7 @@ end
 Base.@kwdef mutable struct PrecipVariables
     precipitation_model::String = "default_precipitation_model"
     mean_rwp::Float64 = 0
+    mean_swp::Float64 = 0
     cutoff_precipitation_rate::Float64 = 0
 end
 function PrecipVariables(namelist, grid::Grid)