From 6826e55e5285230c5bb813ded060085cea0bceb6 Mon Sep 17 00:00:00 2001 From: yairchn Date: Fri, 22 Oct 2021 14:15:48 -0700 Subject: [PATCH] switch to conservative variables in the updrafts --- integration_tests/utils/initial_conditions.jl | 40 +-- integration_tests/utils/main.jl | 8 +- integration_tests/utils/mse_tables.jl | 234 ++++++------- src/EDMF_Updrafts.jl | 18 +- src/Operators.jl | 10 +- src/Turbulence_PrognosticTKE.jl | 311 ++++++++++-------- src/types.jl | 3 - src/update_aux.jl | 59 ++-- 8 files changed, 361 insertions(+), 322 deletions(-) diff --git a/integration_tests/utils/initial_conditions.jl b/integration_tests/utils/initial_conditions.jl index d354ad82a7..73d5fb965f 100644 --- a/integration_tests/utils/initial_conditions.jl +++ b/integration_tests/utils/initial_conditions.jl @@ -33,33 +33,34 @@ end function initialize_updrafts(edmf, grid, state, up::TC.UpdraftVariables, gm::TC.GridMeanVariables) kc_surf = TC.kc_surface(grid) - - prog_up = TC.center_prog_updrafts(state) + aux_up = TC.center_aux_updrafts(state) prog_gm = TC.center_prog_grid_mean(state) aux_up = TC.center_aux_updrafts(state) - prog_up_f = TC.face_prog_updrafts(state) + aux_up_f = TC.face_aux_updrafts(state) aux_gm = TC.center_aux_grid_mean(state) + prog_up = TC.center_prog_updrafts(state) + prog_up_f = TC.face_prog_updrafts(state) @inbounds for i in 1:(up.n_updrafts) @inbounds for k in TC.real_face_indices(grid) - prog_up_f[i].w[k] = 0 + aux_up_f[i].w[k] = 0 + prog_up_f[i].ρaw[k] = 0 end @inbounds for k in TC.real_center_indices(grid) aux_up[i].buoy[k] = 0 # Simple treatment for now, revise when multiple updraft closures # become more well defined - if up.prognostic - prog_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts - else - prog_up[i].area[k] = up.updraft_fraction / up.n_updrafts - end - prog_up[i].q_tot[k] = prog_gm.q_tot[k] - prog_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] + aux_up[i].area[k] = 0 + aux_up[i].q_tot[k] = prog_gm.q_tot[k] + aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] aux_up[i].q_liq[k] = aux_gm.q_liq[k] aux_up[i].T[k] = aux_gm.T[k] + prog_up[i].ρarea[k] = 0 + prog_up[i].ρaq_tot[k] = 0 + prog_up[i].ρaθ_liq_ice[k] = 0 end - prog_up[i].area[kc_surf] = up.updraft_fraction / up.n_updrafts + aux_up[i].area[kc_surf] = up.updraft_fraction / up.n_updrafts end return end @@ -134,9 +135,8 @@ function initialize_updrafts_DryBubble(edmf, grid, state, up::TC.UpdraftVariable 264.1574, 263.6518, 263.1461, 262.6451, 262.1476, 261.6524] #! format: on - prog_up = TC.center_prog_updrafts(state) aux_up = TC.center_aux_updrafts(state) - prog_up_f = TC.face_prog_updrafts(state) + aux_up_f = TC.face_aux_updrafts(state) aux_gm = TC.center_aux_grid_mean(state) prog_gm = TC.center_prog_grid_mean(state) Area_in = TC.pyinterp(grid.zc, z_in, Area_in) @@ -145,22 +145,22 @@ function initialize_updrafts_DryBubble(edmf, grid, state, up::TC.UpdraftVariable @inbounds for i in 1:(up.n_updrafts) @inbounds for k in TC.real_face_indices(grid) if minimum(z_in) <= grid.zf[k] <= maximum(z_in) - prog_up_f[i].w[k] = 0.0 + aux_up_f[i].w[k] = 0.0 end end @inbounds for k in TC.real_center_indices(grid) if minimum(z_in) <= grid.zc[k] <= maximum(z_in) - prog_up[i].area[k] = Area_in[k] #up.updraft_fraction/up.n_updrafts - prog_up[i].θ_liq_ice[k] = θ_liq_in[k] - prog_up[i].q_tot[k] = 0.0 + aux_up[i].area[k] = Area_in[k] #up.updraft_fraction/up.n_updrafts + aux_up[i].θ_liq_ice[k] = θ_liq_in[k] + aux_up[i].q_tot[k] = 0.0 aux_up[i].q_liq[k] = 0.0 # for now temperature is provided as diagnostics from LES aux_up[i].T[k] = T_in[k] else - prog_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts - prog_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] + aux_up[i].area[k] = 0.0 #up.updraft_fraction/up.n_updrafts + aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] aux_up[i].T[k] = aux_gm.T[k] end end diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl index 8b8d434663..9bf85503ae 100644 --- a/integration_tests/utils/main.jl +++ b/integration_tests/utils/main.jl @@ -62,7 +62,7 @@ cent_aux_vars_en_2m(FT) = (; rain_src = FT(0), ) cent_aux_vars_up(FT) = - (; q_liq = FT(0), T = FT(0), RH = FT(0), buoy = FT(0), area_new = FT(0), q_tot_new = FT(0), θ_liq_ice_new = FT(0)) + (; q_liq = FT(0), T = FT(0), RH = FT(0), buoy = FT(0), area = FT(0), q_tot = FT(0), θ_liq_ice = 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)), @@ -91,7 +91,7 @@ cent_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT)..., # Face only face_aux_vars_gm(FT) = () -face_aux_vars_up(FT) = (; w_new = FT(0)) +face_aux_vars_up(FT) = (; w = FT(0)) face_aux_vars_edmf(FT, n_up) = (; turbconv = (; @@ -125,7 +125,7 @@ face_diagnostic_vars(FT, n_up) = (; face_diagnostic_vars_gm(FT)..., face_diagnos # Center only 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_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_edmf(FT, n_up) = (; turbconv = (; @@ -136,7 +136,7 @@ cent_prognostic_vars_edmf(FT, n_up) = (; # Face only face_prognostic_vars(FT, n_up) = (; w = FT(0), face_prognostic_vars_edmf(FT, n_up)...) -face_prognostic_vars_up(FT) = (; w = FT(0)) +face_prognostic_vars_up(FT) = (; ρaw = FT(0)) face_prognostic_vars_edmf(FT, n_up) = (; turbconv = (; up = ntuple(i -> face_prognostic_vars_up(FT), n_up))) # face_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model diff --git a/integration_tests/utils/mse_tables.jl b/integration_tests/utils/mse_tables.jl index 69e9326d5c..d89db90b67 100644 --- a/integration_tests/utils/mse_tables.jl +++ b/integration_tests/utils/mse_tables.jl @@ -5,145 +5,145 @@ all_best_mse = OrderedCollections.OrderedDict() # all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict() -all_best_mse["ARM_SGP"]["qt_mean"] = 2.6809473222721925e-01 -all_best_mse["ARM_SGP"]["updraft_area"] = 342.05657568721784 -all_best_mse["ARM_SGP"]["updraft_w"] = 1.4798746522124213e+02 -all_best_mse["ARM_SGP"]["updraft_qt"] = 28.294790221283453 -all_best_mse["ARM_SGP"]["updraft_thetal"] = 170.9801640138966 +all_best_mse["ARM_SGP"]["qt_mean"] = 0.26734824703687876 +all_best_mse["ARM_SGP"]["updraft_area"] = 342.0195830020112 +all_best_mse["ARM_SGP"]["updraft_w"] = 147.92306653653108 +all_best_mse["ARM_SGP"]["updraft_qt"] = 28.619520341027982 +all_best_mse["ARM_SGP"]["updraft_thetal"] = 170.98037684720282 all_best_mse["ARM_SGP"]["u_mean"] = 5.762470982857257e8 -all_best_mse["ARM_SGP"]["tke_mean"] = 986.3675515161622 -all_best_mse["ARM_SGP"]["temperature_mean"] = 1.2131358476290643e-04 -all_best_mse["ARM_SGP"]["ql_mean"] = 2.0606978346729218e+02 -all_best_mse["ARM_SGP"]["thetal_mean"] = 1.1201123866175818e-04 -all_best_mse["ARM_SGP"]["Hvar_mean"] = 2.9194184175128917e+03 -all_best_mse["ARM_SGP"]["QTvar_mean"] = 1.7213808112900319e+03 +all_best_mse["ARM_SGP"]["tke_mean"] = 982.5851594699658 +all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00012114509170440479 +all_best_mse["ARM_SGP"]["ql_mean"] = 205.7362434059359 +all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00011183644773004053 +all_best_mse["ARM_SGP"]["Hvar_mean"] = 2954.8565861580414 +all_best_mse["ARM_SGP"]["QTvar_mean"] = 1758.2118726364097 # all_best_mse["Bomex"] = OrderedCollections.OrderedDict() -all_best_mse["Bomex"]["qt_mean"] = 1.1110729358738076e-01 -all_best_mse["Bomex"]["updraft_area"] = 1.3209194605797390e+02 -all_best_mse["Bomex"]["updraft_w"] = 1.8757855416502004e+01 -all_best_mse["Bomex"]["updraft_qt"] = 6.5583158094016953e+00 -all_best_mse["Bomex"]["updraft_thetal"] = 6.9463116782644903e+01 -all_best_mse["Bomex"]["v_mean"] = 6.6112855769976591e+01 -all_best_mse["Bomex"]["u_mean"] = 2.2443513442723633e+03 -all_best_mse["Bomex"]["tke_mean"] = 5.5866550617966134e+01 -all_best_mse["Bomex"]["temperature_mean"] = 4.4211958191944796e-05 -all_best_mse["Bomex"]["ql_mean"] = 8.8506606580011589e+00 -all_best_mse["Bomex"]["thetal_mean"] = 4.4860305896182761e-05 -all_best_mse["Bomex"]["Hvar_mean"] = 3.8533215085304478e+03 -all_best_mse["Bomex"]["QTvar_mean"] = 3.8533215085304478e+03 +all_best_mse["Bomex"]["qt_mean"] = 0.11110189419742686 +all_best_mse["Bomex"]["updraft_area"] = 132.08275019675446 +all_best_mse["Bomex"]["updraft_w"] = 18.6888286448069 +all_best_mse["Bomex"]["updraft_qt"] = 6.561125866595388 +all_best_mse["Bomex"]["updraft_thetal"] = 69.46250316275639 +all_best_mse["Bomex"]["v_mean"] = 65.8201562226044 +all_best_mse["Bomex"]["u_mean"] = 2244.2916613239363 +all_best_mse["Bomex"]["tke_mean"] = 55.89170113065913 +all_best_mse["Bomex"]["temperature_mean"] = 4.4233235317380045e-5 +all_best_mse["Bomex"]["ql_mean"] = 8.67853811855733 +all_best_mse["Bomex"]["thetal_mean"] = 4.4890898584002073e-5 +all_best_mse["Bomex"]["Hvar_mean"] = 3961.7043325928034 +all_best_mse["Bomex"]["QTvar_mean"] = 1484.8229281612032 # all_best_mse["DryBubble"] = OrderedCollections.OrderedDict() -all_best_mse["DryBubble"]["updraft_area"] = 1.0700681545005437e-02 -all_best_mse["DryBubble"]["updraft_w"] = 5.1150944866783517e-02 -all_best_mse["DryBubble"]["updraft_thetal"] = 1.2223166344487256e-10 +all_best_mse["DryBubble"]["updraft_area"] = 1015.7706875956121 +all_best_mse["DryBubble"]["updraft_w"] = 634.5754919305659 +all_best_mse["DryBubble"]["updraft_thetal"] = 6.701724935855835e-5 all_best_mse["DryBubble"]["u_mean"] = 0.0 -all_best_mse["DryBubble"]["tke_mean"] = 1.7102419613097049e+00 -all_best_mse["DryBubble"]["temperature_mean"] = 9.1562819359185015e-06 -all_best_mse["DryBubble"]["thetal_mean"] = 9.7916743702040589e-07 -all_best_mse["DryBubble"]["Hvar_mean"] = 5.8429697829269229e+04 +all_best_mse["DryBubble"]["tke_mean"] = 2214.9965297933595 +all_best_mse["DryBubble"]["temperature_mean"] = 2.9080763228912202e-5 +all_best_mse["DryBubble"]["thetal_mean"] = 1.850965399866035e-5 +all_best_mse["DryBubble"]["Hvar_mean"] = 40011.62474788041 # all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF01"]["qt_mean"] = 3.0860725102523578e-02 -all_best_mse["DYCOMS_RF01"]["ql_mean"] = 1.0436898916957983e+01 -all_best_mse["DYCOMS_RF01"]["updraft_area"] = 3.0990562047789613e+01 -all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.1407525231791347e+00 -all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.1583200264873770e+00 -all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 4.6186921370083461e+01 -all_best_mse["DYCOMS_RF01"]["v_mean"] = 1.9100567571407046e+04 -all_best_mse["DYCOMS_RF01"]["u_mean"] = 7.5180517589525669e+04 -all_best_mse["DYCOMS_RF01"]["tke_mean"] = 2.4168930314928438e+01 -all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 1.5183396753484175e-04 -all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 1.5269039216317862e-04 -all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1.2556461967314308e+03 -all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 5.1800594121869506e+02 +all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.02302043313301315 +all_best_mse["DYCOMS_RF01"]["ql_mean"] = 9.997053817264641 +all_best_mse["DYCOMS_RF01"]["updraft_area"] = 30.996336098633655 +all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.143254766653856 +all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.1568446070984995 +all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18691985678771 +all_best_mse["DYCOMS_RF01"]["v_mean"] = 19075.550345532716 +all_best_mse["DYCOMS_RF01"]["u_mean"] = 75082.39928069651 +all_best_mse["DYCOMS_RF01"]["tke_mean"] = 20.93918292858011 +all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 8.031992871110814e-5 +all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 8.105708155388748e-5 +all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1243.0838962717608 +all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 484.1027376489701 # all_best_mse["GABLS"] = OrderedCollections.OrderedDict() -all_best_mse["GABLS"]["updraft_thetal"] = 5.584864486190616e-9 -all_best_mse["GABLS"]["v_mean"] = 0.0007865603279729375 -all_best_mse["GABLS"]["u_mean"] = 6.285670403414794e-5 -all_best_mse["GABLS"]["tke_mean"] = 0.0010760513032549105 -all_best_mse["GABLS"]["temperature_mean"] = 5.535007975017039e-9 -all_best_mse["GABLS"]["thetal_mean"] = 5.565410639688741e-9 -all_best_mse["GABLS"]["Hvar_mean"] = 0.001721862347778764 -all_best_mse["GABLS"]["QTvar_mean"] = 0.017318287479694217 -all_best_mse["GABLS"]["qt_mean"] = 0.21092607675396896 +all_best_mse["GABLS"]["updraft_thetal"] = 2.298193915042689e-12 +all_best_mse["GABLS"]["v_mean"] = 2.7616644378557842e-17 +all_best_mse["GABLS"]["u_mean"] = 1.1285381396515641e-18 +all_best_mse["GABLS"]["tke_mean"] = 3.545473808583206e-17 +all_best_mse["GABLS"]["temperature_mean"] = 2.0068841908688136e-22 +all_best_mse["GABLS"]["thetal_mean"] = 1.9938479687095057e-22 +all_best_mse["GABLS"]["Hvar_mean"] = 5.4556612252050644e-17 +all_best_mse["GABLS"]["QTvar_mean"] = 1.665849985115634e-17 +all_best_mse["GABLS"]["qt_mean"] = 5.835381313381845e-18 # all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict() -all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 8.1384024363088774e-07 -all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 4.3078794816675181e-03 -all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 1.8680603897032191e-03 -all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 2.6607652511720569e-03 -all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 9.1741668270320462e-03 -all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 5.4555871248086177e-06 -all_best_mse["life_cycle_Tan2018"]["v_mean"] = 1.6737318001164233e-03 -all_best_mse["life_cycle_Tan2018"]["u_mean"] = 2.8106431761120461e-06 -all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 8.1161256118218120e-04 -all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 4.6974217695746801e-12 -all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 4.6149603558343825e-12 -all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 2.0410104346414277e+00 -all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 1.7707992854546521e+00 +all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.00018210055205276554 +all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.04670196125126055 +all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.1529913251325874 +all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 1.696570760350363 +all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.029580741673237928 +all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 1.0497730342759248e-5 +all_best_mse["life_cycle_Tan2018"]["v_mean"] = 0.0018955431809035488 +all_best_mse["life_cycle_Tan2018"]["u_mean"] = 4.355490292733853e-6 +all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.01288090795259078 +all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 7.816226544606813e-8 +all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 7.699403927836914e-8 +all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 321.4513636804809 +all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 287.1019803918004 # all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict() -all_best_mse["Nieuwstadt"]["updraft_area"] = 110.65702527331422 -all_best_mse["Nieuwstadt"]["updraft_w"] = 1.2315877017147590e+01 -all_best_mse["Nieuwstadt"]["updraft_thetal"] = 1.1736704568882135e+02 -all_best_mse["Nieuwstadt"]["u_mean"] = 59.73642032443851 -all_best_mse["Nieuwstadt"]["tke_mean"] = 2.2229785052178343e+02 -all_best_mse["Nieuwstadt"]["temperature_mean"] = 9.78812700698016e-6 -all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.0097201478188357e-5 -all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1.2769135065447590e+03 +all_best_mse["Nieuwstadt"]["updraft_area"] = 110.58583393056665 +all_best_mse["Nieuwstadt"]["updraft_w"] = 12.310615938272223 +all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.36704524621437 +all_best_mse["Nieuwstadt"]["u_mean"] = 59.73582873995032 +all_best_mse["Nieuwstadt"]["tke_mean"] = 222.24937865750465 +all_best_mse["Nieuwstadt"]["temperature_mean"] = 9.793124899429011e-6 +all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.010229875340399e-5 +all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1276.0207133521765 # all_best_mse["Rico"] = OrderedCollections.OrderedDict() -all_best_mse["Rico"]["qt_mean"] = 1.3679878955046347e+00 -all_best_mse["Rico"]["updraft_area"] = 483.4022083273521 -all_best_mse["Rico"]["updraft_w"] = 1.0605770172574749e+02 -all_best_mse["Rico"]["updraft_qt"] = 1.1259437743694338e+01 -all_best_mse["Rico"]["updraft_thetal"] = 1.3280902833793590e+02 -all_best_mse["Rico"]["v_mean"] = 2.3916040021134300e+04 -all_best_mse["Rico"]["u_mean"] = 1.6723557502214760e+02 -all_best_mse["Rico"]["tke_mean"] = 8.3753305486816103e+01 -all_best_mse["Rico"]["temperature_mean"] = 6.1968125138192360e-04 -all_best_mse["Rico"]["ql_mean"] = 6.7537841730251046e+01 -all_best_mse["Rico"]["qr_mean"] = 7.5806983602072614e+02 -all_best_mse["Rico"]["thetal_mean"] = 6.1071856044129136e-04 -all_best_mse["Rico"]["Hvar_mean"] = 18654.16088756812 -all_best_mse["Rico"]["QTvar_mean"] = 4036.079013508425 +all_best_mse["Rico"]["qt_mean"] = 1.35174081002589 +all_best_mse["Rico"]["updraft_area"] = 483.46704720888056 +all_best_mse["Rico"]["updraft_w"] = 106.56500016772834 +all_best_mse["Rico"]["updraft_qt"] = 11.35362559214538 +all_best_mse["Rico"]["updraft_thetal"] = 132.80066662552824 +all_best_mse["Rico"]["v_mean"] = 23875.16466072175 +all_best_mse["Rico"]["u_mean"] = 166.8480977418146 +all_best_mse["Rico"]["tke_mean"] = 82.9024763171918 +all_best_mse["Rico"]["temperature_mean"] = 0.0006214599958554891 +all_best_mse["Rico"]["ql_mean"] = 68.58664944698768 +all_best_mse["Rico"]["qr_mean"] = 757.9905334607064 +all_best_mse["Rico"]["thetal_mean"] = 0.0006124442671740663 +all_best_mse["Rico"]["Hvar_mean"] = 21095.14897577657 +all_best_mse["Rico"]["QTvar_mean"] = 4343.2755082370695 # all_best_mse["Soares"] = OrderedCollections.OrderedDict() -all_best_mse["Soares"]["qt_mean"] = 1.2687206447566182e-01 -all_best_mse["Soares"]["updraft_area"] = 112.7273316705158 -all_best_mse["Soares"]["updraft_w"] = 1.1277580224270022e+01 -all_best_mse["Soares"]["updraft_qt"] = 2.3121619450049359e+01 -all_best_mse["Soares"]["updraft_thetal"] = 6.5253349774496741e+01 -all_best_mse["Soares"]["u_mean"] = 1.0182505005403843e+02 -all_best_mse["Soares"]["tke_mean"] = 1.8490467617329168e+02 -all_best_mse["Soares"]["temperature_mean"] = 1.1047172971106984e-05 -all_best_mse["Soares"]["thetal_mean"] = 1.0469297223913385e-05 -all_best_mse["Soares"]["Hvar_mean"] = 1145.5445941618634 +all_best_mse["Soares"]["qt_mean"] = 0.12696028025195605 +all_best_mse["Soares"]["updraft_area"] = 112.25410206505505 +all_best_mse["Soares"]["updraft_w"] = 11.285702848380712 +all_best_mse["Soares"]["updraft_qt"] = 23.122901619149626 +all_best_mse["Soares"]["updraft_thetal"] = 65.25334887109032 +all_best_mse["Soares"]["u_mean"] = 101.8233814645249 +all_best_mse["Soares"]["tke_mean"] = 184.92707805037557 +all_best_mse["Soares"]["temperature_mean"] = 1.1059282606650482e-5 +all_best_mse["Soares"]["thetal_mean"] = 1.0481247345363695e-5 +all_best_mse["Soares"]["Hvar_mean"] = 1144.684750257877 # all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict() -all_best_mse["TRMM_LBA"]["qt_mean"] = 3.6996721397006631e+00 -all_best_mse["TRMM_LBA"]["updraft_area"] = 7.5616268005851007e+03 -all_best_mse["TRMM_LBA"]["updraft_w"] = 1.1169494311940936e+04 -all_best_mse["TRMM_LBA"]["updraft_qt"] = 1.7876926807805918e+02 -all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1.9725084341766897e+03 -all_best_mse["TRMM_LBA"]["v_mean"] = 816.9725417750085 -all_best_mse["TRMM_LBA"]["u_mean"] = 996.5804518730336 -all_best_mse["TRMM_LBA"]["tke_mean"] = 8.2693137314226078e+03 -all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0011323111668659397 -all_best_mse["TRMM_LBA"]["ql_mean"] = 2.6235652901417766e+03 -all_best_mse["TRMM_LBA"]["thetal_mean"] = 1.0479371399644480e-02 -all_best_mse["TRMM_LBA"]["Hvar_mean"] = 2.5385154554204787e+03 -all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2460.6541168787 +all_best_mse["TRMM_LBA"]["qt_mean"] = 3.500120182877692 +all_best_mse["TRMM_LBA"]["updraft_area"] = 7558.43014159832 +all_best_mse["TRMM_LBA"]["updraft_w"] = 10721.449681720187 +all_best_mse["TRMM_LBA"]["updraft_qt"] = 178.88915354774463 +all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1972.3937641987523 +all_best_mse["TRMM_LBA"]["v_mean"] = 817.1073208856625 +all_best_mse["TRMM_LBA"]["u_mean"] = 996.6088819379991 +all_best_mse["TRMM_LBA"]["tke_mean"] = 8741.324999151844 +all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0010537676216305777 +all_best_mse["TRMM_LBA"]["ql_mean"] = 2580.6352613746744 +all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.01045190259488596 +all_best_mse["TRMM_LBA"]["Hvar_mean"] = 2593.6325128090607 +all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2437.826779308088 # all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict() -all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.9052591445508149e+00 -all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7162360331971493e+00 -all_best_mse["LES_driven_SCM"]["u_mean"] = 1.2399488159068501e+00 -all_best_mse["LES_driven_SCM"]["temperature_mean"] = 2.9793961623327896e-03 -all_best_mse["LES_driven_SCM"]["ql_mean"] = 3.9355081949644468e+03 -all_best_mse["LES_driven_SCM"]["thetal_mean"] = 3.2989036089844273e-03 +all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.875464922449264 +all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7100833085198954 +all_best_mse["LES_driven_SCM"]["u_mean"] = 1.2383713052536691 +all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.002969430923084693 +all_best_mse["LES_driven_SCM"]["ql_mean"] = 3459.8779017256375 +all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0032862276838900133 # ################################# ################################# diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl index 79291f0b6e..f4979ea351 100644 --- a/src/EDMF_Updrafts.jl +++ b/src/EDMF_Updrafts.jl @@ -1,4 +1,3 @@ - function initialize_io(up::UpdraftVariables, Stats::NetCDFIO_Stats) add_profile(Stats, "updraft_cloud_fraction") @@ -26,7 +25,7 @@ end function upd_cloud_diagnostics(up::UpdraftVariables, grid, state) up.lwp = 0.0 - prog_up = center_prog_updrafts(state) + aux_up = center_aux_updrafts(state) aux_up = center_aux_updrafts(state) ρ0_c = center_ref_state(state).ρ0 @inbounds for i in 1:(up.n_updrafts) @@ -36,14 +35,14 @@ function upd_cloud_diagnostics(up::UpdraftVariables, grid, state) up.cloud_cover[i] = 0.0 @inbounds for k in real_center_indices(grid) - if prog_up[i].area[k] > 1e-3 + if aux_up[i].area[k] > 1e-3 up.updraft_top[i] = max(up.updraft_top[i], grid.zc[k]) - up.lwp += ρ0_c[k] * aux_up[i].q_liq[k] * prog_up[i].area[k] * grid.Δz + up.lwp += ρ0_c[k] * aux_up[i].q_liq[k] * aux_up[i].area[k] * grid.Δz if aux_up[i].q_liq[k] > 1e-8 up.cloud_base[i] = min(up.cloud_base[i], grid.zc[k]) up.cloud_top[i] = max(up.cloud_top[i], grid.zc[k]) - up.cloud_cover[i] = max(up.cloud_cover[i], prog_up[i].area[k]) + up.cloud_cover[i] = max(up.cloud_cover[i], aux_up[i].area[k]) end end end @@ -65,7 +64,6 @@ function compute_rain_formation_tendencies( ) p0_c = center_ref_state(state).p0 ρ0_c = center_ref_state(state).ρ0 - prog_up = center_prog_updrafts(state) aux_up = center_aux_updrafts(state) prog_ra = center_prog_rain(state) tendencies_ra = center_tendencies_rain(state) @@ -73,7 +71,7 @@ function compute_rain_formation_tendencies( @inbounds for i in 1:(up.n_updrafts) @inbounds for k in real_center_indices(grid) T_up = aux_up[i].T[k] - q_tot_up = prog_up[i].q_tot[k] + q_tot_up = aux_up[i].q_tot[k] ts_up = TD.PhaseEquil_pTq(param_set, p0_c[k], T_up, q_tot_up) # autoconversion and accretion @@ -81,13 +79,13 @@ function compute_rain_formation_tendencies( param_set, rain.rain_model, prog_ra.qr[k], - prog_up[i].area[k], + aux_up[i].area[k], ρ0_c[k], dt, ts_up, ) - up_thermo.qt_tendency_rain_formation[i, k] = mph.qt_tendency * prog_up[i].area[k] - up_thermo.θ_liq_ice_tendency_rain_formation[i, k] = mph.θ_liq_ice_tendency * prog_up[i].area[k] + up_thermo.qt_tendency_rain_formation[i, k] = mph.qt_tendency * aux_up[i].area[k] + up_thermo.θ_liq_ice_tendency_rain_formation[i, k] = mph.θ_liq_ice_tendency * aux_up[i].area[k] end end # TODO - to be deleted once we sum all tendencies elsewhere diff --git a/src/Operators.jl b/src/Operators.jl index eb7f90d904..444df4e3fd 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -576,9 +576,9 @@ function construct_tridiag_diffusion_en( ρ0_f = face_ref_state(state).ρ0 aux_tc = center_aux_tc(state) w_en = face_aux_environment(state).w - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) + aux_up_f = face_aux_updrafts(state) prog_en = center_prog_environment(state) + aux_up = center_aux_updrafts(state) ae = 1 .- aux_tc.bulk.area rho_ae_K_m = face_field(grid) @@ -602,11 +602,11 @@ function construct_tridiag_diffusion_en( D_env = 0.0 @inbounds for i in 1:n_updrafts - if prog_up[i].area[k] > minimum_area + if aux_up[i].area[k] > minimum_area turb_entr = frac_turb_entr[i, k] R_up = pressure_plume_spacing[i] - w_up_c = interpf2c(prog_up_f[i].w, grid, k) - D_env += ρ0_c[k] * prog_up[i].area[k] * w_up_c * (entr_sc[i, k] + turb_entr) + w_up_c = interpf2c(aux_up_f[i].w, grid, k) + D_env += ρ0_c[k] * aux_up[i].area[k] * w_up_c * (entr_sc[i, k] + turb_entr) else D_env = 0.0 end diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index d10e3db3fa..d546853225 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -64,8 +64,8 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti io(edmf.EnvVar, grid, state, Stats) io(edmf.Rain, grid, state, Stats, edmf.UpdThermo, edmf.EnvThermo, TS) - prog_up = center_prog_updrafts(state) aux_tc = center_aux_tc(state) + aux_up = center_aux_updrafts(state) a_up_bulk = aux_tc.bulk.area write_ts(Stats, "rd", StatsBase.mean(edmf.pressure_plume_spacing)) @@ -74,13 +74,13 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti if a_up_bulk[k] > 0.0 @inbounds for i in 1:(edmf.n_updrafts) massflux[k] += interpf2c(edmf.m, grid, k, i) - mean_entr_sc[k] += prog_up[i].area[k] * edmf.entr_sc[i, k] / a_up_bulk[k] - mean_detr_sc[k] += prog_up[i].area[k] * edmf.detr_sc[i, k] / a_up_bulk[k] - mean_asp_ratio[k] += prog_up[i].area[k] * edmf.asp_ratio[i, k] / a_up_bulk[k] - mean_frac_turb_entr[k] += prog_up[i].area[k] * edmf.frac_turb_entr[i, k] / a_up_bulk[k] - mean_horiz_K_eddy[k] += prog_up[i].area[k] * edmf.horiz_K_eddy[i, k] / a_up_bulk[k] - mean_sorting_function[k] += prog_up[i].area[k] * edmf.sorting_function[i, k] / a_up_bulk[k] - mean_b_mix[k] += prog_up[i].area[k] * edmf.b_mix[i, k] / a_up_bulk[k] + mean_entr_sc[k] += aux_up[i].area[k] * edmf.entr_sc[i, k] / a_up_bulk[k] + mean_detr_sc[k] += aux_up[i].area[k] * edmf.detr_sc[i, k] / a_up_bulk[k] + mean_asp_ratio[k] += aux_up[i].area[k] * edmf.asp_ratio[i, k] / a_up_bulk[k] + mean_frac_turb_entr[k] += aux_up[i].area[k] * edmf.frac_turb_entr[i, k] / a_up_bulk[k] + mean_horiz_K_eddy[k] += aux_up[i].area[k] * edmf.horiz_K_eddy[i, k] / a_up_bulk[k] + mean_sorting_function[k] += aux_up[i].area[k] * edmf.sorting_function[i, k] / a_up_bulk[k] + mean_b_mix[k] += aux_up[i].area[k] * edmf.b_mix[i, k] / a_up_bulk[k] end end end @@ -91,7 +91,7 @@ function io(edmf::EDMF_PrognosticTKE, grid, state, Stats::NetCDFIO_Stats, TS::Ti if a_up_bulk_f > 0.0 @inbounds for i in 1:(edmf.n_updrafts) a_up_f = interpc2f( - prog_up[i].area, + aux_up[i].area, grid, k; bottom = SetValue(edmf.area_surface_bc[i]), @@ -155,7 +155,6 @@ function update_radiation end function update_cloud_frac(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables) # update grid-mean cloud fraction and cloud cover - prog_up = center_prog_updrafts(state) aux_tc = center_aux_tc(state) aux_gm = center_aux_grid_mean(state) aux_en = center_aux_environment(state) @@ -178,6 +177,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm, prog_gm = center_prog_grid_mean(state) aux_en = center_aux_environment(state) aux_en_f = face_aux_environment(state) + aux_up = center_aux_updrafts(state) ρ0_f = face_ref_state(state).ρ0 p0_c = center_ref_state(state).p0 α0_c = center_ref_state(state).α0 @@ -254,8 +254,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm, edmf.RainPhys.θ_liq_ice_tendency_rain_evap[k] end - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) + aux_up_f = face_aux_updrafts(state) edmf.massflux_h .= 0.0 edmf.massflux_qt .= 0.0 # Compute the mass flux and associated scalar fluxes @@ -263,9 +262,9 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm, edmf.m[i, kf_surf] = 0.0 a_up_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetZeroGradient()) @inbounds for k in real_face_indices(grid) - a_up = interpc2f(prog_up[i].area, grid, k; a_up_bcs...) + a_up = interpc2f(aux_up[i].area, grid, k; a_up_bcs...) a_en = interpc2f(ae, grid, k; a_up_bcs...) - edmf.m[i, k] = ρ0_f[k] * a_up * a_en * (prog_up_f[i].w[k] - aux_en_f.w[k]) + edmf.m[i, k] = ρ0_f[k] * a_up * a_en * (aux_up_f[i].w[k] - aux_en_f.w[k]) end end @@ -278,8 +277,8 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm, h_en_f = interpc2f(aux_en.θ_liq_ice, grid, k; m_bcs...) qt_en_f = interpc2f(aux_en.q_tot, grid, k; m_bcs...) @inbounds for i in 1:(up.n_updrafts) - h_up_f = interpc2f(prog_up[i].θ_liq_ice, grid, k; m_bcs...) - qt_up_f = interpc2f(prog_up[i].q_tot, grid, k; m_bcs...) + h_up_f = interpc2f(aux_up[i].θ_liq_ice, grid, k; m_bcs...) + qt_up_f = interpc2f(aux_up[i].q_tot, grid, k; m_bcs...) edmf.massflux_h[k] += edmf.m[i, k] * (h_up_f - h_en_f) edmf.massflux_qt[k] += edmf.m[i, k] * (qt_up_f - qt_en_f) end @@ -394,7 +393,6 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca n_updrafts = up.n_updrafts prog_gm = center_prog_grid_mean(state) prog_en = center_prog_environment(state) - prog_up_f = face_prog_updrafts(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 @@ -423,7 +421,6 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca # ----------- TODO: move to compute_tendencies implicit_eqs = edmf.implicit_eqs # 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) common_args = ( grid, @@ -453,6 +450,7 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca ### update ### update_updraft(edmf, grid, state, gm, TS) + filter_updraft_vars(edmf, grid, state, gm) if edmf.Rain.rain_model == "clima_1m" update_rain(edmf.Rain, grid, state, up_thermo, en_thermo, edmf.RainPhys, TS) end @@ -471,6 +469,7 @@ function update(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables, Ca ### ### set values ### + compute_upd_primitives(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) @@ -548,11 +547,11 @@ function get_GMV_CoVar(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, 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) + aux_up_f = face_aux_updrafts(state) aux_en_c = center_aux_environment(state) aux_en_f = face_aux_environment(state) 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) prog_gm = is_tke ? prog_gm_f : prog_gm_c @@ -574,8 +573,8 @@ function get_GMV_CoVar(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, 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) + ϕ_up_var = getproperty(aux_up_f[i], ϕ_sym) + ψ_up_var = getproperty(aux_up_f[i], ψ_sym) ϕ_up_dual = dual_faces(ϕ_up_var, grid, k) ϕ_gm_dual = dual_faces(ϕ_gm, grid, k) ψ_up_dual = dual_faces(ψ_up_var, grid, k) @@ -584,7 +583,7 @@ function get_GMV_CoVar(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, Δψ_dual = ψ_up_dual .- ψ_gm_dual Δϕ = interpf2c(Δϕ_dual, grid, k) Δψ = interpf2c(Δψ_dual, grid, k) - gmv_covar[k] += tke_factor * prog_up[i].area[k] * Δϕ * Δψ + gmv_covar[k] += tke_factor * aux_up[i].area[k] * Δϕ * Δψ end end else @@ -595,11 +594,11 @@ function get_GMV_CoVar(edmf::EDMF_PrognosticTKE, grid, state, covar_sym::Symbol, 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) + ϕ_up_var = getproperty(aux_up[i], ϕ_sym) + ψ_up_var = getproperty(aux_up[i], ψ_sym) Δϕ = ϕ_up_var[k] - ϕ_gm[k] Δψ = ψ_up_var[k] - ψ_gm[k] - gmv_covar[k] += tke_factor * prog_up[i].area[k] * Δϕ * Δψ + gmv_covar[k] += tke_factor * aux_up[i].area[k] * Δϕ * Δψ end end end @@ -626,11 +625,10 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G up = edmf.UpdVar up_thermo = edmf.UpdThermo en = edmf.EnvVar - prog_up = center_prog_updrafts(state) aux_up = center_aux_updrafts(state) aux_en = center_aux_environment(state) aux_en_f = face_aux_environment(state) - prog_up_f = face_prog_updrafts(state) + aux_up_f = face_aux_updrafts(state) tendencies_up = center_tendencies_updrafts(state) tendencies_up_f = face_tendencies_updrafts(state) ρ_0_c = center_ref_state(state).ρ0 @@ -639,12 +637,12 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G @inbounds for i in 1:(up.n_updrafts) @inbounds for k in real_center_indices(grid) - tendencies_up[i].area[k] = 0 - tendencies_up[i].θ_liq_ice[k] = 0 - tendencies_up[i].q_tot[k] = 0 + tendencies_up[i].ρarea[k] = 0 + tendencies_up[i].ρaθ_liq_ice[k] = 0 + tendencies_up[i].ρaq_tot[k] = 0 end @inbounds for k in real_face_indices(grid) - tendencies_up_f[i].w[k] = 0 + tendencies_up_f[i].ρaw[k] = 0 end end @@ -657,19 +655,20 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G @inbounds for k in real_center_indices(grid) @inbounds for i in 1:(up.n_updrafts) is_surface_center(grid, k) && continue - w_up_c = interpf2c(prog_up_f[i].w, grid, k) - adv = upwind_advection_area(ρ_0_c, prog_up[i].area, prog_up_f[i].w, grid, k) - - a_up_c = prog_up[i].area[k] - entr_term = a_up_c * w_up_c * (edmf.entr_sc[i, k]) - detr_term = a_up_c * w_up_c * (-edmf.detr_sc[i, k]) - tendencies_up[i].area[k] = adv + entr_term + detr_term - a_up_candidate = a_up_c + Δt * tendencies_up[i].area[k] - if a_up_candidate > au_lim - a_up_div = a_up_c > 0.0 ? a_up_c : au_lim - a_up_candidate = au_lim - edmf.detr_sc[i, k] = (((au_lim - a_up_c) * dti_ - adv - entr_term) / (-a_up_div * w_up_c)) - tendencies_up[i].area[k] = (a_up_candidate - a_up_c) * dti_ + w_up_c = interpf2c(aux_up_f[i].w, grid, k) + adv = upwind_advection_area(ρ_0_c, aux_up[i].area, aux_up_f[i].w, grid, k) + + ρa_up_c = ρ_0_c[k] * aux_up[i].area[k] + entr_term = ρa_up_c * w_up_c * (edmf.entr_sc[i, k]) + detr_term = ρa_up_c * w_up_c * (-edmf.detr_sc[i, k]) + tendencies_up[i].ρarea[k] = (ρ_0_c[k] * adv + entr_term + detr_term) + ρa_up_candidate = ρa_up_c + Δt * tendencies_up[i].ρarea[k] + ρau_lim = ρ_0_c[k] * au_lim + if ρa_up_candidate > ρau_lim + ρa_up_div = ρa_up_c > 0.0 ? ρa_up_c : ρau_lim + ρa_up_candidate = ρau_lim + edmf.detr_sc[i, k] = (((ρau_lim - ρa_up_c) * dti_ - ρ_0_c[k] * adv - entr_term) / (-ρa_up_div * w_up_c)) + tendencies_up[i].ρarea[k] = (ρa_up_candidate - ρa_up_c) * dti_ end end end @@ -679,20 +678,20 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G @inbounds for k in real_center_indices(grid) @inbounds for i in 1:(up.n_updrafts) - w_up_c = interpf2c(prog_up_f[i].w, grid, k) - m_k = (ρ_0_c[k] * prog_up[i].area[k] * w_up_c) + w_up_c = interpf2c(aux_up_f[i].w, grid, k) + m_k = (ρ_0_c[k] * aux_up[i].area[k] * w_up_c) - adv = upwind_advection_scalar(ρ_0_c, prog_up[i].area, prog_up_f[i].w, prog_up[i].θ_liq_ice, grid, k) + adv = upwind_advection_scalar(ρ_0_c, aux_up[i].area, aux_up_f[i].w, aux_up[i].θ_liq_ice, grid, k) entr = entr_w_c[i, k] * aux_en.θ_liq_ice[k] - detr = detr_w_c[i, k] * prog_up[i].θ_liq_ice[k] + detr = detr_w_c[i, k] * aux_up[i].θ_liq_ice[k] rain = ρ_0_c[k] * up_thermo.θ_liq_ice_tendency_rain_formation[i, k] - tendencies_up[i].θ_liq_ice[k] = -adv + m_k * (entr - detr) + rain + tendencies_up[i].ρaθ_liq_ice[k] = -adv + m_k * (entr - detr) + rain - adv = upwind_advection_scalar(ρ_0_c, prog_up[i].area, prog_up_f[i].w, prog_up[i].q_tot, grid, k) + adv = upwind_advection_scalar(ρ_0_c, aux_up[i].area, aux_up_f[i].w, aux_up[i].q_tot, grid, k) entr = entr_w_c[i, k] * aux_en.q_tot[k] - detr = detr_w_c[i, k] * prog_up[i].q_tot[k] + detr = detr_w_c[i, k] * aux_up[i].q_tot[k] rain = ρ_0_c[k] * up_thermo.qt_tendency_rain_formation[i, k] - tendencies_up[i].q_tot[k] = -adv + m_k * (entr - detr) + rain + tendencies_up[i].ρaq_tot[k] = -adv + m_k * (entr - detr) + rain end end @@ -701,17 +700,17 @@ function compute_updraft_tendencies(edmf::EDMF_PrognosticTKE, grid, state, gm::G is_surface_face(grid, k) && continue @inbounds for i in 1:(up.n_updrafts) a_up_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetZeroGradient()) - a_k = interpc2f(prog_up[i].area, grid, k; a_up_bcs...) + a_k = interpc2f(aux_up[i].area, grid, k; a_up_bcs...) # We know that, since W = 0 at z = 0, these BCs should # not matter in the end: entr_w = interpc2f(entr_w_c, grid, k, i; bottom = SetValue(0), top = SetValue(0)) detr_w = interpc2f(detr_w_c, grid, k, i; bottom = SetValue(0), top = SetValue(0)) B_k = interpc2f(aux_up[i].buoy, grid, k; bottom = SetValue(0), top = SetValue(0)) - adv = upwind_advection_velocity(ρ_0_f, prog_up[i].area, prog_up_f[i].w, grid, k; a_up_bcs) - exch = (ρ_0_f[k] * a_k * prog_up_f[i].w[k] * (entr_w * aux_en_f.w[k] - detr_w * prog_up_f[i].w[k])) + adv = upwind_advection_velocity(ρ_0_f, aux_up[i].area, aux_up_f[i].w, grid, k; a_up_bcs) + exch = (ρ_0_f[k] * a_k * aux_up_f[i].w[k] * (entr_w * aux_en_f.w[k] - detr_w * aux_up_f[i].w[k])) buoy = ρ_0_f[k] * a_k * B_k - tendencies_up_f[i].w[k] = -adv + exch + buoy + edmf.nh_pressure[i, k] + tendencies_up_f[i].ρaw[k] = -adv + exch + buoy + edmf.nh_pressure[i, k] end end return @@ -728,97 +727,144 @@ function update_updraft(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVaria up = edmf.UpdVar up_thermo = edmf.UpdThermo en = edmf.EnvVar + aux_up = center_aux_updrafts(state) tendencies_up = center_tendencies_updrafts(state) tendencies_up_f = face_tendencies_updrafts(state) prog_up = center_prog_updrafts(state) prog_gm = center_prog_grid_mean(state) prog_up_f = face_prog_updrafts(state) - aux_up = center_aux_updrafts(state) aux_up_f = face_aux_updrafts(state) ρ_0_c = center_ref_state(state).ρ0 ρ_0_f = face_ref_state(state).ρ0 @inbounds for i in 1:(up.n_updrafts) - aux_up_f[i].w_new[kf_surf] = edmf.w_surface_bc[i] - aux_up[i].area_new[kc_surf] = edmf.area_surface_bc[i] + prog_up[i].ρarea[kc_surf] = ρ_0_c[kc_surf] * edmf.area_surface_bc[i] + aux_up[i].area[kc_surf] = edmf.area_surface_bc[i] + prog_up_f[i].ρaw[kf_surf] = ρ_0_f[kf_surf] * edmf.w_surface_bc[i] + aux_up_f[i].w[kf_surf] = edmf.w_surface_bc[i] end @inbounds for k in real_center_indices(grid) @inbounds for i in 1:(up.n_updrafts) is_surface_center(grid, k) && continue - a_up_c = prog_up[i].area[k] - a_up_candidate = max(a_up_c + Δt * tendencies_up[i].area[k], 0) - aux_up[i].area_new[k] = a_up_candidate + prog_up[i].ρarea[k] += Δt * tendencies_up[i].ρarea[k] + prog_up[i].ρaθ_liq_ice[k] += Δt * tendencies_up[i].ρaθ_liq_ice[k] + prog_up[i].ρaq_tot[k] += Δt * tendencies_up[i].ρaq_tot[k] end end @inbounds for k in real_face_indices(grid) is_surface_face(grid, k) && continue @inbounds for i in 1:(up.n_updrafts) + prog_up_f[i].ρaw[k] += Δt * tendencies_up_f[i].ρaw[k] + end + end + return +end + +function filter_updraft_vars(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables) + param_set = parameter_set(gm) + kc_surf = kc_surface(grid) + kf_surf = kf_surface(grid) + + up = edmf.UpdVar + prog_up = center_prog_updrafts(state) + prog_gm = center_prog_grid_mean(state) + aux_up = center_aux_updrafts(state) + aux_up_f = face_aux_updrafts(state) + prog_up_f = face_prog_updrafts(state) + ρ_0_c = center_ref_state(state).ρ0 + ρ_0_f = face_ref_state(state).ρ0 + + @inbounds for i in 1:(up.n_updrafts) + prog_up[i].ρarea .= max.(prog_up[i].ρarea, 0) + prog_up[i].ρaθ_liq_ice .= max.(prog_up[i].ρaθ_liq_ice, 0) + prog_up[i].ρaq_tot .= max.(prog_up[i].ρaq_tot, 0) + end + + @inbounds for k in real_face_indices(grid) + is_surface_face(grid, k) && continue + @inbounds for i in 1:(up.n_updrafts) + prog_up_f[i].ρaw[k] = max(prog_up_f[i].ρaw[k], 0) a_up_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetZeroGradient()) - anew_k = interpc2f(aux_up[i].area_new, grid, k; a_up_bcs...) - if anew_k >= edmf.minimum_area - a_k = interpc2f(prog_up[i].area, grid, k; a_up_bcs...) - aux_up_f[i].w_new[k] = - (ρ_0_f[k] * a_k * prog_up_f[i].w[k] + Δt * tendencies_up_f[i].w[k]) / (ρ_0_f[k] * anew_k) - - aux_up_f[i].w_new[k] = max(aux_up_f[i].w_new[k], 0) - # TODO: remove aux_up[i].area_new from this loop. - if aux_up_f[i].w_new[k] <= 0.0 - if !(k.i > length(vec(aux_up[i].area_new))) - aux_up[i].area_new[Cent(k.i)] = 0 - end - end - else - aux_up_f[i].w_new[k] = 0 - if !(k.i > length(vec(aux_up[i].area_new))) - aux_up[i].area_new[Cent(k.i)] = 0 - end + ρa_f = interpc2f(prog_up[i].ρarea, grid, k; a_up_bcs...) + if ρa_f < ρ_0_f[k] * edmf.minimum_area + prog_up_f[i].ρaw[k] = 0 end end end + @inbounds for k in real_center_indices(grid) + @inbounds for i in 1:(up.n_updrafts) + is_surface_center(grid, k) && continue + prog_up[i].ρaq_tot[k] = max(prog_up[i].ρaq_tot[k], 0) + ρaw_up_c = interpf2c(prog_up_f[i].ρaw, grid, k) + if ρaw_up_c <= 0 + prog_up[i].ρarea[k] = 0 + prog_up[i].ρaθ_liq_ice[k] = 0 + prog_up[i].ρaq_tot[k] = 0 + end + # this is needed to make sure Rico is unchanged. + # TODO : look into it further to see why + # a similar filtering of ρaθ_liq_ice breaks the simulation + if prog_up[i].ρarea[k] / ρ_0_c[k] < edmf.minimum_area + prog_up[i].ρaq_tot[k] = 0 + end + end + end + return +end + +function compute_upd_primitives(edmf::EDMF_PrognosticTKE, grid, state, gm::GridMeanVariables) + param_set = parameter_set(gm) + kc_surf = kc_surface(grid) + kf_surf = kf_surface(grid) + + up = edmf.UpdVar + prog_up = center_prog_updrafts(state) + prog_gm = center_prog_grid_mean(state) + aux_up = center_aux_updrafts(state) + aux_up_f = face_aux_updrafts(state) + prog_up_f = face_prog_updrafts(state) + ρ_0_c = center_ref_state(state).ρ0 + ρ_0_f = face_ref_state(state).ρ0 + + # set values @inbounds for k in real_center_indices(grid) @inbounds for i in 1:(up.n_updrafts) if is_surface_center(grid, k) # at the surface - if aux_up[i].area_new[k] >= edmf.minimum_area - aux_up[i].θ_liq_ice_new[k] = edmf.h_surface_bc[i] - aux_up[i].q_tot_new[k] = edmf.qt_surface_bc[i] + if prog_up[i].ρarea[k] / ρ_0_c[k] >= edmf.minimum_area + aux_up[i].θ_liq_ice[k] = edmf.h_surface_bc[i] + aux_up[i].q_tot[k] = edmf.qt_surface_bc[i] else - aux_up[i].θ_liq_ice_new[k] = prog_gm.θ_liq_ice[k] - aux_up[i].q_tot_new[k] = prog_gm.q_tot[k] + aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] + aux_up[i].q_tot[k] = prog_gm.q_tot[k] end continue end - - if aux_up[i].area_new[k] >= edmf.minimum_area - aux_up[i].θ_liq_ice_new[k] = - (ρ_0_c[k] * prog_up[i].area[k] * prog_up[i].θ_liq_ice[k] + Δt * tendencies_up[i].θ_liq_ice[k]) / - (ρ_0_c[k] * aux_up[i].area_new[k]) - - aux_up[i].q_tot_new[k] = max( - (ρ_0_c[k] * prog_up[i].area[k] * prog_up[i].q_tot[k] + Δt * tendencies_up[i].q_tot[k]) / - (ρ_0_c[k] * aux_up[i].area_new[k]), - 0.0, - ) - + if prog_up[i].ρarea[k] / ρ_0_c[k] >= edmf.minimum_area + aux_up[i].θ_liq_ice[k] = prog_up[i].ρaθ_liq_ice[k] / prog_up[i].ρarea[k] + aux_up[i].q_tot[k] = prog_up[i].ρaq_tot[k] / prog_up[i].ρarea[k] + aux_up[i].area[k] = prog_up[i].ρarea[k] / ρ_0_c[k] else - aux_up[i].θ_liq_ice_new[k] = prog_gm.θ_liq_ice[k] - aux_up[i].q_tot_new[k] = prog_gm.q_tot[k] + aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] + aux_up[i].q_tot[k] = prog_gm.q_tot[k] + aux_up[i].area[k] = 0 end end end - # set values - @inbounds for i in 1:(up.n_updrafts) - @inbounds for k in real_face_indices(grid) - prog_up_f[i].w[k] = aux_up_f[i].w_new[k] - end - @inbounds for k in real_center_indices(grid) - prog_up[i].area[k] = aux_up[i].area_new[k] - prog_up[i].q_tot[k] = aux_up[i].q_tot_new[k] - prog_up[i].θ_liq_ice[k] = aux_up[i].θ_liq_ice_new[k] + @inbounds for k in real_face_indices(grid) + is_surface_face(grid, k) && continue + @inbounds for i in 1:(up.n_updrafts) + a_up_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetZeroGradient()) + anew_k = interpc2f(aux_up[i].area, grid, k; a_up_bcs...) + if anew_k >= edmf.minimum_area + aux_up_f[i].w[k] = max(prog_up_f[i].ρaw[k] / (ρ_0_f[k] * anew_k), 0) + else + aux_up_f[i].w[k] = 0 + end end end return @@ -893,11 +939,11 @@ function compute_covariance_interdomain_src( 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_up = center_aux_updrafts(state) + aux_up_f = face_aux_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 + prog_up = is_tke ? aux_up_f : aux_up aux_en_c = center_aux_environment(state) aux_en_f = face_aux_environment(state) aux_en = is_tke ? aux_en_f : aux_en_c @@ -913,7 +959,7 @@ function compute_covariance_interdomain_src( Δϕ = interpf2c(ϕ_up, grid, k) - interpf2c(ϕ_en, grid, k) Δψ = interpf2c(ψ_up, grid, k) - interpf2c(ψ_en, grid, k) - aux_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 * aux_up[i].area[k] * (1.0 - aux_up[i].area[k]) * Δϕ * Δψ end end else @@ -924,7 +970,7 @@ function compute_covariance_interdomain_src( ψ_up = getproperty(prog_up[i], ψ_var) Δϕ = ϕ_up[k] - ϕ_en[k] Δψ = ψ_up[k] - ψ_en[k] - aux_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 * aux_up[i].area[k] * (1.0 - aux_up[i].area[k]) * Δϕ * Δψ end end end @@ -944,12 +990,12 @@ function compute_covariance_entr( 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_up = center_aux_updrafts(state) + aux_up_f = face_aux_updrafts(state) prog_gm_c = center_prog_grid_mean(state) prog_gm_f = face_prog_grid_mean(state) prog_gm = is_tke ? prog_gm_f : prog_gm_c - prog_up = is_tke ? prog_up_f : prog_up_c + prog_up = is_tke ? aux_up_f : aux_up GmvVar1 = getproperty(prog_gm, var1) GmvVar2 = getproperty(prog_gm, var2) aux_en_2m = center_aux_environment_2m(state) @@ -966,7 +1012,7 @@ function compute_covariance_entr( 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] + a_up = aux_up[i].area[k] if a_up > edmf.minimum_area R_up = edmf.pressure_plume_spacing[i] up_var_1 = getproperty(prog_up[i], var1) @@ -980,7 +1026,7 @@ function compute_covariance_entr( eps_turb = edmf.frac_turb_entr[i, k] - w_u = interpf2c(prog_up_f[i].w, grid, k) + w_u = interpf2c(aux_up_f[i].w, grid, k) dynamic_entr = tke_factor * ρ_0_c[k] * @@ -1009,19 +1055,18 @@ end 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_up_f = face_aux_updrafts(state) 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) @inbounds for k in real_center_indices(grid) 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) - aux_covar.detr_loss[k] += prog_up[i].area[k] * abs(w_up_c) * edmf.entr_sc[i, k] + 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] end @@ -1051,16 +1096,16 @@ function en_diffusion_tendencies(grid::Grid, state, TS, covar_sym::Symbol, n_upd 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) + aux_up = center_aux_updrafts(state) ae = center_field(grid) @inbounds for k in real_center_indices(grid) - ae[k] = 1 .- sum(ntuple(i -> prog_up[i].area[k], n_updrafts)) + ae[k] = 1 .- sum(ntuple(i -> aux_up[i].area[k], n_updrafts)) end kc_surf = kc_surface(grid) @@ -1092,22 +1137,22 @@ function GMV_third_m(edmf::EDMF_PrognosticTKE, grid, state, covar_en_sym::Symbol en = edmf.EnvVar aux_tc = center_aux_tc(state) ae = 1 .- aux_tc.bulk.area - prog_up = center_prog_updrafts(state) - prog_up_f = face_prog_updrafts(state) + aux_up_f = face_aux_updrafts(state) is_tke = covar_en_sym == :tke covar_en = getproperty(center_prog_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 + aux_up = center_aux_updrafts(state) var_en = getproperty(aux_en, var) @inbounds for k in real_center_indices(grid) mean_en = is_tke ? interpf2c(var_en, grid, k) : var_en[k] GMVv_ = ae[k] * mean_en @inbounds for i in 1:(up.n_updrafts) - var_up = is_tke ? getproperty(prog_up_f[i], var) : getproperty(prog_up[i], var) + var_up = is_tke ? getproperty(aux_up_f[i], var) : getproperty(aux_up[i], var) mean_up = is_tke ? interpf2c(var_up, grid, k) : var_up[k] - GMVv_ += prog_up[i].area[k] * mean_up + GMVv_ += aux_up[i].area[k] * mean_up end # TODO: report bug: i used outside of scope. @@ -1126,10 +1171,10 @@ function GMV_third_m(edmf::EDMF_PrognosticTKE, grid, state, covar_en_sym::Symbol Upd_cubed = 0.0 GMVcov_ = ae[k] * (Envcov_ + (mean_en - GMVv_)^2) @inbounds for i in 1:(up.n_updrafts) - var_up = is_tke ? getproperty(prog_up_f[i], var) : getproperty(prog_up[i], var) + var_up = is_tke ? getproperty(aux_up_f[i], var) : getproperty(aux_up[i], var) mean_up = is_tke ? interpf2c(var_up, grid, k) : var_up[k] - GMVcov_ += prog_up[i].area[k] * (mean_up - GMVv_)^2 - Upd_cubed += prog_up[i].area[k] * mean_up^3 + GMVcov_ += aux_up[i].area[k] * (mean_up - GMVv_)^2 + Upd_cubed += aux_up[i].area[k] * mean_up^3 end if is_surface_center(grid, k) diff --git a/src/types.jl b/src/types.jl index c5045ec378..f3620db305 100644 --- a/src/types.jl +++ b/src/types.jl @@ -219,7 +219,6 @@ end mutable struct UpdraftVariables{A1} n_updrafts::Int - prognostic::Bool updraft_fraction::Float64 cloud_fraction::A1 cloud_base::A1 @@ -230,7 +229,6 @@ mutable struct UpdraftVariables{A1} function UpdraftVariables(nu, namelist, grid::Grid) n_updrafts = nu - prognostic = true updraft_fraction = namelist["turbulence"]["EDMF_PrognosticTKE"]["surface_area"] # cloud and rain diagnostics for output @@ -245,7 +243,6 @@ mutable struct UpdraftVariables{A1} A1 = typeof(cloud_fraction) return new{A1}( n_updrafts, - prognostic, updraft_fraction, cloud_fraction, cloud_base, diff --git a/src/update_aux.jl b/src/update_aux.jl index 48f4b3420d..aa0be23553 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -18,21 +18,20 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) surface = Case.Sur obukhov_length = surface.obukhov_length FT = eltype(grid) - prog_up = center_prog_updrafts(state) prog_gm = center_prog_grid_mean(state) prog_gm_f = face_prog_grid_mean(state) aux_up = center_aux_updrafts(state) + aux_up_f = face_aux_updrafts(state) aux_en = center_aux_environment(state) aux_en_f = face_aux_environment(state) aux_gm = center_aux_grid_mean(state) - prog_up_f = face_prog_updrafts(state) aux_tc_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)) + aux_tc.bulk.area[k] = sum(ntuple(i -> aux_up[i].area[k], up.n_updrafts)) end ##### @@ -60,8 +59,8 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) if a_bulk_f > 1.0e-20 @inbounds for i in 1:(up.n_updrafts) a_up_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetZeroGradient()) - a_up_f = interpc2f(prog_up[i].area, grid, k; a_up_bcs...) - aux_tc_f.bulk.w[k] += a_up_f * prog_up_f[i].w[k] / a_bulk_f + a_up_f = interpc2f(aux_up[i].area, grid, k; a_up_bcs...) + aux_tc_f.bulk.w[k] += a_up_f * aux_up_f[i].w[k] / a_bulk_f end end # Assuming gm.W = 0! @@ -78,12 +77,12 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) aux_tc.bulk.buoy[k] = 0 if a_bulk_c > 1.0e-20 @inbounds for i in 1:(up.n_updrafts) - aux_tc.bulk.q_tot[k] += prog_up[i].area[k] * prog_up[i].q_tot[k] / a_bulk_c - aux_tc.bulk.q_liq[k] += prog_up[i].area[k] * aux_up[i].q_liq[k] / a_bulk_c - aux_tc.bulk.θ_liq_ice[k] += prog_up[i].area[k] * prog_up[i].θ_liq_ice[k] / a_bulk_c - aux_tc.bulk.T[k] += prog_up[i].area[k] * aux_up[i].T[k] / a_bulk_c - aux_tc.bulk.RH[k] += prog_up[i].area[k] * aux_up[i].RH[k] / a_bulk_c - aux_tc.bulk.buoy[k] += prog_up[i].area[k] * aux_up[i].buoy[k] / a_bulk_c + aux_tc.bulk.q_tot[k] += aux_up[i].area[k] * aux_up[i].q_tot[k] / a_bulk_c + aux_tc.bulk.q_liq[k] += aux_up[i].area[k] * aux_up[i].q_liq[k] / a_bulk_c + aux_tc.bulk.θ_liq_ice[k] += aux_up[i].area[k] * aux_up[i].θ_liq_ice[k] / a_bulk_c + aux_tc.bulk.T[k] += aux_up[i].area[k] * aux_up[i].T[k] / a_bulk_c + aux_tc.bulk.RH[k] += aux_up[i].area[k] * aux_up[i].RH[k] / a_bulk_c + aux_tc.bulk.buoy[k] += aux_up[i].area[k] * aux_up[i].buoy[k] / a_bulk_c end else aux_tc.bulk.q_tot[k] = prog_gm.q_tot[k] @@ -123,17 +122,17 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ##### @inbounds for i in 1:(up.n_updrafts) - if prog_up[i].area[k] > 0.0 - ts_up = thermo_state_pθq(param_set, p0_c[k], prog_up[i].θ_liq_ice[k], prog_up[i].q_tot[k]) + if aux_up[i].area[k] > 0.0 + ts_up = thermo_state_pθq(param_set, p0_c[k], aux_up[i].θ_liq_ice[k], aux_up[i].q_tot[k]) aux_up[i].q_liq[k] = TD.liquid_specific_humidity(ts_up) aux_up[i].T[k] = TD.air_temperature(ts_up) ρ = TD.air_density(ts_up) aux_up[i].buoy[k] = buoyancy_c(param_set, ρ0_c[k], ρ) aux_up[i].RH[k] = TD.relative_humidity(ts_up) elseif k > kc_surf - if prog_up[i].area[k - 1] > 0.0 && edmf.extrapolate_buoyancy - qt = prog_up[i].q_tot[k - 1] - h = prog_up[i].θ_liq_ice[k - 1] + if aux_up[i].area[k - 1] > 0.0 && edmf.extrapolate_buoyancy + qt = aux_up[i].q_tot[k - 1] + h = aux_up[i].θ_liq_ice[k - 1] ts_up = thermo_state_pθq(param_set, p0_c[k], h, qt) ρ = TD.air_density(ts_up) aux_up[i].buoy[k] = buoyancy_c(param_set, ρ0_c[k], ρ) @@ -150,7 +149,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) aux_gm.buoy[k] = (1.0 - aux_tc.bulk.area[k]) * aux_en.buoy[k] @inbounds for i in 1:(up.n_updrafts) - aux_gm.buoy[k] += prog_up[i].area[k] * aux_up[i].buoy[k] + aux_gm.buoy[k] += aux_up[i].area[k] * aux_up[i].buoy[k] end @inbounds for i in 1:(up.n_updrafts) aux_up[i].buoy[k] -= aux_gm.buoy[k] @@ -188,14 +187,14 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) @inbounds for k in real_center_indices(grid) @inbounds for i in 1:(up.n_updrafts) # entrainment - if prog_up[i].area[k] > 0.0 + if aux_up[i].area[k] > 0.0 # compute ∇m at cell centers - a_up_c = prog_up[i].area[k] - w_up_c = interpf2c(prog_up_f[i].w, grid, k) + a_up_c = aux_up[i].area[k] + w_up_c = interpf2c(aux_up_f[i].w, grid, k) w_gm_c = interpf2c(prog_gm_f.w, grid, k) m = a_up_c * (w_up_c - w_gm_c) - a_up_cut = ccut_upwind(prog_up[i].area, grid, k) - w_up_cut = daul_f2c_upwind(prog_up_f[i].w, grid, k) + a_up_cut = ccut_upwind(aux_up[i].area, grid, k) + w_up_cut = daul_f2c_upwind(aux_up_f[i].w, grid, k) w_gm_cut = daul_f2c_upwind(prog_gm_f.w, grid, k) m_cut = a_up_cut .* (w_up_cut .- w_gm_cut) ∇m = FT(c∇_upwind(m_cut, grid, k; bottom = SetValue(0), top = FreeBoundary())) @@ -205,14 +204,14 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) εδ_model = MoistureDeficitEntr(; q_liq_up = aux_up[i].q_liq[k], q_liq_en = aux_en.q_liq[k], - w_up = interpf2c(prog_up_f[i].w, grid, k), + w_up = interpf2c(aux_up_f[i].w, grid, k), 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], dMdz = ∇m, M = m, - a_up = prog_up[i].area[k], + a_up = aux_up[i].area[k], a_en = aux_en.area[k], R_up = edmf.pressure_plume_spacing[i], RH_up = aux_up[i].RH[k], @@ -245,12 +244,12 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) # pressure a_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetValue(0)) - a_kfull = interpc2f(prog_up[i].area, grid, k; a_bcs...) + a_kfull = interpc2f(aux_up[i].area, grid, k; a_bcs...) if a_kfull > 0.0 B = aux_up[i].buoy b_bcs = (; bottom = SetValue(B[kc_surf]), top = SetValue(B[kc_toa])) b_kfull = interpc2f(aux_up[i].buoy, grid, k; b_bcs...) - w_cut = fcut(prog_up_f[i].w, grid, k) + w_cut = fcut(aux_up_f[i].w, grid, k) ∇w_up = f∇(w_cut, grid, k; bottom = SetValue(0), top = SetGradient(0)) asp_ratio = 1.0 nh_pressure_b, nh_pressure_adv, nh_pressure_drag = perturbation_pressure( @@ -259,7 +258,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) a_kfull, b_kfull, ρ0_f[k], - prog_up_f[i].w[k], + aux_up_f[i].w[k], ∇w_up, aux_en_f.w[k], asp_ratio, @@ -287,7 +286,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) V_cut = ccut(prog_gm.v, grid, k) wc_en = interpf2c(aux_en_f.w, grid, k) wc_up = ntuple(up.n_updrafts) do i - interpf2c(prog_up_f[i].w, grid, k) + interpf2c(aux_up_f[i].w, grid, k) end w_dual = dual_faces(aux_en_f.w, grid, k) @@ -401,7 +400,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) a_en = (1 - aux_tc.bulk.area[k]), wc_en = wc_en, wc_up = Tuple(wc_up), - a_up = ntuple(i -> prog_up[i].area[k], up.n_updrafts), + a_up = ntuple(i -> aux_up[i].area[k], up.n_updrafts), ε_turb = ntuple(i -> edmf.frac_turb_entr[i, k], up.n_updrafts), δ_dyn = ntuple(i -> edmf.detr_sc[i, k], up.n_updrafts), N_up = up.n_updrafts, @@ -428,7 +427,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) @inbounds for k in real_center_indices(grid) 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_up_c = interpf2c(aux_up_f[i].w, grid, k) w_en_c = interpf2c(aux_en_f.w, grid, k) press_c = interpf2c(edmf.nh_pressure, grid, k, i) aux_en_2m.tke.press[k] += (w_en_c - w_up_c) * press_c