Use ClimaCore vertical integrals for lwp etc

diagnostics/compute_diagnostics.jl

@@ -52,10 +52,9 @@ diagnostics and still run, at which point we'll be able to export
 the state, auxiliary fields (which the state does depend on), and
 tendencies.
 =#
-function compute_diagnostics!(edmf, gm, grid, state, diagnostics, Case)
+function compute_diagnostics!(edmf, gm, grid, state, diagnostics, Case, Stats)
     FT = eltype(grid)
-    gm.lwp = 0.0
-    gm.iwp = 0.0
+    N_up = TC.n_updrafts(edmf)
     ρ0_c = TC.center_ref_state(state).ρ0
     p0_c = TC.center_ref_state(state).p0
     aux_gm = TC.center_aux_grid_mean(state)
@@ -81,8 +80,6 @@ function compute_diagnostics!(edmf, gm, grid, state, diagnostics, Case)
     diag_tc_f = face_diagnostics_turbconv(diagnostics)
 
     @inbounds for k in TC.real_center_indices(grid)
-        gm.lwp += ρ0_c[k] * aux_gm.q_liq[k] * grid.Δz
-        gm.iwp += ρ0_c[k] * aux_gm.q_ice[k] * grid.Δz
         if TD.has_condensate(aux_gm.q_liq[k] + aux_gm.q_ice[k])
             gm.cloud_base = min(gm.cloud_base, grid.zc[k])
             gm.cloud_top = max(gm.cloud_top, grid.zc[k])
@@ -113,19 +110,13 @@ function compute_diagnostics!(edmf, gm, grid, state, diagnostics, Case)
         @. massflux_s += aux_up_f[i].massflux * (If(aux_up[i].s) - If(aux_en.s))
     end
 
-    up.lwp = 0.0
-    up.iwp = 0.0
-
     @inbounds for i in 1:(up.n_updrafts)
         up.cloud_base[i] = TC.zc_toa(grid)
         up.cloud_top[i] = 0.0
         up.cloud_cover[i] = 0.0
 
         @inbounds for k in TC.real_center_indices(grid)
             if aux_up[i].area[k] > 1e-3
-                up.lwp += ρ0_c[k] * aux_up[i].q_liq[k] * aux_up[i].area[k] * grid.Δz
-                up.iwp += ρ0_c[k] * aux_up[i].q_ice[k] * aux_up[i].area[k] * grid.Δz
-
                 if TD.has_condensate(aux_up[i].q_liq[k] + aux_up[i].q_ice[k])
                     up.cloud_base[i] = min(up.cloud_base[i], grid.zc[k])
                     up.cloud_top[i] = max(up.cloud_top[i], grid.zc[k])
@@ -138,40 +129,15 @@ function compute_diagnostics!(edmf, gm, grid, state, diagnostics, Case)
     en.cloud_top = 0.0
     en.cloud_base = TC.zc_toa(grid)
     en.cloud_cover = 0.0
-    en.lwp = 0.0
-    en.iwp = 0.0
 
     @inbounds for k in TC.real_center_indices(grid)
-        en.lwp += ρ0_c[k] * aux_en.q_liq[k] * aux_en.area[k] * grid.Δz
-        en.iwp += ρ0_c[k] * aux_en.q_ice[k] * aux_en.area[k] * grid.Δz
-
         if TD.has_condensate(aux_en.q_liq[k] + aux_en.q_ice[k]) && aux_en.area[k] > 1e-6
             en.cloud_base = min(en.cloud_base, grid.zc[k])
             en.cloud_top = max(en.cloud_top, grid.zc[k])
             en.cloud_cover = max(en.cloud_cover, aux_en.area[k] * aux_en.cloud_fraction[k])
         end
     end
 
-    precip.mean_rwp = 0.0
-    precip.mean_swp = 0.0
-    precip.cutoff_precipitation_rate = 0.0
-
-    @inbounds for k in TC.real_center_indices(grid)
-        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]
-        if (precip.precipitation_model == "cutoff")
-            precip.cutoff_precipitation_rate -=
-                (aux_en.qt_tendency_precip_formation[k] + aux_bulk.qt_tendency_precip_formation[k]) *
-                ρ0_c[k] *
-                grid.Δz / TC.rho_cloud_liq *
-                3.6 *
-                1e6
-        end
-    end
-
     If = CCO.InterpolateF2C()
     parent(diag_tc.massflux) .= 0
     @inbounds for i in 1:(edmf.n_updrafts)
