From 09f4778429f1eb672b21473080a929f23004433f Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Sun, 31 Oct 2021 20:00:02 -0700 Subject: [PATCH] Solve for conserv, not prim, prog env vars --- integration_tests/utils/initial_conditions.jl | 13 +- integration_tests/utils/main.jl | 6 +- integration_tests/utils/mse_tables.jl | 228 +++++++++--------- src/EDMF_Environment.jl | 27 +-- src/Turbulence_PrognosticTKE.jl | 73 +++--- src/diagnostics.jl | 9 +- src/update_aux.jl | 17 +- 7 files changed, 195 insertions(+), 178 deletions(-) diff --git a/integration_tests/utils/initial_conditions.jl b/integration_tests/utils/initial_conditions.jl index ccaf85d12..06fb75f3a 100644 --- a/integration_tests/utils/initial_conditions.jl +++ b/integration_tests/utils/initial_conditions.jl @@ -22,17 +22,22 @@ function initialize_covariance(edmf::TC.EDMF_PrognosticTKE, grid, state, gm, Cas en = edmf.EnvVar aux_gm = TC.center_aux_grid_mean(state) prog_en = TC.center_prog_environment(state) + aux_en = TC.center_aux_environment(state) + ρ0_c = TC.center_ref_state(state).ρ0 + aux_tc = TC.center_aux_turbconv(state) + ae = 1 .- aux_tc.bulk.area # area of environment - prog_en.tke .= aux_gm.tke + aux_en.tke .= aux_gm.tke + prog_en.ρatke .= aux_en.tke .* ρ0_c .* ae TC.get_GMV_CoVar(edmf, grid, state, :tke, :w) 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 - prog_en.Hvar .= aux_gm.Hvar - prog_en.QTvar .= aux_gm.QTvar - prog_en.HQTcov .= aux_gm.HQTcov + prog_en.ρaHvar .= aux_gm.Hvar .* ρ0_c .* ae + prog_en.ρaQTvar .= aux_gm.QTvar .* ρ0_c .* ae + prog_en.ρaHQTcov .= aux_gm.HQTcov .* ρ0_c .* ae return end diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl index d651f008c..0e4f14239 100644 --- a/integration_tests/utils/main.jl +++ b/integration_tests/utils/main.jl @@ -107,6 +107,10 @@ cent_aux_vars_edmf(FT, n_up) = (; T = FT(0), buoy = FT(0), cloud_fraction = FT(0), + tke = FT(0), + Hvar = FT(0), + QTvar = FT(0), + HQTcov = FT(0), ), en_2m = (; tke = cent_aux_vars_en_2m(FT), @@ -157,7 +161,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), ρaθ_liq_ice = FT(0), ρaq_tot = FT(0)) -cent_prognostic_vars_en(FT) = (; tke = FT(0), Hvar = FT(0), QTvar = FT(0), HQTcov = FT(0)) +cent_prognostic_vars_en(FT) = (; ρatke = FT(0), ρaHvar = FT(0), ρaQTvar = FT(0), ρaHQTcov = 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), pr = (; qr = FT(0)), diff --git a/integration_tests/utils/mse_tables.jl b/integration_tests/utils/mse_tables.jl index 92e73f869..2fc5e3bb5 100644 --- a/integration_tests/utils/mse_tables.jl +++ b/integration_tests/utils/mse_tables.jl @@ -5,149 +5,149 @@ all_best_mse = OrderedCollections.OrderedDict() # all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict() -all_best_mse["ARM_SGP"]["qt_mean"] = 0.22143401926585182 -all_best_mse["ARM_SGP"]["updraft_area"] = 335.16785513690655 -all_best_mse["ARM_SGP"]["updraft_w"] = 141.97803346792176 -all_best_mse["ARM_SGP"]["updraft_qt"] = 27.57483064595048 -all_best_mse["ARM_SGP"]["updraft_thetal"] = 170.98103520062713 +all_best_mse["ARM_SGP"]["qt_mean"] = 0.22336972794853463 +all_best_mse["ARM_SGP"]["updraft_area"] = 335.39395304218874 +all_best_mse["ARM_SGP"]["updraft_w"] = 142.04595107815257 +all_best_mse["ARM_SGP"]["updraft_qt"] = 27.608269061958264 +all_best_mse["ARM_SGP"]["updraft_thetal"] = 170.98104312978447 all_best_mse["ARM_SGP"]["u_mean"] = 1.3375737467153984e-5 -all_best_mse["ARM_SGP"]["tke_mean"] = 1317.678548063657 -all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00010869927305570782 -all_best_mse["ARM_SGP"]["ql_mean"] = 176.03297501190372 +all_best_mse["ARM_SGP"]["tke_mean"] = 1318.6602293296348 +all_best_mse["ARM_SGP"]["temperature_mean"] = 0.000109470808468741 +all_best_mse["ARM_SGP"]["ql_mean"] = 174.19144829957503 all_best_mse["ARM_SGP"]["qi_mean"] = "NA" -all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00010012167033603005 -all_best_mse["ARM_SGP"]["Hvar_mean"] = 10019.605203648243 -all_best_mse["ARM_SGP"]["QTvar_mean"] = 6380.233899055468 +all_best_mse["ARM_SGP"]["thetal_mean"] = 0.0001008500529941707 +all_best_mse["ARM_SGP"]["Hvar_mean"] = 10272.271366887158 +all_best_mse["ARM_SGP"]["QTvar_mean"] = 6520.8612887442105 # all_best_mse["Bomex"] = OrderedCollections.OrderedDict() -all_best_mse["Bomex"]["qt_mean"] = 0.11224921120454945 -all_best_mse["Bomex"]["updraft_area"] = 129.4906967797644 -all_best_mse["Bomex"]["updraft_w"] = 18.517052531724794 -all_best_mse["Bomex"]["updraft_qt"] = 6.542125775887841 -all_best_mse["Bomex"]["updraft_thetal"] = 69.46310949586267 -all_best_mse["Bomex"]["v_mean"] = 67.05698543673745 -all_best_mse["Bomex"]["u_mean"] = 0.3235803034219502 -all_best_mse["Bomex"]["tke_mean"] = 74.27500076671177 -all_best_mse["Bomex"]["temperature_mean"] = 4.4589394073851936e-5 -all_best_mse["Bomex"]["ql_mean"] = 8.84060550534244 +all_best_mse["Bomex"]["qt_mean"] = 0.11089388922127488 +all_best_mse["Bomex"]["updraft_area"] = 129.55617165620623 +all_best_mse["Bomex"]["updraft_w"] = 18.114902960628168 +all_best_mse["Bomex"]["updraft_qt"] = 6.271258519434218 +all_best_mse["Bomex"]["updraft_thetal"] = 69.45693156497397 +all_best_mse["Bomex"]["v_mean"] = 66.90416426051955 +all_best_mse["Bomex"]["u_mean"] = 0.3230366720242709 +all_best_mse["Bomex"]["tke_mean"] = 74.24050507191417 +all_best_mse["Bomex"]["temperature_mean"] = 4.3987566397998385e-5 +all_best_mse["Bomex"]["ql_mean"] = 8.158894286617398 all_best_mse["Bomex"]["qi_mean"] = "NA" -all_best_mse["Bomex"]["thetal_mean"] = 4.528961998613632e-5 -all_best_mse["Bomex"]["Hvar_mean"] = 3999.7718503825167 -all_best_mse["Bomex"]["QTvar_mean"] = 1496.3117112659884 +all_best_mse["Bomex"]["thetal_mean"] = 4.4643405606181006e-5 +all_best_mse["Bomex"]["Hvar_mean"] = 4145.7703137013905 +all_best_mse["Bomex"]["QTvar_mean"] = 1546.822975591833 # all_best_mse["DryBubble"] = OrderedCollections.OrderedDict() -all_best_mse["DryBubble"]["updraft_area"] = 1.5600923378908404e-20 -all_best_mse["DryBubble"]["updraft_w"] = 3.125970197997314e-21 -all_best_mse["DryBubble"]["updraft_thetal"] = 2.2290997257401877e-28 +all_best_mse["DryBubble"]["updraft_area"] = 0.020656115332751076 +all_best_mse["DryBubble"]["updraft_w"] = 0.004892565983035499 +all_best_mse["DryBubble"]["updraft_thetal"] = 1.612480342572087e-10 all_best_mse["DryBubble"]["u_mean"] = 0.0 -all_best_mse["DryBubble"]["tke_mean"] = 8.711612946005953e-21 -all_best_mse["DryBubble"]["temperature_mean"] = 1.900994131825817e-28 -all_best_mse["DryBubble"]["thetal_mean"] = 1.4035572383204598e-28 -all_best_mse["DryBubble"]["Hvar_mean"] = 1.803251133793496e-30 +all_best_mse["DryBubble"]["tke_mean"] = 6.63871246690561 +all_best_mse["DryBubble"]["temperature_mean"] = 1.3094278804958128e-10 +all_best_mse["DryBubble"]["thetal_mean"] = 1.1428633151394506e-10 +all_best_mse["DryBubble"]["Hvar_mean"] = 5.425326274366618e-10 # all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.022674763548033078 -all_best_mse["DYCOMS_RF01"]["ql_mean"] = 10.28174651514448 -all_best_mse["DYCOMS_RF01"]["updraft_area"] = 30.54618421803614 -all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.138810210443974 -all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.1514466616728734 -all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18691260060387 -all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.003016440460640667 -all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.07896949921830682 -all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.338502376891398 -all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 8.08641879994649e-5 -all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 8.159410304007782e-5 -all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1243.0730597273637 -all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 483.9194356419018 +all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.022666046258183616 +all_best_mse["DYCOMS_RF01"]["ql_mean"] = 10.277522414893868 +all_best_mse["DYCOMS_RF01"]["updraft_area"] = 30.54643429538314 +all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.13857029949206 +all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.151426443301224 +all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.1869124644636 +all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.0030074149955970834 +all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.07909006815787396 +all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.33951353772402 +all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 8.087666801443115e-5 +all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 8.160730677044462e-5 +all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1243.0726225476485 +all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 483.91530785298096 # all_best_mse["GABLS"] = OrderedCollections.OrderedDict() all_best_mse["GABLS"]["updraft_thetal"] = 0.0 -all_best_mse["GABLS"]["v_mean"] = 3.0258956903165754e-29 -all_best_mse["GABLS"]["u_mean"] = 5.596072484918322e-31 -all_best_mse["GABLS"]["tke_mean"] = 4.93825088596756e-29 +all_best_mse["GABLS"]["v_mean"] = 2.6346373663372235e-30 +all_best_mse["GABLS"]["u_mean"] = 1.6459036720348012e-31 +all_best_mse["GABLS"]["tke_mean"] = 1.1831209659516639e-29 all_best_mse["GABLS"]["temperature_mean"] = 0.0 all_best_mse["GABLS"]["thetal_mean"] = 0.0 -all_best_mse["GABLS"]["Hvar_mean"] = 3.431124337761387e-28 -all_best_mse["GABLS"]["QTvar_mean"] = 6.2713381408194e-31 -all_best_mse["GABLS"]["qt_mean"] = 2.4026804723790953e-30 +all_best_mse["GABLS"]["Hvar_mean"] = 2.8629794631475323e-29 +all_best_mse["GABLS"]["QTvar_mean"] = 2.0665906915071658e-30 +all_best_mse["GABLS"]["qt_mean"] = 1.5002563914797642e-31 # all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict() -all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 4.2456838042947345e-11 -all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 2.039432562960598e-8 -all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 2.9943754373867077e-10 -all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 2.4372912126737867e-9 -all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 8.847211719680312e-12 -all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 5.0322704745824535e-15 -all_best_mse["life_cycle_Tan2018"]["v_mean"] = 1.744391825500784e-10 -all_best_mse["life_cycle_Tan2018"]["u_mean"] = 4.807055831555187e-13 -all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 1.4390518237145982e-10 -all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 2.1354141352137547e-14 -all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 2.1044103615447052e-14 -all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 0.00563700994258815 -all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 0.004595094880303624 +all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 1.0591223452490358e-5 +all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.00686761168249568 +all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.001312196185899228 +all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.00692406740571909 +all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.0004946712908129363 +all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 2.860414784731618e-7 +all_best_mse["life_cycle_Tan2018"]["v_mean"] = 0.00023425927300220657 +all_best_mse["life_cycle_Tan2018"]["u_mean"] = 1.348741580839515e-6 +all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.00047324877523489896 +all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 4.878384140002433e-9 +all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 4.869462330658505e-9 +all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 924.9168060470322 +all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 744.9248364265027 # all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict() -all_best_mse["Nieuwstadt"]["updraft_area"] = 98.6564977935179 -all_best_mse["Nieuwstadt"]["updraft_w"] = 12.474571100534247 -all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.36697441802559 -all_best_mse["Nieuwstadt"]["u_mean"] = 14.602576944221799 -all_best_mse["Nieuwstadt"]["tke_mean"] = 313.7437656356189 -all_best_mse["Nieuwstadt"]["temperature_mean"] = 9.708439348851666e-6 -all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.0020176199907293e-5 -all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1278.9608116518727 +all_best_mse["Nieuwstadt"]["updraft_area"] = 98.97436274691259 +all_best_mse["Nieuwstadt"]["updraft_w"] = 12.498502801664266 +all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.3669724244686 +all_best_mse["Nieuwstadt"]["u_mean"] = 14.602458009247906 +all_best_mse["Nieuwstadt"]["tke_mean"] = 313.21384305260466 +all_best_mse["Nieuwstadt"]["temperature_mean"] = 9.727469470672468e-6 +all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.0038730380773907e-5 +all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1277.5917192776851 # all_best_mse["Rico"] = OrderedCollections.OrderedDict() -all_best_mse["Rico"]["qt_mean"] = 1.240577610897265 -all_best_mse["Rico"]["updraft_area"] = 477.4723365082512 -all_best_mse["Rico"]["updraft_w"] = 107.93641059929224 -all_best_mse["Rico"]["updraft_qt"] = 12.27330510190275 -all_best_mse["Rico"]["updraft_thetal"] = 133.36345921192816 -all_best_mse["Rico"]["v_mean"] = 0.6113313125463751 -all_best_mse["Rico"]["u_mean"] = 0.6948989327737745 -all_best_mse["Rico"]["tke_mean"] = 82.26914808072146 -all_best_mse["Rico"]["temperature_mean"] = 0.0005737249256216108 -all_best_mse["Rico"]["ql_mean"] = 63.483025341334695 -all_best_mse["Rico"]["qr_mean"] = 760.8610598283299 +all_best_mse["Rico"]["qt_mean"] = 1.2353506326619745 +all_best_mse["Rico"]["updraft_area"] = 477.46844663446234 +all_best_mse["Rico"]["updraft_w"] = 108.20703071711189 +all_best_mse["Rico"]["updraft_qt"] = 12.201832089557566 +all_best_mse["Rico"]["updraft_thetal"] = 133.36308629887355 +all_best_mse["Rico"]["v_mean"] = 0.6109671977037571 +all_best_mse["Rico"]["u_mean"] = 0.6949148036003943 +all_best_mse["Rico"]["tke_mean"] = 82.26337247291202 +all_best_mse["Rico"]["temperature_mean"] = 0.0005722457194843458 +all_best_mse["Rico"]["ql_mean"] = 63.36526466328347 +all_best_mse["Rico"]["qr_mean"] = 760.8106588327406 all_best_mse["Rico"]["qi_mean"] = "NA" -all_best_mse["Rico"]["thetal_mean"] = 0.0005661696874453951 -all_best_mse["Rico"]["Hvar_mean"] = 187519.5898306167 -all_best_mse["Rico"]["QTvar_mean"] = 41933.89043513789 +all_best_mse["Rico"]["thetal_mean"] = 0.000564682405873345 +all_best_mse["Rico"]["Hvar_mean"] = 165606.9909017002 +all_best_mse["Rico"]["QTvar_mean"] = 37237.78655217467 # all_best_mse["Soares"] = OrderedCollections.OrderedDict() -all_best_mse["Soares"]["qt_mean"] = 0.125492505099034 -all_best_mse["Soares"]["updraft_area"] = 99.07613100101109 -all_best_mse["Soares"]["updraft_w"] = 11.354211142559235 -all_best_mse["Soares"]["updraft_qt"] = 23.094442446856103 -all_best_mse["Soares"]["updraft_thetal"] = 65.25330225554521 -all_best_mse["Soares"]["u_mean"] = 96.34890749728717 -all_best_mse["Soares"]["tke_mean"] = 241.84351768806397 -all_best_mse["Soares"]["temperature_mean"] = 1.0933047642371797e-5 -all_best_mse["Soares"]["thetal_mean"] = 1.0360692590042686e-5 -all_best_mse["Soares"]["Hvar_mean"] = 1126.3241374374188 +all_best_mse["Soares"]["qt_mean"] = 0.12623520584762954 +all_best_mse["Soares"]["updraft_area"] = 96.34014526986842 +all_best_mse["Soares"]["updraft_w"] = 11.390128916096858 +all_best_mse["Soares"]["updraft_qt"] = 23.080877966322717 +all_best_mse["Soares"]["updraft_thetal"] = 65.25330119519995 +all_best_mse["Soares"]["u_mean"] = 96.34451750277604 +all_best_mse["Soares"]["tke_mean"] = 241.27289775855488 +all_best_mse["Soares"]["temperature_mean"] = 1.0977617337130468e-5 +all_best_mse["Soares"]["thetal_mean"] = 1.0405848721480761e-5 +all_best_mse["Soares"]["Hvar_mean"] = 1121.2781636769805 # all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict() -all_best_mse["TRMM_LBA"]["qt_mean"] = 3.96721142480849 -all_best_mse["TRMM_LBA"]["updraft_area"] = 7590.326979740429 -all_best_mse["TRMM_LBA"]["updraft_w"] = 29122.667102286345 -all_best_mse["TRMM_LBA"]["updraft_qt"] = 264.28246394543794 -all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1962.213457812691 -all_best_mse["TRMM_LBA"]["v_mean"] = 286.62659583240327 -all_best_mse["TRMM_LBA"]["u_mean"] = 113.79375208017473 -all_best_mse["TRMM_LBA"]["tke_mean"] = 25964.486606169423 -all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.001127051101306071 -all_best_mse["TRMM_LBA"]["ql_mean"] = 11333.082496938974 +all_best_mse["TRMM_LBA"]["qt_mean"] = 3.9639773897846315 +all_best_mse["TRMM_LBA"]["updraft_area"] = 7590.5124280358195 +all_best_mse["TRMM_LBA"]["updraft_w"] = 29117.789630260653 +all_best_mse["TRMM_LBA"]["updraft_qt"] = 264.3743946842157 +all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1962.1733128842725 +all_best_mse["TRMM_LBA"]["v_mean"] = 286.64250267320244 +all_best_mse["TRMM_LBA"]["u_mean"] = 113.79784954567054 +all_best_mse["TRMM_LBA"]["tke_mean"] = 25954.332707353075 +all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0011257656568842474 +all_best_mse["TRMM_LBA"]["ql_mean"] = 11335.71595959313 all_best_mse["TRMM_LBA"]["qi_mean"] = "NA" -all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.02111934332007892 -all_best_mse["TRMM_LBA"]["Hvar_mean"] = 4344.199567761157 -all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2836.009508446295 +all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.021118882573307525 +all_best_mse["TRMM_LBA"]["Hvar_mean"] = 4344.07790745562 +all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2835.250376882707 # all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict() -all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.695140418219589 -all_best_mse["LES_driven_SCM"]["v_mean"] = 3.747051509524756 -all_best_mse["LES_driven_SCM"]["u_mean"] = 1.2506833904476498 -all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0029125238780234835 -all_best_mse["LES_driven_SCM"]["ql_mean"] = 1738.0999059121439 -all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.00322061632952989 +all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.6848941301202665 +all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7461500362773097 +all_best_mse["LES_driven_SCM"]["u_mean"] = 1.2502583540896357 +all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0029087535232275095 +all_best_mse["LES_driven_SCM"]["ql_mean"] = 1640.9315832454179 +all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0032160005045022413 # ################################# ################################# diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index 55930c29f..9bcdb608c 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -62,7 +62,6 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p 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) aux_en = center_aux_environment(state) prog_pr = center_prog_precipitation(state) @@ -93,32 +92,32 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, p @inbounds for k in real_center_indices(grid) if ( - prog_en.QTvar[k] > epsilon && - prog_en.Hvar[k] > epsilon && - abs(prog_en.HQTcov[k]) > epsilon && + aux_en.QTvar[k] > epsilon && + aux_en.Hvar[k] > epsilon && + abs(aux_en.HQTcov[k]) > epsilon && aux_en.q_tot[k] > epsilon && - sqrt(prog_en.QTvar[k]) < aux_en.q_tot[k] + sqrt(aux_en.QTvar[k]) < aux_en.q_tot[k] ) if en_thermo.quadrature_type == "log-normal" # Lognormal parameters (mu, sd) from mean and variance - sd_q = sqrt(log(prog_en.QTvar[k] / aux_en.q_tot[k] / aux_en.q_tot[k] + 1.0)) - sd_h = sqrt(log(prog_en.Hvar[k] / aux_en.θ_liq_ice[k] / aux_en.θ_liq_ice[k] + 1.0)) + sd_q = sqrt(log(aux_en.QTvar[k] / aux_en.q_tot[k] / aux_en.q_tot[k] + 1.0)) + sd_h = sqrt(log(aux_en.Hvar[k] / aux_en.θ_liq_ice[k] / aux_en.θ_liq_ice[k] + 1.0)) # Enforce Schwarz"s inequality - corr = max(min(prog_en.HQTcov[k] / sqrt(prog_en.Hvar[k] * prog_en.QTvar[k]), 1.0), -1.0) + corr = max(min(aux_en.HQTcov[k] / sqrt(aux_en.Hvar[k] * aux_en.QTvar[k]), 1.0), -1.0) sd2_hq = - log(corr * sqrt(prog_en.Hvar[k] * prog_en.QTvar[k]) / aux_en.θ_liq_ice[k] / aux_en.q_tot[k] + 1.0) + log(corr * sqrt(aux_en.Hvar[k] * aux_en.QTvar[k]) / aux_en.θ_liq_ice[k] / aux_en.q_tot[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(aux_en.q_tot[k] * aux_en.q_tot[k] / sqrt(aux_en.q_tot[k] * aux_en.q_tot[k] + prog_en.QTvar[k])) + log(aux_en.q_tot[k] * aux_en.q_tot[k] / sqrt(aux_en.q_tot[k] * aux_en.q_tot[k] + aux_en.QTvar[k])) mu_h = log( aux_en.θ_liq_ice[k] * aux_en.θ_liq_ice[k] / - sqrt(aux_en.θ_liq_ice[k] * aux_en.θ_liq_ice[k] + prog_en.Hvar[k]), + sqrt(aux_en.θ_liq_ice[k] * aux_en.θ_liq_ice[k] + aux_en.Hvar[k]), ) else - 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) + sd_q = sqrt(aux_en.QTvar[k]) + sd_h = sqrt(aux_en.Hvar[k]) + corr = max(min(aux_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 - aux_en.q_tot[k]) / (sqrt2 * abscissas[1]) diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index a7cfbe75a..5235f3171 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -278,10 +278,10 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca compute_gm_tendencies!(edmf, grid, state, Case, gm, TS) compute_updraft_tendencies(edmf, grid, state, gm) - parent(tendencies_en.tke) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :tke, n_updrafts) - parent(tendencies_en.Hvar) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :Hvar, n_updrafts) - parent(tendencies_en.QTvar) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :QTvar, n_updrafts) - parent(tendencies_en.HQTcov) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :HQTcov, n_updrafts) + parent(tendencies_en.ρatke) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :tke, n_updrafts) + parent(tendencies_en.ρaHvar) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :Hvar, n_updrafts) + parent(tendencies_en.ρaQTvar) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :QTvar, n_updrafts) + parent(tendencies_en.ρaHQTcov) .= en_diffusion_tendencies(edmf, grid, state, param_set, TS, :HQTcov, n_updrafts) ### ### update @@ -292,10 +292,10 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca end ae = 1 .- aux_tc.bulk.area for k in real_center_indices(grid) - prog_en.tke[k] += TS.dt * tendencies_en.tke[k] / (ρ0_c[k] * ae[k]) - prog_en.Hvar[k] += TS.dt * tendencies_en.Hvar[k] / (ρ0_c[k] * ae[k]) - prog_en.QTvar[k] += TS.dt * tendencies_en.QTvar[k] / (ρ0_c[k] * ae[k]) - prog_en.HQTcov[k] += TS.dt * tendencies_en.HQTcov[k] / (ρ0_c[k] * ae[k]) + prog_en.ρatke[k] += TS.dt * tendencies_en.ρatke[k] + prog_en.ρaHvar[k] += TS.dt * tendencies_en.ρaHvar[k] + prog_en.ρaQTvar[k] += TS.dt * tendencies_en.ρaQTvar[k] + prog_en.ρaHQTcov[k] += TS.dt * tendencies_en.ρaHQTcov[k] end @inbounds for k in real_center_indices(grid) prog_gm.u[k] += tendencies_gm.u[k] * TS.dt @@ -310,11 +310,11 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca set_edmf_surface_bc(edmf, grid, state, up, Case.Sur) filter_updraft_vars(edmf, grid, state, gm) @inbounds for k in real_center_indices(grid) - 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])) + prog_en.ρatke[k] = max(prog_en.ρatke[k], 0.0) + prog_en.ρaHvar[k] = max(prog_en.ρaHvar[k], 0.0) + prog_en.ρaQTvar[k] = max(prog_en.ρaQTvar[k], 0.0) + prog_en.ρaHQTcov[k] = max(prog_en.ρaHQTcov[k], -sqrt(prog_en.ρaHvar[k] * prog_en.ρaQTvar[k])) + prog_en.ρaHQTcov[k] = min(prog_en.ρaHQTcov[k], sqrt(prog_en.ρaHvar[k] * prog_en.ρaQTvar[k])) end return @@ -329,6 +329,7 @@ function set_edmf_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, up, surface) prog_up = center_prog_updrafts(state) prog_en = center_prog_environment(state) prog_up_f = face_prog_updrafts(state) + aux_tc = center_aux_turbconv(state) @inbounds for i in 1:(up.n_updrafts) prog_up[i].ρarea[kc_surf] = ρ_0_c[kc_surf] * edmf.area_surface_bc[i] prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * edmf.h_surface_bc[i] @@ -342,11 +343,14 @@ function set_edmf_surface_bc(edmf::EDMF_PrognosticTKE, grid, state, up, surface) ustar = surface.ustar oblength = surface.obukhov_length α0LL = center_ref_state(state).α0[kc_surf] + ae = 1 .- aux_tc.bulk.area # area of environment + + ρ0_ae = ρ_0_c[kc_surf] * ae[kc_surf] - prog_en.tke[kc_surf] = get_surface_tke(surface.ustar, grid.zc[kc_surf], surface.obukhov_length) - 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) + prog_en.ρatke[kc_surf] = ρ0_ae * get_surface_tke(surface.ustar, grid.zc[kc_surf], surface.obukhov_length) + prog_en.ρaHvar[kc_surf] = ρ0_ae * get_surface_variance(flux1 * α0LL, flux1 * α0LL, ustar, zLL, oblength) + prog_en.ρaQTvar[kc_surf] = ρ0_ae * get_surface_variance(flux2 * α0LL, flux2 * α0LL, ustar, zLL, oblength) + prog_en.ρaHQTcov[kc_surf] = ρ0_ae * get_surface_variance(flux1 * α0LL, flux2 * α0LL, ustar, zLL, oblength) end @@ -405,7 +409,7 @@ function get_GMV_CoVar(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, aux_en = is_tke ? aux_en_f : aux_en_c aux_up = center_aux_updrafts(state) gmv_covar = getproperty(center_aux_grid_mean(state), covar_sym) - covar_e = getproperty(center_prog_environment(state), covar_sym) + covar_e = getproperty(center_aux_environment(state), covar_sym) prog_gm = is_tke ? prog_gm_f : prog_gm_c ϕ_gm = getproperty(prog_gm, ϕ_sym) ψ_gm = getproperty(prog_gm, ψ_sym) @@ -759,9 +763,8 @@ function compute_covariance_entr( 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) aux_en_c = center_aux_environment(state) + covar = getproperty(aux_en_c, covar_sym) aux_en_f = face_aux_environment(state) aux_en = is_tke ? aux_en_f : aux_en_c EnvVar1 = getproperty(aux_en, var1) @@ -803,7 +806,7 @@ function compute_covariance_entr( ((envvar1 - gmvvar1) * (updvar2 - envvar2) + (envvar2 - gmvvar2) * (updvar1 - envvar1)) 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] + tke_factor * ρ_0_c[k] * a_up * abs(w_u) * (edmf.entr_sc[i, k] + eps_turb) * covar[k] end end end @@ -818,8 +821,8 @@ function compute_covariance_detr(edmf::EDMF_PrognosticTKE, grid, state, covar_sy aux_en_2m = center_aux_environment_2m(state) aux_covar = getproperty(aux_en_2m, covar_sym) aux_up = center_aux_updrafts(state) - prog_en = center_prog_environment(state) - prog_covar = getproperty(prog_en, covar_sym) + aux_en = center_aux_environment(state) + covar = getproperty(aux_en, covar_sym) @inbounds for k in real_center_indices(grid) aux_covar.detr_loss[k] = 0.0 @@ -827,7 +830,7 @@ function compute_covariance_detr(edmf::EDMF_PrognosticTKE, grid, state, covar_sy w_up_c = interpf2c(aux_up_f[i].w, grid, k) aux_covar.detr_loss[k] += aux_up[i].area[k] * abs(w_up_c) * edmf.entr_sc[i, k] end - aux_covar.detr_loss[k] *= ρ0_c[k] * prog_covar[k] + aux_covar.detr_loss[k] *= ρ0_c[k] * covar[k] end return end @@ -839,14 +842,14 @@ function compute_covariance_dissipation(edmf::EDMF_PrognosticTKE, grid, state, c ae = 1 .- aux_tc.bulk.area ρ0_c = center_ref_state(state).ρ0 prog_en = center_prog_environment(state) + aux_en = center_aux_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) + covar = getproperty(aux_en, covar_sym) @inbounds for k in real_center_indices(grid) 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) + (ρ0_c[k] * ae[k] * covar[k] * max(aux_en.tke[k], 0)^0.5 / max(edmf.mixing_length[k], 1.0e-3) * c_d) end return end @@ -861,7 +864,8 @@ function en_diffusion_tendencies(edmf, grid::Grid, state, param_set, TS, covar_s prog_en = center_prog_environment(state) aux_en_2m = center_aux_environment_2m(state) aux_up_f = face_aux_updrafts(state) - prog_covar = getproperty(prog_en, covar_sym) + aux_en = center_aux_environment(state) + covar = getproperty(aux_en, covar_sym) aux_covar = getproperty(aux_en_2m, covar_sym) aux_up = center_aux_updrafts(state) w_en = face_aux_environment(state).w @@ -883,7 +887,7 @@ function en_diffusion_tendencies(edmf, grid::Grid, state, param_set, TS, covar_s @inbounds for k in real_face_indices(grid) ρ_ae_K[k] = interpc2f(aeK, grid, k; aeK_bcs...) * ρ0_f[k] - ϕ_dual = dual_centers(prog_covar, grid, k) + ϕ_dual = dual_centers(covar, grid, k) ρ_ae_K∇ϕ[k] = ρ_ae_K[k] * ∇c2f(ϕ_dual, grid, k; prog_bcs...) end @inbounds for k in real_center_indices(grid) @@ -907,11 +911,11 @@ function en_diffusion_tendencies(edmf, grid::Grid, state, param_set, TS, covar_s @inbounds for k in real_center_indices(grid) w_en_c[k] = interpf2c(w_en, grid, k) - ρaew_en_ϕ[k] = ρ0_c[k] * ae[k] * w_en_c[k] * prog_covar[k] + ρaew_en_ϕ[k] = ρ0_c[k] * ae[k] * w_en_c[k] * covar[k] end kc_surf = kc_surface(grid) - covar_surf = prog_covar[kc_surf] + covar_surf = covar[kc_surf] @inbounds for k in real_center_indices(grid) D_env = sum(1:n_updrafts) do i @@ -925,7 +929,7 @@ function en_diffusion_tendencies(edmf, grid::Grid, state, param_set, TS, covar_s end D_env_i end - dissipation = ρ0_c[k] * ae[k] * c_d * sqrt(max(prog_en.tke[k], 0)) / max(mixing_length[k], 1) + dissipation = ρ0_c[k] * ae[k] * c_d * sqrt(max(aux_en.tke[k], 0)) / max(mixing_length[k], 1) ρaew_en_ϕ_cut = ccut_downwind(ρaew_en_ϕ, grid, k) ∇ρaew_en_ϕ = c∇_downwind(ρaew_en_ϕ_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) if is_surface_center(grid, k) @@ -936,8 +940,7 @@ function en_diffusion_tendencies(edmf, grid::Grid, state, param_set, TS, covar_s aux_covar.buoy[k] + aux_covar.shear[k] + aux_covar.entr_gain[k] + - aux_covar.rain_src[k] - D_env * prog_covar[k] - dissipation * prog_covar[k] - ∇ρaew_en_ϕ + - ∇ρ_ae_K∇ϕ[k] + aux_covar.rain_src[k] - D_env * covar[k] - dissipation * covar[k] - ∇ρaew_en_ϕ + ∇ρ_ae_K∇ϕ[k] ) end end @@ -955,7 +958,7 @@ function GMV_third_m(edmf::EDMF_PrognosticTKE, grid, state, covar_en_sym::Symbol ae = 1 .- aux_tc.bulk.area aux_up_f = face_aux_updrafts(state) is_tke = covar_en_sym == :tke - covar_en = getproperty(center_prog_environment(state), covar_en_sym) + covar_en = getproperty(center_aux_environment(state), covar_en_sym) aux_en_c = center_aux_environment(state) aux_en_f = face_aux_environment(state) aux_en = is_tke ? aux_en_f : aux_en_c diff --git a/src/diagnostics.jl b/src/diagnostics.jl index 8b627bf3b..8cc9e6b9b 100644 --- a/src/diagnostics.jl +++ b/src/diagnostics.jl @@ -64,11 +64,10 @@ 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_turbconv(state).KM), "eddy_diffusivity" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(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), + "env_tke" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment(state).tke), + "env_Hvar" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment(state).Hvar), + "env_QTvar" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_environment(state).QTvar), + "env_HQTcov" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_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), diff --git a/src/update_aux.jl b/src/update_aux.jl index 70dd57cb3..f2b01d277 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -80,6 +80,14 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) for k in real_center_indices(grid) aux_tc.bulk.area[k] = sum(ntuple(i -> aux_up[i].area[k], up.n_updrafts)) end + ae = 1 .- aux_tc.bulk.area # area of environment + + for k in real_center_indices(grid) + aux_en.tke[k] = prog_en.ρatke[k] / (ρ0_c[k] * ae[k]) + aux_en.Hvar[k] = prog_en.ρaHvar[k] / (ρ0_c[k] * ae[k]) + aux_en.QTvar[k] = prog_en.ρaQTvar[k] / (ρ0_c[k] * ae[k]) + aux_en.HQTcov[k] = prog_en.ρaHQTcov[k] / (ρ0_c[k] * ae[k]) + end ##### ##### diagnose_GMV_moments @@ -261,7 +269,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) w_en = interpf2c(aux_en_f.w, grid, k), b_up = aux_up[i].buoy[k], b_en = aux_en.buoy[k], - tke = prog_en.tke[k], + tke = aux_en.tke[k], dMdz = ∇m, M = m, a_up = aux_up[i].area[k], @@ -442,13 +450,13 @@ 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 = prog_en.tke[kc_surf], + tke_surf = aux_en.tke[kc_surf], ustar = surface.ustar, Pr = edmf.prandtl_nvec[k], p0 = p0_c[k], ∇b = bg, Shear² = Shear², - tke = prog_en.tke[k], + tke = aux_en.tke[k], a_en = (1 - aux_tc.bulk.area[k]), wc_en = wc_en, wc_up = Tuple(wc_up), @@ -463,7 +471,7 @@ 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(prog_en.tke[k], 0.0)) + KM[k] = c_m * edmf.mixing_length[k] * sqrt(max(aux_en.tke[k], 0.0)) KH[k] = KM[k] / edmf.prandtl_nvec[k] aux_en_2m.tke.buoy[k] = -ml_model.a_en * ρ0_c[k] * KH[k] * bg.∂b∂z @@ -501,7 +509,6 @@ 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) 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]