diff --git a/integration_tests/utils/mse_tables.jl b/integration_tests/utils/mse_tables.jl index ea5eac7df3..06cfccb4a3 100644 --- a/integration_tests/utils/mse_tables.jl +++ b/integration_tests/utils/mse_tables.jl @@ -16,8 +16,8 @@ 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.0001008500529941707 -all_best_mse["ARM_SGP"]["Hvar_mean"] = 10272.271366887158 -all_best_mse["ARM_SGP"]["QTvar_mean"] = 6520.8612887442105 +all_best_mse["ARM_SGP"]["Hvar_mean"] = 10272.025247957436 +all_best_mse["ARM_SGP"]["QTvar_mean"] = 6520.655602593681 # all_best_mse["Bomex"] = OrderedCollections.OrderedDict() all_best_mse["Bomex"]["qt_mean"] = 0.11089388922127488 @@ -32,44 +32,44 @@ 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.4643405606181006e-5 -all_best_mse["Bomex"]["Hvar_mean"] = 4145.770313701388 -all_best_mse["Bomex"]["QTvar_mean"] = 1546.822975591833 +all_best_mse["Bomex"]["Hvar_mean"] = 4145.578933520064 +all_best_mse["Bomex"]["QTvar_mean"] = 1546.7494339261843 # all_best_mse["DryBubble"] = OrderedCollections.OrderedDict() -all_best_mse["DryBubble"]["updraft_area"] = 5.342733081510248e-21 -all_best_mse["DryBubble"]["updraft_w"] = 1.504306725196539e-21 -all_best_mse["DryBubble"]["updraft_thetal"] = 2.085357197131444e-28 +all_best_mse["DryBubble"]["updraft_area"] = 0.0 +all_best_mse["DryBubble"]["updraft_w"] = 0.0 +all_best_mse["DryBubble"]["updraft_thetal"] = 0.0 all_best_mse["DryBubble"]["u_mean"] = 0.0 -all_best_mse["DryBubble"]["tke_mean"] = 1.447582969419808e-21 -all_best_mse["DryBubble"]["temperature_mean"] = 5.051124771997088e-29 -all_best_mse["DryBubble"]["thetal_mean"] = 3.8122194338679874e-29 -all_best_mse["DryBubble"]["Hvar_mean"] = 8.062691041323753e-31 +all_best_mse["DryBubble"]["tke_mean"] = 0.0 +all_best_mse["DryBubble"]["temperature_mean"] = 0.0 +all_best_mse["DryBubble"]["thetal_mean"] = 0.0 +all_best_mse["DryBubble"]["Hvar_mean"] = 3.6376577831233795e-14 # all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict() -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["DYCOMS_RF01"]["qt_mean"] = 0.02266768949492931 +all_best_mse["DYCOMS_RF01"]["ql_mean"] = 10.279853226025203 +all_best_mse["DYCOMS_RF01"]["updraft_area"] = 30.546343419259646 +all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.13860507148924 +all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.151491944788959 +all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18691270440001 +all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.003007505042897622 +all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.07909085313193075 +all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.33883567622849 +all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 8.085863436032554e-5 +all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 8.15893797083448e-5 +all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1243.0781411983992 +all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 483.9134520366122 # all_best_mse["GABLS"] = OrderedCollections.OrderedDict() all_best_mse["GABLS"]["updraft_thetal"] = 0.0 -all_best_mse["GABLS"]["v_mean"] = 3.714897400978881e-61 +all_best_mse["GABLS"]["v_mean"] = 0.0 all_best_mse["GABLS"]["u_mean"] = 0.0 -all_best_mse["GABLS"]["tke_mean"] = 1.4676586011780295e-78 +all_best_mse["GABLS"]["tke_mean"] = 0.0 all_best_mse["GABLS"]["temperature_mean"] = 0.0 all_best_mse["GABLS"]["thetal_mean"] = 0.0 -all_best_mse["GABLS"]["Hvar_mean"] = 1.0682670526253067e-29 -all_best_mse["GABLS"]["QTvar_mean"] = 1.638999269581565e-63 -all_best_mse["GABLS"]["qt_mean"] = 1.4111973896803769e-62 +all_best_mse["GABLS"]["Hvar_mean"] = 6.622785430136218e-9 +all_best_mse["GABLS"]["QTvar_mean"] = 2.0949520033688362e-11 +all_best_mse["GABLS"]["qt_mean"] = 0.0 # all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict() all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0 @@ -83,8 +83,8 @@ all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.0 all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.0 all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 0.0 all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 8.398066020086122e-38 -all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 1.7777704758595357e-40 +all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 5.122332880474674e-7 +all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 5.204566313335327e-7 # all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict() all_best_mse["Nieuwstadt"]["updraft_area"] = 98.97436274691259 @@ -94,7 +94,7 @@ 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.591719277685 +all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1277.483163122067 # all_best_mse["Rico"] = OrderedCollections.OrderedDict() all_best_mse["Rico"]["qt_mean"] = 1.2392247687753224 @@ -110,8 +110,8 @@ all_best_mse["Rico"]["ql_mean"] = 63.36370009749748 all_best_mse["Rico"]["qi_mean"] = "NA" all_best_mse["Rico"]["qr_mean"] = 760.8095424564954 all_best_mse["Rico"]["thetal_mean"] = 0.0005657435436849961 -all_best_mse["Rico"]["Hvar_mean"] = 179284.00789323114 -all_best_mse["Rico"]["QTvar_mean"] = 40234.39819335395 +all_best_mse["Rico"]["Hvar_mean"] = 179281.32561909297 +all_best_mse["Rico"]["QTvar_mean"] = 40233.81857158096 # all_best_mse["Soares"] = OrderedCollections.OrderedDict() all_best_mse["Soares"]["qt_mean"] = 0.12623520584762954 @@ -123,7 +123,7 @@ 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["Soares"]["Hvar_mean"] = 1121.220961389257 # all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict() all_best_mse["TRMM_LBA"]["qt_mean"] = 3.963977389784696 @@ -138,8 +138,8 @@ all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.00112576565688425 all_best_mse["TRMM_LBA"]["ql_mean"] = 11335.71595959235 all_best_mse["TRMM_LBA"]["qi_mean"] = "NA" all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.021118882573307667 -all_best_mse["TRMM_LBA"]["Hvar_mean"] = 4344.077907455965 -all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2835.250376882502 +all_best_mse["TRMM_LBA"]["Hvar_mean"] = 4343.288145343326 +all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2834.894441259608 # all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict() all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.6848941301202665 diff --git a/src/update_aux.jl b/src/update_aux.jl index 2388cf5697..00527addcd 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -34,102 +34,87 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) prog_up_f = face_prog_updrafts(state) ##### - ##### Set primitive variables + ##### center variables ##### - @inbounds for i in 1:(up.n_updrafts) - - # at the surface - if prog_up[i].ρarea[kc_surf] / ρ0_c[kc_surf] >= edmf.minimum_area - aux_up[i].θ_liq_ice[kc_surf] = edmf.h_surface_bc[i] - aux_up[i].q_tot[kc_surf] = edmf.qt_surface_bc[i] - aux_up[i].area[kc_surf] = edmf.area_surface_bc[i] - aux_up_f[i].w[kf_surf] = edmf.w_surface_bc[i] - else - aux_up[i].θ_liq_ice[kc_surf] = prog_gm.θ_liq_ice[kc_surf] - aux_up[i].q_tot[kc_surf] = prog_gm.q_tot[kc_surf] - end - - @inbounds for k in real_center_indices(grid) - is_surface_center(grid, k) && continue - 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[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 - - @inbounds for k in real_face_indices(grid) - is_surface_face(grid, k) && continue + @inbounds for k in real_center_indices(grid) + ##### + ##### Set primitive variables + ##### @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) + if is_surface_center(grid, k) + 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] + aux_up[i].area[k] = edmf.area_surface_bc[i] + else + aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k] + aux_up[i].q_tot[k] = prog_gm.q_tot[k] + end else - aux_up_f[i].w[k] = 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[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 - end - @inbounds for k in real_center_indices(grid) + ##### + ##### compute bulk + ##### + aux_tc.bulk.q_tot[k] = 0 + aux_tc.bulk.θ_liq_ice[k] = 0 aux_tc.bulk.area[k] = sum(ntuple(i -> aux_up[i].area[k], up.n_updrafts)) - end - aux_en.area .= 1 .- aux_tc.bulk.area # area of environment - - @inbounds for k in real_center_indices(grid) + if aux_tc.bulk.area[k] > 0 + @inbounds for i in 1:(up.n_updrafts) + aux_tc.bulk.q_tot[k] += aux_up[i].area[k] * aux_up[i].q_tot[k] / aux_tc.bulk.area[k] + aux_tc.bulk.θ_liq_ice[k] += aux_up[i].area[k] * aux_up[i].θ_liq_ice[k] / aux_tc.bulk.area[k] + end + else + aux_tc.bulk.q_tot[k] = prog_gm.q_tot[k] + aux_tc.bulk.θ_liq_ice[k] = prog_gm.θ_liq_ice[k] + end + aux_en.area[k] = 1 - aux_tc.bulk.area[k] aux_en.tke[k] = prog_en.ρatke[k] / (ρ0_c[k] * aux_en.area[k]) aux_en.Hvar[k] = prog_en.ρaHvar[k] / (ρ0_c[k] * aux_en.area[k]) aux_en.QTvar[k] = prog_en.ρaQTvar[k] / (ρ0_c[k] * aux_en.area[k]) aux_en.HQTcov[k] = prog_en.ρaHQTcov[k] / (ρ0_c[k] * aux_en.area[k]) - end - - ##### - ##### diagnose_GMV_moments - ##### - #! format: off - get_GMV_CoVar(edmf, grid, state, :Hvar, :θ_liq_ice) - get_GMV_CoVar(edmf, grid, state, :QTvar, :q_tot) - get_GMV_CoVar(edmf, grid, state, :HQTcov, :θ_liq_ice, :q_tot) - GMV_third_m(edmf, grid, state, :Hvar, :θ_liq_ice, :H_third_m) - GMV_third_m(edmf, grid, state, :QTvar, :q_tot, :QT_third_m) - GMV_third_m(edmf, grid, state, :tke, :w, :W_third_m) - #! format: on - ##### - ##### decompose_environment - ##### - # (Find values of environmental variables by subtracting updraft values from grid mean values. - # Make sure the "bulkvalues" of the updraft variables are updated first) - # velocity (face indicies) - @inbounds for k in real_face_indices(grid) - aux_tc_f.bulk.w[k] = 0 - a_bulk_bcs = (; bottom = SetValue(sum(edmf.area_surface_bc)), top = SetZeroGradient()) - a_bulk_f = interpc2f(aux_tc.bulk.area, grid, k; a_bulk_bcs...) - 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(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! - aux_en_f.w[k] = -a_bulk_f / (1 - a_bulk_f) * aux_tc_f.bulk.w[k] - end - @inbounds for k in real_center_indices(grid) + ##### + ##### decompose_environment + ##### + # (Find values of environmental variables by subtracting updraft values from grid mean values. + # Make sure the "bulkvalues" of the updraft variables are updated first) + # velocity (face indicies) a_bulk_c = aux_tc.bulk.area[k] - aux_tc.bulk.q_tot[k] = 0 - aux_tc.bulk.q_liq[k] = 0 - aux_tc.bulk.q_ice[k] = 0 - aux_tc.bulk.θ_liq_ice[k] = 0 - aux_tc.bulk.T[k] = 0 - aux_tc.bulk.RH[k] = 0 - aux_tc.bulk.buoy[k] = 0 + val1 = 1 / (1 - a_bulk_c) + val2 = a_bulk_c * val1 + + aux_en.area[k] = 1 - a_bulk_c + aux_en.q_tot[k] = max(val1 * prog_gm.q_tot[k] - val2 * aux_tc.bulk.q_tot[k], 0) #Yair - this is here to prevent negative QT + aux_en.θ_liq_ice[k] = val1 * prog_gm.θ_liq_ice[k] - val2 * aux_tc.bulk.θ_liq_ice[k] + + ##### + ##### saturation_adjustment and buoyancy + ##### + ts = thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k]) + aux_tc.θ_virt[k] = TD.virtual_pottemp(ts) + + ts_en = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k]) + + aux_en.T[k] = TD.air_temperature(ts_en) + aux_en.q_liq[k] = TD.liquid_specific_humidity(ts_en) + aux_en.q_ice[k] = TD.ice_specific_humidity(ts_en) + rho = TD.air_density(ts_en) + aux_en.buoy[k] = buoyancy_c(param_set, ρ0_c[k], rho) + + update_sat_unsat(en_thermo, state, k, ts_en) + aux_en.RH[k] = TD.relative_humidity(ts_en) @inbounds for i in 1:(up.n_updrafts) if aux_up[i].area[k] < edmf.minimum_area && k > kc_surf && aux_up[i].area[k - 1] > 0.0 @@ -146,84 +131,105 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) aux_up[i].buoy[k] = buoyancy_c(param_set, ρ0_c[k], ρ) aux_up[i].RH[k] = TD.relative_humidity(ts_up) end + 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] += 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] + end + aux_en.buoy[k] -= aux_gm.buoy[k] - if a_bulk_c > 1.0e-20 + + ##### + ##### compute bulk thermodynamics + ##### + aux_tc.bulk.q_liq[k] = 0 + aux_tc.bulk.q_ice[k] = 0 + aux_tc.bulk.T[k] = 0 + aux_tc.bulk.RH[k] = 0 + aux_tc.bulk.buoy[k] = 0 + if a_bulk_c > 0 @inbounds for i in 1:(up.n_updrafts) - 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.q_ice[k] += aux_up[i].area[k] * aux_up[i].q_ice[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] - aux_tc.bulk.θ_liq_ice[k] = prog_gm.θ_liq_ice[k] - aux_tc.bulk.RH[k] = aux_gm.RH[k] # TODO - here we are using previous timestep values - aux_tc.bulk.T[k] = aux_gm.T[k] # TODO - here we are using previous timestep values + aux_tc.bulk.RH[k] = aux_en.RH[k] + aux_tc.bulk.T[k] = aux_en.T[k] end + + ##### + ##### update_GMV_diagnostics + ##### + aux_gm.q_liq[k] = (aux_tc.bulk.area[k] * aux_tc.bulk.q_liq[k] + (1 - aux_tc.bulk.area[k]) * aux_en.q_liq[k]) + aux_gm.q_ice[k] = (aux_tc.bulk.area[k] * aux_tc.bulk.q_ice[k] + (1 - aux_tc.bulk.area[k]) * aux_en.q_ice[k]) + aux_gm.T[k] = (aux_tc.bulk.area[k] * aux_tc.bulk.T[k] + (1 - aux_tc.bulk.area[k]) * aux_en.T[k]) + aux_gm.buoy[k] = (aux_tc.bulk.area[k] * aux_tc.bulk.buoy[k] + (1 - aux_tc.bulk.area[k]) * aux_en.buoy[k]) + has_condensate = TD.has_condensate(aux_tc.bulk.q_liq[k] + aux_tc.bulk.q_ice[k]) aux_tc.bulk.cloud_fraction[k] = if has_condensate && a_bulk_c > 1e-3 1 else 0 end + end - val1 = 1 / (1 - a_bulk_c) - val2 = a_bulk_c * val1 - - aux_en.area[k] = 1 - a_bulk_c - aux_en.q_tot[k] = max(val1 * prog_gm.q_tot[k] - val2 * aux_tc.bulk.q_tot[k], 0) #Yair - this is here to prevent negative QT - aux_en.θ_liq_ice[k] = val1 * prog_gm.θ_liq_ice[k] - val2 * aux_tc.bulk.θ_liq_ice[k] - - ##### - ##### saturation_adjustment - ##### - ts_en = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k]) - - aux_en.T[k] = TD.air_temperature(ts_en) - aux_en.q_liq[k] = TD.liquid_specific_humidity(ts_en) - aux_en.q_ice[k] = TD.ice_specific_humidity(ts_en) - rho = TD.air_density(ts_en) - aux_en.buoy[k] = buoyancy_c(param_set, ρ0_c[k], rho) - - update_sat_unsat(en_thermo, state, k, ts_en) - aux_en.RH[k] = TD.relative_humidity(ts_en) - - ##### - ##### buoyancy - ##### - 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] += aux_up[i].area[k] * aux_up[i].buoy[k] + ##### + ##### face variables: diagnose primitive, diagnose env and compute bulk + ##### + @inbounds for k in real_face_indices(grid) + if is_surface_face(grid, k) + @inbounds for i in 1:(up.n_updrafts) + aux_up_f[i].w[k] = edmf.w_surface_bc[i] + end + else + @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 - @inbounds for i in 1:(up.n_updrafts) - aux_up[i].buoy[k] -= aux_gm.buoy[k] + aux_tc_f.bulk.w[k] = 0 + a_bulk_bcs = (; bottom = SetValue(sum(edmf.area_surface_bc)), top = SetZeroGradient()) + a_bulk_f = interpc2f(aux_tc.bulk.area, grid, k; a_bulk_bcs...) + if a_bulk_f > 0 + @inbounds for i in 1:(up.n_updrafts) + a_up_bcs = (; bottom = SetValue(edmf.area_surface_bc[i]), top = SetZeroGradient()) + 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 - aux_en.buoy[k] -= aux_gm.buoy[k] + # Assuming gm.W = 0! + aux_en_f.w[k] = -a_bulk_f / (1 - a_bulk_f) * aux_tc_f.bulk.w[k] end + + ##### + ##### diagnose_GMV_moments + ##### + #! format: off + get_GMV_CoVar(edmf, grid, state, :Hvar, :θ_liq_ice) + get_GMV_CoVar(edmf, grid, state, :QTvar, :q_tot) + get_GMV_CoVar(edmf, grid, state, :HQTcov, :θ_liq_ice, :q_tot) + GMV_third_m(edmf, grid, state, :Hvar, :θ_liq_ice, :H_third_m) + GMV_third_m(edmf, grid, state, :QTvar, :q_tot, :QT_third_m) + GMV_third_m(edmf, grid, state, :tke, :w, :W_third_m) + #! format: on + # TODO - use this inversion in free_convection_windspeed and not compute zi twice - @inbounds for k in real_center_indices(grid) - ts = thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k]) - aux_tc.θ_virt[k] = TD.virtual_pottemp(ts) - end edmf.zi = get_inversion(param_set, aux_tc.θ_virt, prog_gm.u, prog_gm.v, grid, surface.Ri_bulk_crit) update_surface(Case, grid, state, gm, TS, param_set) update_forcing(Case, grid, state, gm, TS, param_set) update_radiation(Case, grid, state, gm, TS, param_set) - ##### - ##### update_GMV_diagnostics - ##### - a_up_bulk = aux_tc.bulk.area - @inbounds for k in real_center_indices(grid) - aux_gm.q_liq[k] = (a_up_bulk[k] * aux_tc.bulk.q_liq[k] + (1 - a_up_bulk[k]) * aux_en.q_liq[k]) - aux_gm.q_ice[k] = (a_up_bulk[k] * aux_tc.bulk.q_ice[k] + (1 - a_up_bulk[k]) * aux_en.q_ice[k]) - aux_gm.T[k] = (a_up_bulk[k] * aux_tc.bulk.T[k] + (1 - a_up_bulk[k]) * aux_en.T[k]) - aux_gm.buoy[k] = (a_up_bulk[k] * aux_tc.bulk.buoy[k] + (1 - a_up_bulk[k]) * aux_en.buoy[k]) - end compute_pressure_plume_spacing(edmf, param_set) #####