From ce8b17739c372442a546232c2eb90ccd0b747902 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Sun, 17 Oct 2021 10:37:50 -0700 Subject: [PATCH 1/2] Make some edmf variables fields --- integration_tests/utils/main.jl | 5 +++- src/Operators.jl | 6 ++-- src/Turbulence_PrognosticTKE.jl | 53 +++++++++++++++------------------ src/diagnostics.jl | 2 ++ src/types.jl | 16 ---------- src/update_aux.jl | 4 +-- 6 files changed, 35 insertions(+), 51 deletions(-) diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl index b05cd235b..6c619a04b 100644 --- a/integration_tests/utils/main.jl +++ b/integration_tests/utils/main.jl @@ -55,13 +55,16 @@ cent_aux_vars_edmf(FT, n_up) = (; turbconv = (; bulk = (; area = FT(0), θ_liq_ice = FT(0), RH = FT(0), buoy = FT(0), q_tot = FT(0), q_liq = FT(0), T = FT(0)), up = ntuple(i -> cent_aux_vars_up(FT), n_up), + en = (; area = FT(0)), + KM = FT(0), + KH = FT(0), ), ) cent_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT)..., cent_aux_vars_edmf(FT, n_up)...) # Face only face_aux_vars_gm(FT) = () -face_aux_vars_edmf(FT, n_up) = (; turbconv = (; bulk = (; w = FT(0)))) +face_aux_vars_edmf(FT, n_up) = (; turbconv = (; bulk = (; w = FT(0)), ρ_ae_KM = FT(0), ρ_ae_KH = FT(0))) face_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., face_aux_vars_gm(FT)..., face_aux_vars_edmf(FT, n_up)...) ##### diff --git a/src/Operators.jl b/src/Operators.jl index d174f840e..c1a6c7e8c 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -619,8 +619,6 @@ function construct_tridiag_diffusion_en( param_set::APS, state, TS, - KM, - KH, w_en, tke_en, n_updrafts::Int, @@ -646,10 +644,12 @@ function construct_tridiag_diffusion_en( prog_up = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) - ae = vec(1 .- aux_tc.bulk.area) + ae = 1 .- aux_tc.bulk.area rho_ae_K_m = face_field(grid) w_en_c = center_field(grid) D_env = 0.0 + KM = center_aux_tc(state).KM + KH = center_aux_tc(state).KH aeK = is_tke ? ae .* KM : ae .* KH aeK_bcs = (; bottom = SetValue(aeK[kc_surf]), top = SetValue(aeK[kc_toa])) diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index ecfe95aa0..59b8c73d4 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -16,8 +16,6 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats) initialize_io(edmf.EnvVar, Stats) initialize_io(edmf.Rain, Stats) - add_profile(Stats, "eddy_viscosity") - add_profile(Stats, "eddy_diffusivity") add_profile(Stats, "entrainment_sc") add_profile(Stats, "detrainment_sc") add_profile(Stats, "nh_pressure") @@ -107,8 +105,6 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti aux_tc = center_aux_tc(state) a_up_bulk = aux_tc.bulk.area - write_profile(Stats, "eddy_viscosity", diffusivity_m(edmf).values) - write_profile(Stats, "eddy_diffusivity", diffusivity_h(edmf).values) write_ts(Stats, "rd", StatsBase.mean(edmf.pressure_plume_spacing)) @inbounds for k in real_center_indices(grid) @@ -242,6 +238,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm, parent(tendencies_gm.θ_liq_ice) .= 0 param_set = parameter_set(gm) prog_gm = center_prog_grid_mean(state) + aux_en = center_aux_environment(state) ρ0_f = face_ref_state(state).ρ0 p0_c = center_ref_state(state).p0 α0_c = center_ref_state(state).α0 @@ -368,10 +365,10 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm, tendencies_gm.q_tot[k] += mf_tend_qt end - aeKHq_tot_bc = Case.Sur.rho_qtflux / edmf.ae[kc_surf] - aeKHθ_liq_ice_bc = Case.Sur.rho_hflux / edmf.ae[kc_surf] - aeKHu_bc = Case.Sur.rho_uflux / edmf.ae[kc_surf] - aeKHv_bc = Case.Sur.rho_vflux / edmf.ae[kc_surf] + aeKHq_tot_bc = Case.Sur.rho_qtflux / aux_en.area[kc_surf] + aeKHθ_liq_ice_bc = Case.Sur.rho_hflux / aux_en.area[kc_surf] + aeKHu_bc = Case.Sur.rho_uflux / aux_en.area[kc_surf] + aeKHv_bc = Case.Sur.rho_vflux / aux_en.area[kc_surf] @inbounds for k in real_center_indices(grid) aeKH_q_tot_cut = dual_faces(edmf.diffusive_flux_qt, grid, k) ∇aeKH_q_tot = ∇f2c(aeKH_q_tot_cut, grid, k; bottom = SetValue(aeKHq_tot_bc), top = SetValue(0)) @@ -402,11 +399,13 @@ function compute_diffusive_fluxes( ) ρ0_f = face_ref_state(state).ρ0 aux_tc = center_aux_tc(state) - edmf.ae .= 1 .- vec(aux_tc.bulk.area) # area of environment - KM = diffusivity_m(edmf).values - KH = diffusivity_h(edmf).values - aeKM = edmf.ae .* KM - aeKH = edmf.ae .* KH + aux_tc_f = face_aux_tc(state) + aux_en = center_aux_environment(state) + aux_en.area .= 1 .- aux_tc.bulk.area # area of environment + KM = center_aux_tc(state).KM + KH = center_aux_tc(state).KH + aeKM = aux_en.area .* KM + aeKH = aux_en.area .* KH kc_surf = kc_surface(grid) kc_toa = kc_top_of_atmos(grid) kf_surf = kf_surface(grid) @@ -415,31 +414,31 @@ function compute_diffusive_fluxes( aeKH_bcs = (; bottom = SetValue(aeKH[kc_surf]), top = SetValue(aeKH[kc_toa])) @inbounds for k in real_face_indices(grid) - edmf.rho_ae_KH[k] = interpc2f(aeKH, grid, k; aeKH_bcs...) * ρ0_f[k] - edmf.rho_ae_KM[k] = interpc2f(aeKM, grid, k; aeKM_bcs...) * ρ0_f[k] + aux_tc_f.ρ_ae_KH[k] = interpc2f(aeKH, grid, k; aeKH_bcs...) * ρ0_f[k] + aux_tc_f.ρ_ae_KM[k] = interpc2f(aeKM, grid, k; aeKM_bcs...) * ρ0_f[k] end - aeKHq_tot_bc = -Case.Sur.rho_qtflux / edmf.ae[kc_surf] / edmf.rho_ae_KH[kc_surf] - aeKHθ_liq_ice_bc = -Case.Sur.rho_hflux / edmf.ae[kc_surf] / edmf.rho_ae_KH[kc_surf] - aeKHu_bc = -Case.Sur.rho_uflux / edmf.ae[kc_surf] / edmf.rho_ae_KM[kc_surf] - aeKHv_bc = -Case.Sur.rho_vflux / edmf.ae[kc_surf] / edmf.rho_ae_KM[kc_surf] + aeKHq_tot_bc = -Case.Sur.rho_qtflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf] + aeKHθ_liq_ice_bc = -Case.Sur.rho_hflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf] + aeKHu_bc = -Case.Sur.rho_uflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf] + aeKHv_bc = -Case.Sur.rho_vflux / aux_en.area[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf] @inbounds for k in real_face_indices(grid) q_dual = dual_centers(edmf.EnvVar.QT.values, grid, k) ∇q_tot_f = ∇c2f(q_dual, grid, k; bottom = SetGradient(aeKHq_tot_bc), top = SetGradient(0)) - edmf.diffusive_flux_qt[k] = -edmf.rho_ae_KH[k] * ∇q_tot_f + edmf.diffusive_flux_qt[k] = -aux_tc_f.ρ_ae_KH[k] * ∇q_tot_f θ_liq_ice_dual = dual_centers(edmf.EnvVar.H.values, grid, k) ∇θ_liq_ice_f = ∇c2f(θ_liq_ice_dual, grid, k; bottom = SetGradient(aeKHθ_liq_ice_bc), top = SetGradient(0)) - edmf.diffusive_flux_h[k] = -edmf.rho_ae_KH[k] * ∇θ_liq_ice_f + edmf.diffusive_flux_h[k] = -aux_tc_f.ρ_ae_KH[k] * ∇θ_liq_ice_f u_dual = dual_centers(prog_gm.u, grid, k) ∇u_f = ∇c2f(u_dual, grid, k; bottom = SetGradient(aeKHu_bc), top = SetGradient(0)) - edmf.diffusive_flux_u[k] = -edmf.rho_ae_KM[k] * ∇u_f + edmf.diffusive_flux_u[k] = -aux_tc_f.ρ_ae_KM[k] * ∇u_f v_dual = dual_centers(prog_gm.v, grid, k) ∇v_f = ∇c2f(v_dual, grid, k; bottom = SetGradient(aeKHv_bc), top = SetGradient(0)) - edmf.diffusive_flux_v[k] = -edmf.rho_ae_KM[k] * ∇v_f + edmf.diffusive_flux_v[k] = -aux_tc_f.ρ_ae_KM[k] * ∇v_f end return end @@ -478,15 +477,11 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca # Matrix is the same for all variables that use the same eddy diffusivity, we can construct once and reuse prog_up = center_prog_updrafts(state) - KM = diffusivity_m(edmf).values - KH = diffusivity_h(edmf).values common_args = ( grid, param_set, state, TS, - KM, - KH, en.W.values, en.TKE.values, up.n_updrafts, @@ -907,12 +902,12 @@ function compute_covariance_shear( aux_tc = center_aux_tc(state) ae = 1 .- aux_tc.bulk.area # area of environment - KH = diffusivity_h(edmf).values + aux_tc = center_aux_tc(state) ρ0_c = center_ref_state(state).ρ0 prog_gm = center_prog_grid_mean(state) is_tke = Covar.name == "tke" tke_factor = is_tke ? 0.5 : 1 - k_eddy = is_tke ? diffusivity_m(edmf).values : diffusivity_h(edmf).values + k_eddy = is_tke ? aux_tc.KM : aux_tc.KH if is_tke @inbounds for k in real_center_indices(grid) diff --git a/src/diagnostics.jl b/src/diagnostics.jl index d553c65b2..d3b89904e 100644 --- a/src/diagnostics.jl +++ b/src/diagnostics.jl @@ -59,6 +59,8 @@ function io_dictionary_aux(state) "w_mean" => (; dims = ("zf", "t"), group = "profiles", field = face_prog_grid_mean(state).w), "qt_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_grid_mean(state).q_tot), "thetal_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_grid_mean(state).θ_liq_ice), + "eddy_viscosity" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_tc(state).KM), + "eddy_diffusivity" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_tc(state).KH), ) return io_dict end diff --git a/src/types.jl b/src/types.jl index 4ccf620ca..9ae938430 100644 --- a/src/types.jl +++ b/src/types.jl @@ -645,8 +645,6 @@ function center_field_tridiagonal_matrix(grid::Grid) end mutable struct EDMF_PrognosticTKE{A1, A2, IE} - KM::VariableDiagnostic - KH::VariableDiagnostic Ri_bulk_crit::Float64 zi::Float64 n_updrafts::Int @@ -681,9 +679,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} m::A2 mixing_length::A1 implicit_eqs::IE - ae::A1 - rho_ae_KM::A1 - rho_ae_KH::A1 horiz_K_eddy::A2 area_surface_bc::A1 w_surface_bc::A1 @@ -712,8 +707,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} detr_surface_bc::Float64 sde_model::sde_struct function EDMF_PrognosticTKE(namelist, grid::Grid, param_set::PS) where {PS} - KM = VariableDiagnostic(grid, "half", "diffusivity", "m^2/s") # eddy viscosity - KH = VariableDiagnostic(grid, "half", "viscosity", "m^2/s") # eddy diffusivity # get values from namelist prandtl_number = namelist["turbulence"]["prandtl_number_0"] Ri_bulk_crit = namelist["turbulence"]["Ri_bulk_crit"] @@ -822,10 +815,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} mixing_length = center_field(grid) horiz_K_eddy = center_field(grid, n_updrafts) - ae = center_field(grid) - rho_ae_KM = face_field(grid) - rho_ae_KH = face_field(grid) - # Near-surface BC of updraft area fraction area_surface_bc = zeros(n_updrafts) w_surface_bc = zeros(n_updrafts) @@ -887,8 +876,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} A2 = typeof(horiz_K_eddy) IE = typeof(implicit_eqs) return new{A1, A2, IE}( - KM, - KH, Ri_bulk_crit, zi, n_updrafts, @@ -923,9 +910,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} m, mixing_length, implicit_eqs, - ae, - rho_ae_KM, - rho_ae_KH, horiz_K_eddy, area_surface_bc, w_surface_bc, diff --git a/src/update_aux.jl b/src/update_aux.jl index 88e6d4b3d..4a8d5c6f7 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -13,8 +13,8 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) α0_c = center_ref_state(state).α0 g = CPP.grav(param_set) c_m = CPEDMF.c_m(param_set) - KM = diffusivity_m(edmf).values - KH = diffusivity_h(edmf).values + KM = center_aux_tc(state).KM + KH = center_aux_tc(state).KH surface = Case.Sur obukhov_length = surface.obukhov_length FT = eltype(grid) From 82f577b2a3d0f12eb9f52724a4d64ea5a4a9b5cd Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Sun, 17 Oct 2021 13:10:49 -0700 Subject: [PATCH 2/2] Make 2m env variables fields --- integration_tests/utils/main.jl | 18 ++- src/EDMF_Environment.jl | 39 +++--- src/Operators.jl | 4 +- src/Turbulence_PrognosticTKE.jl | 226 ++++++++++++++------------------ src/Variables.jl | 6 +- src/diagnostics.jl | 36 +++++ src/dycore_api.jl | 1 + src/types.jl | 68 +--------- src/update_aux.jl | 60 +++++---- 9 files changed, 207 insertions(+), 251 deletions(-) diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl index 6c619a04b..a693b2b4f 100644 --- a/integration_tests/utils/main.jl +++ b/integration_tests/utils/main.jl @@ -50,12 +50,28 @@ cent_aux_vars_gm(FT) = (; W_third_m = FT(0), QT_third_m = FT(0), ) +cent_aux_vars_en_2m(FT) = (; + dissipation = FT(0), + shear = FT(0), + entr_gain = FT(0), + detr_loss = FT(0), + press = FT(0), + buoy = FT(0), + interdomain = FT(0), + rain_src = FT(0), +) cent_aux_vars_up(FT) = (; q_liq = FT(0), T = FT(0), RH = FT(0), buoy = FT(0)) cent_aux_vars_edmf(FT, n_up) = (; turbconv = (; bulk = (; area = FT(0), θ_liq_ice = FT(0), RH = FT(0), buoy = FT(0), q_tot = FT(0), q_liq = FT(0), T = FT(0)), up = ntuple(i -> cent_aux_vars_up(FT), n_up), en = (; area = FT(0)), + en_2m = (; + tke = cent_aux_vars_en_2m(FT), + Hvar = cent_aux_vars_en_2m(FT), + QTvar = cent_aux_vars_en_2m(FT), + HQTcov = cent_aux_vars_en_2m(FT), + ), KM = FT(0), KH = FT(0), ), @@ -89,7 +105,7 @@ face_diagnostic_vars(FT, n_up) = (; face_diagnostic_vars_gm(FT)..., face_diagnos cent_prognostic_vars(FT, n_up) = (; cent_prognostic_vars_gm(FT)..., cent_prognostic_vars_edmf(FT, n_up)...) cent_prognostic_vars_gm(FT) = (; u = FT(0), v = FT(0), θ_liq_ice = FT(0), q_tot = FT(0)) cent_prognostic_vars_up(FT) = (; area = FT(0), θ_liq_ice = FT(0), q_tot = FT(0)) -cent_prognostic_vars_en(FT) = () +cent_prognostic_vars_en(FT) = (; tke = FT(0), Hvar = FT(0), QTvar = FT(0), HQTcov = FT(0)) cent_prognostic_vars_edmf(FT, n_up) = (; turbconv = (; en = cent_prognostic_vars_en(FT), up = ntuple(i -> cent_prognostic_vars_up(FT), n_up))) # cent_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index a34c50222..7911ecb57 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -6,10 +6,6 @@ function initialize_io(en::EnvironmentVariables, Stats::NetCDFIO_Stats) add_profile(Stats, "env_temperature") add_profile(Stats, "env_RH") add_profile(Stats, "env_thetal") - add_profile(Stats, "env_tke") - add_profile(Stats, "env_Hvar") - add_profile(Stats, "env_QTvar") - add_profile(Stats, "env_HQTcov") add_profile(Stats, "env_cloud_fraction") @@ -29,10 +25,6 @@ function io(en::EnvironmentVariables, grid, state, Stats::NetCDFIO_Stats) write_profile(Stats, "env_temperature", en.T.values) write_profile(Stats, "env_RH", en.RH.values) write_profile(Stats, "env_thetal", en.H.values) - write_profile(Stats, "env_tke", en.TKE.values) - write_profile(Stats, "env_Hvar", en.Hvar.values) - write_profile(Stats, "env_QTvar", en.QTvar.values) - write_profile(Stats, "env_HQTcov", en.HQTcov.values) write_profile(Stats, "env_cloud_fraction", en.cloud_fraction.values) @@ -111,6 +103,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r a, w = FastGaussQuadrature.gausshermite(en_thermo.quadrature_order) p0_c = center_ref_state(state).p0 ρ0_c = center_ref_state(state).ρ0 + prog_en = center_prog_environment(state) #TODO - remember you output source terms multipierd by dt (bec. of instanteneous autoconcv) #TODO - add tendencies for gm H, QT and QR due to rain @@ -139,30 +132,28 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r @inbounds for k in real_center_indices(grid) if ( - en.QTvar.values[k] > epsilon && - en.Hvar.values[k] > epsilon && - abs(en.HQTcov.values[k]) > epsilon && + prog_en.QTvar[k] > epsilon && + prog_en.Hvar[k] > epsilon && + abs(prog_en.HQTcov[k]) > epsilon && en.QT.values[k] > epsilon && - sqrt(en.QTvar.values[k]) < en.QT.values[k] + sqrt(prog_en.QTvar[k]) < en.QT.values[k] ) if en_thermo.quadrature_type == "log-normal" # Lognormal parameters (mu, sd) from mean and variance - sd_q = sqrt(log(en.QTvar.values[k] / en.QT.values[k] / en.QT.values[k] + 1.0)) - sd_h = sqrt(log(en.Hvar.values[k] / en.H.values[k] / en.H.values[k] + 1.0)) + sd_q = sqrt(log(prog_en.QTvar[k] / en.QT.values[k] / en.QT.values[k] + 1.0)) + sd_h = sqrt(log(prog_en.Hvar[k] / en.H.values[k] / en.H.values[k] + 1.0)) # Enforce Schwarz"s inequality - corr = max(min(en.HQTcov.values[k] / sqrt(en.Hvar.values[k] * en.QTvar.values[k]), 1.0), -1.0) - sd2_hq = - log(corr * sqrt(en.Hvar.values[k] * en.QTvar.values[k]) / en.H.values[k] / en.QT.values[k] + 1.0) + corr = max(min(prog_en.HQTcov[k] / sqrt(prog_en.Hvar[k] * prog_en.QTvar[k]), 1.0), -1.0) + sd2_hq = log(corr * sqrt(prog_en.Hvar[k] * prog_en.QTvar[k]) / en.H.values[k] / en.QT.values[k] + 1.0) sd_cond_h_q = sqrt(max(sd_h * sd_h - sd2_hq * sd2_hq / sd_q / sd_q, 0.0)) - mu_q = log( - en.QT.values[k] * en.QT.values[k] / sqrt(en.QT.values[k] * en.QT.values[k] + en.QTvar.values[k]), - ) - mu_h = log(en.H.values[k] * en.H.values[k] / sqrt(en.H.values[k] * en.H.values[k] + en.Hvar.values[k])) + mu_q = + log(en.QT.values[k] * en.QT.values[k] / sqrt(en.QT.values[k] * en.QT.values[k] + prog_en.QTvar[k])) + mu_h = log(en.H.values[k] * en.H.values[k] / sqrt(en.H.values[k] * en.H.values[k] + prog_en.Hvar[k])) else - sd_q = sqrt(en.QTvar.values[k]) - sd_h = sqrt(en.Hvar.values[k]) - corr = max(min(en.HQTcov.values[k] / max(sd_h * sd_q, 1e-13), 1.0), -1.0) + sd_q = sqrt(prog_en.QTvar[k]) + sd_h = sqrt(prog_en.Hvar[k]) + corr = max(min(prog_en.HQTcov[k] / max(sd_h * sd_q, 1e-13), 1.0), -1.0) # limit sd_q to prevent negative qt_hat sd_q_lim = (1e-10 - en.QT.values[k]) / (sqrt2 * abscissas[1]) diff --git a/src/Operators.jl b/src/Operators.jl index c1a6c7e8c..7244e2286 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -620,7 +620,6 @@ function construct_tridiag_diffusion_en( state, TS, w_en, - tke_en, n_updrafts::Int, minimum_area::Float64, pressure_plume_spacing::Vector, @@ -643,6 +642,7 @@ function construct_tridiag_diffusion_en( aux_tc = center_aux_tc(state) prog_up = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) + prog_en = center_prog_environment(state) ae = 1 .- aux_tc.bulk.area rho_ae_K_m = face_field(grid) @@ -690,7 +690,7 @@ function construct_tridiag_diffusion_en( rho_ae_K_m[k + 1] * Δzi * Δzi + rho_ae_K_m[k] * Δzi * Δzi + D_env + - ρ0_c[k] * ae[k] * c_d * sqrt(max(tke_en[k], 0)) / max(mixing_length[k], 1) + ρ0_c[k] * ae[k] * c_d * sqrt(max(prog_en.tke[k], 0)) / max(mixing_length[k], 1) ) if is_toa_center(grid, k) c[k] = 0.0 diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 59b8c73d4..b3f002b25 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -52,32 +52,6 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats) add_profile(Stats, "ed_length_scheme") add_profile(Stats, "mixing_length_ratio") add_profile(Stats, "entdet_balance_length") - add_profile(Stats, "interdomain_tke_t") - add_profile(Stats, "tke_buoy") - add_profile(Stats, "tke_dissipation") - add_profile(Stats, "tke_entr_gain") - add_profile(Stats, "tke_detr_loss") - add_profile(Stats, "tke_shear") - add_profile(Stats, "tke_pressure") - add_profile(Stats, "tke_interdomain") - add_profile(Stats, "Hvar_dissipation") - add_profile(Stats, "QTvar_dissipation") - add_profile(Stats, "HQTcov_dissipation") - add_profile(Stats, "Hvar_entr_gain") - add_profile(Stats, "QTvar_entr_gain") - add_profile(Stats, "Hvar_detr_loss") - add_profile(Stats, "QTvar_detr_loss") - add_profile(Stats, "HQTcov_detr_loss") - add_profile(Stats, "HQTcov_entr_gain") - add_profile(Stats, "Hvar_shear") - add_profile(Stats, "QTvar_shear") - add_profile(Stats, "HQTcov_shear") - add_profile(Stats, "Hvar_rain") - add_profile(Stats, "QTvar_rain") - add_profile(Stats, "HQTcov_rain") - add_profile(Stats, "Hvar_interdomain") - add_profile(Stats, "QTvar_interdomain") - add_profile(Stats, "HQTcov_interdomain") return end @@ -172,41 +146,16 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti write_profile(Stats, "ed_length_scheme", edmf.mls) write_profile(Stats, "mixing_length_ratio", edmf.ml_ratio) write_profile(Stats, "entdet_balance_length", edmf.l_entdet) - write_profile(Stats, "interdomain_tke_t", edmf.b) - compute_covariance_dissipation(edmf, grid, state, edmf.EnvVar.TKE, param_set) - write_profile(Stats, "tke_dissipation", edmf.EnvVar.TKE.dissipation) - write_profile(Stats, "tke_entr_gain", edmf.EnvVar.TKE.entr_gain) - compute_covariance_detr(edmf, grid, state, edmf.EnvVar.TKE) - write_profile(Stats, "tke_detr_loss", edmf.EnvVar.TKE.detr_loss) - write_profile(Stats, "tke_shear", edmf.EnvVar.TKE.shear) - write_profile(Stats, "tke_buoy", edmf.EnvVar.TKE.buoy) - write_profile(Stats, "tke_pressure", edmf.EnvVar.TKE.press) - write_profile(Stats, "tke_interdomain", edmf.EnvVar.TKE.interdomain) - - compute_covariance_dissipation(edmf, grid, state, edmf.EnvVar.Hvar, param_set) - write_profile(Stats, "Hvar_dissipation", edmf.EnvVar.Hvar.dissipation) - compute_covariance_dissipation(edmf, grid, state, edmf.EnvVar.QTvar, param_set) - write_profile(Stats, "QTvar_dissipation", edmf.EnvVar.QTvar.dissipation) - compute_covariance_dissipation(edmf, grid, state, edmf.EnvVar.HQTcov, param_set) - write_profile(Stats, "HQTcov_dissipation", edmf.EnvVar.HQTcov.dissipation) - write_profile(Stats, "Hvar_entr_gain", edmf.EnvVar.Hvar.entr_gain) - write_profile(Stats, "QTvar_entr_gain", edmf.EnvVar.QTvar.entr_gain) - write_profile(Stats, "HQTcov_entr_gain", edmf.EnvVar.HQTcov.entr_gain) - compute_covariance_detr(edmf, grid, state, edmf.EnvVar.Hvar) - compute_covariance_detr(edmf, grid, state, edmf.EnvVar.QTvar) - compute_covariance_detr(edmf, grid, state, edmf.EnvVar.HQTcov) - write_profile(Stats, "Hvar_detr_loss", edmf.EnvVar.Hvar.detr_loss) - write_profile(Stats, "QTvar_detr_loss", edmf.EnvVar.QTvar.detr_loss) - write_profile(Stats, "HQTcov_detr_loss", edmf.EnvVar.HQTcov.detr_loss) - write_profile(Stats, "Hvar_shear", edmf.EnvVar.Hvar.shear) - write_profile(Stats, "QTvar_shear", edmf.EnvVar.QTvar.shear) - write_profile(Stats, "HQTcov_shear", edmf.EnvVar.HQTcov.shear) - write_profile(Stats, "Hvar_rain", edmf.EnvVar.Hvar.rain_src) - write_profile(Stats, "QTvar_rain", edmf.EnvVar.QTvar.rain_src) - write_profile(Stats, "HQTcov_rain", edmf.EnvVar.HQTcov.rain_src) - write_profile(Stats, "Hvar_interdomain", edmf.EnvVar.Hvar.interdomain) - write_profile(Stats, "QTvar_interdomain", edmf.EnvVar.QTvar.interdomain) - write_profile(Stats, "HQTcov_interdomain", edmf.EnvVar.HQTcov.interdomain) + + compute_covariance_dissipation(edmf, grid, state, :tke, param_set) + compute_covariance_detr(edmf, grid, state, :tke) + + compute_covariance_dissipation(edmf, grid, state, :Hvar, param_set) + compute_covariance_dissipation(edmf, grid, state, :QTvar, param_set) + compute_covariance_dissipation(edmf, grid, state, :HQTcov, param_set) + compute_covariance_detr(edmf, grid, state, :Hvar) + compute_covariance_detr(edmf, grid, state, :QTvar) + compute_covariance_detr(edmf, grid, state, :HQTcov) return end @@ -455,6 +404,7 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca n_updrafts = up.n_updrafts prog_gm = center_prog_grid_mean(state) tendencies_gm = center_tendencies_grid_mean(state) + prog_en = center_prog_environment(state) # 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 @@ -483,7 +433,6 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca state, TS, en.W.values, - en.TKE.values, up.n_updrafts, edmf.minimum_area, edmf.pressure_plume_spacing, @@ -497,10 +446,10 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca implicit_eqs.A_QTvar .= construct_tridiag_diffusion_en(common_args..., false) implicit_eqs.A_HQTcov .= construct_tridiag_diffusion_en(common_args..., false) - implicit_eqs.b_TKE .= en_diffusion_tendencies(grid, state, TS, en.TKE, n_updrafts) - implicit_eqs.b_Hvar .= en_diffusion_tendencies(grid, state, TS, en.Hvar, n_updrafts) - implicit_eqs.b_QTvar .= en_diffusion_tendencies(grid, state, TS, en.QTvar, n_updrafts) - implicit_eqs.b_HQTcov .= en_diffusion_tendencies(grid, state, TS, en.HQTcov, n_updrafts) + implicit_eqs.b_TKE .= en_diffusion_tendencies(grid, state, TS, :tke, n_updrafts) + implicit_eqs.b_Hvar .= en_diffusion_tendencies(grid, state, TS, :Hvar, n_updrafts) + implicit_eqs.b_QTvar .= en_diffusion_tendencies(grid, state, TS, :QTvar, n_updrafts) + implicit_eqs.b_HQTcov .= en_diffusion_tendencies(grid, state, TS, :HQTcov, n_updrafts) # ----------- ### @@ -511,10 +460,10 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca update_rain(edmf.Rain, grid, state, up_thermo, en_thermo, edmf.RainPhys, TS) end - en.TKE.values .= tridiag_solve(implicit_eqs.b_TKE, implicit_eqs.A_TKE) - en.Hvar.values .= tridiag_solve(implicit_eqs.b_Hvar, implicit_eqs.A_Hvar) - en.QTvar.values .= tridiag_solve(implicit_eqs.b_QTvar, implicit_eqs.A_QTvar) - en.HQTcov.values .= tridiag_solve(implicit_eqs.b_HQTcov, implicit_eqs.A_HQTcov) + parent(prog_en.tke) .= tridiag_solve(implicit_eqs.b_TKE, implicit_eqs.A_TKE) + parent(prog_en.Hvar) .= tridiag_solve(implicit_eqs.b_Hvar, implicit_eqs.A_Hvar) + parent(prog_en.QTvar) .= tridiag_solve(implicit_eqs.b_QTvar, implicit_eqs.A_QTvar) + parent(prog_en.HQTcov) .= tridiag_solve(implicit_eqs.b_HQTcov, implicit_eqs.A_HQTcov) @inbounds for k in real_center_indices(grid) prog_gm.u[k] += tendencies_gm.u[k] * TS.dt prog_gm.v[k] += tendencies_gm.v[k] * TS.dt @@ -526,11 +475,11 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca ### set values ### @inbounds for k in real_center_indices(grid) - en.TKE.values[k] = max(en.TKE.values[k], 0.0) - en.Hvar.values[k] = max(en.Hvar.values[k], 0.0) - en.QTvar.values[k] = max(en.QTvar.values[k], 0.0) - en.HQTcov.values[k] = max(en.HQTcov.values[k], -sqrt(en.Hvar.values[k] * en.QTvar.values[k])) - en.HQTcov.values[k] = min(en.HQTcov.values[k], sqrt(en.Hvar.values[k] * en.QTvar.values[k])) + prog_en.tke[k] = max(prog_en.tke[k], 0.0) + prog_en.Hvar[k] = max(prog_en.Hvar[k], 0.0) + prog_en.QTvar[k] = max(prog_en.QTvar[k], 0.0) + prog_en.HQTcov[k] = max(prog_en.HQTcov[k], -sqrt(prog_en.Hvar[k] * prog_en.QTvar[k])) + prog_en.HQTcov[k] = min(prog_en.HQTcov[k], sqrt(prog_en.Hvar[k] * prog_en.QTvar[k])) end # set values @@ -584,13 +533,14 @@ function reset_surface_covariance(edmf::EDMF_PrognosticTKE, grid, state, gm, Cas gm = gm up = edmf.UpdVar en = edmf.EnvVar + prog_en = center_prog_environment(state) - en.TKE.values[kc_surf] = get_surface_tke(Case.Sur.ustar, grid.zc[kc_surf], Case.Sur.obukhov_length) - get_GMV_CoVar(edmf, grid, state, en.W, en.W, en.TKE, :tke, :w) + prog_en.tke[kc_surf] = get_surface_tke(Case.Sur.ustar, grid.zc[kc_surf], Case.Sur.obukhov_length) + get_GMV_CoVar(edmf, grid, state, en.W, en.W, :tke, :w) - en.Hvar.values[kc_surf] = get_surface_variance(flux1 * α0LL, flux1 * α0LL, ustar, zLL, oblength) - en.QTvar.values[kc_surf] = get_surface_variance(flux2 * α0LL, flux2 * α0LL, ustar, zLL, oblength) - en.HQTcov.values[kc_surf] = get_surface_variance(flux1 * α0LL, flux2 * α0LL, ustar, zLL, oblength) + prog_en.Hvar[kc_surf] = get_surface_variance(flux1 * α0LL, flux1 * α0LL, ustar, zLL, oblength) + prog_en.QTvar[kc_surf] = get_surface_variance(flux2 * α0LL, flux2 * α0LL, ustar, zLL, oblength) + prog_en.HQTcov[kc_surf] = get_surface_variance(flux1 * α0LL, flux2 * α0LL, ustar, zLL, oblength) return end @@ -602,21 +552,21 @@ function get_GMV_CoVar( state, ϕ_en::EnvironmentVariable, ψ_en::EnvironmentVariable, - covar_e::EnvironmentVariable_2m, - gm_covar_sym::Symbol, + covar_sym::Symbol, ϕ_sym::Symbol, ψ_sym::Symbol = ϕ_sym, ) aux_tc = center_aux_tc(state) ae = 1 .- aux_tc.bulk.area - is_tke = covar_e.name == "tke" + is_tke = covar_sym == :tke tke_factor = is_tke ? 0.5 : 1 prog_gm_c = center_prog_grid_mean(state) prog_gm_f = face_prog_grid_mean(state) prog_up = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) - gmv_covar = getproperty(center_aux_grid_mean(state), gm_covar_sym) + gmv_covar = getproperty(center_aux_grid_mean(state), covar_sym) + covar_e = getproperty(center_prog_environment(state), covar_sym) prog_gm = is_tke ? prog_gm_f : prog_gm_c ϕ_gm = getproperty(prog_gm, ϕ_sym) ψ_gm = getproperty(prog_gm, ψ_sym) @@ -632,7 +582,7 @@ function get_GMV_CoVar( Δϕ = interpf2c(Δϕ_dual, grid, k) Δψ = interpf2c(Δψ_dual, grid, k) - gmv_covar[k] = tke_factor * ae[k] * Δϕ * Δψ + ae[k] * covar_e.values[k] + gmv_covar[k] = tke_factor * ae[k] * Δϕ * Δψ + ae[k] * covar_e[k] @inbounds for i in 1:(edmf.n_updrafts) ϕ_up_var = getproperty(prog_up_f[i], ϕ_sym) ψ_up_var = getproperty(prog_up_f[i], ψ_sym) @@ -653,7 +603,7 @@ function get_GMV_CoVar( Δϕ = ϕ_en.values[k] - ϕ_gm[k] Δψ = ψ_en.values[k] - ψ_gm[k] - gmv_covar[k] = tke_factor * ae[k] * Δϕ * Δψ + ae[k] * covar_e.values[k] + gmv_covar[k] = tke_factor * ae[k] * Δϕ * Δψ + ae[k] * covar_e[k] @inbounds for i in 1:(edmf.n_updrafts) ϕ_up_var = getproperty(prog_up[i], ϕ_sym) ψ_up_var = getproperty(prog_up[i], ψ_sym) @@ -876,17 +826,18 @@ function initialize_covariance(edmf::EDMF_PrognosticTKE, grid, state, gm, Case:: kc_surf = kc_surface(grid) en = edmf.EnvVar aux_gm = center_aux_grid_mean(state) + prog_en = center_prog_environment(state) - en.TKE.values .= vec(aux_gm.tke) + prog_en.tke .= aux_gm.tke reset_surface_covariance(edmf, grid, state, gm, Case) aux_gm.Hvar .= aux_gm.Hvar[kc_surf] .* aux_gm.tke aux_gm.QTvar .= aux_gm.QTvar[kc_surf] .* aux_gm.tke aux_gm.HQTcov .= aux_gm.HQTcov[kc_surf] .* aux_gm.tke - en.Hvar.values .= vec(aux_gm.Hvar) - en.QTvar.values .= vec(aux_gm.QTvar) - en.HQTcov.values .= vec(aux_gm.HQTcov) + prog_en.Hvar .= aux_gm.Hvar + prog_en.QTvar .= aux_gm.QTvar + prog_en.HQTcov .= aux_gm.HQTcov return end @@ -895,7 +846,7 @@ function compute_covariance_shear( grid, state, gm::GridMeanVariables, - Covar::EnvironmentVariable_2m, + covar_sym::Symbol, EnvVar1, EnvVar2, ) @@ -905,9 +856,11 @@ function compute_covariance_shear( aux_tc = center_aux_tc(state) ρ0_c = center_ref_state(state).ρ0 prog_gm = center_prog_grid_mean(state) - is_tke = Covar.name == "tke" + is_tke = covar_sym == :tke tke_factor = is_tke ? 0.5 : 1 k_eddy = is_tke ? aux_tc.KM : aux_tc.KH + aux_en_2m = center_aux_environment_2m(state) + aux_covar = getproperty(aux_en_2m, covar_sym) if is_tke @inbounds for k in real_center_indices(grid) @@ -923,7 +876,7 @@ function compute_covariance_shear( ∇var2 = ∇f2c(var2_dual, grid, k; bottom = SetValue(0), top = SetGradient(0)) ∇var1 = ∇f2c(var1_dual, grid, k; bottom = SetValue(0), top = SetGradient(0)) - Covar.shear[k] = tke_factor * 2 * (ρ0_c[k] * ae[k] * k_eddy[k] * (∇var1 * ∇var2 + ∇u^2 + ∇v^2)) + aux_covar.shear[k] = tke_factor * 2 * (ρ0_c[k] * ae[k] * k_eddy[k] * (∇var1 * ∇var2 + ∇u^2 + ∇v^2)) end else @inbounds for k in real_center_indices(grid) @@ -934,7 +887,7 @@ function compute_covariance_shear( var2_cut = ccut(EnvVar2, grid, k) ∇var2 = c∇(var2_cut, grid, k; bottom = Extrapolate(), top = SetGradient(0)) - Covar.shear[k] = tke_factor * 2 * (ρ0_c[k] * ae[k] * k_eddy[k] * (∇var1 * ∇var2)) + aux_covar.shear[k] = tke_factor * 2 * (ρ0_c[k] * ae[k] * k_eddy[k] * (∇var1 * ∇var2)) end end return @@ -946,37 +899,39 @@ function compute_covariance_interdomain_src( state, ϕ_en::EnvironmentVariable, ψ_en::EnvironmentVariable, - Covar::EnvironmentVariable_2m, + Covar_sym::Symbol, ϕ_var::Symbol, ψ_var::Symbol = ϕ_var, ) - is_tke = Covar.name == "tke" + is_tke = Covar_sym == :tke tke_factor = is_tke ? 0.5 : 1 prog_up_c = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) + aux_en_2m = center_aux_environment_2m(state) + aux_covar = getproperty(aux_en_2m, Covar_sym) prog_up = is_tke ? prog_up_f : prog_up_c if is_tke @inbounds for k in real_center_indices(grid) - Covar.interdomain[k] = 0.0 + aux_covar.interdomain[k] = 0.0 @inbounds for i in 1:(edmf.n_updrafts) ϕ_up = getproperty(prog_up[i], ϕ_var) ψ_up = getproperty(prog_up[i], ψ_var) Δϕ = interpf2c(ϕ_up, grid, k) - interpf2c(ϕ_en.values, grid, k) Δψ = interpf2c(ψ_up, grid, k) - interpf2c(ψ_en.values, grid, k) - Covar.interdomain[k] += tke_factor * prog_up_c[i].area[k] * (1.0 - prog_up_c[i].area[k]) * Δϕ * Δψ + aux_covar.interdomain[k] += tke_factor * prog_up_c[i].area[k] * (1.0 - prog_up_c[i].area[k]) * Δϕ * Δψ end end else @inbounds for k in real_center_indices(grid) - Covar.interdomain[k] = 0.0 + aux_covar.interdomain[k] = 0.0 @inbounds for i in 1:(edmf.n_updrafts) ϕ_up = getproperty(prog_up[i], ϕ_var) ψ_up = getproperty(prog_up[i], ψ_var) Δϕ = ϕ_up[k] - ϕ_en.values[k] Δψ = ψ_up[k] - ψ_en.values[k] - Covar.interdomain[k] += tke_factor * prog_up_c[i].area[k] * (1.0 - prog_up_c[i].area[k]) * Δϕ * Δψ + aux_covar.interdomain[k] += tke_factor * prog_up_c[i].area[k] * (1.0 - prog_up_c[i].area[k]) * Δϕ * Δψ end end end @@ -987,7 +942,7 @@ function compute_covariance_entr( edmf::EDMF_PrognosticTKE, grid, state, - Covar::EnvironmentVariable_2m, + covar_sym::Symbol, EnvVar1::EnvironmentVariable, EnvVar2::EnvironmentVariable, var1::Symbol, @@ -996,7 +951,7 @@ function compute_covariance_entr( ρ_0_c = center_ref_state(state).ρ0 - is_tke = Covar.name == "tke" + is_tke = covar_sym == :tke tke_factor = is_tke ? 0.5 : 1 prog_up_c = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) @@ -1006,10 +961,14 @@ function compute_covariance_entr( prog_up = is_tke ? prog_up_f : prog_up_c GmvVar1 = getproperty(prog_gm, var1) GmvVar2 = getproperty(prog_gm, var2) + aux_en_2m = center_aux_environment_2m(state) + aux_covar = getproperty(aux_en_2m, covar_sym) + prog_en = center_prog_environment(state) + prog_covar = getproperty(prog_en, covar_sym) @inbounds for k in real_center_indices(grid) - Covar.entr_gain[k] = 0.0 - Covar.detr_loss[k] = 0.0 + aux_covar.entr_gain[k] = 0.0 + aux_covar.detr_loss[k] = 0.0 @inbounds for i in 1:(edmf.n_updrafts) a_up = prog_up_c[i].area[k] if a_up > edmf.minimum_area @@ -1041,9 +1000,9 @@ function compute_covariance_entr( abs(w_u) * eps_turb * ((envvar1 - gmvvar1) * (updvar2 - envvar2) + (envvar2 - gmvvar2) * (updvar1 - envvar1)) - Covar.entr_gain[k] += dynamic_entr + turbulent_entr - Covar.detr_loss[k] += - tke_factor * ρ_0_c[k] * a_up * abs(w_u) * (edmf.entr_sc[i, k] + eps_turb) * Covar.values[k] + aux_covar.entr_gain[k] += dynamic_entr + turbulent_entr + aux_covar.detr_loss[k] += + tke_factor * ρ_0_c[k] * a_up * abs(w_u) * (edmf.entr_sc[i, k] + eps_turb) * prog_covar[k] end end end @@ -1051,42 +1010,56 @@ function compute_covariance_entr( return end -function compute_covariance_detr(edmf::EDMF_PrognosticTKE, grid, state, Covar::EnvironmentVariable_2m) +function compute_covariance_detr(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol) up = edmf.UpdVar ρ0_c = center_ref_state(state).ρ0 prog_up = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) + + aux_en_2m = center_aux_environment_2m(state) + aux_covar = getproperty(aux_en_2m, covar_sym) + prog_en = center_prog_environment(state) + prog_covar = getproperty(prog_en, covar_sym) + @inbounds for k in real_center_indices(grid) - Covar.detr_loss[k] = 0.0 + aux_covar.detr_loss[k] = 0.0 @inbounds for i in 1:(up.n_updrafts) w_up_c = interpf2c(prog_up_f[i].w, grid, k) - Covar.detr_loss[k] += prog_up[i].area[k] * abs(w_up_c) * edmf.entr_sc[i, k] + aux_covar.detr_loss[k] += prog_up[i].area[k] * abs(w_up_c) * edmf.entr_sc[i, k] end - Covar.detr_loss[k] *= ρ0_c[k] * Covar.values[k] + aux_covar.detr_loss[k] *= ρ0_c[k] * prog_covar[k] end return end -function compute_covariance_dissipation(edmf::EDMF_PrognosticTKE, grid, state, Covar::EnvironmentVariable_2m, param_set) +function compute_covariance_dissipation(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, param_set) en = edmf.EnvVar c_d = CPEDMF.c_d(param_set) aux_tc = center_aux_tc(state) ae = 1 .- aux_tc.bulk.area ρ0_c = center_ref_state(state).ρ0 + prog_en = center_prog_environment(state) + aux_en_2m = center_aux_environment_2m(state) + aux_covar = getproperty(aux_en_2m, covar_sym) + prog_covar = getproperty(prog_en, covar_sym) + prog_en = center_prog_environment(state) @inbounds for k in real_center_indices(grid) - Covar.dissipation[k] = ( - ρ0_c[k] * ae[k] * Covar.values[k] * max(en.TKE.values[k], 0)^0.5 / max(edmf.mixing_length[k], 1.0e-3) * c_d - ) + aux_covar.dissipation[k] = + (ρ0_c[k] * ae[k] * prog_covar[k] * max(prog_en.tke[k], 0)^0.5 / max(edmf.mixing_length[k], 1.0e-3) * c_d) end return end -function en_diffusion_tendencies(grid::Grid, state, TS, covar, n_updrafts) +function en_diffusion_tendencies(grid::Grid, state, TS, covar_sym::Symbol, n_updrafts) dti = TS.dti b = center_field(grid) ρ0_c = center_ref_state(state).ρ0 prog_up = center_prog_updrafts(state) + prog_en = center_prog_environment(state) + aux_en_2m = center_aux_environment_2m(state) + prog_covar = getproperty(prog_en, covar_sym) + aux_covar = getproperty(aux_en_2m, covar_sym) ae = center_field(grid) @@ -1095,19 +1068,19 @@ function en_diffusion_tendencies(grid::Grid, state, TS, covar, n_updrafts) end kc_surf = kc_surface(grid) - covar_surf = covar.values[kc_surf] + covar_surf = prog_covar[kc_surf] @inbounds for k in real_center_indices(grid) if is_surface_center(grid, k) b[k] = covar_surf else b[k] = ( - ρ0_c[k] * ae[k] * covar.values[k] * dti + - covar.press[k] + - covar.buoy[k] + - covar.shear[k] + - covar.entr_gain[k] + - covar.rain_src[k] + ρ0_c[k] * ae[k] * prog_covar[k] * dti + + aux_covar.press[k] + + aux_covar.buoy[k] + + aux_covar.shear[k] + + aux_covar.entr_gain[k] + + aux_covar.rain_src[k] ) end end @@ -1119,7 +1092,7 @@ function GMV_third_m( edmf::EDMF_PrognosticTKE, grid, state, - env_covar::EnvironmentVariable_2m, + covar_en_sym::Symbol, env_mean::EnvironmentVariable, up_var::Symbol, gm_third_m_sym::Symbol, @@ -1133,7 +1106,8 @@ function GMV_third_m( ae = 1 .- aux_tc.bulk.area prog_up = center_prog_updrafts(state) prog_up_f = face_prog_updrafts(state) - is_tke = env_covar.name == "tke" + is_tke = covar_en_sym == :tke + covar_en = getproperty(center_prog_environment(state), covar_en_sym) @inbounds for k in real_center_indices(grid) mean_en = is_tke ? interpf2c(env_mean.values, grid, k) : env_mean.values[k] @@ -1154,7 +1128,7 @@ function GMV_third_m( ∇w_en = ∇f2c(w_en_dual, grid, k; w_bcs...) Envcov_ = -edmf.horiz_K_eddy[i_last, k] * ∇w_en else - Envcov_ = env_covar.values[k] + Envcov_ = covar_en[k] end Upd_cubed = 0.0 diff --git a/src/Variables.jl b/src/Variables.jl index 591e113b1..fc3d26cb9 100644 --- a/src/Variables.jl +++ b/src/Variables.jl @@ -24,9 +24,9 @@ function satadjust(gm::GridMeanVariables, grid, state) prog_gm = center_prog_grid_mean(state) param_set = parameter_set(gm) @inbounds for k in real_center_indices(grid) - h = prog_gm.θ_liq_ice[k] - qt = prog_gm.q_tot[k] - ts = TD.PhaseEquil_pθq(param_set, p0_c[k], h, qt) + θ_liq_ice = prog_gm.θ_liq_ice[k] + q_tot = prog_gm.q_tot[k] + ts = TD.PhaseEquil_pθq(param_set, p0_c[k], θ_liq_ice, q_tot) aux_gm.q_liq[k] = TD.liquid_specific_humidity(ts) aux_gm.T[k] = TD.air_temperature(ts) ρ = TD.air_density(ts) diff --git a/src/diagnostics.jl b/src/diagnostics.jl index d3b89904e..b27e57c11 100644 --- a/src/diagnostics.jl +++ b/src/diagnostics.jl @@ -31,6 +31,7 @@ function io_dictionary_ref_state(state) end #! format: off +# TODO: We probably don't need to split the aux/prog dictionaries. Only static vs dynamic. function io_dictionary_aux(state) DT = NamedTuple{(:dims, :group, :field), Tuple{Tuple{String, String}, String, Any}} io_dict = Dict{String, DT}( @@ -61,6 +62,41 @@ function io_dictionary_aux(state) "thetal_mean" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_grid_mean(state).θ_liq_ice), "eddy_viscosity" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_tc(state).KM), "eddy_diffusivity" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_tc(state).KH), + "env_tke" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_environment(state).tke), + "env_Hvar" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_environment(state).Hvar), + "env_QTvar" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_environment(state).QTvar), + + "env_HQTcov" => (; dims = ("zc", "t"), group = "profiles", field = center_prog_environment(state).HQTcov), + + "tke_buoy" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.buoy), + "tke_pressure" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.press), + "tke_dissipation" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.dissipation), + "tke_entr_gain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.entr_gain), + "tke_detr_loss" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.detr_loss), + "tke_shear" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.shear), + "tke_interdomain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).tke.interdomain), + + "Hvar_dissipation" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).Hvar.dissipation), + "Hvar_entr_gain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).Hvar.entr_gain), + "Hvar_detr_loss" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).Hvar.detr_loss), + "Hvar_interdomain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).Hvar.interdomain), + "Hvar_rain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).Hvar.rain_src), + "Hvar_shear" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).Hvar.shear), + + "QTvar_dissipation" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).QTvar.dissipation), + "QTvar_entr_gain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).QTvar.entr_gain), + "QTvar_detr_loss" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).QTvar.detr_loss), + "QTvar_shear" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).QTvar.shear), + "QTvar_rain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).QTvar.rain_src), + "QTvar_interdomain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).QTvar.interdomain), + + "HQTcov_rain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).HQTcov.rain_src), + "HQTcov_dissipation" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).HQTcov.dissipation), + "HQTcov_entr_gain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).HQTcov.entr_gain), + "HQTcov_detr_loss" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).HQTcov.detr_loss), + "HQTcov_shear" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).HQTcov.shear), + "HQTcov_interdomain" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment_2m(state).HQTcov.interdomain), + ) return io_dict end diff --git a/src/dycore_api.jl b/src/dycore_api.jl index c2ff2a5e6..82cfc382f 100644 --- a/src/dycore_api.jl +++ b/src/dycore_api.jl @@ -55,6 +55,7 @@ face_aux_tc(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_environment_2m(state) = aux_turbconv(state, CentField()).en_2m face_aux_environment(state) = aux_turbconv(state, FaceField()).en #= Tendency fields for TurbulenceConvection =# diff --git a/src/types.jl b/src/types.jl index 9ae938430..73e2eb15f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -373,44 +373,6 @@ struct EnvironmentVariable{T} end end -struct EnvironmentVariable_2m{A1} - values::A1 - dissipation::A1 - shear::A1 - entr_gain::A1 - detr_loss::A1 - press::A1 - buoy::A1 - interdomain::A1 - rain_src::A1 - name::String - units::String - function EnvironmentVariable_2m(grid, name, units) - values = center_field(grid) - dissipation = center_field(grid) - entr_gain = center_field(grid) - detr_loss = center_field(grid) - buoy = center_field(grid) - press = center_field(grid) - shear = center_field(grid) - interdomain = center_field(grid) - rain_src = center_field(grid) - return new{typeof(values)}( - values, - dissipation, - shear, - entr_gain, - detr_loss, - press, - buoy, - interdomain, - rain_src, - name, - units, - ) - end -end - Base.@kwdef mutable struct EnvironmentVariables W::EnvironmentVariable Area::EnvironmentVariable @@ -421,10 +383,6 @@ Base.@kwdef mutable struct EnvironmentVariables T::EnvironmentVariable B::EnvironmentVariable cloud_fraction::EnvironmentVariable - TKE::EnvironmentVariable_2m - Hvar::EnvironmentVariable_2m - QTvar::EnvironmentVariable_2m - HQTcov::EnvironmentVariable_2m cloud_base::Float64 = 0 cloud_top::Float64 = 0 cloud_cover::Float64 = 0 @@ -444,27 +402,8 @@ function EnvironmentVariables(namelist, grid::Grid) # TODO - the flag setting is repeated from Variables.pyx logic EnvThermo_scheme = parse_namelist(namelist, "thermodynamics", "sgs"; default = "mean") - TKE = EnvironmentVariable_2m(grid, "tke", "m^2/s^2") - QTvar = EnvironmentVariable_2m(grid, "qt_var", "kg^2/kg^2") - Hvar = EnvironmentVariable_2m(grid, "thetal_var", "K^2") - HQTcov = EnvironmentVariable_2m(grid, "thetal_qt_covar", "K(kg/kg)") - - return EnvironmentVariables(; - W, - Area, - QT, - QL, - H, - RH, - T, - B, - cloud_fraction, - TKE, - Hvar, - QTvar, - HQTcov, - EnvThermo_scheme, - ) + + return EnvironmentVariables(; W, Area, QT, QL, H, RH, T, B, cloud_fraction, EnvThermo_scheme) end struct EnvironmentThermodynamics{A1} @@ -701,7 +640,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} mls::A1 ml_ratio::A1 l_entdet::A1 - b::A1 wstar::Float64 entr_surface_bc::Float64 detr_surface_bc::Float64 @@ -868,7 +806,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} mls = center_field(grid) ml_ratio = center_field(grid) l_entdet = center_field(grid) - b = center_field(grid) wstar = 0 entr_surface_bc = 0 detr_surface_bc = 0 @@ -932,7 +869,6 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} mls, ml_ratio, l_entdet, - b, wstar, entr_surface_bc, detr_surface_bc, diff --git a/src/update_aux.jl b/src/update_aux.jl index 4a8d5c6f7..c66e3e451 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -26,6 +26,8 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) prog_up_f = face_prog_updrafts(state) aux_up_f = face_aux_tc(state) aux_tc = center_aux_tc(state) + prog_en = center_prog_environment(state) + aux_en_2m = center_aux_environment_2m(state) for k in real_center_indices(grid) aux_tc.bulk.area[k] = sum(ntuple(i -> prog_up[i].area[k], up.n_updrafts)) @@ -35,12 +37,12 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ##### diagnose_GMV_moments ##### #! format: off - get_GMV_CoVar(edmf, grid, state, en.H, en.H, en.Hvar, :Hvar, :θ_liq_ice) - get_GMV_CoVar(edmf, grid, state, en.QT, en.QT, en.QTvar, :QTvar, :q_tot) - get_GMV_CoVar(edmf, grid, state, en.H, en.QT, en.HQTcov, :HQTcov, :θ_liq_ice, :q_tot) - GMV_third_m(edmf, grid, state, en.Hvar, en.H, :θ_liq_ice, :H_third_m) - GMV_third_m(edmf, grid, state, en.QTvar, en.QT, :q_tot, :QT_third_m) - GMV_third_m(edmf, grid, state, en.TKE, en.W, :w, :W_third_m) + get_GMV_CoVar(edmf, grid, state, en.H, en.H, :Hvar, :θ_liq_ice) + get_GMV_CoVar(edmf, grid, state, en.QT, en.QT, :QTvar, :q_tot) + get_GMV_CoVar(edmf, grid, state, en.H, en.QT, :HQTcov, :θ_liq_ice, :q_tot) + GMV_third_m(edmf, grid, state, :Hvar, en.H, :θ_liq_ice, :H_third_m) + GMV_third_m(edmf, grid, state, :QTvar, en.QT, :q_tot, :QT_third_m) + GMV_third_m(edmf, grid, state, :tke, en.W, :w, :W_third_m) #! format: on ##### @@ -205,7 +207,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) w_en = interpf2c(en.W.values, grid, k), b_up = aux_up[i].buoy[k], b_en = en.B.values[k], - tke = en.TKE.values[k], + tke = prog_en.tke[k], dMdz = ∇m, M = m, a_up = prog_up[i].area[k], @@ -330,7 +332,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ml_model = MinDisspLen(; z = grid.zc[k].z, obukhov_length = obukhov_length, - tke_surf = en.TKE.values[kc_surf], + tke_surf = prog_en.tke[kc_surf], ustar = surface.ustar, Pr = edmf.prandtl_nvec[k], p0 = p0_c[k], @@ -341,7 +343,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ∂qt∂z = ∂qt∂z, ∂θl∂z = ∂θl∂z, θv = θv, - tke = en.TKE.values[k], + tke = prog_en.tke[k], a_en = (1 - aux_tc.bulk.area[k]), wc_en = wc_en, wc_up = Tuple(wc_up), @@ -361,38 +363,38 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) edmf.mixing_length[k] = ml.mixing_length edmf.ml_ratio[k] = ml.ml_ratio - KM[k] = c_m * edmf.mixing_length[k] * sqrt(max(en.TKE.values[k], 0.0)) + KM[k] = c_m * edmf.mixing_length[k] * sqrt(max(prog_en.tke[k], 0.0)) KH[k] = KM[k] / edmf.prandtl_nvec[k] - en.TKE.buoy[k] = -ml_model.a_en * ρ0_c[k] * KH[k] * (bg.∂b∂z_θl + bg.∂b∂z_qt) + aux_en_2m.tke.buoy[k] = -ml_model.a_en * ρ0_c[k] * KH[k] * (bg.∂b∂z_θl + bg.∂b∂z_qt) end - compute_covariance_entr(edmf, grid, state, en.TKE, en.W, en.W, :w) - compute_covariance_shear(edmf, grid, state, gm, en.TKE, en.W.values, en.W.values) - compute_covariance_interdomain_src(edmf, grid, state, en.W, en.W, en.TKE, :w) + compute_covariance_entr(edmf, grid, state, :tke, en.W, en.W, :w) + compute_covariance_shear(edmf, grid, state, gm, :tke, en.W.values, en.W.values) + compute_covariance_interdomain_src(edmf, grid, state, en.W, en.W, :tke, :w) ##### ##### compute_tke_pressure ##### @inbounds for k in real_center_indices(grid) - en.TKE.press[k] = 0.0 + aux_en_2m.tke.press[k] = 0.0 @inbounds for i in 1:(up.n_updrafts) w_up_c = interpf2c(prog_up_f[i].w, grid, k) w_en_c = interpf2c(en.W.values, grid, k) press_c = interpf2c(edmf.nh_pressure, grid, k, i) - en.TKE.press[k] += (w_en_c - w_up_c) * press_c + aux_en_2m.tke.press[k] += (w_en_c - w_up_c) * press_c end end - compute_covariance_entr(edmf, grid, state, en.Hvar, en.H, en.H, :θ_liq_ice) - compute_covariance_entr(edmf, grid, state, en.QTvar, en.QT, en.QT, :q_tot) - compute_covariance_entr(edmf, grid, state, en.HQTcov, en.H, en.QT, :θ_liq_ice, :q_tot) - compute_covariance_shear(edmf, grid, state, gm, en.Hvar, en.H.values, en.H.values) - compute_covariance_shear(edmf, grid, state, gm, en.QTvar, en.QT.values, en.QT.values) - compute_covariance_shear(edmf, grid, state, gm, en.HQTcov, en.H.values, en.QT.values) - compute_covariance_interdomain_src(edmf, grid, state, en.H, en.H, en.Hvar, :θ_liq_ice) - compute_covariance_interdomain_src(edmf, grid, state, en.QT, en.QT, en.QTvar, :q_tot) - compute_covariance_interdomain_src(edmf, grid, state, en.H, en.QT, en.HQTcov, :θ_liq_ice, :q_tot) + compute_covariance_entr(edmf, grid, state, :Hvar, en.H, en.H, :θ_liq_ice) + compute_covariance_entr(edmf, grid, state, :QTvar, en.QT, en.QT, :q_tot) + compute_covariance_entr(edmf, grid, state, :HQTcov, en.H, en.QT, :θ_liq_ice, :q_tot) + compute_covariance_shear(edmf, grid, state, gm, :Hvar, en.H.values, en.H.values) + compute_covariance_shear(edmf, grid, state, gm, :QTvar, en.QT.values, en.QT.values) + compute_covariance_shear(edmf, grid, state, gm, :HQTcov, en.H.values, en.QT.values) + compute_covariance_interdomain_src(edmf, grid, state, en.H, en.H, :Hvar, :θ_liq_ice) + compute_covariance_interdomain_src(edmf, grid, state, en.QT, en.QT, :QTvar, :q_tot) + compute_covariance_interdomain_src(edmf, grid, state, en.H, en.QT, :HQTcov, :θ_liq_ice, :q_tot) ##### ##### compute_covariance_rain # need to update this one @@ -401,10 +403,10 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) # TODO defined again in compute_covariance_shear and compute_covaraince ae = 1 .- aux_tc.bulk.area # area of environment @inbounds for k in real_center_indices(grid) - en.TKE.rain_src[k] = 0 - en.Hvar.rain_src[k] = ρ0_c[k] * ae[k] * 2 * en_thermo.Hvar_rain_dt[k] - en.QTvar.rain_src[k] = ρ0_c[k] * ae[k] * 2 * en_thermo.QTvar_rain_dt[k] - en.HQTcov.rain_src[k] = ρ0_c[k] * ae[k] * en_thermo.HQTcov_rain_dt[k] + aux_en_2m.tke.rain_src[k] = 0 + aux_en_2m.Hvar.rain_src[k] = ρ0_c[k] * ae[k] * 2 * en_thermo.Hvar_rain_dt[k] + aux_en_2m.QTvar.rain_src[k] = ρ0_c[k] * ae[k] * 2 * en_thermo.QTvar_rain_dt[k] + aux_en_2m.HQTcov.rain_src[k] = ρ0_c[k] * ae[k] * en_thermo.HQTcov_rain_dt[k] end reset_surface_covariance(edmf, grid, state, gm, Case)