@@ -228,5 +194,27 @@ function compute_diagnostics!(edmf, gm, grid, state, diagnostics, Case)
     TC.update_cloud_frac(edmf, grid, state, gm)
 
 
+    TC.write_ts(Stats, "lwp_mean", sum(ρ0_c .* aux_gm.q_liq))
+    TC.write_ts(Stats, "iwp_mean", sum(ρ0_c .* aux_gm.q_ice))
+    TC.write_ts(Stats, "env_lwp", sum(ρ0_c .* aux_en.q_liq .* aux_en.area))
+    TC.write_ts(Stats, "env_iwp", sum(ρ0_c .* aux_en.q_ice .* aux_en.area))
+
+    TC.write_ts(Stats, "rwp_mean", sum(ρ0_c .* prog_pr.q_rai))
+    TC.write_ts(Stats, "swp_mean", sum(ρ0_c .* prog_pr.q_sno))
+    #TODO - change to rain rate that depends on rain model choice
+
+    # TODO: Move TC.rho_cloud_liq to CLIMAParameters
+    if (precip.precipitation_model == "cutoff")
+        f =
+            (aux_en.qt_tendency_precip_formation .+ aux_bulk.qt_tendency_precip_formation) .* ρ0_c ./
+            TC.rho_cloud_liq .* 3.6 .* 1e6
+        TC.write_ts(Stats, "cutoff_precipitation_rate", sum(f))
+    end
+
+    lwp = sum(i -> sum(ρ0_c .* aux_up[i].q_liq .* aux_up[i].area .* (aux_up[i].area .> 1e-3)), 1:N_up)
+    iwp = sum(i -> sum(ρ0_c .* aux_up[i].q_ice .* aux_up[i].area .* (aux_up[i].area .> 1e-3)), 1:N_up)
+    TC.write_ts(Stats, "updraft_lwp", lwp)
+    TC.write_ts(Stats, "updraft_iwp", iwp)
+
     return
 end

driver/callbacks.jl

@@ -19,7 +19,7 @@ function affect_io!(integrator)
 
     param_set = TC.parameter_set(gm)
     # TODO: is this the best location to call diagnostics?
-    compute_diagnostics!(edmf, gm, grid, state, diagnostics, case)
+    compute_diagnostics!(edmf, gm, grid, state, diagnostics, case, Stats)
 
     # TODO: remove `vars` hack that avoids
     # https://github.com/Alexander-Barth/NCDatasets.jl/issues/135

driver/initial_conditions.jl

@@ -89,7 +89,6 @@ function initialize_updrafts(edmf, grid, state, up::TC.UpdraftVariables, gm::TC.
 end
 
 function initialize_updrafts_DryBubble(edmf, grid, state, up::TC.UpdraftVariables, gm::TC.GridMeanVariables)
-    dz = grid.Δz
 
     # criterion 2: b>1e-4
     #! format: off

src/Turbulence_PrognosticTKE.jl

@@ -345,7 +345,6 @@ end
 function compute_updraft_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, Case::CasesBase)
     kc_surf = kc_surface(grid)
 
-    Δzi = grid.Δzi
     zLL = grid.zc[kc_surf]
     ustar = Case.Sur.ustar
     oblength = Case.Sur.obukhov_length
@@ -356,7 +355,7 @@ function compute_updraft_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, Case:
 
     if Case.Sur.bflux > 0.0
         a_total = edmf.surface_area
-        edmf.entr_surface_bc = 2.0 * Δzi
+        edmf.entr_surface_bc = 1 / zLL
         edmf.detr_surface_bc = 0.0
         a_ = a_total / edmf.n_updrafts
         @inbounds for i in 1:(edmf.n_updrafts)
@@ -369,7 +368,7 @@ function compute_updraft_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, Case:
         end
     else
         edmf.entr_surface_bc = 0.0
-        edmf.detr_surface_bc = 2.0 * Δzi
+        edmf.detr_surface_bc = 1 / zLL
         @inbounds for i in 1:(edmf.n_updrafts)
             edmf.area_surface_bc[i] = 0
             edmf.w_surface_bc[i] = 0.0

src/io.jl

