From 13154bd9af174570418e4208eb003ab664003b74 Mon Sep 17 00:00:00 2001 From: yairchn Date: Wed, 27 Apr 2022 17:22:18 -0700 Subject: [PATCH] =?UTF-8?q?change=20prognostic=20variables=20from=20(?= =?UTF-8?q?=CF=81=CE=B8=5Fliq=5Fice,=CF=81q=5Ftot)=20to=20(=CF=81e=5Ftot,?= =?UTF-8?q?=CF=81q=5Ftot)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit WIP --- .buildkite/pipeline.yml | 10 +- driver/Forcing.jl | 6 +- driver/Surface.jl | 8 +- driver/compute_diagnostics.jl | 2 - driver/dycore.jl | 133 +++++++++----- driver/generate_namelist.jl | 3 +- driver/initial_conditions.jl | 15 +- driver/main.jl | 4 +- integration_tests/3dBomex.jl | 4 +- integration_tests/cli_options.jl | 3 - integration_tests/driver.jl | 1 - post_processing/mse_tables.jl | 290 +++++++++++++++---------------- src/EDMF_Environment.jl | 42 +++-- src/EDMF_Precipitation.jl | 13 +- src/EDMF_Updrafts.jl | 6 +- src/EDMF_functions.jl | 156 +++++++++++++---- src/TurbulenceConvection.jl | 4 +- src/diagnostics.jl | 2 +- src/microphysics_coupling.jl | 25 ++- src/thermodynamics.jl | 52 ++++++ src/types.jl | 25 +-- src/update_aux.jl | 13 +- src/variables.jl | 9 +- 23 files changed, 523 insertions(+), 303 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index aa14faed29..26f6ed49ad 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -154,12 +154,12 @@ steps: command: "julia --color=yes --project=integration_tests integration_tests/driver.jl --case Rico --thermo_covariance_model prognostic --skip_tests true --suffix _prognostic_covarinaces" artifact_paths: "Output.Rico.01_prognostic_covarinaces/stats/comparison/*" - - group: "3D Experiments" - steps: + # - group: "3D Experiments" + # steps: - - label: ":package: 3D Bomex" - command: "julia --color=yes --project=integration_tests integration_tests/3dBomex.jl" - artifact_paths: "integration_tests/output/3dBomex/*" + # - label: ":package: 3D Bomex" + # command: "julia --color=yes --project=integration_tests integration_tests/3dBomex.jl" + # artifact_paths: "integration_tests/output/3dBomex/*" - group: "Performance monitoring" steps: diff --git a/driver/Forcing.jl b/driver/Forcing.jl index 64c5379cc9..edc25fac7e 100644 --- a/driver/Forcing.jl +++ b/driver/Forcing.jl @@ -9,7 +9,7 @@ function initialize(self::TC.ForcingBase{TC.ForcingLES}, grid, state, LESDat::TC zc_les = Array(TC.get_nc_data(data, "zc")) dTdt_hadv = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dtdt_hadv", imin, imax)) - H_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "thetali_mean", imin, imax)) + T_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "temperature_mean", imin, imax)) dTdt_fluc = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dtdt_fluc", imin, imax)) dqtdt_hadv = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "dqtdt_hadv", imin, imax)) qt_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "qt_mean", imin, imax)) @@ -17,11 +17,11 @@ function initialize(self::TC.ForcingBase{TC.ForcingLES}, grid, state, LESDat::TC subsidence = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "ls_subsidence", imin, imax)) u_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "u_mean", imin, imax)) v_nudge = TC.pyinterp(grid.zc, zc_les, TC.mean_nc_data(data, "profiles", "v_mean", imin, imax)) - (; dTdt_hadv, H_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge) + (; dTdt_hadv, T_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge) end for k in TC.real_center_indices(grid) aux_gm.dTdt_hadv[k] = nt.dTdt_hadv[k] - aux_gm.H_nudge[k] = nt.H_nudge[k] + aux_gm.T_nudge[k] = nt.T_nudge[k] aux_gm.dTdt_fluc[k] = nt.dTdt_fluc[k] aux_gm.dqtdt_hadv[k] = nt.dqtdt_hadv[k] aux_gm.qt_nudge[k] = nt.qt_nudge[k] diff --git a/driver/Surface.jl b/driver/Surface.jl index e277e13a9c..343a91bbc7 100644 --- a/driver/Surface.jl +++ b/driver/Surface.jl @@ -68,7 +68,7 @@ function get_surface( ch = result.Ch, ρu_flux = surf_params.zero_uv_fluxes ? FT(0) : result.ρτxz, ρv_flux = surf_params.zero_uv_fluxes ? FT(0) : result.ρτyz, - ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in), + ρe_tot_flux = shf + lhf, ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in), wstar = convective_vel, ρq_liq_flux = FT(0), @@ -126,7 +126,7 @@ function get_surface( ustar = result.ustar, ρu_flux = result.ρτxz, ρv_flux = result.ρτyz, - ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in), + ρe_tot_flux = shf + lhf, ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in), bflux = result.buoy_flux, wstar = convective_vel, @@ -183,7 +183,7 @@ function get_surface( ustar = result.ustar, ρu_flux = result.ρτxz, ρv_flux = result.ρτyz, - ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in), + ρe_tot_flux = shf + lhf, ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in), bflux = result.buoy_flux, wstar = convective_vel, @@ -245,7 +245,7 @@ function get_surface( ustar = result.ustar, ρu_flux = result.ρτxz, ρv_flux = result.ρτyz, - ρθ_liq_ice_flux = shf / TD.cp_m(param_set, ts_in), + ρe_tot_flux = shf + lhf, ρq_tot_flux = lhf / TD.latent_heat_vapor(param_set, ts_in), bflux = result.buoy_flux, wstar = convective_vel, diff --git a/driver/compute_diagnostics.jl b/driver/compute_diagnostics.jl index cdffc2cc37..bdb2453309 100644 --- a/driver/compute_diagnostics.jl +++ b/driver/compute_diagnostics.jl @@ -139,8 +139,6 @@ function compute_diagnostics!( @inbounds for k in TC.real_center_indices(grid) aux_gm.s[k] = TD.specific_entropy(param_set, ts_gm[k]) - aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k] - aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k] aux_en.s[k] = TD.specific_entropy(param_set, ts_en[k]) end @inbounds for k in TC.real_center_indices(grid) diff --git a/driver/dycore.jl b/driver/dycore.jl index 7e7bb6bece..bb41edf3ff 100644 --- a/driver/dycore.jl +++ b/driver/dycore.jl @@ -62,7 +62,7 @@ cent_aux_vars_gm(FT, edmf) = (; subsidence = FT(0), #Large-scale subsidence dTdt_hadv = FT(0), #Horizontal advection of temperature dqtdt_hadv = FT(0), #Horizontal advection of moisture - H_nudge = FT(0), #Reference H profile for relaxation tendency + T_nudge = FT(0), #Reference T profile for relaxation tendency qt_nudge = FT(0), #Reference qt profile for relaxation tendency dTdt_fluc = FT(0), #Vertical turbulent advection of temperature dqtdt_fluc = FT(0), #Vertical turbulent advection of moisture @@ -70,7 +70,7 @@ cent_aux_vars_gm(FT, edmf) = (; v_nudge = FT(0), #Reference v profile for relaxation tendency ug = FT(0), #Geostrophic u velocity vg = FT(0), #Geostrophic v velocity - ∇θ_liq_ice_gm = FT(0), + ∇MSE_gm = FT(0), ∇q_tot_gm = FT(0), cent_aux_vars_gm_moisture(FT, edmf.moisture_model)..., θ_virt = FT(0), @@ -89,7 +89,7 @@ face_aux_vars_gm(FT, edmf) = (; diffusive_flux_s = FT(0), total_flux_s = FT(0), f_rad = FT(0), - sgs_flux_θ_liq_ice = FT(0), + sgs_flux_h_tot = FT(0), sgs_flux_q_tot = FT(0), face_aux_vars_gm_moisture(FT, edmf.moisture_model)..., sgs_flux_u = FT(0), @@ -139,7 +139,7 @@ cent_prognostic_vars_gm(::Type{FT}, local_geometry, edmf) where {FT} = (; ρ = FT(0), u = FT(0), v = FT(0), - ρθ_liq_ice = FT(0), + ρe_tot = FT(0), ρq_tot = FT(0), # TODO: Change to: # uₕ = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry), @@ -221,7 +221,6 @@ function compute_ref_state!( @info "z_span = $z_span" prob = ODE.ODEProblem(rhs, logp, z_span) sol = ODE.solve(prob, ODE.Tsit5(), reltol = 1e-12, abstol = 1e-12) - parent(p_f) .= sol.(vec(grid.zf)) parent(p_c) .= sol.(vec(grid.zc)) @@ -241,23 +240,43 @@ function compute_ref_state!( return nothing end -function set_prog_from_aux!(state, ::TC.ρθVar) + +function set_thermo_state_peq!(state, grid, moisture_model, param_set) + Ic = CCO.InterpolateF2C() + FT = eltype(grid) + ts_gm = TC.center_aux_grid_mean(state).ts prog_gm = TC.center_prog_grid_mean(state) + prog_gm_f = TC.face_prog_grid_mean(state) aux_gm = TC.center_aux_grid_mean(state) - @. prog_gm.ρθ_liq_ice = prog_gm.ρ * aux_gm.θ_liq_ice - @. prog_gm.ρq_tot .= prog_gm.ρ * aux_gm.q_tot + p_c = aux_gm.p + ρ_c = prog_gm.ρ + w_c = copy(prog_gm.ρe_tot) + @. w_c = Ic(FT(0) + prog_gm_f.w) + @inbounds for k in TC.real_center_indices(grid) + thermo_args = if moisture_model isa TC.EquilibriumMoisture + () + elseif moisture_model isa TC.NonEquilibriumMoisture + (prog_gm.q_liq[k], prog_gm.q_ice[k]) + else + error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium") + end + e_kin = TC.kinetic_energy(prog_gm.u[k], prog_gm.v[k], w_c[k]) + e_pot = TC.geopotential(param_set, grid.zc.z[k]) + e_int = prog_gm.ρe_tot[k] / ρ_c[k] - e_kin - e_pot + ts_gm[k] = TC.thermo_state_peq(param_set, p_c[k], e_int, aux_gm.q_tot[k], thermo_args...) + aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp(param_set, ts_gm[k]) + aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k] + end return nothing end -function set_thermo_state!(state, grid, moisture_model, param_set) +function set_thermo_state_pθq!(state, grid, moisture_model, param_set) + Ic = CCO.InterpolateF2C() ts_gm = TC.center_aux_grid_mean(state).ts prog_gm = TC.center_prog_grid_mean(state) aux_gm = TC.center_aux_grid_mean(state) p_c = aux_gm.p - ρ_c = prog_gm.ρ @inbounds for k in TC.real_center_indices(grid) - aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k] - aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k] thermo_args = if moisture_model isa TC.EquilibriumMoisture () elseif moisture_model isa TC.NonEquilibriumMoisture @@ -270,14 +289,31 @@ function set_thermo_state!(state, grid, moisture_model, param_set) return nothing end +function set_grid_mean_from_thermo_state!(param_set, state, grid) + Ic = CCO.InterpolateF2C() + ts_gm = TC.center_aux_grid_mean(state).ts + prog_gm = TC.center_prog_grid_mean(state) + prog_gm_f = TC.face_prog_grid_mean(state) + aux_gm = TC.center_aux_grid_mean(state) + ρ_c = prog_gm.ρ + FT = eltype(grid) + w_c = copy(prog_gm.ρe_tot) + @. w_c = Ic(FT(0) + prog_gm_f.w) + @inbounds for k in TC.real_center_indices(grid) + e_kin = TC.kinetic_energy(prog_gm.u[k], prog_gm.v[k], w_c[k]) + e_pot = TC.geopotential(param_set, grid.zc.z[k]) + prog_gm.ρe_tot[k] = ρ_c[k] * TD.total_energy(param_set, ts_gm[k], e_kin, e_pot) + prog_gm.ρq_tot[k] = ρ_c[k] * aux_gm.q_tot[k] + end + return nothing +end + function assign_thermo_aux!(state, grid, moisture_model, param_set) aux_gm = TC.center_aux_grid_mean(state) prog_gm = TC.center_prog_grid_mean(state) ts_gm = TC.center_aux_grid_mean(state).ts ρ_c = prog_gm.ρ @inbounds for k in TC.real_center_indices(grid) - aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k] - aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k] ts = ts_gm[k] aux_gm.q_liq[k] = TD.liquid_specific_humidity(param_set, ts) aux_gm.q_ice[k] = TD.ice_specific_humidity(param_set, ts) @@ -311,7 +347,7 @@ function ∑tendencies!(tendencies::FV, prog::FV, params::NT, t::Real) where {NT state = TC.State(prog, aux, tendencies) - set_thermo_state!(state, grid, edmf.moisture_model, param_set) + set_thermo_state_peq!(state, grid, edmf.moisture_model, param_set) # TODO: where should this live? aux_gm = TC.center_aux_grid_mean(state) @@ -365,14 +401,19 @@ function compute_gm_tendencies!( force::TC.ForcingBase, param_set::APS, ) + Ic = CCO.InterpolateF2C() + R_d = CPP.R_d(param_set) + T_0 = CPP.T_0(param_set) + Lv_0 = CPP.LH_v0(param_set) tendencies_gm = TC.center_tendencies_grid_mean(state) kc_toa = TC.kc_top_of_atmos(grid) kf_surf = TC.kf_surface(grid) FT = eltype(grid) prog_gm = TC.center_prog_grid_mean(state) + prog_gm_f = TC.face_prog_grid_mean(state) aux_gm = TC.center_aux_grid_mean(state) aux_gm_f = TC.face_aux_grid_mean(state) - ∇θ_liq_ice_gm = TC.center_aux_grid_mean(state).∇θ_liq_ice_gm + ∇MSE_gm = TC.center_aux_grid_mean(state).∇MSE_gm ∇q_tot_gm = TC.center_aux_grid_mean(state).∇q_tot_gm aux_en = TC.center_aux_environment(state) aux_en_f = TC.face_aux_environment(state) @@ -384,14 +425,20 @@ function compute_gm_tendencies!( aux_tc = TC.center_aux_turbconv(state) ts_gm = TC.center_aux_grid_mean(state).ts - θ_liq_ice_gm_toa = aux_gm.θ_liq_ice[kc_toa] - q_tot_gm_toa = aux_gm.q_tot[kc_toa] - RBθ = CCO.RightBiasedC2F(; top = CCO.SetValue(θ_liq_ice_gm_toa)) + e_kin = copy(prog_gm.ρe_tot) + w_c = copy(prog_gm.ρe_tot) + h_tot_gm = copy(prog_gm.ρe_tot) + @. w_c = Ic(FT(0) + prog_gm_f.w) + @. h_tot_gm = TC.anelastic_total_enthalpy(param_set, prog_gm.ρe_tot / ρ_c, ts_gm) + @. e_kin = TC.kinetic_energy(prog_gm.u, prog_gm.v, w_c) + MSE_gm_toa = h_tot_gm[kc_toa] - e_kin[kc_toa] + q_tot_gm_toa = prog_gm.ρq_tot[kc_toa] / ρ_c[kc_toa] + RBe = CCO.RightBiasedC2F(; top = CCO.SetValue(MSE_gm_toa)) RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa)) wvec = CC.Geometry.WVector ∇c = CCO.DivergenceF2C() - @. ∇θ_liq_ice_gm = ∇c(wvec(RBθ(aux_gm.θ_liq_ice))) - @. ∇q_tot_gm = ∇c(wvec(RBq(aux_gm.q_tot))) + @. ∇MSE_gm = ∇c(wvec(RBe(h_tot_gm - e_kin))) + @. ∇q_tot_gm = ∇c(wvec(RBq(prog_gm.ρq_tot / ρ_c))) if edmf.moisture_model isa TC.NonEquilibriumMoisture ∇q_liq_gm = TC.center_aux_grid_mean(state).∇q_liq_gm @@ -406,7 +453,12 @@ function compute_gm_tendencies!( # Apply forcing and radiation @inbounds for k in TC.real_center_indices(grid) - Π = TD.exner(param_set, ts_gm[k]) + cp_m = TD.cp_m(param_set, ts_gm[k]) + cp_v = CPP.cp_v(param_set) + cv_v = CPP.cv_v(param_set) + R_v = CPP.R_v(param_set) + cv_m = TD.cv_m(param_set, ts_gm[k]) + h_v = cv_v * (aux_gm.T[k] - T_0) + Lv_0 - R_v * T_0 # Coriolis if force.apply_coriolis @@ -414,20 +466,20 @@ function compute_gm_tendencies!( tendencies_gm.v[k] += force.coriolis_param * (aux_gm.ug[k] - prog_gm.u[k]) end # LS Subsidence - tendencies_gm.ρθ_liq_ice[k] -= ρ_c[k] * ∇θ_liq_ice_gm[k] * aux_gm.subsidence[k] - tendencies_gm.ρq_tot[k] -= ρ_c[k] * ∇q_tot_gm[k] * aux_gm.subsidence[k] + tendencies_gm.ρe_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇MSE_gm[k] + tendencies_gm.ρq_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇q_tot_gm[k] if edmf.moisture_model isa TC.NonEquilibriumMoisture tendencies_gm.q_liq[k] -= ∇q_liq_gm[k] * aux_gm.subsidence[k] tendencies_gm.q_ice[k] -= ∇q_ice_gm[k] * aux_gm.subsidence[k] end # Radiation if TC.rad_type(radiation) <: Union{TC.RadiationDYCOMS_RF01, TC.RadiationLES, TC.RadiationTRMM_LBA} - tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_rad[k] / Π + tendencies_gm.ρe_tot[k] += ρ_c[k] * cv_m * aux_gm.dTdt_rad[k] end # LS advection tendencies_gm.ρq_tot[k] += ρ_c[k] * aux_gm.dqtdt_hadv[k] if !(TC.force_type(force) <: TC.ForcingDYCOMS_RF01) - tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_hadv[k] / Π + tendencies_gm.ρe_tot[k] += ρ_c[k] * (cp_m * aux_gm.dTdt_hadv[k] + h_v * aux_gm.dqtdt_hadv[k]) end if edmf.moisture_model isa TC.NonEquilibriumMoisture tendencies_gm.q_liq[k] += aux_gm.dqldt[k] @@ -436,21 +488,22 @@ function compute_gm_tendencies!( # LES specific forcings if TC.force_type(force) <: TC.ForcingLES - H_fluc = aux_gm.dTdt_fluc[k] / Π + T_fluc = aux_gm.dTdt_fluc[k] gm_U_nudge_k = (aux_gm.u_nudge[k] - prog_gm.u[k]) / force.wind_nudge_τᵣ gm_V_nudge_k = (aux_gm.v_nudge[k] - prog_gm.v[k]) / force.wind_nudge_τᵣ Γᵣ = TC.compute_les_Γᵣ(grid.zc[k], force.scalar_nudge_τᵣ, force.scalar_nudge_zᵢ, force.scalar_nudge_zᵣ) - gm_H_nudge_k = Γᵣ * (aux_gm.H_nudge[k] - aux_gm.θ_liq_ice[k]) + gm_T_nudge_k = Γᵣ * (aux_gm.T_nudge[k] - aux_gm.T[k]) gm_q_tot_nudge_k = Γᵣ * (aux_gm.qt_nudge[k] - aux_gm.q_tot[k]) if edmf.moisture_model isa TC.NonEquilibriumMoisture gm_q_liq_nudge_k = Γᵣ * (aux_gm.ql_nudge[k] - prog_gm.q_liq[k]) gm_q_ice_nudge_k = Γᵣ * (aux_gm.qi_nudge[k] - prog_gm.q_ice[k]) end - tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * (gm_H_nudge_k + H_fluc) tendencies_gm.ρq_tot[k] += ρ_c[k] * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k]) + tendencies_gm.ρe_tot[k] += + ρ_c[k] * (cv_m * (gm_T_nudge_k + T_fluc) + h_v * (gm_q_tot_nudge_k + aux_gm.dqtdt_fluc[k])) tendencies_gm.u[k] += gm_U_nudge_k tendencies_gm.v[k] += gm_V_nudge_k if edmf.moisture_model isa TC.NonEquilibriumMoisture @@ -466,38 +519,38 @@ function compute_gm_tendencies!( aux_en.qt_tendency_precip_formation[k] + aux_tc.qt_tendency_precip_sinks[k] ) - tendencies_gm.ρθ_liq_ice[k] += + + tendencies_gm.ρe_tot[k] += ρ_c[k] * ( - aux_bulk.θ_liq_ice_tendency_precip_formation[k] + - aux_en.θ_liq_ice_tendency_precip_formation[k] + - aux_tc.θ_liq_ice_tendency_precip_sinks[k] + aux_bulk.e_tot_tendency_precip_formation[k] + + aux_en.e_tot_tendency_precip_formation[k] + + aux_tc.e_tot_tendency_precip_sinks[k] ) + if edmf.moisture_model isa TC.NonEquilibriumMoisture tendencies_gm.q_liq[k] += aux_bulk.ql_tendency_precip_formation[k] + aux_en.ql_tendency_precip_formation[k] tendencies_gm.q_ice[k] += aux_bulk.qi_tendency_precip_formation[k] + aux_en.qi_tendency_precip_formation[k] end end - TC.compute_sgs_flux!(edmf, grid, state, surf) - - sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice + TC.compute_sgs_flux!(edmf, grid, state, surf, param_set) + sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot sgs_flux_u = aux_gm_f.sgs_flux_u sgs_flux_v = aux_gm_f.sgs_flux_v # apply surface BC as SGS flux at lowest level - sgs_flux_θ_liq_ice[kf_surf] = surf.ρθ_liq_ice_flux + sgs_flux_h_tot[kf_surf] = surf.ρe_tot_flux sgs_flux_q_tot[kf_surf] = surf.ρq_tot_flux sgs_flux_u[kf_surf] = surf.ρu_flux sgs_flux_v[kf_surf] = surf.ρv_flux - tends_ρθ_liq_ice = tendencies_gm.ρθ_liq_ice + tends_ρe_tot = tendencies_gm.ρe_tot tends_ρq_tot = tendencies_gm.ρq_tot tends_u = tendencies_gm.u tends_v = tendencies_gm.v ∇sgs = CCO.DivergenceF2C() - - @. tends_ρθ_liq_ice += -∇sgs(wvec(sgs_flux_θ_liq_ice)) + @. tends_ρe_tot += -∇sgs(wvec(sgs_flux_h_tot)) @. tends_ρq_tot += -∇sgs(wvec(sgs_flux_q_tot)) @. tends_u += -∇sgs(wvec(sgs_flux_u)) / ρ_c @. tends_v += -∇sgs(wvec(sgs_flux_v)) / ρ_c diff --git a/driver/generate_namelist.jl b/driver/generate_namelist.jl index 15f0d35716..01264055a1 100644 --- a/driver/generate_namelist.jl +++ b/driver/generate_namelist.jl @@ -79,7 +79,6 @@ function default_namelist( namelist_defaults["logging"] = Dict() namelist_defaults["logging"]["truncate_stack_trace"] = truncate_stack_trace - namelist_defaults["energy_var"] = "rhotheta" namelist_defaults["turbulence"] = Dict() namelist_defaults["turbulence"]["EDMF_PrognosticTKE"] = Dict() @@ -351,7 +350,7 @@ function TRMM_LBA(namelist_defaults) namelist["meta"]["casename"] = "TRMM_LBA" - namelist["grid"]["nz"] = 80 + namelist["grid"]["nz"] = 82 namelist["grid"]["dz"] = 200 namelist["time_stepping"]["adapt_dt"] = true diff --git a/driver/initial_conditions.jl b/driver/initial_conditions.jl index 6eaf5aa7c2..7b0c5e15d1 100644 --- a/driver/initial_conditions.jl +++ b/driver/initial_conditions.jl @@ -83,7 +83,7 @@ function initialize_updrafts(edmf, grid, state, surf) 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 + prog_up[i].ρae_tot[k] = 0 end if edmf.entr_closure isa TC.PrognosticNoisyRelaxationProcess @. prog_up[i].ε_nondim = 0 @@ -137,14 +137,23 @@ function initialize_updrafts_DryBubble(edmf, grid, state) # for now temperature is provided as diagnostics from LES aux_up[i].T[k] = prof_T(z) prog_up[i].ρarea[k] = ρ_0_c[k] * aux_up[i].area[k] - prog_up[i].ρaθ_liq_ice[k] = prog_up[i].ρarea[k] * aux_up[i].θ_liq_ice[k] prog_up[i].ρaq_tot[k] = prog_up[i].ρarea[k] * aux_up[i].q_tot[k] + # thermostate here YAIR-YAIR + e_tot = total_energy(param_set, + grid.zc[k].z, + prog_gm.u[k], + prog_gm.v[k], + prog_gm.w[k], + aux_gm.p[k], + aux_up[i].θ_liq_ice[k], + aux_up[i].q_tot[k]) + prog_up[i].ρae_tot[k] = prog_up[i].ρarea[k] * e_tot else aux_up[i].area[k] = 0.0 aux_up[i].θ_liq_ice[k] = aux_gm.θ_liq_ice[k] aux_up[i].T[k] = aux_gm.T[k] prog_up[i].ρarea[k] = 0.0 - prog_up[i].ρaθ_liq_ice[k] = 0.0 + prog_up[i].ρae_tot[k] = 0.0 prog_up[i].ρaq_tot[k] = 0.0 end end diff --git a/driver/main.jl b/driver/main.jl index 87bc5b8b21..aa15b2bff0 100644 --- a/driver/main.jl +++ b/driver/main.jl @@ -185,8 +185,8 @@ function initialize(sim::Simulation1d) FT = eltype(sim.grid) t = FT(0) Cases.initialize_profiles(sim.case, sim.grid, sim.param_set, state) - set_prog_from_aux!(state, sim.edmf.evar) - set_thermo_state!(state, sim.grid, sim.edmf.moisture_model, sim.param_set) + set_thermo_state_pθq!(state, sim.grid, sim.edmf.moisture_model, sim.param_set) + set_grid_mean_from_thermo_state!(sim.param_set, sim.state, sim.grid) assign_thermo_aux!(state, sim.grid, sim.edmf.moisture_model, sim.param_set) Cases.initialize_forcing(sim.case, sim.grid, state, sim.param_set) diff --git a/integration_tests/3dBomex.jl b/integration_tests/3dBomex.jl index 8eb8026b4c..ac333395bb 100644 --- a/integration_tests/3dBomex.jl +++ b/integration_tests/3dBomex.jl @@ -129,8 +129,8 @@ function init_state(edmf, coords, face_coords, hv_center_space, hv_face_space) ρuₕ = CCG.UVVector(FT(0), FT(0)), u = FT(0), v = FT(0), - θ_liq_ice = FT(0), - q_tot = FT(0), + ρe_tot = FT(0), + ρq_tot = FT(0), TC.cent_prognostic_vars_edmf(FT, edmf)..., ), ) diff --git a/integration_tests/cli_options.jl b/integration_tests/cli_options.jl index 57a99f64bb..6419974ebb 100644 --- a/integration_tests/cli_options.jl +++ b/integration_tests/cli_options.jl @@ -41,9 +41,6 @@ function parse_commandline() arg_type = Int "--moisture_model" # Moisture model (equilibrium or non-equilibrium) arg_type = String - "--energy_var" # Energy variable [rho-theta] - arg_type = String - default = "rhotheta" "--precipitation_model" # Precipitation model (None, cutoff or clima_1m) arg_type = String "--precip_fraction_model" # Precipitation model (prescribed or cloud_cover) diff --git a/integration_tests/driver.jl b/integration_tests/driver.jl index cb98737be3..365bc84636 100644 --- a/integration_tests/driver.jl +++ b/integration_tests/driver.jl @@ -43,7 +43,6 @@ overwrite_namelist_map = Dict( "prescribed_precip_frac_value" => (nl, pa, key) -> (nl["microphysics"]["prescribed_precip_frac_value"] = pa[key]), "precip_fraction_limiter" => (nl, pa, key) -> (nl["microphysics"]["precip_fraction_limiter"] = pa[key]), "thermo_covariance_model" => (nl, pa, key) -> (nl["thermodynamics"]["thermo_covariance_model"] = pa[key]), -"energy_var" => (nl, pa, key) -> (nl["energy_var"] = pa[key]), ) no_overwrites = ( "case", # default_namelist already overwrites namelist["meta"]["casename"] diff --git a/post_processing/mse_tables.jl b/post_processing/mse_tables.jl index a3215064ea..408d04455a 100644 --- a/post_processing/mse_tables.jl +++ b/post_processing/mse_tables.jl @@ -5,182 +5,182 @@ all_best_mse = OrderedCollections.OrderedDict() # all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict() -all_best_mse["ARM_SGP"]["qt_mean"] = 0.2521234151375062 -all_best_mse["ARM_SGP"]["updraft_area"] = 327.08864676538116 -all_best_mse["ARM_SGP"]["updraft_w"] = 153.60644490005336 -all_best_mse["ARM_SGP"]["updraft_qt"] = 30.63670608026404 -all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.01251579142894 +all_best_mse["ARM_SGP"]["qt_mean"] = 0.12006325559150749 +all_best_mse["ARM_SGP"]["updraft_area"] = 314.32846258747117 +all_best_mse["ARM_SGP"]["updraft_w"] = 91.47951373482611 +all_best_mse["ARM_SGP"]["updraft_qt"] = 23.569101319280463 +all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.03760524692865 all_best_mse["ARM_SGP"]["u_mean"] = 1.3235797273549681e-5 -all_best_mse["ARM_SGP"]["tke_mean"] = 1101.575341290362 -all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00012552650010396978 -all_best_mse["ARM_SGP"]["ql_mean"] = 211.93337253648957 +all_best_mse["ARM_SGP"]["tke_mean"] = 359.93566013342814 +all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00021704618840129086 +all_best_mse["ARM_SGP"]["ql_mean"] = 187.54092611860383 all_best_mse["ARM_SGP"]["qi_mean"] = "NA" -all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00011819973647039193 -all_best_mse["ARM_SGP"]["Hvar_mean"] = 2195.6201522931256 -all_best_mse["ARM_SGP"]["QTvar_mean"] = 1309.931007883647 +all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00018325900360641501 +all_best_mse["ARM_SGP"]["Hvar_mean"] = 717.6682880462903 +all_best_mse["ARM_SGP"]["QTvar_mean"] = 525.2935827948501 # all_best_mse["Bomex"] = OrderedCollections.OrderedDict() -all_best_mse["Bomex"]["qt_mean"] = 1.0153399034723752e-01 -all_best_mse["Bomex"]["updraft_area"] = 1.2758624591266978e+02 -all_best_mse["Bomex"]["updraft_w"] = 1.8472338989587136e+01 -all_best_mse["Bomex"]["updraft_qt"] = 6.9599379424859746e+00 -all_best_mse["Bomex"]["updraft_thetal"] = 6.9794466362390068e+01 -all_best_mse["Bomex"]["v_mean"] = 6.2743623073176813e+01 -all_best_mse["Bomex"]["u_mean"] = 2.5579199210944609e-01 -all_best_mse["Bomex"]["tke_mean"] = 7.2838020933193349e+01 -all_best_mse["Bomex"]["temperature_mean"] = 4.1584662439429652e-05 -all_best_mse["Bomex"]["ql_mean"] = 1.0188759556276752e+01 +all_best_mse["Bomex"]["qt_mean"] = 0.09405733846274923 +all_best_mse["Bomex"]["updraft_area"] = 130.63846811594675 +all_best_mse["Bomex"]["updraft_w"] = 15.13887832382455 +all_best_mse["Bomex"]["updraft_qt"] = 6.807142999113192 +all_best_mse["Bomex"]["updraft_thetal"] = 69.76068508427221 +all_best_mse["Bomex"]["v_mean"] = 60.30123223847712 +all_best_mse["Bomex"]["u_mean"] = 0.24949204675388187 +all_best_mse["Bomex"]["tke_mean"] = 85.37822528752113 +all_best_mse["Bomex"]["temperature_mean"] = 4.086274893991192e-5 +all_best_mse["Bomex"]["ql_mean"] = 8.564592215956198 all_best_mse["Bomex"]["qi_mean"] = "NA" -all_best_mse["Bomex"]["thetal_mean"] = 4.2183773454122099e-05 -all_best_mse["Bomex"]["Hvar_mean"] = 2.7841209080224403e+03 -all_best_mse["Bomex"]["QTvar_mean"] = 1.0220020765444779e+03 +all_best_mse["Bomex"]["thetal_mean"] = 4.239544214918139e-5 +all_best_mse["Bomex"]["Hvar_mean"] = 249.02558404984347 +all_best_mse["Bomex"]["QTvar_mean"] = 91.07364422285154 # all_best_mse["DryBubble"] = OrderedCollections.OrderedDict() -all_best_mse["DryBubble"]["updraft_area"] = 1.8868758686558362e-22 -all_best_mse["DryBubble"]["updraft_w"] = 8.03826701999901e-22 -all_best_mse["DryBubble"]["updraft_thetal"] = 5.318251222730423e-29 -all_best_mse["DryBubble"]["u_mean"] = 2.7072625790928257e-21 -all_best_mse["DryBubble"]["tke_mean"] = 1.0152784228781806e-21 -all_best_mse["DryBubble"]["temperature_mean"] = 1.9514588075231292e-29 -all_best_mse["DryBubble"]["thetal_mean"] = 1.3784349509741648e-29 -all_best_mse["DryBubble"]["Hvar_mean"] = 1.287385151761912e-21 +all_best_mse["DryBubble"]["updraft_area"] = 0.5316757604544952 +all_best_mse["DryBubble"]["updraft_w"] = 5.981525336075205 +all_best_mse["DryBubble"]["updraft_thetal"] = 8.536255338456069e-8 +all_best_mse["DryBubble"]["u_mean"] = 0.0035498203156191194 +all_best_mse["DryBubble"]["tke_mean"] = 91.21039171740239 +all_best_mse["DryBubble"]["temperature_mean"] = 7.958413884972341e-7 +all_best_mse["DryBubble"]["thetal_mean"] = 6.597678094123517e-7 +all_best_mse["DryBubble"]["Hvar_mean"] = 48.74500884217395 # all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.018639933430779777 -all_best_mse["DYCOMS_RF01"]["ql_mean"] = 5.334566062625586 -all_best_mse["DYCOMS_RF01"]["updraft_area"] = 29.95179703088047 -all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.289899088761302 -all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.141607789809973 -all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.188242922404484 -all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.0027609048314120445 -all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.08875414330146213 -all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.124943747027338 -all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 7.135714470964393e-5 -all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 7.174069491523632e-5 -all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1287.4864740380342 -all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 513.0844806409915 +all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.04486808349954573 +all_best_mse["DYCOMS_RF01"]["ql_mean"] = 3.190495182574239 +all_best_mse["DYCOMS_RF01"]["updraft_area"] = 29.078417892777647 +all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.298568952382179 +all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.412251484518239 +all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.21605772582681 +all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.006908825603208253 +all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.09678517985175016 +all_best_mse["DYCOMS_RF01"]["tke_mean"] = 24.050865361974907 +all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 1.6702057064174004e-5 +all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 1.7006767937397086e-5 +all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1289.6890825161463 +all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 514.8654445474223 # all_best_mse["DYCOMS_RF02"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF02"]["qt_mean"] = 0.06464355324743704 -all_best_mse["DYCOMS_RF02"]["ql_mean"] = 3.649974401217585 -all_best_mse["DYCOMS_RF02"]["qr_mean"] = 7.6271118150093065 -all_best_mse["DYCOMS_RF02"]["updraft_area"] = 27.753115487782903 -all_best_mse["DYCOMS_RF02"]["updraft_w"] = 6.710130525850158 -all_best_mse["DYCOMS_RF02"]["updraft_qt"] = 5.0625970497522275 -all_best_mse["DYCOMS_RF02"]["updraft_thetal"] = 40.549664544151085 -all_best_mse["DYCOMS_RF02"]["v_mean"] = 42.870288538109754 -all_best_mse["DYCOMS_RF02"]["u_mean"] = 20.17393420658208 -all_best_mse["DYCOMS_RF02"]["tke_mean"] = 21.409525048284575 -all_best_mse["DYCOMS_RF02"]["temperature_mean"] = 2.7353259967934458e-5 -all_best_mse["DYCOMS_RF02"]["thetal_mean"] = 3.009403591053178e-5 -all_best_mse["DYCOMS_RF02"]["Hvar_mean"] = 1465.0320663327811 -all_best_mse["DYCOMS_RF02"]["QTvar_mean"] = 514.653226995155 +all_best_mse["DYCOMS_RF02"]["qt_mean"] = 0.08553109432458665 +all_best_mse["DYCOMS_RF02"]["ql_mean"] = 15.125322376624721 +all_best_mse["DYCOMS_RF02"]["qr_mean"] = 19.237402292640738 +all_best_mse["DYCOMS_RF02"]["updraft_area"] = 26.218032303143875 +all_best_mse["DYCOMS_RF02"]["updraft_w"] = 4.0948700125828585 +all_best_mse["DYCOMS_RF02"]["updraft_qt"] = 5.130703137571791 +all_best_mse["DYCOMS_RF02"]["updraft_thetal"] = 40.580341324888224 +all_best_mse["DYCOMS_RF02"]["v_mean"] = 42.81530578847082 +all_best_mse["DYCOMS_RF02"]["u_mean"] = 20.255222400273368 +all_best_mse["DYCOMS_RF02"]["tke_mean"] = 25.96213432675999 +all_best_mse["DYCOMS_RF02"]["temperature_mean"] = 8.623113138084161e-5 +all_best_mse["DYCOMS_RF02"]["thetal_mean"] = 9.839308804421931e-5 +all_best_mse["DYCOMS_RF02"]["Hvar_mean"] = 1465.2952200901275 +all_best_mse["DYCOMS_RF02"]["QTvar_mean"] = 568.0513098636278 # all_best_mse["GABLS"] = OrderedCollections.OrderedDict() -all_best_mse["GABLS"]["updraft_thetal"] = 0.0 -all_best_mse["GABLS"]["v_mean"] = 1.2389165526931998e-30 -all_best_mse["GABLS"]["u_mean"] = 1.735310521433882e-32 -all_best_mse["GABLS"]["tke_mean"] = 1.3423807154074417e-30 -all_best_mse["GABLS"]["temperature_mean"] = 0.0 -all_best_mse["GABLS"]["thetal_mean"] = 0.0 -all_best_mse["GABLS"]["Hvar_mean"] = 1.857573961661239e-30 +all_best_mse["GABLS"]["updraft_thetal"] = 7.408095686075098e-6 +all_best_mse["GABLS"]["v_mean"] = 0.1818306066384794 +all_best_mse["GABLS"]["u_mean"] = 0.03450365338040005 +all_best_mse["GABLS"]["tke_mean"] = 1.3091231702072166 +all_best_mse["GABLS"]["temperature_mean"] = 7.400927512973925e-6 +all_best_mse["GABLS"]["thetal_mean"] = 7.378921499105522e-6 +all_best_mse["GABLS"]["Hvar_mean"] = 4.4207336171941725 # all_best_mse["GATE_III"] = OrderedCollections.OrderedDict() -all_best_mse["GATE_III"]["qt_mean"] = 3.014277474867634e-16 -all_best_mse["GATE_III"]["updraft_area"] = 1.1253246655431564e-11 -all_best_mse["GATE_III"]["updraft_w"] = 1.9487118831294562e-14 -all_best_mse["GATE_III"]["updraft_qt"] = 2.230160680670979e-16 -all_best_mse["GATE_III"]["updraft_thetal"] = 5.411051140073406e-19 -all_best_mse["GATE_III"]["u_mean"] = 1.746776715939859e-18 -all_best_mse["GATE_III"]["tke_mean"] = 2.2067260004431243e-15 -all_best_mse["GATE_III"]["temperature_mean"] = 5.591331696450916e-19 -all_best_mse["GATE_III"]["ql_mean"] = 1.873345898184791e-13 -all_best_mse["GATE_III"]["qi_mean"] = 6.354265932615595e-17 -all_best_mse["GATE_III"]["qr_mean"] = 5.299542121227325e-12 -all_best_mse["GATE_III"]["qs_mean"] = 4.898284583069764e-14 -all_best_mse["GATE_III"]["thetal_mean"] = 3.0588913651700388e-19 -all_best_mse["GATE_III"]["Hvar_mean"] = 4.114722666257766e-14 -all_best_mse["GATE_III"]["QTvar_mean"] = 1.7233402926254475e-12 +all_best_mse["GATE_III"]["qt_mean"] = 0.04419790332160301 +all_best_mse["GATE_III"]["updraft_area"] = 272.0097701640771 +all_best_mse["GATE_III"]["updraft_w"] = 59.16741846405505 +all_best_mse["GATE_III"]["updraft_qt"] = 13.266501806585472 +all_best_mse["GATE_III"]["updraft_thetal"] = 0.009902680342795536 +all_best_mse["GATE_III"]["u_mean"] = 0.0007375563064614529 +all_best_mse["GATE_III"]["tke_mean"] = 0.3294752952920766 +all_best_mse["GATE_III"]["temperature_mean"] = 0.00035662247063435903 +all_best_mse["GATE_III"]["ql_mean"] = 4.9269634519113215 +all_best_mse["GATE_III"]["qi_mean"] = 0.3109701769106098 +all_best_mse["GATE_III"]["qr_mean"] = 0.0007233473703977575 +all_best_mse["GATE_III"]["qs_mean"] = 0.0008387348095542681 +all_best_mse["GATE_III"]["thetal_mean"] = 0.0002911936656412958 +all_best_mse["GATE_III"]["Hvar_mean"] = 1.924837141665807 +all_best_mse["GATE_III"]["QTvar_mean"] = 13.050497732527838 # all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict() -all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 1.1882499692883095e-6 -all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.0023452730980281673 -all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.0005815686348367982 -all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.0008199589201687151 -all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.0005255788935356413 -all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 3.024541418498564e-7 -all_best_mse["life_cycle_Tan2018"]["v_mean"] = 1.7665923350882435e-5 -all_best_mse["life_cycle_Tan2018"]["u_mean"] = 5.042443328446725e-8 -all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 5.74141900039723e-5 -all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 5.768060410772945e-10 -all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 5.692815552156728e-10 -all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 55.57118312521745 -all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 30.651261631568822 +all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0017639570580040651 +all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.46015254217329166 +all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 2.603159777994327 +all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.4096760216983021 +all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.006109920991945191 +all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 4.7456194496877675e-6 +all_best_mse["life_cycle_Tan2018"]["v_mean"] = 1.1080045784688166 +all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.008112379084763118 +all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 11.406586948376232 +all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 2.282718590221987e-6 +all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 2.3239118547946528e-6 +all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 1133.91225246113 +all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 461.82364949012646 # all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict() -all_best_mse["Nieuwstadt"]["updraft_area"] = 99.95518948417036 -all_best_mse["Nieuwstadt"]["updraft_w"] = 14.234210108108666 -all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.6059354414369 -all_best_mse["Nieuwstadt"]["u_mean"] = 13.553312452894682 -all_best_mse["Nieuwstadt"]["tke_mean"] = 283.45279336077976 -all_best_mse["Nieuwstadt"]["temperature_mean"] = 1.1069166228395703e-5 -all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.1163599173153909e-5 -all_best_mse["Nieuwstadt"]["Hvar_mean"] = 719.387813253075 +all_best_mse["Nieuwstadt"]["updraft_area"] = 115.036365375561 +all_best_mse["Nieuwstadt"]["updraft_w"] = 12.107567820871138 +all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.61124792810575 +all_best_mse["Nieuwstadt"]["u_mean"] = 13.701404730849868 +all_best_mse["Nieuwstadt"]["tke_mean"] = 108.84673651124068 +all_best_mse["Nieuwstadt"]["temperature_mean"] = 5.599701169886529e-5 +all_best_mse["Nieuwstadt"]["thetal_mean"] = 5.441653507544631e-5 +all_best_mse["Nieuwstadt"]["Hvar_mean"] = 662.5243482385383 # all_best_mse["Rico"] = OrderedCollections.OrderedDict() -all_best_mse["Rico"]["qt_mean"] = 0.8732617812127805 -all_best_mse["Rico"]["updraft_area"] = 478.6103900655811 -all_best_mse["Rico"]["updraft_w"] = 75.0403670264857 -all_best_mse["Rico"]["updraft_qt"] = 16.277480831993625 -all_best_mse["Rico"]["updraft_thetal"] = 133.87298592651348 -all_best_mse["Rico"]["v_mean"] = 0.5566971415416723 -all_best_mse["Rico"]["u_mean"] = 0.37707238163166606 -all_best_mse["Rico"]["tke_mean"] = 144.6261226338786 -all_best_mse["Rico"]["temperature_mean"] = 0.0005406092778061022 -all_best_mse["Rico"]["ql_mean"] = 81.78345675815379 +all_best_mse["Rico"]["qt_mean"] = 2.152843322568303 +all_best_mse["Rico"]["updraft_area"] = 489.56441878463653 +all_best_mse["Rico"]["updraft_w"] = 100.76481261815032 +all_best_mse["Rico"]["updraft_qt"] = 11.066192478467672 +all_best_mse["Rico"]["updraft_thetal"] = 133.4380515938658 +all_best_mse["Rico"]["v_mean"] = 0.4398379732285477 +all_best_mse["Rico"]["u_mean"] = 0.38908260369328856 +all_best_mse["Rico"]["tke_mean"] = 346.7111244905397 +all_best_mse["Rico"]["temperature_mean"] = 0.0013714577808486821 +all_best_mse["Rico"]["ql_mean"] = 130.0201032299897 all_best_mse["Rico"]["qi_mean"] = "NA" -all_best_mse["Rico"]["qr_mean"] = 751.7749016991701 -all_best_mse["Rico"]["thetal_mean"] = 0.0005282588047823552 -all_best_mse["Rico"]["Hvar_mean"] = 41772.53419410204 -all_best_mse["Rico"]["QTvar_mean"] = 9689.617792098537 +all_best_mse["Rico"]["qr_mean"] = 760.4335012745742 +all_best_mse["Rico"]["thetal_mean"] = 0.0013538937214399528 +all_best_mse["Rico"]["Hvar_mean"] = 88386.47582661136 +all_best_mse["Rico"]["QTvar_mean"] = 14382.447523416013 # all_best_mse["Soares"] = OrderedCollections.OrderedDict() -all_best_mse["Soares"]["qt_mean"] = 0.1427328325123903 -all_best_mse["Soares"]["updraft_area"] = 97.54831238309664 -all_best_mse["Soares"]["updraft_w"] = 13.105743423554053 -all_best_mse["Soares"]["updraft_qt"] = 23.64110329423397 -all_best_mse["Soares"]["updraft_thetal"] = 65.72138701636423 -all_best_mse["Soares"]["u_mean"] = 93.89558704063909 -all_best_mse["Soares"]["tke_mean"] = 216.45184688981107 -all_best_mse["Soares"]["temperature_mean"] = 1.2381752395765417e-5 -all_best_mse["Soares"]["thetal_mean"] = 1.206493134678075e-5 -all_best_mse["Soares"]["Hvar_mean"] = 680.944052028659 +all_best_mse["Soares"]["qt_mean"] = 0.01484511949623237 +all_best_mse["Soares"]["updraft_area"] = 123.23470646440379 +all_best_mse["Soares"]["updraft_w"] = 8.915136618537609 +all_best_mse["Soares"]["updraft_qt"] = 22.908486551242916 +all_best_mse["Soares"]["updraft_thetal"] = 65.72524659018117 +all_best_mse["Soares"]["u_mean"] = 93.7653486916643 +all_best_mse["Soares"]["tke_mean"] = 72.36680557316008 +all_best_mse["Soares"]["temperature_mean"] = 6.252681841322685e-5 +all_best_mse["Soares"]["thetal_mean"] = 5.6552441218968035e-5 +all_best_mse["Soares"]["Hvar_mean"] = 553.8263398401125 # all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict() -all_best_mse["TRMM_LBA"]["qt_mean"] = 2.0998380118645743 -all_best_mse["TRMM_LBA"]["updraft_area"] = 1268.38083109634 -all_best_mse["TRMM_LBA"]["updraft_w"] = 15532.47743794566 -all_best_mse["TRMM_LBA"]["updraft_qt"] = 285.106942667409 -all_best_mse["TRMM_LBA"]["updraft_thetal"] = 490.0926293239174 -all_best_mse["TRMM_LBA"]["v_mean"] = 67.45554061935081 -all_best_mse["TRMM_LBA"]["u_mean"] = 27.227766500062224 -all_best_mse["TRMM_LBA"]["tke_mean"] = 135664.71364864224 -all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0006728449387002054 -all_best_mse["TRMM_LBA"]["ql_mean"] = 235473.9870620743 +all_best_mse["TRMM_LBA"]["qt_mean"] = 2.8755703383637883 +all_best_mse["TRMM_LBA"]["updraft_area"] = 1316.2743930819443 +all_best_mse["TRMM_LBA"]["updraft_w"] = 16743.16021753859 +all_best_mse["TRMM_LBA"]["updraft_qt"] = 289.89413560750694 +all_best_mse["TRMM_LBA"]["updraft_thetal"] = 534.4330404476216 +all_best_mse["TRMM_LBA"]["v_mean"] = 64.62771685999486 +all_best_mse["TRMM_LBA"]["u_mean"] = 28.84351327921325 +all_best_mse["TRMM_LBA"]["tke_mean"] = 1.2912970947081042e6 +all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0011264127513910369 +all_best_mse["TRMM_LBA"]["ql_mean"] = 215664.5497128382 all_best_mse["TRMM_LBA"]["qi_mean"] = "NA" all_best_mse["TRMM_LBA"]["qr_mean"] = "NA" all_best_mse["TRMM_LBA"]["qs_mean"] = "NA" -all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.004140237903002215 -all_best_mse["TRMM_LBA"]["Hvar_mean"] = 430415.19992669823 -all_best_mse["TRMM_LBA"]["QTvar_mean"] = 6659.921911738702 +all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.0056394176808409165 +all_best_mse["TRMM_LBA"]["Hvar_mean"] = 525094.9836790846 +all_best_mse["TRMM_LBA"]["QTvar_mean"] = 6430.568392018743 # all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict() -all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.5656551128069336 -all_best_mse["LES_driven_SCM"]["v_mean"] = 1.3786065803844898 -all_best_mse["LES_driven_SCM"]["u_mean"] = 0.47643802037799565 -all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.001283074413616987 -all_best_mse["LES_driven_SCM"]["ql_mean"] = 51725.39021346529 -all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0015164076617609533 +all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.5085239201936997 +all_best_mse["LES_driven_SCM"]["v_mean"] = 1.37648422221553 +all_best_mse["LES_driven_SCM"]["u_mean"] = 0.47500950600491726 +all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0014752952109381328 +all_best_mse["LES_driven_SCM"]["ql_mean"] = 32773.53878280255 +all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.001717394645557076 # ################################# ################################# diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index e0f5c10b44..450d56311a 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -30,6 +30,7 @@ function microphysics( prog_pr.q_sno[k], aux_en.area[k], ρ_c[k], + grid.zc[k].z, Δt, ts, precip_fraction, @@ -51,10 +52,9 @@ function microphysics( end # update_env_precip_tendencies - # TODO: move qt_tendency_precip_formation and θ_liq_ice_tendency_precip_formation - # to diagnostics + # TODO: move ..._tendency_precip_formation to diagnostics aux_en.qt_tendency_precip_formation[k] = mph.qt_tendency * aux_en.area[k] - aux_en.θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_tendency * aux_en.area[k] + aux_en.e_tot_tendency_precip_formation[k] = mph.e_tot_tendency * aux_en.area[k] if edmf.moisture_model isa NonEquilibriumMoisture aux_en.ql_tendency_precip_formation[k] = mph.ql_tendency * aux_en.area[k] aux_en.qi_tendency_precip_formation[k] = mph.qi_tendency * aux_en.area[k] @@ -68,9 +68,9 @@ end function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt::Real) env_len = 8 - src_len = 8 + src_len = 9 i_ql, i_qi, i_T, i_cf, i_qt_sat, i_qt_unsat, i_T_sat, i_T_unsat = 1:env_len - i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH, i_Sqr, i_Sqs = 1:src_len + i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH, i_Sqr, i_Sqs, i_Se_tot = 1:src_len quadrature_type = en_thermo.quadrature_type quad_order = quadrature_order(en_thermo) @@ -81,7 +81,8 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt: # θl - liquid ice potential temperature # _mean and ′ - subdomain mean and (co)variances # q_rai, q_sno - grid mean precipitation - UnPack.@unpack qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, subdomain_area, q_rai, q_sno, ρ_c, p_c, precip_frac = vars + UnPack.@unpack qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, subdomain_area, q_rai, q_sno, ρ_c, p_c, zc, precip_frac = + vars FT = eltype(ρ_c) @@ -160,8 +161,18 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt: q_ice_en = TD.ice_specific_humidity(param_set, ts) T = TD.air_temperature(param_set, ts) # autoconversion and accretion - mph = - precipitation_formation(param_set, precip_model, q_rai, q_sno, subdomain_area, ρ_c, Δt, ts, precip_frac) + mph = precipitation_formation( + param_set, + precip_model, + q_rai, + q_sno, + subdomain_area, + ρ_c, + zc, + Δt, + ts, + precip_frac, + ) # environmental variables inner_env[i_ql] += q_liq_en * weights[m_h] * sqpi_inv @@ -181,6 +192,7 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt: inner_src[i_Sqr] += mph.qr_tendency * weights[m_h] * sqpi_inv inner_src[i_Sqs] += mph.qs_tendency * weights[m_h] * sqpi_inv inner_src[i_SH] += mph.θ_liq_ice_tendency * weights[m_h] * sqpi_inv + inner_src[i_Se_tot] += mph.e_tot_tendency * weights[m_h] * sqpi_inv inner_src[i_Sqt_H] += mph.qt_tendency * h_hat * weights[m_h] * sqpi_inv inner_src[i_Sqt_qt] += mph.qt_tendency * qt_hat * weights[m_h] * sqpi_inv inner_src[i_SH_H] += mph.θ_liq_ice_tendency * h_hat * weights[m_h] * sqpi_inv @@ -201,6 +213,7 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt: Sqt_qt = outer_src[i_Sqt_qt], Sqt = outer_src[i_Sqt], SH = outer_src[i_SH], + Se_tot = outer_src[i_Se_tot], Sqr = outer_src[i_Sqr], Sqs = outer_src[i_Sqs], ) @@ -272,6 +285,7 @@ function microphysics( ρ_c = ρ_c[k], p_c = p_c[k], precip_frac = precip_fraction, + zc = grid.zc[k].z, ) outer_env, outer_src = quad_loop(en_thermo, precip_model, vars, param_set, Δt) @@ -280,12 +294,12 @@ function microphysics( # update_env_precip_tendencies qt_tendency = outer_src.Sqt θ_liq_ice_tendency = outer_src.SH + e_tot_tendency = outer_src.Se_tot qr_tendency = outer_src.Sqr qs_tendency = outer_src.Sqs - # TODO: move qt_tendency_precip_formation and θ_liq_ice_tendency_precip_formation - # to diagnostics + # TODO: move ..._tendency_precip_formation to diagnostics aux_en.qt_tendency_precip_formation[k] = qt_tendency * aux_en.area[k] - aux_en.θ_liq_ice_tendency_precip_formation[k] = θ_liq_ice_tendency * aux_en.area[k] + aux_en.e_tot_tendency_precip_formation[k] = e_tot_tendency * aux_en.area[k] tendencies_pr.q_rai[k] += qr_tendency * aux_en.area[k] tendencies_pr.q_sno[k] += qs_tendency * aux_en.area[k] @@ -333,16 +347,16 @@ function microphysics( prog_pr.q_sno[k], aux_en.area[k], ρ_c[k], + grid.zc[k].z, Δt, ts, precip_fraction, ) # update_env_precip_tendencies - # TODO: move qt_tendency_precip_formation and θ_liq_ice_tendency_precip_formation - # to diagnostics + # TODO: move ..._tendency_precip_formation to diagnostics aux_en.qt_tendency_precip_formation[k] = mph.qt_tendency * aux_en.area[k] - aux_en.θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_tendency * aux_en.area[k] + aux_en.e_tot_tendency_precip_formation[k] = mph.e_tot_tendency * aux_en.area[k] tendencies_pr.q_rai[k] += mph.qr_tendency * aux_en.area[k] tendencies_pr.q_sno[k] += mph.qs_tendency * aux_en.area[k] diff --git a/src/EDMF_Precipitation.jl b/src/EDMF_Precipitation.jl index aaef422538..a32b83df0d 100644 --- a/src/EDMF_Precipitation.jl +++ b/src/EDMF_Precipitation.jl @@ -121,6 +121,11 @@ function compute_precipitation_sink_tendencies( L_s = TD.latent_heat_sublim(param_set, ts) L_f = TD.latent_heat_fusion(param_set, ts) + I_l = TD.internal_energy_liquid(param_set, ts) + I_i = TD.internal_energy_ice(param_set, ts) + I = TD.internal_energy(param_set, ts) + Φ = geopotential(param_set, grid.zc.z[k]) + α_evp = ICP.microph_scaling(param_set) α_dep_sub = ICP.microph_scaling_dep_sub(param_set) α_melt = ICP.microph_scaling_melt(param_set) @@ -146,12 +151,8 @@ function compute_precipitation_sink_tendencies( tendencies_pr.q_sno[k] += S_qs_sub_dep + S_qs_melt aux_tc.qt_tendency_precip_sinks[k] = -S_qr_evap - S_qs_sub_dep - aux_tc.θ_liq_ice_tendency_precip_sinks[k] = - 1 / Π_m / c_pm * ( - S_qr_evap * (L_v - R_v * T_gm) * (1 + R_m / c_vm) + - S_qs_sub_dep * (L_s - R_v * T_gm) * (1 + R_m / c_vm) + - S_qs_melt * L_f * (1 + R_m / c_vm) - ) + aux_tc.e_tot_tendency_precip_sinks[k] = + -S_qr_evap * (I_l + Φ - I) - S_qs_sub_dep * (I_i + Φ - I) + S_qs_melt * L_f end return nothing end diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl index 46fc47568b..8cd7a6bb95 100644 --- a/src/EDMF_Updrafts.jl +++ b/src/EDMF_Updrafts.jl @@ -88,12 +88,14 @@ function compute_precipitation_formation_tendencies( prog_pr.q_sno[k], aux_up[i].area[k], ρ_c[k], + grid.zc[k].z, Δt, ts_up, precip_fraction, ) aux_up[i].qt_tendency_precip_formation[k] = mph.qt_tendency * aux_up[i].area[k] aux_up[i].θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_tendency * aux_up[i].area[k] + aux_up[i].e_tot_tendency_precip_formation[k] = mph.e_tot_tendency * aux_up[i].area[k] if edmf.moisture_model isa NonEquilibriumMoisture aux_up[i].ql_tendency_precip_formation[k] = mph.ql_tendency * aux_up[i].area[k] aux_up[i].qi_tendency_precip_formation[k] = mph.qi_tendency * aux_up[i].area[k] @@ -104,10 +106,10 @@ function compute_precipitation_formation_tendencies( end # TODO - to be deleted once we sum all tendencies elsewhere @inbounds for k in real_center_indices(grid) - aux_bulk.θ_liq_ice_tendency_precip_formation[k] = 0 + aux_bulk.e_tot_tendency_precip_formation[k] = 0 aux_bulk.qt_tendency_precip_formation[k] = 0 @inbounds for i in 1:N_up - aux_bulk.θ_liq_ice_tendency_precip_formation[k] += aux_up[i].θ_liq_ice_tendency_precip_formation[k] + aux_bulk.e_tot_tendency_precip_formation[k] += aux_up[i].e_tot_tendency_precip_formation[k] aux_bulk.qt_tendency_precip_formation[k] += aux_up[i].qt_tendency_precip_formation[k] end if edmf.moisture_model isa NonEquilibriumMoisture diff --git a/src/EDMF_functions.jl b/src/EDMF_functions.jl index 0ebfd72a13..81922d2c5b 100755 --- a/src/EDMF_functions.jl +++ b/src/EDMF_functions.jl @@ -47,7 +47,7 @@ function compute_turbconv_tendencies!( return nothing end -function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase) +function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase, param_set::APS) N_up = n_updrafts(edmf) tendencies_gm = center_tendencies_grid_mean(state) FT = eltype(grid) @@ -77,17 +77,32 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf a_en = aux_en.area w_en = aux_en_f.w w_gm = prog_gm_f.w - θ_liq_ice_en = aux_en.θ_liq_ice - θ_liq_ice_gm = aux_gm.θ_liq_ice + h_tot_en = copy(a_en) + h_tot_gm = copy(a_en) q_tot_gm = aux_gm.q_tot q_tot_en = aux_en.q_tot a_en_bcs = a_en_boundary_conditions(surf, edmf) Ifae = CCO.InterpolateC2F(; a_en_bcs...) If = CCO.InterpolateC2F(; bottom = CCO.SetValue(FT(0)), top = CCO.SetValue(FT(0))) + Ic = CCO.InterpolateF2C() + # compute total enthalpies + ts_en = center_aux_environment(state).ts + ts_gm = center_aux_grid_mean(state).ts + @. h_tot_gm = anelastic_total_enthalpy(param_set, prog_gm.ρe_tot / ρ_c, ts_gm) + @. h_tot_en = anelastic_total_enthalpy( + param_set, + TD.total_energy( + param_set, + ts_en, + kinetic_energy(prog_gm.u, prog_gm.v, Ic(FT(0) + w_en)), + geopotential(param_set, grid.zc.z), + ), + ts_en, + ) # Compute the mass flux and associated scalar fluxes @. massflux = ρ_f * Ifae(a_en) * (w_en - w_gm) - @. massflux_h = ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm)) + @. massflux_h = ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(h_tot_en) - If(h_tot_gm)) @. massflux_qt = ρ_f * Ifae(a_en) * (w_en - w_gm) * (If(q_tot_en) - If(q_tot_gm)) @inbounds for i in 1:N_up aux_up_f_i = aux_up_f[i] @@ -95,17 +110,33 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf a_up_bcs = a_up_boundary_conditions(surf, edmf, i) Ifau = CCO.InterpolateC2F(; a_up_bcs...) a_up = aux_up[i].area - w_up = aux_up_f[i].w + w_up_i = aux_up_f[i].w θ_liq_ice_up = aux_up_i.θ_liq_ice q_tot_up = aux_up_i.q_tot - @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up - w_gm) + @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up_i - w_gm) # We know that, since W = 0 at z = 0, m = 0 also, and - # therefore θ_liq_ice / q_tot values do not matter - @. massflux_h += ρ_f * (Ifau(a_up) * (w_up - w_gm) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm))) - @. massflux_qt += ρ_f * (Ifau(a_up) * (w_up - w_gm) * (If(q_tot_up) - If(q_tot_gm))) + # therefore h_tot / q_tot values do not matter + ts_up_i = copy(ts_en) + h_tot_up_i = copy(q_tot_up) + thermo_args = if edmf.moisture_model isa EquilibriumMoisture + () + elseif edmf.moisture_model isa NonEquilibriumMoisture + (aux_up[i].q_liq, aux_up[i].q_ice) + end + @. ts_up_i = thermo_state_pθq(param_set, p_c, θ_liq_ice_up, q_tot_up, thermo_args...) + @. h_tot_up_i = anelastic_total_enthalpy( + param_set, + TD.total_energy( + param_set, + ts_up_i, + kinetic_energy(prog_gm.u, prog_gm.v, Ic(FT(0) + w_up_i)), + geopotential(param_set, grid.zc.z), + ), + ts_up_i, + ) + @. massflux_h += ρ_f * (Ifau(a_up) * (w_up_i - w_gm) * (If(h_tot_up_i) - If(h_tot_gm))) + @. massflux_qt += ρ_f * (Ifau(a_up) * (w_up_i - w_gm) * (If(q_tot_up) - If(q_tot_gm))) end - massflux[kf_surf] = 0 - massflux_h[kf_surf] = 0 if edmf.moisture_model isa NonEquilibriumMoisture massflux_en = aux_tc_f.massflux_en @@ -131,6 +162,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf massflux_qi[kf_surf] = 0 end + massflux_h[kf_surf] = 0 + massflux_qt[kf_surf] = 0 + massflux_tendency_h = aux_tc.massflux_tendency_h massflux_tendency_qt = aux_tc.massflux_tendency_qt # Compute the mass flux tendencies @@ -144,12 +178,12 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf diffusive_flux_u = aux_tc_f.diffusive_flux_u diffusive_flux_v = aux_tc_f.diffusive_flux_v - sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice + sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot sgs_flux_u = aux_gm_f.sgs_flux_u sgs_flux_v = aux_gm_f.sgs_flux_v - @. sgs_flux_θ_liq_ice = diffusive_flux_h + massflux_h + @. sgs_flux_h_tot = diffusive_flux_h + massflux_h @. sgs_flux_q_tot = diffusive_flux_qt + massflux_qt @. sgs_flux_u = diffusive_flux_u # + massflux_u @. sgs_flux_v = diffusive_flux_v # + massflux_v @@ -178,7 +212,9 @@ function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, sur FT = eltype(grid) aux_bulk = center_aux_bulk(state) aux_tc_f = face_aux_turbconv(state) + aux_en_f = face_aux_environment(state) aux_en = center_aux_environment(state) + aux_gm = center_aux_grid_mean(state) aux_gm_f = face_aux_grid_mean(state) KM = center_aux_turbconv(state).KM KH = center_aux_turbconv(state).KH @@ -188,29 +224,45 @@ function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, sur a_en = aux_en.area @. aeKM = a_en * KM @. aeKH = a_en * KH + w_en_c = copy(a_en) kc_surf = kc_surface(grid) kc_toa = kc_top_of_atmos(grid) kf_surf = kf_surface(grid) prog_gm = center_prog_grid_mean(state) IfKH = CCO.InterpolateC2F(; bottom = CCO.SetValue(aeKH[kc_surf]), top = CCO.SetValue(aeKH[kc_toa])) IfKM = CCO.InterpolateC2F(; bottom = CCO.SetValue(aeKM[kc_surf]), top = CCO.SetValue(aeKM[kc_toa])) + Ic = CCO.InterpolateF2C() + + h_tot_en = copy(a_en) + ts_gm = center_aux_grid_mean(state).ts + ts_en = center_aux_environment(state).ts + @. h_tot_en = anelastic_total_enthalpy( + param_set, + TD.total_energy( + param_set, + ts_en, + kinetic_energy(prog_gm.u, prog_gm.v, Ic(FT(0) + aux_en_f.w)), + geopotential(param_set, grid.zc.z), + ), + ts_en, + ) @. aux_tc_f.ρ_ae_KH = IfKH(aeKH) * ρ_f @. aux_tc_f.ρ_ae_KM = IfKM(aeKM) * ρ_f aeKHq_tot_bc = -surf.ρq_tot_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf] - aeKHθ_liq_ice_bc = -surf.ρθ_liq_ice_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf] + aeKHh_tot_bc = -surf.ρe_tot_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KH[kf_surf] aeKMu_bc = -surf.ρu_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf] aeKMv_bc = -surf.ρv_flux / a_en[kc_surf] / aux_tc_f.ρ_ae_KM[kf_surf] ∇q_tot_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHq_tot_bc), top = CCO.SetDivergence(FT(0))) - ∇θ_liq_ice_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHθ_liq_ice_bc), top = CCO.SetDivergence(FT(0))) + ∇h_tot_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHh_tot_bc), top = CCO.SetDivergence(FT(0))) ∇u_gm = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKMu_bc), top = CCO.SetDivergence(FT(0))) ∇v_gm = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKMv_bc), top = CCO.SetDivergence(FT(0))) wvec = CC.Geometry.WVector @. aux_tc_f.diffusive_flux_qt = -aux_tc_f.ρ_ae_KH * ∇q_tot_en(wvec(aux_en.q_tot)) - @. aux_tc_f.diffusive_flux_h = -aux_tc_f.ρ_ae_KH * ∇θ_liq_ice_en(wvec(aux_en.θ_liq_ice)) + @. aux_tc_f.diffusive_flux_h = -aux_tc_f.ρ_ae_KH * ∇h_tot_en(wvec(h_tot_en)) @. aux_tc_f.diffusive_flux_u = -aux_tc_f.ρ_ae_KM * ∇u_gm(wvec(prog_gm.u)) @. aux_tc_f.diffusive_flux_v = -aux_tc_f.ρ_ae_KM * ∇v_gm(wvec(prog_gm.v)) @@ -247,7 +299,7 @@ function affect_filter!( ### Filters ### set_edmf_surface_bc(edmf, grid, state, surf, param_set) - filter_updraft_vars(edmf, grid, state, surf) + filter_updraft_vars(edmf, grid, state, surf, param_set) @inbounds for k in real_center_indices(grid) prog_en.ρatke[k] = max(prog_en.ρatke[k], 0.0) @@ -266,21 +318,32 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su N_up = n_updrafts(edmf) kc_surf = kc_surface(grid) kf_surf = kf_surface(grid) + aux_gm = center_aux_grid_mean(state) + aux_gm_f = face_aux_grid_mean(state) + prog_gm = center_prog_grid_mean(state) prog_up = center_prog_updrafts(state) prog_en = center_prog_environment(state) prog_up_f = face_prog_updrafts(state) aux_bulk = center_aux_bulk(state) - aux_en = center_aux_environment(state) - prog_gm = center_prog_grid_mean(state) - aux_gm_f = face_aux_grid_mean(state) + ts_gm = aux_gm.ts + cp = TD.cp_m(param_set, ts_gm[kc_surf]) ρ_c = prog_gm.ρ ρ_f = aux_gm_f.ρ @inbounds for i in 1:N_up - θ_surf = θ_surface_bc(surf, grid, state, edmf, i) + θ_surf = θ_surface_bc(surf, grid, state, edmf, i, param_set) q_surf = q_surface_bc(surf, grid, state, edmf, i) a_surf = area_surface_bc(surf, edmf, i) prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf - prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * θ_surf + e_tot_surf = total_energy(param_set, + grid.zc[kc_surf].z, + prog_gm.u[kc_surf], + prog_gm.v[kc_surf], + prog_gm.w[kc_surf], + aux_gm.p[kc_surf], + θ_surf, + q_surf) + + prog_up[i].ρae_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * e_tot_surf prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf if edmf.moisture_model isa NonEquilibriumMoisture q_liq_surf = FT(0) @@ -291,7 +354,7 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su prog_up_f[i].ρaw[kf_surf] = ρ_f[kf_surf] * w_surface_bc(surf) end - flux1 = surf.ρθ_liq_ice_flux + flux1 = surf.shf / cp flux2 = surf.ρq_tot_flux zLL = grid.zc[kc_surf].z ustar = surf.ustar @@ -346,22 +409,34 @@ end function w_surface_bc(::SurfaceBase{FT})::FT where {FT} return FT(0) end -function θ_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::EDMFModel, i::Int)::FT where {FT} - prog_gm = center_prog_grid_mean(state) +function θ_surface_bc( + surf::SurfaceBase{FT}, + grid::Grid, + state::State, + edmf::EDMFModel, + i::Int, + param_set::APS, +)::FT where {FT} aux_gm = center_aux_grid_mean(state) + prog_gm = center_prog_grid_mean(state) + ρ_c = prog_gm.ρ kc_surf = kc_surface(grid) + ts_gm = aux_gm.ts + c_p = TD.cp_m(param_set, ts_gm[kc_surf]) + UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state) + surf.bflux > 0 || return FT(0) a_total = edmf.surface_area a_ = area_surface_bc(surf, edmf, i) - UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state) - ρθ_liq_ice_flux = surf.ρθ_liq_ice_flux + ρθ_liq_ice_flux = surf.shf / c_p # assuming no ql,qi flux h_var = get_surface_variance(ρθ_liq_ice_flux / ρLL, ρθ_liq_ice_flux / ρLL, ustar, zLL, oblength) surface_scalar_coeff = percentile_bounds_mean_norm(1 - a_total + (i - 1) * a_, 1 - a_total + i * a_, 1000) return aux_gm.θ_liq_ice[kc_surf] + surface_scalar_coeff * sqrt(h_var) end function q_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::EDMFModel, i::Int)::FT where {FT} - prog_gm = center_prog_grid_mean(state) aux_gm = center_aux_grid_mean(state) + prog_gm = center_prog_grid_mean(state) + ρ_c = prog_gm.ρ kc_surf = kc_surface(grid) surf.bflux > 0 || return aux_gm.q_tot[kc_surf] a_total = edmf.surface_area @@ -502,14 +577,15 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param qt_tendency_precip_formation = aux_up_i.qt_tendency_precip_formation tends_ρarea = tendencies_up[i].ρarea - tends_ρaθ_liq_ice = tendencies_up[i].ρaθ_liq_ice + tends_ρae_tot = tendencies_up[i].ρae_tot tends_ρaq_tot = tendencies_up[i].ρaq_tot @. tends_ρarea = -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up))) + (ρ_c * a_up * Ic(w_up) * entr_turb_dyn) - (ρ_c * a_up * Ic(w_up) * detr_turb_dyn) - @. tends_ρaθ_liq_ice = + # enthalpy here YAIR-YAIR + @. tends_ρae_tot = -∇c(wvec(LBF(Ic(w_up) * ρ_c * a_up * θ_liq_ice_up))) + (ρ_c * a_up * Ic(w_up) * entr_turb_dyn * θ_liq_ice_en) - (ρ_c * a_up * Ic(w_up) * detr_turb_dyn * θ_liq_ice_up) + (ρ_c * θ_liq_ice_tendency_precip_formation) @@ -563,7 +639,7 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param end tends_ρarea[kc_surf] = 0 - tends_ρaθ_liq_ice[kc_surf] = 0 + tends_ρae_tot[kc_surf] = 0 tends_ρaq_tot[kc_surf] = 0 end @@ -600,7 +676,7 @@ function compute_up_tendencies!(edmf::EDMFModel, grid::Grid, state::State, param return nothing end -function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase) +function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase, param_set::APS) N_up = n_updrafts(edmf) kc_surf = kc_surface(grid) kf_surf = kf_surface(grid) @@ -620,7 +696,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su @inbounds for i in 1:N_up 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].ρae_tot .= max.(prog_up[i].ρae_tot, 0) prog_up[i].ρaq_tot .= max.(prog_up[i].ρaq_tot, 0) if edmf.entr_closure isa PrognosticNoisyRelaxationProcess @. prog_up[i].ε_nondim = max(prog_up[i].ε_nondim, 0) @@ -676,14 +752,22 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su Ic = CCO.InterpolateF2C() @inbounds for i in 1:N_up @. prog_up[i].ρarea = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρarea) - @. prog_up[i].ρaθ_liq_ice = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρaθ_liq_ice) + @. prog_up[i].ρae_tot = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρae_tot) @. prog_up[i].ρaq_tot = ifelse(Ic(prog_up_f[i].ρaw) <= 0, FT(0), prog_up[i].ρaq_tot) - θ_surf = θ_surface_bc(surf, grid, state, edmf, i) + θ_surf = θ_surface_bc(surf, grid, state, edmf, i, param_set) q_surf = q_surface_bc(surf, grid, state, edmf, i) a_surf = area_surface_bc(surf, edmf, i) prog_up[i].ρarea[kc_surf] = ρ_c[kc_surf] * a_surf - prog_up[i].ρaθ_liq_ice[kc_surf] = prog_up[i].ρarea[kc_surf] * θ_surf prog_up[i].ρaq_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * q_surf + e_tot_surf = total_energy(param_set, + grid.zc[kc_surf].z, + prog_gm.u[kc_surf], + prog_gm.v[kc_surf], + prog_gm.w[kc_surf], + aux_gm.p[kc_surf], + θ_surf, + q_surf) + prog_up[i].ρae_tot[kc_surf] = prog_up[i].ρarea[kc_surf] * e_tot_surf end if edmf.moisture_model isa NonEquilibriumMoisture @inbounds for i in 1:N_up diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl index c801ad8f97..0f2ffe7746 100644 --- a/src/TurbulenceConvection.jl +++ b/src/TurbulenceConvection.jl @@ -98,10 +98,10 @@ function debug_state(state, code_location::String) ###### vars_positive = [ - vec(prog_gm.ρθ_liq_ice), + vec(prog_gm.ρe_tot), vec(prog_gm_f.w), vec(prog_up[1].ρarea), - vec(prog_up[1].ρaθ_liq_ice), + vec(prog_up[1].ρae_tot), vec(prog_up_f[1].ρaw), vec(aux_en.area), vec(aux_en.θ_liq_ice), diff --git a/src/diagnostics.jl b/src/diagnostics.jl index d6fe3b49c0..7807217bcc 100644 --- a/src/diagnostics.jl +++ b/src/diagnostics.jl @@ -103,7 +103,7 @@ function io_dictionary_aux() "updraft_cloud_fraction" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).cloud_fraction), "updraft_qt_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).qt_tendency_precip_formation), - "updraft_thetal_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).θ_liq_ice_tendency_precip_formation), + "updraft_e_tot_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).e_tot_tendency_precip_formation), "massflux_tendency_h" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).massflux_tendency_h), "massflux_tendency_qt" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_turbconv(state).massflux_tendency_qt), diff --git a/src/microphysics_coupling.jl b/src/microphysics_coupling.jl index 1a493bf0bd..77c7fd142a 100644 --- a/src/microphysics_coupling.jl +++ b/src/microphysics_coupling.jl @@ -50,6 +50,7 @@ function precipitation_formation( qs::FT, area::FT, ρ::FT, + z::FT, Δt::Real, ts, precip_fraction, @@ -63,6 +64,7 @@ function precipitation_formation( qr_tendency = FT(0) qs_tendency = FT(0) θ_liq_ice_tendency = FT(0) + e_tot_tendency = FT(0) if area > 0 @@ -72,6 +74,10 @@ function precipitation_formation( c_pm = TD.cp_m(param_set, ts) L_v0 = ICP.LH_v0(param_set) L_s0 = ICP.LH_s0(param_set) + I_l = TD.internal_energy_liquid(param_set, ts) + I_i = TD.internal_energy_ice(param_set, ts) + I = TD.internal_energy(param_set, ts) + Φ = geopotential(param_set, z) if precip_model isa Clima0M qsat = TD.q_vap_saturation(param_set, ts) @@ -85,6 +91,7 @@ function precipitation_formation( ql_tendency += S_qt * λ qi_tendency += S_qt * (1 - λ) θ_liq_ice_tendency -= S_qt / Π_m / c_pm * (L_v0 * λ + L_s0 * (1 - λ)) + e_tot_tendency += (λ * I_l + (1 - λ) * I_i + Φ - I) * S_qt end if precip_model isa Clima1M @@ -110,6 +117,7 @@ function precipitation_formation( ql_tendency += S_qt_rain qi_tendency += S_qt_snow θ_liq_ice_tendency -= 1 / Π_m / c_pm * (L_v0 * S_qt_rain + L_s0 * S_qt_snow) + e_tot_tendency += S_qt_rain * (I_l + Φ - I) + S_qt_snow * (I_i + Φ - I) # accretion cloud water + rain S_qr = min(q.liq / Δt, CM1.accretion(param_set, liq_type, rain_type, q.liq, qr, ρ)) * precip_fraction @@ -117,6 +125,7 @@ function precipitation_formation( qt_tendency -= S_qr ql_tendency -= S_qr θ_liq_ice_tendency += S_qr / Π_m / c_pm * L_v0 + e_tot_tendency -= S_qr * (I_l + Φ - I) # accretion cloud ice + snow S_qs = min(q.ice / Δt, CM1.accretion(param_set, ice_type, snow_type, q.ice, qs, ρ)) * precip_fraction @@ -124,6 +133,7 @@ function precipitation_formation( qt_tendency -= S_qs qi_tendency -= S_qs θ_liq_ice_tendency += S_qs / Π_m / c_pm * L_s0 + e_tot_tendency -= S_qs * (I_i + Φ - I) # sink of cloud water via accretion cloud water + snow S_qt = -min(q.liq / Δt, CM1.accretion(param_set, liq_type, snow_type, q.liq, qs, ρ)) * precip_fraction @@ -132,6 +142,7 @@ function precipitation_formation( qt_tendency += S_qt ql_tendency += S_qt θ_liq_ice_tendency -= S_qt / Π_m / c_pm * Lf * (1 + Rm / c_vm) + e_tot_tendency += S_qt * (I_i + Φ - I) else # snow melts, both cloud water and snow become rain α::FT = c_vl / Lf * (T - T_fr) qt_tendency += S_qt @@ -139,6 +150,7 @@ function precipitation_formation( qs_tendency += S_qt * α qr_tendency -= S_qt * (1 + α) θ_liq_ice_tendency += S_qt / Π_m / c_pm * (Lf * (1 + Rm / c_vm) * α - L_v0) + e_tot_tendency += S_qt * ((1 + α) * I_l - α * I_i + Φ - I) end # sink of cloud ice via accretion cloud ice - rain @@ -150,6 +162,8 @@ function precipitation_formation( qr_tendency += S_qr qs_tendency += -(S_qt + S_qr) θ_liq_ice_tendency -= 1 / Π_m / c_pm * (S_qr * Lf * (1 + Rm / c_vm) + S_qt * L_s0) + e_tot_tendency += S_qt * (I_i + Φ - I) + e_tot_tendency -= S_qr * Lf # accretion rain - snow if T < T_fr @@ -166,7 +180,16 @@ function precipitation_formation( qs_tendency += S_qs qr_tendency -= S_qs θ_liq_ice_tendency += S_qs * Lf / Π_m / c_vm + e_tot_tendency += S_qs * Lf end end - return PrecipFormation{FT}(θ_liq_ice_tendency, qt_tendency, ql_tendency, qi_tendency, qr_tendency, qs_tendency) + return PrecipFormation{FT}( + θ_liq_ice_tendency, + e_tot_tendency, + qt_tendency, + ql_tendency, + qi_tendency, + qr_tendency, + qs_tendency, + ) end diff --git a/src/thermodynamics.jl b/src/thermodynamics.jl index 67c78e2d9d..1a9b13697c 100644 --- a/src/thermodynamics.jl +++ b/src/thermodynamics.jl @@ -9,3 +9,55 @@ function thermo_state_pθq(param_set::APS, p::FT, θ_liq_ice::FT, q_tot::FT, q_l q = TD.PhasePartition(q_tot, q_liq, q_ice) return TD.PhaseNonEquil_pθq(param_set, p, θ_liq_ice, q, config...) end + + +function thermo_state_peq(param_set::APS, p::FT, e_int::FT, q_tot::FT) where {FT} + # config = (50, 1e-3, RootSolvers.RegulaFalsiMethod) + # config = (50, 1e-3, RootSolvers.NewtonMethodAD) + config = () + return TD.PhaseEquil_peq(param_set, p, e_int, q_tot, config...) +end + +function thermo_state_peq(param_set::APS, p::FT, e_int::FT, q_tot::FT, q_liq::FT, q_ice::FT) where {FT} + config = () + q = TD.PhasePartition(q_tot, q_liq, q_ice) + return TD.PhaseNonEquil_peq(param_set, p, e_int, q, config...) +end + +function geopotential(param_set, z::FT) where {FT} + grav = FT(ICP.grav(param_set)) + return grav * z +end + +function anelastic_total_enthalpy(param_set::APS, e_tot::FT, ts) where {FT} + Rm = TD.gas_constant_air(param_set, ts) + T = TD.air_temperature(param_set, ts) + return e_tot + Rm * T +end + +function anelastic_total_enthalpy(param_set::APS, e_tot::FT, p::FT, ρ::FT) where {FT} + return e_tot + p / ρ +end + +function kinetic_energy(u::FT, v::FT, w::FT) where {FT} + return FT(0.5) * (u^2 + v^2 + w^2) +end + +function total_energy(param_set::APS, z::FT, u::FT, v::FT, w::FT, p::FT, θ_liq_ice::FT, q_tot::FT) where {FT} + FT = eltype(grid) + config = () + e_kin = kinetic_energy(u::FT, v::FT, w::FT) + e_pot = geopotential(param_set, z::FT) + e_int = TD.internal_energy(TD.PhaseEquil_pθq(param_set, p, θ_liq_ice, q_tot, config...)) + return e_kin + e_pot + e_int +end + +function total_energy(param_set::APS, z::FT, u::FT, v::FT, w::FT, p::FT, θ_liq_ice::FT, q_tot::FT, q_liq::FT, q_ice::FT) where {FT} + FT = eltype(grid) + config = () + q = TD.PhasePartition(q_tot, q_liq, q_ice) + e_kin = kinetic_energy(u::FT, v::FT, w::FT) + e_pot = geopotential(param_set, z::FT) + e_int = TD.internal_energy(TD.PhaseEquil_pθq(param_set, p, θ_liq_ice, q, config...)) + return e_kin + e_pot + e_int +end diff --git a/src/types.jl b/src/types.jl index 33d60bb196..e6683f92f7 100644 --- a/src/types.jl +++ b/src/types.jl @@ -7,6 +7,7 @@ $(DocStringExtensions.FIELDS) """ Base.@kwdef struct PrecipFormation{FT} θ_liq_ice_tendency::FT + e_tot_tendency::FT qt_tendency::FT ql_tendency::FT qi_tendency::FT @@ -349,7 +350,7 @@ Base.@kwdef struct SurfaceBase{FT} ρq_tot_flux::FT = 0 ρq_liq_flux::FT = 0 ρq_ice_flux::FT = 0 - ρθ_liq_ice_flux::FT = 0 + ρe_tot_flux::FT = 0 ρu_flux::FT = 0 ρv_flux::FT = 0 obukhov_length::FT = 0 @@ -450,11 +451,7 @@ function CasesBase(case::T; inversion_type, surf_params, Fo, Rad, LESDat = nothi ) end -abstract type AbstractEnergyVariable end -struct ρeVar <: AbstractEnergyVariable end -struct ρθVar <: AbstractEnergyVariable end - -struct EDMFModel{N_up, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG, EV} +struct EDMFModel{N_up, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG} surface_area::FT max_area::FT minimum_area::FT @@ -469,7 +466,6 @@ struct EDMFModel{N_up, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG, EV} entr_dim_scale::EDS detr_dim_scale::DDS entr_pi_subset::EPG - evar::EV function EDMFModel(namelist, precip_model) where {PS} # TODO: move this into arg list FT = Float64 @@ -502,16 +498,6 @@ struct EDMFModel{N_up, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG, EV} error("Something went wrong. Invalid moisture model: '$moisture_model_name'") end - evar_name = parse_namelist(namelist, "energy_var"; default = "rhotheta") - - evar = if evar_name == "rhotheta" - ρθVar() - elseif evar_name == "rhoe" - ρeVar() - else - error("Bad energy variable given") - end - thermo_covariance_model_name = parse_namelist(namelist, "thermodynamics", "thermo_covariance_model"; default = "prognostic") @@ -669,8 +655,7 @@ struct EDMFModel{N_up, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG, EV} EBGC = typeof(bg_closure) ENT = typeof(en_thermo) EPG = typeof(entr_pi_subset) - EV = typeof(evar) - return new{n_updrafts, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG, EV}( + return new{n_updrafts, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG}( surface_area, max_area, minimum_area, @@ -685,7 +670,6 @@ struct EDMFModel{N_up, FT, MM, TCM, PM, PFM, ENT, EBGC, EC, EDS, DDS, EPG, EV} entr_dim_scale, detr_dim_scale, entr_pi_subset, - evar, ) end end @@ -694,7 +678,6 @@ n_updrafts(::EDMFModel{N_up}) where {N_up} = N_up Base.eltype(::EDMFModel{N_up, FT}) where {N_up, FT} = FT n_Π_groups(m::EDMFModel) = length(m.entr_pi_subset) entrainment_Π_subset(m::EDMFModel) = m.entr_pi_subset -energy_var(m::EDMFModel) = m.evar Base.broadcastable(edmf::EDMFModel) = Ref(edmf) diff --git a/src/update_aux.jl b/src/update_aux.jl index 59b0ea9834..1c832cb881 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -40,7 +40,6 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas ##### ##### center variables ##### - @inbounds for k in real_center_indices(grid) ##### ##### Set primitive variables @@ -48,7 +47,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas @inbounds for i in 1:N_up if is_surface_center(grid, k) if prog_up[i].ρarea[k] / ρ_c[k] >= edmf.minimum_area - θ_surf = θ_surface_bc(surf, grid, state, edmf, i) + θ_surf = θ_surface_bc(surf, grid, state, edmf, i, param_set) q_surf = q_surface_bc(surf, grid, state, edmf, i) a_surf = area_surface_bc(surf, edmf, i) aux_up[i].θ_liq_ice[k] = θ_surf @@ -60,8 +59,13 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas end else if prog_up[i].ρarea[k] / ρ_c[k] >= edmf.minimum_area - aux_up[i].θ_liq_ice[k] = prog_up[i].ρaθ_liq_ice[k] / prog_up[i].ρarea[k] + e_tot = prog_up[i].ρae_tot[k] / prog_up[i].ρarea[k] aux_up[i].q_tot[k] = prog_up[i].ρaq_tot[k] / prog_up[i].ρarea[k] + e_kin = kinetic_energy(prog_gm.u[k], prog_gm.v[k], prog_up[i].w[k]) + e_pot = geopotential(param_set, grid.zc[k].z) + e_int = e_tot - e_kin - e_pot + ts_up_i = TD.PhaseEquil_peq(p_c[k], e_int, aux_up[i].q_tot[k]) + aux_up[i].θ_liq_ice[k] = TD.liquid_ice_pottemp(ts_up_i) aux_up[i].area[k] = prog_up[i].ρarea[k] / ρ_c[k] else aux_up[i].θ_liq_ice[k] = aux_gm.θ_liq_ice[k] @@ -266,7 +270,6 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas 0 end end - ##### ##### face variables: diagnose primitive, diagnose env and compute bulk ##### @@ -504,7 +507,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas ### Diagnostic thermodynamiccovariances if edmf.thermo_covariance_model isa DiagnosticThermoCovariances - flux1 = surf.ρθ_liq_ice_flux + flux1 = surf.shf / TD.cp_m(param_set, ts_gm[kc_surf]) flux2 = surf.ρq_tot_flux zLL = grid.zc[kc_surf].z ustar = surf.ustar diff --git a/src/variables.jl b/src/variables.jl index ce0bbc5c15..7e0f562a81 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -38,6 +38,7 @@ cent_aux_vars_up(FT, edmf) = (; q_tot = FT(0), θ_liq_ice = FT(0), θ_liq_ice_tendency_precip_formation = FT(0), + e_tot_tendency_precip_formation = FT(0), qt_tendency_precip_formation = FT(0), cent_aux_vars_up_moisture(FT, edmf.moisture_model)..., entr_sc = FT(0), @@ -85,7 +86,7 @@ cent_aux_vars_edmf(FT, edmf) = (; q_ice = FT(0), T = FT(0), cloud_fraction = FT(0), - θ_liq_ice_tendency_precip_formation = FT(0), + e_tot_tendency_precip_formation = FT(0), qt_tendency_precip_formation = FT(0), cent_aux_vars_edmf_bulk_moisture(FT, edmf.moisture_model)..., ), @@ -109,7 +110,8 @@ cent_aux_vars_edmf(FT, edmf) = (; Hvar = FT(0), QTvar = FT(0), HQTcov = FT(0), - θ_liq_ice_tendency_precip_formation = FT(0), + #θ_liq_ice_tendency_precip_formation = FT(0), + e_tot_tendency_precip_formation = FT(0), qt_tendency_precip_formation = FT(0), cent_aux_vars_edmf_en_moisture(FT, edmf.moisture_model)..., unsat = (; q_tot = FT(0), θ_dry = FT(0), θ_virt = FT(0)), @@ -119,6 +121,7 @@ cent_aux_vars_edmf(FT, edmf) = (; HQTcov_rain_dt = FT(0), ), θ_liq_ice_tendency_precip_sinks = FT(0), + e_tot_tendency_precip_sinks = FT(0), qt_tendency_precip_sinks = FT(0), qr_tendency_evap = FT(0), qs_tendency_melt = FT(0), @@ -249,7 +252,7 @@ cent_prognostic_vars_up_moisture(::Type{FT}, ::EquilibriumMoisture) where {FT} = cent_prognostic_vars_up_moisture(::Type{FT}, ::NonEquilibriumMoisture) where {FT} = (; ρaq_liq = FT(0), ρaq_ice = FT(0)) cent_prognostic_vars_up(::Type{FT}, edmf) where {FT} = (; ρarea = FT(0), - ρaθ_liq_ice = FT(0), + ρae_tot = FT(0), ρaq_tot = FT(0), cent_prognostic_vars_up_noisy_relaxation(FT, edmf.entr_closure)..., cent_prognostic_vars_up_moisture(FT, edmf.moisture_model)...,