From 9f4a376b78562d2bb5732f06fb98575ea126be10 Mon Sep 17 00:00:00 2001 From: Ignacio Date: Thu, 14 Oct 2021 22:28:22 -0700 Subject: [PATCH] Add consistent buoyancy gradient for environmental mean and quadratures. --- integration_tests/utils/generate_namelist.jl | 1 + integration_tests/utils/mse_tables.jl | 156 +++++++++---------- src/EDMF_Environment.jl | 7 +- src/Operators.jl | 48 ++++++ src/closures/buoyancy_gradients.jl | 63 ++++++-- src/closures/mixing_length.jl | 3 +- src/turbulence_functions.jl | 4 +- src/types.jl | 73 ++++++--- src/update_aux.jl | 96 +++++++++--- 9 files changed, 313 insertions(+), 138 deletions(-) diff --git a/integration_tests/utils/generate_namelist.jl b/integration_tests/utils/generate_namelist.jl index 707c79346..ad5912af8 100644 --- a/integration_tests/utils/generate_namelist.jl +++ b/integration_tests/utils/generate_namelist.jl @@ -142,6 +142,7 @@ function default_namelist(case_name::String) namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["constant_area"] = false namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["calculate_tke"] = true namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["mixing_length"] = "sbtd_eq" + namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["env_buoy_grad"] = "quadratures" namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_buoy"] = "normalmode" namelist_defaults["turbulence"]["EDMF_PrognosticTKE"]["pressure_closure_drag"] = "normalmode" diff --git a/integration_tests/utils/mse_tables.jl b/integration_tests/utils/mse_tables.jl index f408639d0..0e528ed1c 100644 --- a/integration_tests/utils/mse_tables.jl +++ b/integration_tests/utils/mse_tables.jl @@ -5,7 +5,7 @@ all_best_mse = OrderedCollections.OrderedDict() # all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict() -all_best_mse["ARM_SGP"]["qt_mean"] = 0.26446718607725306 +all_best_mse["ARM_SGP"]["qt_mean"] = 2.6506823775591898e-01 all_best_mse["ARM_SGP"]["updraft_area"] = 342.05657568721784 all_best_mse["ARM_SGP"]["updraft_w"] = 147.98545782669382 all_best_mse["ARM_SGP"]["updraft_qt"] = 28.294790221283453 @@ -13,50 +13,50 @@ all_best_mse["ARM_SGP"]["updraft_thetal"] = 170.9801640138966 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"] = 0.00012053009592043601 -all_best_mse["ARM_SGP"]["ql_mean"] = 204.94997013272453 +all_best_mse["ARM_SGP"]["ql_mean"] = 2.0547918525992608e+02 all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00011122266036009708 -all_best_mse["ARM_SGP"]["Hvar_mean"] = 2633.2221109819784 -all_best_mse["ARM_SGP"]["QTvar_mean"] = 1559.544125199823 +all_best_mse["ARM_SGP"]["Hvar_mean"] = 2.9194184175128917e+03 +all_best_mse["ARM_SGP"]["QTvar_mean"] = 1.7213808112900319e+03 # all_best_mse["Bomex"] = OrderedCollections.OrderedDict() -all_best_mse["Bomex"]["qt_mean"] = 0.11034670891224195 -all_best_mse["Bomex"]["updraft_area"] = 132.0788865649432 -all_best_mse["Bomex"]["updraft_w"] = 18.392509423445038 -all_best_mse["Bomex"]["updraft_qt"] = 6.288835389580369 -all_best_mse["Bomex"]["updraft_thetal"] = 69.45747637501809 -all_best_mse["Bomex"]["v_mean"] = 66.15195917633564 -all_best_mse["Bomex"]["u_mean"] = 2244.359495472965 -all_best_mse["Bomex"]["tke_mean"] = 55.75185915047925 -all_best_mse["Bomex"]["temperature_mean"] = 4.3858868869668704e-5 -all_best_mse["Bomex"]["ql_mean"] = 8.228711881484596 -all_best_mse["Bomex"]["thetal_mean"] = 4.448304470557014e-5 -all_best_mse["Bomex"]["Hvar_mean"] = 4019.5984832795702 -all_best_mse["Bomex"]["QTvar_mean"] = 1505.300250338699 +all_best_mse["Bomex"]["qt_mean"] = 1.0983320350059607e-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.5576868612678965e+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.5835472696934382e+01 +all_best_mse["Bomex"]["temperature_mean"] = 4.3677471984664055e-05 +all_best_mse["Bomex"]["ql_mean"] = 8.8506606580011589e+00 +all_best_mse["Bomex"]["thetal_mean"] = 4.4330671168211216e-05 +all_best_mse["Bomex"]["Hvar_mean"] = 4.2126196703263922e+03 +all_best_mse["Bomex"]["QTvar_mean"] = 1.5794628293458845e+03 # all_best_mse["DryBubble"] = OrderedCollections.OrderedDict() -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"]["updraft_area"] = 8.8530673193084747e-03 +all_best_mse["DryBubble"]["updraft_w"] = 1.2479889474537974e-02 +all_best_mse["DryBubble"]["updraft_thetal"] = 1.2223166344487256e-10 all_best_mse["DryBubble"]["u_mean"] = 0.0 -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"] = 0.0 +all_best_mse["DryBubble"]["tke_mean"] = 2.0135570156946084e-01 +all_best_mse["DryBubble"]["temperature_mean"] = 3.2783493600492222e-07 +all_best_mse["DryBubble"]["thetal_mean"] = 4.5557150098548290e-07 +all_best_mse["DryBubble"]["Hvar_mean"] = 9.6713965877036813e+01 # all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.04290066515065224 -all_best_mse["DYCOMS_RF01"]["ql_mean"] = 0.8976960436098287 -all_best_mse["DYCOMS_RF01"]["updraft_area"] = 31.57502795434024 -all_best_mse["DYCOMS_RF01"]["updraft_w"] = 5.415789637698043 -all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 1.734023675326015 -all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.186016892032455 -all_best_mse["DYCOMS_RF01"]["v_mean"] = 19017.018104720682 -all_best_mse["DYCOMS_RF01"]["u_mean"] = 74833.98050860969 -all_best_mse["DYCOMS_RF01"]["tke_mean"] = 18.837749675831937 -all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 0.00016283538235415605 -all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 0.00016376758990380806 -all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1227.7967628483952 -all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 478.07194167693547 +all_best_mse["DYCOMS_RF01"]["qt_mean"] = 3.0682630810336644e-02 +all_best_mse["DYCOMS_RF01"]["ql_mean"] = 2.7907515016732569e+00 +all_best_mse["DYCOMS_RF01"]["updraft_area"] = 2.9996855083383483e+01 +all_best_mse["DYCOMS_RF01"]["updraft_w"] = 3.6786963992723485e+00 +all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 1.8434661312139911e+00 +all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 4.6186138798587415e+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.4119502493062683e+01 +all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 1.5145561673937695e-04 +all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 1.5230990980765101e-04 +all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1.2554483566036963e+03 +all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 5.1751103217917216e+02 # all_best_mse["GABLS"] = OrderedCollections.OrderedDict() all_best_mse["GABLS"]["updraft_thetal"] = 5.584864486190616e-9 @@ -70,80 +70,80 @@ all_best_mse["GABLS"]["QTvar_mean"] = 0.017318287479694217 all_best_mse["GABLS"]["qt_mean"] = 0.21092607675396896 # all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict() -all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["v_mean"] = 0.0 -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"] = 0.0 -all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 0.0 +all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 3.9000477871154622e-09 +all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 3.2093573944402762e-06 +all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 1.6350375660088284e-07 +all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 1.5548562332198521e-06 +all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 7.0074885089062560e-10 +all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 2.4326279199247991e-13 +all_best_mse["life_cycle_Tan2018"]["v_mean"] = 5.3858193910463953e-07 +all_best_mse["life_cycle_Tan2018"]["u_mean"] = 3.1825523774843913e-09 +all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 1.9880489190175280e-07 +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.7791377140463445e-02 +all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 2.4478780510096252e-02 # all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict() all_best_mse["Nieuwstadt"]["updraft_area"] = 110.65702527331422 all_best_mse["Nieuwstadt"]["updraft_w"] = 12.298314607367429 -all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.3670456600826 +all_best_mse["Nieuwstadt"]["updraft_thetal"] = 1.1736704568882135e+02 all_best_mse["Nieuwstadt"]["u_mean"] = 59.73642032443851 -all_best_mse["Nieuwstadt"]["tke_mean"] = 222.28508474371915 +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"] = 1276.761018651613 +all_best_mse["Nieuwstadt"]["Hvar_mean"] = 1.2769135065447590e+03 # all_best_mse["Rico"] = OrderedCollections.OrderedDict() -all_best_mse["Rico"]["qt_mean"] = 1.3562728436394182 +all_best_mse["Rico"]["qt_mean"] = 1.3576174162823518e+00 all_best_mse["Rico"]["updraft_area"] = 483.4022083273521 -all_best_mse["Rico"]["updraft_w"] = 105.6746441121623 -all_best_mse["Rico"]["updraft_qt"] = 10.859715229316084 +all_best_mse["Rico"]["updraft_w"] = 1.0605422046690978e+02 +all_best_mse["Rico"]["updraft_qt"] = 1.0914863189252404e+01 all_best_mse["Rico"]["updraft_thetal"] = 132.79730672922278 -all_best_mse["Rico"]["v_mean"] = 23893.44384637507 -all_best_mse["Rico"]["u_mean"] = 167.0481771976619 -all_best_mse["Rico"]["tke_mean"] = 82.85087674864991 -all_best_mse["Rico"]["temperature_mean"] = 0.0006129402538716675 -all_best_mse["Rico"]["ql_mean"] = 67.02280413541271 +all_best_mse["Rico"]["v_mean"] = 2.3913389154144726e+04 +all_best_mse["Rico"]["u_mean"] = 1.6720325360173803e+02 +all_best_mse["Rico"]["tke_mean"] = 8.2957963768830496e+01 +all_best_mse["Rico"]["temperature_mean"] = 6.1310148877327848e-04 +all_best_mse["Rico"]["ql_mean"] = 6.7286920396505195e+01 all_best_mse["Rico"]["qr_mean"] = 750.8723638131503 -all_best_mse["Rico"]["thetal_mean"] = 0.0006043563413443778 +all_best_mse["Rico"]["thetal_mean"] = 6.0448834814914290e-04 all_best_mse["Rico"]["Hvar_mean"] = 18654.16088756812 all_best_mse["Rico"]["QTvar_mean"] = 4036.079013508425 # all_best_mse["Soares"] = OrderedCollections.OrderedDict() -all_best_mse["Soares"]["qt_mean"] = 0.12636886303765799 +all_best_mse["Soares"]["qt_mean"] = 1.2659089197243104e-01 all_best_mse["Soares"]["updraft_area"] = 112.7273316705158 -all_best_mse["Soares"]["updraft_w"] = 11.265484888980165 -all_best_mse["Soares"]["updraft_qt"] = 23.119052296484572 -all_best_mse["Soares"]["updraft_thetal"] = 65.25334973003676 -all_best_mse["Soares"]["u_mean"] = 101.82502556650786 -all_best_mse["Soares"]["tke_mean"] = 184.8671644305064 +all_best_mse["Soares"]["updraft_w"] = 1.1266999168717131e+01 +all_best_mse["Soares"]["updraft_qt"] = 2.3120736199996990e+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.8487029454315532e+02 all_best_mse["Soares"]["temperature_mean"] = 1.1031965842549474e-5 all_best_mse["Soares"]["thetal_mean"] = 1.0453842793106344e-5 all_best_mse["Soares"]["Hvar_mean"] = 1145.5445941618634 # all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict() all_best_mse["TRMM_LBA"]["qt_mean"] = 3.696894480356394 -all_best_mse["TRMM_LBA"]["updraft_area"] = 7561.609964448662 -all_best_mse["TRMM_LBA"]["updraft_w"] = 11162.65860693764 +all_best_mse["TRMM_LBA"]["updraft_area"] = 7.5616159743597864e+03 +all_best_mse["TRMM_LBA"]["updraft_w"] = 1.1162662767994134e+04 all_best_mse["TRMM_LBA"]["updraft_qt"] = 178.62874777200045 all_best_mse["TRMM_LBA"]["updraft_thetal"] = 1972.5026068596098 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"] = 8195.78791967079 +all_best_mse["TRMM_LBA"]["tke_mean"] = 8.1959925236028794e+03 all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0011323111668659397 -all_best_mse["TRMM_LBA"]["ql_mean"] = 2623.5510869685795 +all_best_mse["TRMM_LBA"]["ql_mean"] = 2.6235652901417766e+03 all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.010479008728209362 all_best_mse["TRMM_LBA"]["Hvar_mean"] = 2511.796313996477 all_best_mse["TRMM_LBA"]["QTvar_mean"] = 2460.6541168787 # all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict() -all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.491835618015334 -all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7057083566611517 -all_best_mse["LES_driven_SCM"]["u_mean"] = 1.235432758754188 -all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0028509097874652232 -all_best_mse["LES_driven_SCM"]["ql_mean"] = 252.08833547882853 -all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0031436563856298365 +all_best_mse["LES_driven_SCM"]["qt_mean"] = 6.6765992515646833e+00 +all_best_mse["LES_driven_SCM"]["v_mean"] = 3.7142882513530031e+00 +all_best_mse["LES_driven_SCM"]["u_mean"] = 1.2391794278032586e+00 +all_best_mse["LES_driven_SCM"]["temperature_mean"] = 2.9024394724690302e-03 +all_best_mse["LES_driven_SCM"]["ql_mean"] = 1.1397368402146715e+03 +all_best_mse["LES_driven_SCM"]["thetal_mean"] = 3.2049121067961291e-03 # ################################# ################################# diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index a449536b3..df2735b59 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -52,13 +52,15 @@ function update_cloud_dry(en_thermo::EnvironmentThermodynamics, state, k, ts) aux_en = center_aux_environment(state) if q_liq > 0.0 aux_en.cloud_fraction[k] = 1.0 - en_thermo.th_cloudy[k] = TD.liquid_ice_pottemp(ts) + en_thermo.th_cloudy[k] = TD.dry_pottemp(ts) + en_thermo.thl_cloudy[k] = TD.liquid_ice_pottemp(ts) en_thermo.t_cloudy[k] = TD.air_temperature(ts) en_thermo.qt_cloudy[k] = TD.total_specific_humidity(ts) en_thermo.qv_cloudy[k] = TD.vapor_specific_humidity(ts) else aux_en.cloud_fraction[k] = 0.0 - en_thermo.th_dry[k] = TD.liquid_ice_pottemp(ts) + en_thermo.th_dry[k] = TD.dry_pottemp(ts) + en_thermo.thv_dry[k] = TD.virtual_pottemp(ts) en_thermo.qt_dry[k] = TD.total_specific_humidity(ts) end return @@ -256,6 +258,7 @@ function sgs_quadrature(en_thermo::EnvironmentThermodynamics, grid, state, en, r en_thermo.qt_cloudy[k] = outer_env[i_qt_cld] ts_cld = TD.PhaseEquil_pTq(param_set, p0_c[k], en_thermo.t_cloudy[k], en_thermo.qt_cloudy[k]) en_thermo.th_cloudy[k] = TD.dry_pottemp(ts_cld) + en_thermo.thl_cloudy[k] = TD.liquid_ice_pottemp(ts_cld) # update var/covar rain sources en_thermo.Hvar_rain_dt[k] = outer_src[i_SH_H] - outer_src[i_SH] * aux_en.θ_liq_ice[k] diff --git a/src/Operators.jl b/src/Operators.jl index ac0d52a2a..b11a3bbca 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -250,6 +250,54 @@ function c∇(f::SA.SVector, grid::Grid, ::BottomBCTag, ::Extrapolate) return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2 end +# ∇(center data) for possibly vanishing subdomains. + +""" + c∇_vanishing_subdomain + +Used when the vertical gradient of a field must be computed over a conditional subdomain which may vanish +below or above the current cell. If the subdomains vanishes, a default user-defined gradient is returned. + +Inputs: + - f_cut: Slice of field. + - sd_cut: Slice of subdomain volume fraction. + - default∇: Gradient used in vanishing subdomains. +""" +function c∇_vanishing_subdomain( + f_cut::SA.SVector, + sd_cut::SA.SVector, + default∇::FT, + grid::Grid, + k; + bottom = NoBCGivenError(), + top = NoBCGivenError(), +) where {FT <: Real} + if is_surface_center(grid, k) + return c∇(f_cut, grid, BottomBCTag(), bottom) + elseif is_toa_center(grid, k) + return c∇(f_cut, grid, TopBCTag(), top) + else + return c∇_vanishing_subdomain(f_cut, sd_cut, default∇, grid, InteriorTag()) + end +end +function c∇_vanishing_subdomain( + f::SA.SVector, + sd::SA.SVector, + default∇::FT, + grid::Grid, + ::InteriorTag, +) where {FT <: Real} + @assert length(f) == 3 + @assert length(sd) == 3 + if sd[1] * sd[3] ≈ FT(0) + return default∇ + else + f_dual⁺ = SA.SVector(f[2], f[3]) + f_dual⁻ = SA.SVector(f[1], f[2]) + return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2 + end +end + ##### ##### ∇(face data) ##### diff --git a/src/closures/buoyancy_gradients.jl b/src/closures/buoyancy_gradients.jl index ded569e4d..f8d74ca68 100644 --- a/src/closures/buoyancy_gradients.jl +++ b/src/closures/buoyancy_gradients.jl @@ -1,36 +1,71 @@ -function buoyancy_gradients(param_set, bg_model::Tan2018) +""" + buoyancy_gradients(param_set, bg_model::EnvBuoyGrad{FT, EBG}) where {FT <: Real, EBG <: EnvBuoyGradClosure} + +Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions. +The dispatch on EnvBuoyGrad type is performed at the EnvBuoyGrad construction time, and the analytical solutions +used here are consistent for both mean fields and conditional fields obtained from assumed distributions +over the conserved thermodynamic variables. +""" +function buoyancy_gradients(param_set, bg_model::EnvBuoyGrad{FT, EBG}) where {FT <: Real, EBG <: EnvBuoyGradClosure} g = CPP.grav(param_set) molmass_ratio = CPP.molmass_ratio(param_set) R_d = CPP.R_d(param_set) R_v = CPP.R_v(param_set) - # buoyancy_gradients phase_part = TD.PhasePartition(0.0, 0.0, 0.0) # assuming R_d = R_m Π = TD.exner_given_pressure(param_set, bg_model.p0, phase_part) ∂b∂θv = g * (R_d / bg_model.alpha0 / bg_model.p0) * Π - ∂b∂θl_dry = ∂b∂θv * (1 + (molmass_ratio - 1) * bg_model.qt_dry) - ∂b∂qt_dry = ∂b∂θv * bg_model.th_dry * (molmass_ratio - 1) if bg_model.en_cld_frac > 0.0 - ts_cloudy = thermo_state_pθq(param_set, bg_model.p0, bg_model.th_cloudy, bg_model.qt_cloudy) - lh = TD.latent_heat_vapor(param_set, bg_model.t_cloudy) + ts_cloudy = thermo_state_pθq(param_set, bg_model.p0, bg_model.θ_liq_ice_cloudy, bg_model.qt_cloudy) + phase_part = TD.PhasePartition(ts_cloudy) + lh = TD.latent_heat_liq_ice(param_set, phase_part) cp_m = TD.cp_m(ts_cloudy) ∂b∂θl_cld = ( ∂b∂θv * (1 + molmass_ratio * (1 + lh / R_v / bg_model.t_cloudy) * bg_model.qv_cloudy - bg_model.qt_cloudy) / (1 + lh * lh / cp_m / R_v / bg_model.t_cloudy / bg_model.t_cloudy * bg_model.qv_cloudy) ) - ∂b∂qt_cld = (lh / cp_m / bg_model.t_cloudy * ∂b∂θl_cld - ∂b∂θv) * bg_model.th_cloudy + ∂b∂qt_cld = (lh / cp_m / bg_model.t_cloudy * ∂b∂θl_cld - ∂b∂θv) * bg_model.θ_cloudy + else + ∂b∂θl_cld = FT(0) + ∂b∂qt_cld = FT(0) + end + + ∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy = buoyancy_gradient_chain_rule(bg_model, ∂b∂θv, ∂b∂θl_cld, ∂b∂qt_cld) + return GradBuoy(∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy) +end + +""" + buoyancy_gradient_chain_rule( + bg_model::EnvBuoyGrad{FT, EBG}, + ∂b∂θv::FT, + ∂b∂θl_cld::FT, + ∂b∂qt_cld::FT, + ) where {FT <: Real, EBG <: EnvBuoyGradClosure} + +Returns the vertical buoyancy gradients in the environment, as well as in its dry and cloudy volume fractions, +from the partial derivatives with respect to thermodynamic variables in dry and cloudy volumes. +""" +function buoyancy_gradient_chain_rule( + bg_model::EnvBuoyGrad{FT, EBG}, + ∂b∂θv::FT, + ∂b∂θl_cld::FT, + ∂b∂qt_cld::FT, +) where {FT <: Real, EBG <: EnvBuoyGradClosure} + if bg_model.en_cld_frac > FT(0) + ∂b∂z_θl_cld = ∂b∂θl_cld * bg_model.∂θl∂z_cloudy + ∂b∂z_qt_cld = ∂b∂qt_cld * bg_model.∂qt∂z_cloudy else - ∂b∂θl_cld = 0 - ∂b∂qt_cld = 0 + ∂b∂z_θl_cld = FT(0) + ∂b∂z_qt_cld = FT(0) end - ∂b∂θl = (bg_model.en_cld_frac * ∂b∂θl_cld + (1 - bg_model.en_cld_frac) * ∂b∂θl_dry) - ∂b∂qt = (bg_model.en_cld_frac * ∂b∂qt_cld + (1 - bg_model.en_cld_frac) * ∂b∂qt_dry) - ∂b∂z_θl = ∂b∂θl * bg_model.∂θl∂z - ∂b∂z_qt = ∂b∂qt * bg_model.∂qt∂z + ∂b∂z_dry = bg_model.en_cld_frac < FT(1) ? ∂b∂θv * bg_model.∂θv∂z_dry : FT(0) + + ∂b∂z_cloudy = ∂b∂z_θl_cld + ∂b∂z_qt_cld + ∂b∂z = (1 - bg_model.en_cld_frac) * ∂b∂z_dry + bg_model.en_cld_frac * ∂b∂z_cloudy - return GradBuoy(∂b∂θl, ∂b∂qt, ∂b∂z_θl, ∂b∂z_qt) + return ∂b∂z, ∂b∂z_dry, ∂b∂z_cloudy end diff --git a/src/closures/mixing_length.jl b/src/closures/mixing_length.jl index ed4d62a13..98c68d0d9 100644 --- a/src/closures/mixing_length.jl +++ b/src/closures/mixing_length.jl @@ -23,8 +23,7 @@ function mixing_length(param_set, ml_model::MinDisspLen{FT}) where {FT} end # compute l_TKE - the production/destruction term - a_pd = - c_m * (ml_model.Shear² - ml_model.∂b∂z_θl / ml_model.Pr - ml_model.∂b∂z_qt / ml_model.Pr) * sqrt(ml_model.tke) + a_pd = c_m * (ml_model.Shear² - ml_model.∂b∂z / ml_model.Pr) * sqrt(ml_model.tke) # Dissipation term c_neg = c_d * ml_model.tke * sqrt(ml_model.tke) # Subdomain exchange term diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl index 2460f99df..5f06ce192 100644 --- a/src/turbulence_functions.jl +++ b/src/turbulence_functions.jl @@ -71,8 +71,8 @@ function set_cloudbase_flag(ql, current_flag) return new_flag end -function gradient_Richardson_number(∂b∂θ_l, ∂b∂q_tot, Shear², ϵ) - return min(∂b∂θ_l / max(Shear², ϵ) + ∂b∂q_tot / max(Shear², ϵ), 0.25) +function gradient_Richardson_number(∂b∂z, Shear², ϵ) + return min(∂b∂z / max(Shear², ϵ), 0.25) end function turbulent_Prandtl_number(param_set, obukhov_length, ∇Ri) diff --git a/src/types.jl b/src/types.jl index d2e189744..2b6e7f241 100644 --- a/src/types.jl +++ b/src/types.jl @@ -99,12 +99,10 @@ Base.@kwdef struct MinDisspLen{FT, T} Pr::FT "reference pressure" p0::FT - "partial change in buoyancy due to θl changes" - ∂b∂z_θl::FT + "environmental vertical buoyancy gradient" + ∂b∂z::FT "env shear" Shear²::FT - "partial change in buoyancy due to qt changes" - ∂b∂z_qt::FT "virtual potential temperature gradient" ∂θv∂z::FT "total specific humidity gradient" @@ -145,27 +143,31 @@ end """ GradBuoy +Environmental buoyancy gradients. + $(DocStringExtensions.FIELDS) """ Base.@kwdef struct GradBuoy{FT} - ∂b∂θl::FT - ∂b∂qt::FT - ∂b∂z_qt::FT - ∂b∂z_θl::FT + "environmental vertical buoyancy gradient" + ∂b∂z::FT + "vertical buoyancy gradient in the dry part of the environment" + ∂b∂z_dry::FT + "vertical buoyancy gradient in the cloudy part of the environment" + ∂b∂z_cloudy::FT end +abstract type EnvBuoyGradClosure end +struct BuoyGradMean <: EnvBuoyGradClosure end +struct BuoyGradQuadratures <: EnvBuoyGradClosure end + """ - Tan2018 + EnvBuoyGrad -The buoyancy gradient estimating in Eqs. () - () in Tan2018 +Variables used in the environmental buoyancy gradient computation. $(DocStringExtensions.FIELDS) """ -Base.@kwdef struct Tan2018{FT} - "specific humidity in the non cloudy part" - qt_dry::FT - "potential temperature in the non cloudy part" - th_dry::FT +Base.@kwdef struct EnvBuoyGrad{FT, EBC <: EnvBuoyGradClosure} "temperature in the cloudy part" t_cloudy::FT "vapor specific humidity in the cloudy part" @@ -173,11 +175,15 @@ Base.@kwdef struct Tan2018{FT} "total specific humidity in the cloudy part" qt_cloudy::FT "potential temperature in the cloudy part" - th_cloudy::FT - "total specific humidity gradient" - ∂qt∂z::FT - "potential temperature gradient" - ∂θl∂z::FT + θ_cloudy::FT + "liquid ice potential temperature in the cloudy part" + θ_liq_ice_cloudy::FT + "virtual potential temperature gradient in the non cloudy part" + ∂θv∂z_dry::FT + "total specific humidity gradient in the cloudy part" + ∂qt∂z_cloudy::FT + "liquid ice potential temperature gradient in the cloudy part" + ∂θl∂z_cloudy::FT "reference pressure" p0::FT "cloud fraction" @@ -185,6 +191,9 @@ Base.@kwdef struct Tan2018{FT} "specific volume" alpha0::FT end +function EnvBuoyGrad(::EBG; t_cloudy::FT, bg_kwargs...) where {FT <: Real, EBG <: EnvBuoyGradClosure} + return EnvBuoyGrad{FT, EBG}(; t_cloudy, bg_kwargs...) +end Base.@kwdef mutable struct RainVariables rain_model::String = "default_rain_model" @@ -331,10 +340,12 @@ struct EnvironmentThermodynamics{A1} quadrature_type::String qt_dry::A1 th_dry::A1 + thv_dry::A1 t_cloudy::A1 qv_cloudy::A1 qt_cloudy::A1 th_cloudy::A1 + thl_cloudy::A1 Hvar_rain_dt::A1 QTvar_rain_dt::A1 HQTcov_rain_dt::A1 @@ -346,11 +357,13 @@ struct EnvironmentThermodynamics{A1} qt_dry = center_field(grid) th_dry = center_field(grid) + thv_dry = center_field(grid) t_cloudy = center_field(grid) qv_cloudy = center_field(grid) qt_cloudy = center_field(grid) th_cloudy = center_field(grid) + thl_cloudy = center_field(grid) Hvar_rain_dt = center_field(grid) QTvar_rain_dt = center_field(grid) @@ -364,10 +377,12 @@ struct EnvironmentThermodynamics{A1} quadrature_type, qt_dry, th_dry, + thv_dry, t_cloudy, qv_cloudy, qt_cloudy, th_cloudy, + thl_cloudy, Hvar_rain_dt, QTvar_rain_dt, HQTcov_rain_dt, @@ -564,6 +579,7 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} entr_surface_bc::Float64 detr_surface_bc::Float64 sde_model::sde_struct + bg_closure::EnvBuoyGradClosure function EDMF_PrognosticTKE(namelist, grid::Grid, param_set::PS) where {PS} # get values from namelist prandtl_number = namelist["turbulence"]["prandtl_number_0"] @@ -719,6 +735,22 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} end sde_model = sde_struct{closure_type}(u0 = 1, dt = dt) + bg_type = parse_namelist( + namelist, + "turbulence", + "EDMF_PrognosticTKE", + "env_buoy_grad"; + default = "mean", + valid_options = ["mean", "quadratures"], + ) + bg_closure = if bg_type == "mean" + BuoyGradMean() + elseif bg_type == "quadratures" + BuoyGradQuadratures() + else + error("Something went wrong. Invalid environmental buoyancy gradient closure type '$bg_type'") + end + # Added by Ignacio : Length scheme in use (mls), and smooth min effect (ml_ratio) # Variable Prandtl number initialized as neutral value. ones_vec = center_field(grid) @@ -793,6 +825,7 @@ mutable struct EDMF_PrognosticTKE{A1, A2, IE} entr_surface_bc, detr_surface_bc, sde_model, + bg_closure, ) end end diff --git a/src/update_aux.jl b/src/update_aux.jl index 13630e44a..f8fd3fbdc 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -300,23 +300,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ∂qt∂z = c∇(QT_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0)) THL_cut = ccut(aux_en.θ_liq_ice, grid, k) ∂θl∂z = c∇(THL_cut, grid, k; bottom = SetGradient(0), top = SetGradient(0)) - # buoyancy_gradients - bg_model = Tan2018(; - qt_dry = en_thermo.qt_dry[k], - th_dry = en_thermo.th_dry[k], - t_cloudy = en_thermo.t_cloudy[k], - qv_cloudy = en_thermo.qv_cloudy[k], - qt_cloudy = en_thermo.qt_cloudy[k], - th_cloudy = en_thermo.th_cloudy[k], - ∂qt∂z = ∂qt∂z, - ∂θl∂z = ∂θl∂z, - p0 = p0_c[k], - en_cld_frac = aux_en.cloud_fraction[k], - alpha0 = α0_c[k], - ) - bg = buoyancy_gradients(param_set, bg_model) - # Limiting stratification scale (Deardorff, 1976) p0_cut = ccut(p0_c, grid, k) T_cut = ccut(aux_en.T, grid, k) QT_cut = ccut(aux_en.q_tot, grid, k) @@ -325,10 +309,83 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) thv_cut = TD.virtual_pottemp.(ts_cut) ts = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k]) + θ = TD.dry_pottemp(ts) θv = TD.virtual_pottemp(ts) ∂θv∂z = c∇(thv_cut, grid, k; bottom = SetGradient(0), top = Extrapolate()) + + # buoyancy_gradients + if edmf.bg_closure == BuoyGradMean() + # First order approximation: Use environmental mean fields. + bg_kwargs = (; + t_cloudy = aux_en.T[k], + # WARNING: ICE MISSING + qv_cloudy = (aux_en.q_tot - aux_en.q_liq)[k], + qt_cloudy = aux_en.q_tot[k], + θ_cloudy = θ, + θ_liq_ice_cloudy = aux_en.θ_liq_ice[k], + ∂θv∂z_dry = ∂θv∂z, + ∂qt∂z_cloudy = ∂qt∂z, + ∂θl∂z_cloudy = ∂θl∂z, + p0 = p0_c[k], + en_cld_frac = aux_en.cloud_fraction[k], + alpha0 = α0_c[k], + ) + bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...) + + elseif edmf.bg_closure == BuoyGradQuadratures() + # Second order approximation: Use dry and cloudy environmental fields. + cf_cut = ccut(aux_en.cloud_fraction, grid, k) + QT_cloudy_cut = ccut(en_thermo.qt_cloudy, grid, k) + ∂qt∂z_cloudy = c∇_vanishing_subdomain( + QT_cloudy_cut, + cf_cut, + ∂qt∂z, + grid, + k; + bottom = SetGradient(0), + top = SetGradient(0), + ) + THL_cloudy_cut = ccut(en_thermo.thl_cloudy, grid, k) + ∂θl∂z_cloudy = c∇_vanishing_subdomain( + THL_cloudy_cut, + cf_cut, + ∂θl∂z, + grid, + k; + bottom = SetGradient(0), + top = SetGradient(0), + ) + THV_dry_cut = ccut(en_thermo.thv_dry, grid, k) + ∂θv∂z_dry = c∇_vanishing_subdomain( + THV_dry_cut, + cf_cut, + ∂θv∂z, + grid, + k; + bottom = SetGradient(0), + top = SetGradient(0), + ) + + bg_kwargs = (; + t_cloudy = en_thermo.t_cloudy[k], + qv_cloudy = en_thermo.qv_cloudy[k], + qt_cloudy = en_thermo.qt_cloudy[k], + θ_cloudy = en_thermo.th_cloudy[k], + θ_liq_ice_cloudy = en_thermo.thl_cloudy[k], + ∂θv∂z_dry = ∂θv∂z_dry, + ∂qt∂z_cloudy = ∂qt∂z_cloudy, + ∂θl∂z_cloudy = ∂θl∂z_cloudy, + p0 = p0_c[k], + en_cld_frac = aux_en.cloud_fraction[k], + alpha0 = α0_c[k], + ) + bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...) + end + bg = buoyancy_gradients(param_set, bg_model) + + # Limiting stratification scale (Deardorff, 1976) # compute ∇Ri and Pr - ∇_Ri = gradient_Richardson_number(bg.∂b∂z_θl, bg.∂b∂z_qt, Shear², eps(0.0)) + ∇_Ri = gradient_Richardson_number(bg.∂b∂z, Shear², eps(0.0)) edmf.prandtl_nvec[k] = turbulent_Prandtl_number(param_set, obukhov_length, ∇_Ri) ml_model = MinDisspLen(; @@ -338,9 +395,8 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) ustar = surface.ustar, Pr = edmf.prandtl_nvec[k], p0 = p0_c[k], - ∂b∂z_θl = bg.∂b∂z_θl, + ∂b∂z = bg.∂b∂z, Shear² = Shear², - ∂b∂z_qt = bg.∂b∂z_qt, ∂θv∂z = ∂θv∂z, ∂qt∂z = ∂qt∂z, ∂θl∂z = ∂θl∂z, @@ -368,7 +424,7 @@ function update_aux!(edmf, gm, grid, state, Case, param_set, TS) KM[k] = c_m * edmf.mixing_length[k] * sqrt(max(prog_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_θl + bg.∂b∂z_qt) + aux_en_2m.tke.buoy[k] = -ml_model.a_en * ρ0_c[k] * KH[k] * bg.∂b∂z end compute_covariance_entr(edmf, grid, state, :tke, :w)