@@ -19,8 +19,6 @@ function io(en::EnvironmentVariables, grid, state, Stats::NetCDFIO_Stats)
     write_ts(Stats, "env_cloud_cover", en.cloud_cover)
     write_ts(Stats, "env_cloud_base", en.cloud_base)
     write_ts(Stats, "env_cloud_top", en.cloud_top)
-    write_ts(Stats, "env_lwp", en.lwp)
-    write_ts(Stats, "env_iwp", en.iwp)
     return
 end
 
@@ -31,14 +29,7 @@ function initialize_io(precip::PrecipVariables, Stats::NetCDFIO_Stats)
     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)
-    return
-end
+io(precip::PrecipVariables, grid, state, Stats::NetCDFIO_Stats) = nothing
 
 function initialize_io(up::UpdraftVariables, Stats::NetCDFIO_Stats)
 
@@ -58,8 +49,6 @@ function io(up::UpdraftVariables, grid, state, Stats::NetCDFIO_Stats)
     write_ts(Stats, "updraft_cloud_cover", sum(up.cloud_cover))
     write_ts(Stats, "updraft_cloud_base", minimum(abs.(up.cloud_base)))
     write_ts(Stats, "updraft_cloud_top", maximum(abs.(up.cloud_top)))
-    write_ts(Stats, "updraft_lwp", up.lwp)
-    write_ts(Stats, "updraft_iwp", up.iwp)
     return
 end
 
@@ -74,9 +63,6 @@ end
 
 function io(gm::GridMeanVariables, grid, state, Stats::NetCDFIO_Stats)
     write_ts(Stats, "cloud_cover_mean", gm.cloud_cover)
-
-    write_ts(Stats, "lwp_mean", gm.lwp)
-    write_ts(Stats, "iwp_mean", gm.iwp)
     write_ts(Stats, "cloud_base_mean", gm.cloud_base)
     write_ts(Stats, "cloud_top_mean", gm.cloud_top)
     return

src/types.jl

@@ -176,11 +176,8 @@ Base.@kwdef struct MinDisspLen{FT}
     b_exch::FT
 end
 
-Base.@kwdef mutable struct PrecipVariables
+Base.@kwdef 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)
 
@@ -199,14 +196,12 @@ function PrecipVariables(namelist, grid::Grid)
     return PrecipVariables(; precipitation_model)
 end
 
-mutable struct UpdraftVariables{A1}
+struct UpdraftVariables{A1}
     n_updrafts::Int
     cloud_base::A1
     cloud_top::A1
     cloud_cover::A1
     updraft_top::A1
-    lwp::Float64
-    iwp::Float64
     function UpdraftVariables(nu, namelist, grid::Grid)
         n_updrafts = nu
 
@@ -216,42 +211,32 @@ mutable struct UpdraftVariables{A1}
         cloud_cover = zeros(nu)
         updraft_top = zeros(nu)
 
-        lwp = 0.0
-        iwp = 0.0
-
         A1 = typeof(cloud_base)
-        return new{A1}(n_updrafts, cloud_base, cloud_top, cloud_cover, updraft_top, lwp, iwp)
+        return new{A1}(n_updrafts, cloud_base, cloud_top, cloud_cover, updraft_top)
     end
 end
 
 Base.@kwdef mutable struct GridMeanVariables{PS}
     param_set::PS
-    lwp::Float64
-    iwp::Float64
     cloud_base::Float64
     cloud_top::Float64
     cloud_cover::Float64
     EnvThermo_scheme::String
 end
 function GridMeanVariables(namelist, grid::Grid, param_set::PS) where {PS}
-    lwp = 0.0
-    iwp = 0.0
-
     cloud_base = 0.0
     cloud_top = 0.0
     cloud_cover = 0.0
 
     EnvThermo_scheme = parse_namelist(namelist, "thermodynamics", "sgs"; default = "mean")
 
-    return GridMeanVariables(; param_set, lwp, iwp, cloud_base, cloud_top, cloud_cover, EnvThermo_scheme)
+    return GridMeanVariables(; param_set, cloud_base, cloud_top, cloud_cover, EnvThermo_scheme)
 end
 
 Base.@kwdef mutable struct EnvironmentVariables
     cloud_base::Float64 = 0
     cloud_top::Float64 = 0
     cloud_cover::Float64 = 0
-    lwp::Float64 = 0
-    iwp::Float64 = 0
     EnvThermo_scheme::String = "default_EnvThermo_scheme"
 end
 function EnvironmentVariables(namelist, grid::Grid)