From 455326ff3d23e6550eb3b4ee33e65919e5e4f18d Mon Sep 17 00:00:00 2001 From: yairchn Date: Fri, 8 Jul 2022 15:25:08 -0700 Subject: [PATCH] change grid mean variable to theta --- README.md | 2 +- driver/Forcing.jl | 6 +- driver/Surface.jl | 6 +- driver/compute_diagnostics.jl | 2 + driver/dycore.jl | 88 ++++------ driver/dycore_variables.jl | 10 +- driver/main.jl | 2 +- integration_tests/sphere_utils.jl | 2 +- post_processing/mse_tables.jl | 260 +++++++++++++++--------------- src/EDMF_Environment.jl | 32 +--- src/EDMF_Precipitation.jl | 7 +- src/EDMF_Updrafts.jl | 6 +- src/EDMF_functions.jl | 52 +++--- src/TurbulenceConvection.jl | 2 +- src/diagnostics.jl | 2 +- src/microphysics_coupling.jl | 23 +-- src/thermodynamics.jl | 45 ------ src/types.jl | 3 +- src/update_aux.jl | 20 +-- src/variables.jl | 15 +- 20 files changed, 222 insertions(+), 363 deletions(-) diff --git a/README.md b/README.md index 1bd3bf661..57c2cf273 100644 --- a/README.md +++ b/README.md @@ -128,7 +128,7 @@ Upon running a particular experiment (described above), comparison plots (agains | `u` | [✓] | [x] | [x] | | `v` | [✓] | [x] | [x] | | `w` | [-] | [x] | [✓] | -| `e_tot` | [✓] | [x] | [✓] | +| `θ_liq_ice` | [✓] | [x] | [✓] | | `q_tot` | [✓] | [x] | [✓] | | `a` | [-] | [x] | [✓] | | `tke` | [x] | [✓] | [-] | diff --git a/driver/Forcing.jl b/driver/Forcing.jl index d4f6026ae..63274a95f 100644 --- a/driver/Forcing.jl +++ b/driver/Forcing.jl @@ -10,7 +10,7 @@ function initialize(::ForcingBase{ForcingLES}, grid, state, LESDat::LESData) getvar(var) = TC.pyinterp(vec(grid.zc.z), zc_les, TC.mean_nc_data(data, "profiles", var, imin, imax)) dTdt_hadv = getvar("dtdt_hadv") - T_nudge = getvar("temperature_mean") + H_nudge = getvar("thetali_mean") dTdt_fluc = getvar("dtdt_fluc") dqtdt_hadv = getvar("dqtdt_hadv") qt_nudge = getvar("qt_mean") @@ -18,11 +18,11 @@ function initialize(::ForcingBase{ForcingLES}, grid, state, LESDat::LESData) subsidence = getvar("ls_subsidence") u_nudge = getvar("u_mean") v_nudge = getvar("v_mean") - (; dTdt_hadv, T_nudge, dTdt_fluc, dqtdt_hadv, qt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge) + (; dTdt_hadv, H_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.T_nudge[k] = nt.T_nudge[k] + aux_gm.H_nudge[k] = nt.H_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 e12e64184..5bc41f4e0 100644 --- a/driver/Surface.jl +++ b/driver/Surface.jl @@ -66,7 +66,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, - ρe_tot_flux = shf + lhf, + ρθ_liq_ice_flux = shf / TD.cp_m(thermo_params, ts_in), ρq_tot_flux = lhf / TD.latent_heat_vapor(thermo_params, ts_in), wstar = convective_vel, ρq_liq_flux = FT(0), @@ -125,7 +125,7 @@ function get_surface( ustar = result.ustar, ρu_flux = result.ρτxz, ρv_flux = result.ρτyz, - ρe_tot_flux = shf + lhf, + ρθ_liq_ice_flux = shf / TD.cp_m(thermo_params, ts_in), ρq_tot_flux = lhf / TD.latent_heat_vapor(thermo_params, 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, - ρe_tot_flux = shf + lhf, + ρθ_liq_ice_flux = shf / TD.cp_m(thermo_params, ts_in), ρq_tot_flux = lhf / TD.latent_heat_vapor(thermo_params, ts_in), bflux = result.buoy_flux, wstar = convective_vel, diff --git a/driver/compute_diagnostics.jl b/driver/compute_diagnostics.jl index 0a9f0de7f..7aab4fb46 100644 --- a/driver/compute_diagnostics.jl +++ b/driver/compute_diagnostics.jl @@ -139,6 +139,8 @@ function compute_diagnostics!( surf = get_surface(surf_params, grid, state, t, param_set) @. aux_gm.s = TD.specific_entropy(thermo_params, ts_gm) + @. aux_gm.θ_liq_ice = prog_gm.ρθ_liq_ice / ρ_c + @. aux_gm.q_tot = prog_gm.ρq_tot / ρ_c @. aux_en.s = TD.specific_entropy(thermo_params, ts_en) @inbounds for k in TC.real_center_indices(grid) diff --git a/driver/dycore.jl b/driver/dycore.jl index 5fe75b467..74a8bdb5c 100644 --- a/driver/dycore.jl +++ b/driver/dycore.jl @@ -110,18 +110,15 @@ function compute_ref_state!( end -function set_thermo_state_peq!(state, grid, moisture_model, param_set) +function set_thermo_state_from_prog!(state, grid, moisture_model, param_set) Ic = CCO.InterpolateF2C() thermo_params = TCP.thermodynamics_params(param_set) 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_uₕ = TC.grid_mean_uₕ(state) p_c = aux_gm.p ρ_c = prog_gm.ρ - C123 = CCG.Covariant123Vector - @. aux_gm.e_kin = LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(prog_gm_f.w))) / 2 @inbounds for k in TC.real_center_indices(grid) thermo_args = if moisture_model isa TC.EquilibriumMoisture @@ -131,14 +128,18 @@ function set_thermo_state_peq!(state, grid, moisture_model, param_set) else error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium") end - e_pot = TC.geopotential(param_set, grid.zc.z[k]) - e_int = prog_gm.ρe_tot[k] / ρ_c[k] - aux_gm.e_kin[k] - e_pot - ts_gm[k] = TC.thermo_state_peq(param_set, p_c[k], e_int, prog_gm.ρq_tot[k] / ρ_c[k], thermo_args...) + ts_gm[k] = TC.thermo_state_pθq( + param_set, + p_c[k], + prog_gm.ρθ_liq_ice[k] / ρ_c[k], + prog_gm.ρq_tot[k] / ρ_c[k], + thermo_args..., + ) end return nothing end -function set_thermo_state_pθq!(state, grid, moisture_model, param_set) +function set_thermo_state_from_aux!(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) @@ -159,23 +160,13 @@ end function set_grid_mean_from_thermo_state!(param_set, state, grid) thermo_params = TCP.thermodynamics_params(param_set) - 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) - prog_gm_uₕ = TC.grid_mean_uₕ(state) + ts_gm = TC.center_aux_grid_mean(state).ts ρ_c = prog_gm.ρ - C123 = CCG.Covariant123Vector - @. prog_gm.ρe_tot = - ρ_c * TD.total_energy( - thermo_params, - ts_gm, - LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(prog_gm_f.w))) / 2, - TC.geopotential(param_set, grid.zc.z), - ) - @. prog_gm.ρq_tot = ρ_c * aux_gm.q_tot - + @. prog_gm.ρθ_liq_ice = ρ_c * TD.liquid_ice_pottemp(thermo_params, ts_gm) + @. prog_gm.ρq_tot = ρ_c * TD.total_specific_humidity(thermo_params, ts_gm) return nothing end @@ -188,12 +179,11 @@ function assign_thermo_aux!(state, grid, moisture_model, param_set) @inbounds for k in TC.real_center_indices(grid) ts = ts_gm[k] aux_gm.q_tot[k] = prog_gm.ρq_tot[k] / ρ_c[k] + aux_gm.θ_liq_ice[k] = prog_gm.ρθ_liq_ice[k] / ρ_c[k] aux_gm.q_liq[k] = TD.liquid_specific_humidity(thermo_params, ts) aux_gm.q_ice[k] = TD.ice_specific_humidity(thermo_params, ts) aux_gm.T[k] = TD.air_temperature(thermo_params, ts) aux_gm.RH[k] = TD.relative_humidity(thermo_params, ts) - aux_gm.θ_liq_ice[k] = TD.liquid_ice_pottemp(thermo_params, ts) - aux_gm.h_tot[k] = TC.total_enthalpy(param_set, prog_gm.ρe_tot[k] / ρ_c[k], ts) end return end @@ -233,7 +223,7 @@ function ∑tendencies!(tendencies::FV, prog::FV, params::NT, t::Real) where {NT state = TC.column_state(prog, aux, tendencies, inds...) grid = TC.Grid(state) - set_thermo_state_peq!(state, grid, edmf.moisture_model, param_set) + set_thermo_state_from_prog!(state, grid, edmf.moisture_model, param_set) assign_thermo_aux!(state, grid, edmf.moisture_model, param_set) aux_gm = TC.center_aux_grid_mean(state) @@ -288,9 +278,6 @@ function compute_gm_tendencies!( ) thermo_params = TCP.thermodynamics_params(param_set) Ic = CCO.InterpolateF2C() - R_d = TCP.R_d(param_set) - T_0 = TCP.T_0(param_set) - Lv_0 = TCP.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) @@ -299,7 +286,7 @@ function compute_gm_tendencies!( 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) - ∇MSE_gm = TC.center_aux_grid_mean(state).∇MSE_gm + ∇θ_liq_ice_gm = TC.center_aux_grid_mean(state).∇θ_liq_ice_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) @@ -311,13 +298,13 @@ function compute_gm_tendencies!( aux_tc = TC.center_aux_turbconv(state) ts_gm = TC.center_aux_grid_mean(state).ts - MSE_gm_toa = aux_gm.h_tot[kc_toa] - aux_gm.e_kin[kc_toa] + θ_liq_ice_gm_toa = prog_gm.ρθ_liq_ice[kc_toa] / ρ_c[kc_toa] q_tot_gm_toa = prog_gm.ρq_tot[kc_toa] / ρ_c[kc_toa] - RBe = CCO.RightBiasedC2F(; top = CCO.SetValue(MSE_gm_toa)) + RBθ = CCO.RightBiasedC2F(; top = CCO.SetValue(θ_liq_ice_gm_toa)) RBq = CCO.RightBiasedC2F(; top = CCO.SetValue(q_tot_gm_toa)) wvec = CC.Geometry.WVector ∇c = CCO.DivergenceF2C() - @. ∇MSE_gm = ∇c(wvec(RBe(aux_gm.h_tot - aux_gm.e_kin))) + @. ∇θ_liq_ice_gm = ∇c(wvec(RBθ(prog_gm.ρθ_liq_ice / ρ_c))) @. ∇q_tot_gm = ∇c(wvec(RBq(prog_gm.ρq_tot / ρ_c))) if edmf.moisture_model isa TC.NonEquilibriumMoisture @@ -334,9 +321,7 @@ function compute_gm_tendencies!( # Apply forcing and radiation prog_gm_uₕ = TC.grid_mean_uₕ(state) aux_gm_uₕ_g = TC.grid_mean_uₕ_g(state) - # prog_gm_v = TC.grid_mean_v(state) tendencies_gm_uₕ = TC.tendencies_grid_mean_uₕ(state) - # tendencies_gm_v = TC.tendencies_grid_mean_v(state) prog_gm_u = TC.physical_grid_mean_u(state) prog_gm_v = TC.physical_grid_mean_v(state) @@ -355,15 +340,9 @@ function compute_gm_tendencies!( @inbounds for k in TC.real_center_indices(grid) - cp_m = TD.cp_m(thermo_params, ts_gm[k]) - cp_v = TCP.cp_v(param_set) - cv_v = TCP.cv_v(param_set) - R_v = TCP.R_v(param_set) - cv_m = TD.cv_m(thermo_params, ts_gm[k]) - h_v = cv_v * (aux_gm.T[k] - T_0) + Lv_0 - R_v * T_0 - + Π = TD.exner(thermo_params, ts_gm[k]) # LS Subsidence - tendencies_gm.ρe_tot[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇MSE_gm[k] + tendencies_gm.ρθ_liq_ice[k] -= ρ_c[k] * aux_gm.subsidence[k] * ∇θ_liq_ice_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] @@ -371,12 +350,12 @@ function compute_gm_tendencies!( end # Radiation if Cases.rad_type(radiation) <: Union{Cases.RadiationDYCOMS_RF01, Cases.RadiationLES, Cases.RadiationTRMM_LBA} - tendencies_gm.ρe_tot[k] += ρ_c[k] * cv_m * aux_gm.dTdt_rad[k] + tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_rad[k] / Π end # LS advection tendencies_gm.ρq_tot[k] += ρ_c[k] * aux_gm.dqtdt_hadv[k] if !(Cases.force_type(force) <: Cases.ForcingDYCOMS_RF01) - tendencies_gm.ρe_tot[k] += ρ_c[k] * (cp_m * aux_gm.dTdt_hadv[k] + h_v * aux_gm.dqtdt_hadv[k]) + tendencies_gm.ρθ_liq_ice[k] += ρ_c[k] * aux_gm.dTdt_hadv[k] / Π end if edmf.moisture_model isa TC.NonEquilibriumMoisture tendencies_gm.q_liq[k] += aux_gm.dqldt[k] @@ -385,13 +364,13 @@ function compute_gm_tendencies!( # LES specific forcings if Cases.force_type(force) <: Cases.ForcingLES - T_fluc = aux_gm.dTdt_fluc[k] + H_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_τᵣ Γᵣ = compute_les_Γᵣ(grid.zc[k].z, force.scalar_nudge_τᵣ, force.scalar_nudge_zᵢ, force.scalar_nudge_zᵣ) - gm_T_nudge_k = Γᵣ * (aux_gm.T_nudge[k] - aux_gm.T[k]) + gm_H_nudge_k = Γᵣ * (aux_gm.H_nudge[k] - aux_gm.θ_liq_ice[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]) @@ -399,8 +378,7 @@ function compute_gm_tendencies!( end 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.ρθ_liq_ice[k] += ρ_c[k] * (gm_H_nudge_k + H_fluc) tendencies_gm_uₕ[k] += CCG.Covariant12Vector(CCG.UVVector(gm_U_nudge_k, gm_V_nudge_k), lg[k]) if edmf.moisture_model isa TC.NonEquilibriumMoisture tendencies_gm.q_liq[k] += aux_gm.dqldt_hadv[k] + gm_q_liq_nudge_k + aux_gm.dqldt_fluc[k] @@ -416,11 +394,11 @@ function compute_gm_tendencies!( aux_tc.qt_tendency_precip_sinks[k] ) - tendencies_gm.ρe_tot[k] += + tendencies_gm.ρθ_liq_ice[k] += ρ_c[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] + aux_bulk.θ_liq_ice_tendency_precip_formation[k] + + aux_en.θ_liq_ice_tendency_precip_formation[k] + + aux_tc.θ_liq_ice_tendency_precip_sinks[k] ) if edmf.moisture_model isa TC.NonEquilibriumMoisture @@ -428,16 +406,16 @@ function compute_gm_tendencies!( 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, param_set) - sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot + TC.compute_sgs_flux!(edmf, grid, state, surf) + sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot sgs_flux_uₕ = aux_gm_f.sgs_flux_uₕ - tends_ρe_tot = tendencies_gm.ρe_tot + tends_ρθ_liq_ice = tendencies_gm.ρθ_liq_ice tends_ρq_tot = tendencies_gm.ρq_tot tends_uₕ = TC.tendencies_grid_mean_uₕ(state) ∇sgs = CCO.DivergenceF2C() - @. tends_ρe_tot += -∇sgs(wvec(sgs_flux_h_tot)) + @. tends_ρθ_liq_ice += -∇sgs(wvec(sgs_flux_θ_liq_ice)) @. tends_ρq_tot += -∇sgs(wvec(sgs_flux_q_tot)) @. tends_uₕ += -∇sgs(sgs_flux_uₕ) / ρ_c diff --git a/driver/dycore_variables.jl b/driver/dycore_variables.jl index 300badf5c..aa4836344 100644 --- a/driver/dycore_variables.jl +++ b/driver/dycore_variables.jl @@ -48,14 +48,14 @@ cent_aux_vars_gm(FT, local_geometry, edmf) = (; subsidence = FT(0), #Large-scale subsidence dTdt_hadv = FT(0), #Horizontal advection of temperature dqtdt_hadv = FT(0), #Horizontal advection of moisture - T_nudge = FT(0), #Reference T profile for relaxation tendency + H_nudge = FT(0), #Reference thetali 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 u_nudge = FT(0), #Reference u profile for relaxation tendency v_nudge = FT(0), #Reference v profile for relaxation tendency uₕ_g = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry), #Geostrophic u velocity - ∇MSE_gm = FT(0), + ∇θ_liq_ice_gm = FT(0), ∇q_tot_gm = FT(0), cent_aux_vars_gm_moisture(FT, edmf.moisture_model)..., θ_virt = FT(0), @@ -63,8 +63,6 @@ cent_aux_vars_gm(FT, local_geometry, edmf) = (; θ_liq_ice = FT(0), q_tot = FT(0), p = FT(0), - e_kin = FT(0), - h_tot = FT(0), ) cent_aux_vars(FT, local_geometry, edmf) = (; cent_aux_vars_gm(FT, local_geometry, edmf)..., TC.cent_aux_vars_edmf(FT, local_geometry, edmf)...) @@ -77,7 +75,7 @@ face_aux_vars_gm(FT, local_geometry, edmf) = (; diffusive_flux_s = FT(0), total_flux_s = FT(0), f_rad = FT(0), - sgs_flux_h_tot = FT(0), + sgs_flux_θ_liq_ice = FT(0), sgs_flux_q_tot = FT(0), face_aux_vars_gm_moisture(FT, edmf.moisture_model)..., sgs_flux_uₕ = CCG.Covariant3Vector(FT(0)) ⊗ CCG.Covariant12Vector(FT(0), FT(0)), @@ -128,7 +126,7 @@ cent_prognostic_vars_gm_moisture(::Type{FT}, ::TC.EquilibriumMoisture) where {FT cent_prognostic_vars_gm(::Type{FT}, local_geometry, edmf) where {FT} = (; ρ = FT(0), uₕ = CCG.Covariant12Vector(CCG.UVVector(FT(0), FT(0)), local_geometry), - ρe_tot = FT(0), + ρθ_liq_ice = FT(0), ρq_tot = FT(0), cent_prognostic_vars_gm_moisture(FT, edmf.moisture_model)..., ) diff --git a/driver/main.jl b/driver/main.jl index 700ea147a..5bcd0ece8 100644 --- a/driver/main.jl +++ b/driver/main.jl @@ -266,7 +266,7 @@ function initialize(sim::Simulation1d) end end Cases.initialize_profiles(case, grid, param_set, state; les_data_kwarg...) - set_thermo_state_pθq!(state, grid, edmf.moisture_model, param_set) + set_thermo_state_from_aux!(state, grid, edmf.moisture_model, param_set) set_grid_mean_from_thermo_state!(param_set, state, grid) assign_thermo_aux!(state, grid, edmf.moisture_model, param_set) Cases.initialize_forcing(case, forcing, grid, state, param_set; les_data_kwarg...) diff --git a/integration_tests/sphere_utils.jl b/integration_tests/sphere_utils.jl index c95f602ec..0172e601f 100644 --- a/integration_tests/sphere_utils.jl +++ b/integration_tests/sphere_utils.jl @@ -79,7 +79,7 @@ function test_zero_horizontal_variance(Y; filter_prop_chain = pc -> true) (:cent, :ρ) => false, (:cent, :uₕ, :components, :data, 1) => false, (:cent, :uₕ, :components, :data, 2) => false, - (:cent, :ρe_tot) => false, + (:cent, :ρθ_liq_ice) => false, (:cent, :ρq_tot) => false, (:cent, :turbconv, :en, :ρatke) => false, (:cent, :turbconv, :up, 1, :ρarea) => false, diff --git a/post_processing/mse_tables.jl b/post_processing/mse_tables.jl index 5cc311b77..8bbc652d3 100644 --- a/post_processing/mse_tables.jl +++ b/post_processing/mse_tables.jl @@ -5,165 +5,165 @@ all_best_mse = OrderedCollections.OrderedDict() # all_best_mse["ARM_SGP"] = OrderedCollections.OrderedDict() -all_best_mse["ARM_SGP"]["qt_mean"] = 0.0957229254528298 -all_best_mse["ARM_SGP"]["updraft_area"] = 319.4000272588446 -all_best_mse["ARM_SGP"]["updraft_w"] = 97.87630299882935 -all_best_mse["ARM_SGP"]["updraft_qt"] = 24.620317112482173 -all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.04137385462744 +all_best_mse["ARM_SGP"]["qt_mean"] = 0.2547158947572649 +all_best_mse["ARM_SGP"]["updraft_area"] = 327.4648381189871 +all_best_mse["ARM_SGP"]["updraft_w"] = 154.33612466069923 +all_best_mse["ARM_SGP"]["updraft_qt"] = 30.03992043607562 +all_best_mse["ARM_SGP"]["updraft_thetal"] = 172.0117106042987 all_best_mse["ARM_SGP"]["u_mean"] = 1.3235797273549681e-5 -all_best_mse["ARM_SGP"]["tke_mean"] = 545.4969096763037 -all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00021245440433044206 -all_best_mse["ARM_SGP"]["ql_mean"] = 218.56477273777267 +all_best_mse["ARM_SGP"]["tke_mean"] = 1107.1177081770727 +all_best_mse["ARM_SGP"]["temperature_mean"] = 0.00012731778118573874 +all_best_mse["ARM_SGP"]["ql_mean"] = 212.242418342972 all_best_mse["ARM_SGP"]["qi_mean"] = "NA" -all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00017702514712616304 -all_best_mse["ARM_SGP"]["Hvar_mean"] = 3269.1246326339074 -all_best_mse["ARM_SGP"]["QTvar_mean"] = 1844.1903620926616 +all_best_mse["ARM_SGP"]["thetal_mean"] = 0.00011913333708523678 +all_best_mse["ARM_SGP"]["Hvar_mean"] = 1137.8827730362934 +all_best_mse["ARM_SGP"]["QTvar_mean"] = 745.1763763345691 # all_best_mse["Bomex"] = OrderedCollections.OrderedDict() -all_best_mse["Bomex"]["qt_mean"] = 0.14717448881738263 -all_best_mse["Bomex"]["updraft_area"] = 137.47850762307968 -all_best_mse["Bomex"]["updraft_w"] = 17.360636564028198 -all_best_mse["Bomex"]["updraft_qt"] = 6.561881005993562 -all_best_mse["Bomex"]["updraft_thetal"] = 69.76738309609604 -all_best_mse["Bomex"]["v_mean"] = 49.2245869661139 -all_best_mse["Bomex"]["u_mean"] = 0.10356807041387374 -all_best_mse["Bomex"]["tke_mean"] = 236.79221761729715 -all_best_mse["Bomex"]["temperature_mean"] = 5.78644277240463e-5 -all_best_mse["Bomex"]["ql_mean"] = 22.16679293843376 +all_best_mse["Bomex"]["qt_mean"] = 0.10129372376733127 +all_best_mse["Bomex"]["updraft_area"] = 127.66714249325291 +all_best_mse["Bomex"]["updraft_w"] = 18.206140756564217 +all_best_mse["Bomex"]["updraft_qt"] = 7.078682148317164 +all_best_mse["Bomex"]["updraft_thetal"] = 69.79384019189457 +all_best_mse["Bomex"]["v_mean"] = 62.47158884512445 +all_best_mse["Bomex"]["u_mean"] = 0.255429889697546 +all_best_mse["Bomex"]["tke_mean"] = 72.56761522891617 +all_best_mse["Bomex"]["temperature_mean"] = 4.1194985385604214e-5 +all_best_mse["Bomex"]["ql_mean"] = 9.479895025499541 all_best_mse["Bomex"]["qi_mean"] = "NA" -all_best_mse["Bomex"]["thetal_mean"] = 6.01292633415193e-5 -all_best_mse["Bomex"]["Hvar_mean"] = 9655.675978036958 -all_best_mse["Bomex"]["QTvar_mean"] = 3170.6892285167132 +all_best_mse["Bomex"]["thetal_mean"] = 4.1878936713814785e-5 +all_best_mse["Bomex"]["Hvar_mean"] = 2386.664143185325 +all_best_mse["Bomex"]["QTvar_mean"] = 868.8201442465643 # all_best_mse["DryBubble"] = OrderedCollections.OrderedDict() -all_best_mse["DryBubble"]["updraft_area"] = 1.2420615202367218e-5 -all_best_mse["DryBubble"]["updraft_w"] = 0.00034476053312610907 -all_best_mse["DryBubble"]["updraft_thetal"] = 1.3657383140501062e-12 -all_best_mse["DryBubble"]["u_mean"] = 1.9106710242038562e-5 -all_best_mse["DryBubble"]["tke_mean"] = 0.0021417835070300507 -all_best_mse["DryBubble"]["temperature_mean"] = 1.7377049913414865e-11 -all_best_mse["DryBubble"]["thetal_mean"] = 1.5224505591146085e-11 -all_best_mse["DryBubble"]["Hvar_mean"] = 0.009862311144017531 +all_best_mse["DryBubble"]["updraft_area"] = 0.493762139856581 +all_best_mse["DryBubble"]["updraft_w"] = 6.011782176716596 +all_best_mse["DryBubble"]["updraft_thetal"] = 8.017439134017947e-8 +all_best_mse["DryBubble"]["u_mean"] = 0.003349490718428698 +all_best_mse["DryBubble"]["tke_mean"] = 111.05775466853102 +all_best_mse["DryBubble"]["temperature_mean"] = 7.786467452320724e-7 +all_best_mse["DryBubble"]["thetal_mean"] = 6.504125487505511e-7 +all_best_mse["DryBubble"]["Hvar_mean"] = 48.36457777988673 # all_best_mse["DYCOMS_RF01"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.03545748659531928 -all_best_mse["DYCOMS_RF01"]["ql_mean"] = 0.27084019611335747 -all_best_mse["DYCOMS_RF01"]["updraft_area"] = 28.923798054739034 -all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.498010535840437 -all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.3416602463392286 -all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.21582271825862 -all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.002122282580053137 -all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.08521781532913578 -all_best_mse["DYCOMS_RF01"]["tke_mean"] = 20.89501104721309 -all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 2.4430721263798322e-5 -all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 2.402869950855993e-5 -all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1289.247167722031 -all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 513.3492355706912 +all_best_mse["DYCOMS_RF01"]["qt_mean"] = 0.018045420296732124 +all_best_mse["DYCOMS_RF01"]["ql_mean"] = 3.485041337094255 +all_best_mse["DYCOMS_RF01"]["updraft_area"] = 30.03055645110736 +all_best_mse["DYCOMS_RF01"]["updraft_w"] = 4.3200951277825395 +all_best_mse["DYCOMS_RF01"]["updraft_qt"] = 2.130027153656729 +all_best_mse["DYCOMS_RF01"]["updraft_thetal"] = 46.18824334991148 +all_best_mse["DYCOMS_RF01"]["v_mean"] = 0.002887414638915658 +all_best_mse["DYCOMS_RF01"]["u_mean"] = 0.08828127568713307 +all_best_mse["DYCOMS_RF01"]["tke_mean"] = 21.324141485523885 +all_best_mse["DYCOMS_RF01"]["temperature_mean"] = 7.18490422610244e-5 +all_best_mse["DYCOMS_RF01"]["thetal_mean"] = 7.23503697573604e-5 +all_best_mse["DYCOMS_RF01"]["Hvar_mean"] = 1287.456160459314 +all_best_mse["DYCOMS_RF01"]["QTvar_mean"] = 513.1223107144899 # all_best_mse["DYCOMS_RF02"] = OrderedCollections.OrderedDict() -all_best_mse["DYCOMS_RF02"]["qt_mean"] = 0.10731034933376828 -all_best_mse["DYCOMS_RF02"]["ql_mean"] = 21.014747563221093 -all_best_mse["DYCOMS_RF02"]["qr_mean"] = 12.252328654138124 -all_best_mse["DYCOMS_RF02"]["updraft_area"] = 25.858245790929555 -all_best_mse["DYCOMS_RF02"]["updraft_w"] = 3.9691635533260716 -all_best_mse["DYCOMS_RF02"]["updraft_qt"] = 5.151402482152906 -all_best_mse["DYCOMS_RF02"]["updraft_thetal"] = 40.581660628288304 -all_best_mse["DYCOMS_RF02"]["v_mean"] = 42.847475583149496 -all_best_mse["DYCOMS_RF02"]["u_mean"] = 20.198387737817814 -all_best_mse["DYCOMS_RF02"]["tke_mean"] = 27.654214928029784 -all_best_mse["DYCOMS_RF02"]["temperature_mean"] = 0.0001235204286817038 -all_best_mse["DYCOMS_RF02"]["thetal_mean"] = 0.0001466546426767952 -all_best_mse["DYCOMS_RF02"]["Hvar_mean"] = 1488.3578301861028 -all_best_mse["DYCOMS_RF02"]["QTvar_mean"] = 853.5369476921379 +all_best_mse["DYCOMS_RF02"]["qt_mean"] = 0.07954271879125009 +all_best_mse["DYCOMS_RF02"]["ql_mean"] = 8.799228881307583 +all_best_mse["DYCOMS_RF02"]["qr_mean"] = 6.430686369966681 +all_best_mse["DYCOMS_RF02"]["updraft_area"] = 28.085934862282514 +all_best_mse["DYCOMS_RF02"]["updraft_w"] = 8.389588468411938 +all_best_mse["DYCOMS_RF02"]["updraft_qt"] = 5.121426021575978 +all_best_mse["DYCOMS_RF02"]["updraft_thetal"] = 40.551767054185625 +all_best_mse["DYCOMS_RF02"]["v_mean"] = 42.99394133720737 +all_best_mse["DYCOMS_RF02"]["u_mean"] = 20.04098869694078 +all_best_mse["DYCOMS_RF02"]["tke_mean"] = 22.349854113393853 +all_best_mse["DYCOMS_RF02"]["temperature_mean"] = 6.266974778894454e-5 +all_best_mse["DYCOMS_RF02"]["thetal_mean"] = 6.600697249731386e-5 +all_best_mse["DYCOMS_RF02"]["Hvar_mean"] = 1464.6190584343099 +all_best_mse["DYCOMS_RF02"]["QTvar_mean"] = 514.8544151811988 # all_best_mse["GABLS"] = OrderedCollections.OrderedDict() -all_best_mse["GABLS"]["updraft_thetal"] = 0.0 -all_best_mse["GABLS"]["v_mean"] = 0.0 -all_best_mse["GABLS"]["u_mean"] = 0.0 -all_best_mse["GABLS"]["tke_mean"] = 0.0 -all_best_mse["GABLS"]["temperature_mean"] = 0.0 -all_best_mse["GABLS"]["thetal_mean"] = 0.0 -all_best_mse["GABLS"]["Hvar_mean"] = 0.0 +all_best_mse["GABLS"]["updraft_thetal"] = 7.370200362804219e-6 +all_best_mse["GABLS"]["v_mean"] = 0.16957612628684732 +all_best_mse["GABLS"]["u_mean"] = 0.035304094008498016 +all_best_mse["GABLS"]["tke_mean"] = 1.0025596150147267 +all_best_mse["GABLS"]["temperature_mean"] = 7.392188550948941e-6 +all_best_mse["GABLS"]["thetal_mean"] = 7.370200362804219e-6 +all_best_mse["GABLS"]["Hvar_mean"] = 5.788951267721536 # all_best_mse["life_cycle_Tan2018"] = OrderedCollections.OrderedDict() -all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.0006793419796923063 -all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 0.4223103576324509 -all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 0.28501770251347985 -all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 0.16893861214312877 -all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.00434583560575179 -all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 1.4901114839547677e-6 -all_best_mse["life_cycle_Tan2018"]["v_mean"] = 1.9344142199629524 -all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.0024047438848902075 -all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 0.799788533111904 -all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 2.4499461764092775e-7 -all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 2.319750620411135e-7 -all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 26561.574872179564 -all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 14253.619456402466 +all_best_mse["life_cycle_Tan2018"]["qt_mean"] = 0.00917106749132344 +all_best_mse["life_cycle_Tan2018"]["ql_mean"] = 6.493950775313687 +all_best_mse["life_cycle_Tan2018"]["updraft_area"] = 9.18749428417205 +all_best_mse["life_cycle_Tan2018"]["updraft_w"] = 3.6572697607902693 +all_best_mse["life_cycle_Tan2018"]["updraft_qt"] = 0.010843993727610806 +all_best_mse["life_cycle_Tan2018"]["updraft_thetal"] = 4.5654280418170335e-6 +all_best_mse["life_cycle_Tan2018"]["v_mean"] = 31.09321053446677 +all_best_mse["life_cycle_Tan2018"]["u_mean"] = 0.0688779240697305 +all_best_mse["life_cycle_Tan2018"]["tke_mean"] = 32.43511480181172 +all_best_mse["life_cycle_Tan2018"]["temperature_mean"] = 6.142346580014008e-6 +all_best_mse["life_cycle_Tan2018"]["thetal_mean"] = 6.003896494625166e-6 +all_best_mse["life_cycle_Tan2018"]["Hvar_mean"] = 8489.761225775786 +all_best_mse["life_cycle_Tan2018"]["QTvar_mean"] = 4885.725673654171 # all_best_mse["Nieuwstadt"] = OrderedCollections.OrderedDict() -all_best_mse["Nieuwstadt"]["updraft_area"] = 111.0115511709332 -all_best_mse["Nieuwstadt"]["updraft_w"] = 12.02108001495304 -all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.6112481780205 -all_best_mse["Nieuwstadt"]["u_mean"] = 13.700796002736682 -all_best_mse["Nieuwstadt"]["tke_mean"] = 108.98689584274143 -all_best_mse["Nieuwstadt"]["temperature_mean"] = 5.77324171994201e-5 -all_best_mse["Nieuwstadt"]["thetal_mean"] = 5.4307129370935114e-5 -all_best_mse["Nieuwstadt"]["Hvar_mean"] = 664.9274745155124 +all_best_mse["Nieuwstadt"]["updraft_area"] = 98.70082646564556 +all_best_mse["Nieuwstadt"]["updraft_w"] = 14.184262620204711 +all_best_mse["Nieuwstadt"]["updraft_thetal"] = 117.60593470655523 +all_best_mse["Nieuwstadt"]["u_mean"] = 13.553417467513016 +all_best_mse["Nieuwstadt"]["tke_mean"] = 283.6622287396672 +all_best_mse["Nieuwstadt"]["temperature_mean"] = 1.1367531237369026e-5 +all_best_mse["Nieuwstadt"]["thetal_mean"] = 1.1155541949121985e-5 +all_best_mse["Nieuwstadt"]["Hvar_mean"] = 717.3487709472662 # all_best_mse["Rico"] = OrderedCollections.OrderedDict() -all_best_mse["Rico"]["qt_mean"] = 3.3793232431058327 -all_best_mse["Rico"]["updraft_area"] = 491.75586735299726 -all_best_mse["Rico"]["updraft_w"] = 135.85242239121547 -all_best_mse["Rico"]["updraft_qt"] = 9.550461904407854 -all_best_mse["Rico"]["updraft_thetal"] = 133.29492887700837 -all_best_mse["Rico"]["v_mean"] = 0.3531392261411799 -all_best_mse["Rico"]["u_mean"] = 0.29262536503289555 -all_best_mse["Rico"]["tke_mean"] = 541.949305663487 -all_best_mse["Rico"]["temperature_mean"] = 0.0017181411832877024 -all_best_mse["Rico"]["ql_mean"] = 175.94308418250262 +all_best_mse["Rico"]["qt_mean"] = 0.9613962072647245 +all_best_mse["Rico"]["updraft_area"] = 477.5831767593734 +all_best_mse["Rico"]["updraft_w"] = 77.66215281843687 +all_best_mse["Rico"]["updraft_qt"] = 16.85052739385283 +all_best_mse["Rico"]["updraft_thetal"] = 133.8740333117735 +all_best_mse["Rico"]["v_mean"] = 0.5551880545138421 +all_best_mse["Rico"]["u_mean"] = 0.3553124250721766 +all_best_mse["Rico"]["tke_mean"] = 142.97586711076562 +all_best_mse["Rico"]["temperature_mean"] = 0.0005430301553147434 +all_best_mse["Rico"]["ql_mean"] = 130.0849649957846 all_best_mse["Rico"]["qi_mean"] = "NA" -all_best_mse["Rico"]["qr_mean"] = 690.496245724426 -all_best_mse["Rico"]["thetal_mean"] = 0.0016636133937053483 -all_best_mse["Rico"]["Hvar_mean"] = 86813.00588246963 -all_best_mse["Rico"]["QTvar_mean"] = 14277.081614395107 +all_best_mse["Rico"]["qr_mean"] = 678.7643434304638 +all_best_mse["Rico"]["thetal_mean"] = 0.0005188908711384763 +all_best_mse["Rico"]["Hvar_mean"] = 64412.68057258725 +all_best_mse["Rico"]["QTvar_mean"] = 14582.516453015409 # all_best_mse["Soares"] = OrderedCollections.OrderedDict() -all_best_mse["Soares"]["qt_mean"] = 0.016623077511468968 -all_best_mse["Soares"]["updraft_area"] = 120.84628443038324 -all_best_mse["Soares"]["updraft_w"] = 8.848702232570785 -all_best_mse["Soares"]["updraft_qt"] = 22.975618329088768 -all_best_mse["Soares"]["updraft_thetal"] = 65.72527174155819 -all_best_mse["Soares"]["u_mean"] = 93.82461696990283 -all_best_mse["Soares"]["tke_mean"] = 78.85574467887265 -all_best_mse["Soares"]["temperature_mean"] = 6.502830780241758e-5 -all_best_mse["Soares"]["thetal_mean"] = 5.645806724553706e-5 -all_best_mse["Soares"]["Hvar_mean"] = 570.5288921821567 +all_best_mse["Soares"]["qt_mean"] = 0.14271215380491728 +all_best_mse["Soares"]["updraft_area"] = 94.42322359287252 +all_best_mse["Soares"]["updraft_w"] = 13.054667599240542 +all_best_mse["Soares"]["updraft_qt"] = 23.635846787398997 +all_best_mse["Soares"]["updraft_thetal"] = 65.72137426540347 +all_best_mse["Soares"]["u_mean"] = 93.89561563090922 +all_best_mse["Soares"]["tke_mean"] = 216.5798143634326 +all_best_mse["Soares"]["temperature_mean"] = 1.3143615855019983e-5 +all_best_mse["Soares"]["thetal_mean"] = 1.206813295763359e-5 +all_best_mse["Soares"]["Hvar_mean"] = 679.5403978722024 # all_best_mse["TRMM_LBA"] = OrderedCollections.OrderedDict() -all_best_mse["TRMM_LBA"]["qt_mean"] = 3.1675432498041762 -all_best_mse["TRMM_LBA"]["updraft_area"] = 1565.8738087698493 -all_best_mse["TRMM_LBA"]["updraft_w"] = 9326.019529872992 -all_best_mse["TRMM_LBA"]["updraft_qt"] = 227.5109580786232 -all_best_mse["TRMM_LBA"]["updraft_thetal"] = 542.7674494382694 -all_best_mse["TRMM_LBA"]["v_mean"] = 68.86600294731369 -all_best_mse["TRMM_LBA"]["u_mean"] = 29.73490960969987 -all_best_mse["TRMM_LBA"]["tke_mean"] = 81802.70872855624 -all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0006388401486651304 -all_best_mse["TRMM_LBA"]["ql_mean"] = 214176.04627885437 +all_best_mse["TRMM_LBA"]["qt_mean"] = 2.7027676543718586 +all_best_mse["TRMM_LBA"]["updraft_area"] = 1566.030013951543 +all_best_mse["TRMM_LBA"]["updraft_w"] = 9050.830076664302 +all_best_mse["TRMM_LBA"]["updraft_qt"] = 226.65788400374382 +all_best_mse["TRMM_LBA"]["updraft_thetal"] = 542.1440083827489 +all_best_mse["TRMM_LBA"]["v_mean"] = 69.36685373865427 +all_best_mse["TRMM_LBA"]["u_mean"] = 30.0014676696119 +all_best_mse["TRMM_LBA"]["tke_mean"] = 77186.86238635755 +all_best_mse["TRMM_LBA"]["temperature_mean"] = 0.0003917448495922493 +all_best_mse["TRMM_LBA"]["ql_mean"] = 238273.87485678124 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.0004834764358320113 -all_best_mse["TRMM_LBA"]["Hvar_mean"] = 411711.728982357 -all_best_mse["TRMM_LBA"]["QTvar_mean"] = 6219.713466364576 +all_best_mse["TRMM_LBA"]["thetal_mean"] = 0.0003644659154349736 +all_best_mse["TRMM_LBA"]["Hvar_mean"] = 459671.2135908613 +all_best_mse["TRMM_LBA"]["QTvar_mean"] = 6372.959094492408 # all_best_mse["LES_driven_SCM"] = OrderedCollections.OrderedDict() -all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.606074644315858 -all_best_mse["LES_driven_SCM"]["v_mean"] = 1.3477991190958085 -all_best_mse["LES_driven_SCM"]["u_mean"] = 0.4688043856274925 -all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0014871026793565234 -all_best_mse["LES_driven_SCM"]["ql_mean"] = 56341.54777517539 -all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.001757073568632076 +all_best_mse["LES_driven_SCM"]["qt_mean"] = 3.578203745164621 +all_best_mse["LES_driven_SCM"]["v_mean"] = 1.3789743057943122 +all_best_mse["LES_driven_SCM"]["u_mean"] = 0.47628634045978246 +all_best_mse["LES_driven_SCM"]["temperature_mean"] = 0.0012878801731970955 +all_best_mse["LES_driven_SCM"]["ql_mean"] = 54805.55242360169 +all_best_mse["LES_driven_SCM"]["thetal_mean"] = 0.0015208242323557642 # ################################# ################################# diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index fb10b1f06..43ad9f636 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -32,7 +32,6 @@ function microphysics( prog_pr.q_sno[k], aux_en.area[k], ρ_c[k], - FT(grid.zc[k].z), Δt, ts, precip_fraction, @@ -56,7 +55,7 @@ function microphysics( # update_env_precip_tendencies # TODO: move ..._tendency_precip_formation to diagnostics aux_en.qt_tendency_precip_formation[k] = mph.qt_tendency * aux_en.area[k] - aux_en.e_tot_tendency_precip_formation[k] = mph.e_tot_tendency * aux_en.area[k] + aux_en.θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_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] @@ -70,9 +69,9 @@ end function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt::Real) env_len = 8 - src_len = 9 + src_len = 8 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, i_Se_tot = 1:src_len + i_SH_qt, i_Sqt_H, i_SH_H, i_Sqt_qt, i_Sqt, i_SH, i_Sqr, i_Sqs = 1:src_len thermo_params = TCP.thermodynamics_params(param_set) quadrature_type = en_thermo.quadrature_type @@ -84,8 +83,7 @@ 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, zc, precip_frac = - vars + UnPack.@unpack qt′qt′, qt_mean, θl′θl′, θl_mean, θl′qt′, subdomain_area, q_rai, q_sno, ρ_c, p_c, precip_frac = vars FT = eltype(ρ_c) @@ -164,18 +162,8 @@ function quad_loop(en_thermo::SGSQuadrature, precip_model, vars, param_set, Δt: q_ice_en = TD.ice_specific_humidity(thermo_params, ts) T = TD.air_temperature(thermo_params, ts) # autoconversion and accretion - mph = precipitation_formation( - param_set, - precip_model, - q_rai, - q_sno, - subdomain_area, - ρ_c, - zc, - Δt, - ts, - precip_frac, - ) + mph = + precipitation_formation(param_set, precip_model, q_rai, q_sno, subdomain_area, ρ_c, Δt, ts, precip_frac) # environmental variables inner_env[i_ql] += q_liq_en * weights[m_h] * sqpi_inv @@ -195,7 +183,6 @@ 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 @@ -216,7 +203,6 @@ 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], ) @@ -299,12 +285,11 @@ 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 ..._tendency_precip_formation to diagnostics aux_en.qt_tendency_precip_formation[k] = qt_tendency * aux_en.area[k] - aux_en.e_tot_tendency_precip_formation[k] = e_tot_tendency * aux_en.area[k] + aux_en.θ_liq_ice_tendency_precip_formation[k] = θ_liq_ice_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] @@ -352,7 +337,6 @@ function microphysics( prog_pr.q_sno[k], aux_en.area[k], ρ_c[k], - FT(grid.zc[k].z), Δt, ts, precip_fraction, @@ -361,7 +345,7 @@ function microphysics( # update_env_precip_tendencies # TODO: move ..._tendency_precip_formation to diagnostics aux_en.qt_tendency_precip_formation[k] = mph.qt_tendency * aux_en.area[k] - aux_en.e_tot_tendency_precip_formation[k] = mph.e_tot_tendency * aux_en.area[k] + aux_en.θ_liq_ice_tendency_precip_formation[k] = mph.θ_liq_ice_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 81245cdb7..820237bb4 100644 --- a/src/EDMF_Precipitation.jl +++ b/src/EDMF_Precipitation.jl @@ -151,7 +151,12 @@ 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.e_tot_tendency_precip_sinks[k] = -S_qr_evap * (I_l + Φ) - S_qs_sub_dep * (I_i + Φ) + S_qs_melt * L_f + 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) + ) end return nothing end diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl index 1125cd296..e5086e8ae 100644 --- a/src/EDMF_Updrafts.jl +++ b/src/EDMF_Updrafts.jl @@ -91,14 +91,12 @@ function compute_precipitation_formation_tendencies( prog_pr.q_sno[k], aux_up[i].area[k], ρ_c[k], - FT(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] @@ -109,10 +107,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.e_tot_tendency_precip_formation[k] = 0 + aux_bulk.θ_liq_ice_tendency_precip_formation[k] = 0 aux_bulk.qt_tendency_precip_formation[k] = 0 @inbounds for i in 1:N_up - aux_bulk.e_tot_tendency_precip_formation[k] += aux_up[i].e_tot_tendency_precip_formation[k] + aux_bulk.θ_liq_ice_tendency_precip_formation[k] += aux_up[i].θ_liq_ice_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 767e01de9..8b48b2eb3 100755 --- a/src/EDMF_functions.jl +++ b/src/EDMF_functions.jl @@ -29,7 +29,7 @@ function compute_turbconv_tendencies!( return nothing end -function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase, param_set::APS) +function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase) N_up = n_updrafts(edmf) tendencies_gm = center_tendencies_grid_mean(state) FT = float_type(state) @@ -59,21 +59,17 @@ 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 - h_tot_gm = aux_gm.h_tot + θ_liq_ice_en = aux_en.θ_liq_ice + θ_liq_ice_gm = aux_gm.θ_liq_ice 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 = total_enthalpy(param_set, prog_gm.ρe_tot / ρ_c, ts_gm) # Compute the mass flux and associated scalar fluxes @. massflux = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm)) - @. massflux_h = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm)) * (If(aux_en.h_tot) - If(h_tot_gm)) + @. massflux_h = ρ_f * Ifae(a_en) * (w_en - toscalar(w_gm)) * (If(θ_liq_ice_en) - If(θ_liq_ice_gm)) @. massflux_qt = ρ_f * Ifae(a_en) * (w_en - toscalar(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] @@ -83,9 +79,9 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf a_up = aux_up[i].area w_up_i = aux_up_f[i].w q_tot_up = aux_up_i.q_tot - h_tot_up = aux_up_i.h_tot + θ_liq_ice_up = aux_up_i.θ_liq_ice @. aux_up_f[i].massflux = ρ_f * Ifau(a_up) * (w_up_i - toscalar(w_gm)) - @. massflux_h += ρ_f * (Ifau(a_up) * (w_up_i - toscalar(w_gm)) * (If(h_tot_up) - If(h_tot_gm))) + @. massflux_h += ρ_f * (Ifau(a_up) * (w_up_i - toscalar(w_gm)) * (If(θ_liq_ice_up) - If(θ_liq_ice_gm))) @. massflux_qt += ρ_f * (Ifau(a_up) * (w_up_i - toscalar(w_gm)) * (If(q_tot_up) - If(q_tot_gm))) end @@ -128,17 +124,17 @@ function compute_sgs_flux!(edmf::EDMFModel, grid::Grid, state::State, surf::Surf diffusive_flux_qt = aux_tc_f.diffusive_flux_qt diffusive_flux_uₕ = aux_tc_f.diffusive_flux_uₕ - sgs_flux_h_tot = aux_gm_f.sgs_flux_h_tot + sgs_flux_θ_liq_ice = aux_gm_f.sgs_flux_θ_liq_ice sgs_flux_q_tot = aux_gm_f.sgs_flux_q_tot sgs_flux_uₕ = aux_gm_f.sgs_flux_uₕ - @. sgs_flux_h_tot = diffusive_flux_h + massflux_h + @. sgs_flux_θ_liq_ice = diffusive_flux_h + massflux_h @. sgs_flux_q_tot = diffusive_flux_qt + massflux_qt @. sgs_flux_uₕ = diffusive_flux_uₕ # + massflux_u # apply surface BC as SGS flux at lowest level lg_surf = CC.Fields.local_geometry_field(axes(ρ_f))[kf_surf] - sgs_flux_h_tot[kf_surf] = surf.ρe_tot_flux + sgs_flux_θ_liq_ice[kf_surf] = surf.ρθ_liq_ice_flux sgs_flux_q_tot[kf_surf] = surf.ρq_tot_flux sgs_flux_uₕ[kf_surf] = CCG.Covariant3Vector(wvec(FT(1)), lg_surf) ⊗ @@ -191,20 +187,19 @@ function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, sur 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() @. 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] - aeKHh_tot_bc = -surf.ρe_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] 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] aeKMuₕ_bc = CCG.UVVector(aeKMu_bc, aeKMv_bc) ∇q_tot_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHq_tot_bc), top = CCO.SetDivergence(FT(0))) - ∇h_tot_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHh_tot_bc), top = CCO.SetDivergence(FT(0))) + ∇θ_liq_ice_en = CCO.DivergenceC2F(; bottom = CCO.SetDivergence(aeKHθ_liq_ice_bc), top = CCO.SetDivergence(FT(0))) # CCG.Covariant3Vector(FT(1)) ⊗ CCG.Covariant12Vector(FT(aeKMu_bc),FT(aeKMv_bc)) local_geometry_surf = CC.Fields.local_geometry_field(axes(ρ_f))[kf_surf] wvec = CC.Geometry.WVector @@ -219,7 +214,7 @@ function compute_diffusive_fluxes(edmf::EDMFModel, grid::Grid, state::State, sur ) @. 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 * ∇h_tot_en(wvec(aux_en.h_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_uₕ = -aux_tc_f.ρ_ae_KM * ∇uₕ_gm(prog_gm_uₕ) if edmf.moisture_model isa NonEquilibriumMoisture @@ -243,7 +238,7 @@ function affect_filter!(edmf::EDMFModel, grid::Grid, state::State, param_set::AP ### Filters ### set_edmf_surface_bc(edmf, grid, state, surf, param_set) - filter_updraft_vars(edmf, grid, state, surf, param_set) + filter_updraft_vars(edmf, grid, state, surf) @inbounds for k in real_center_indices(grid) prog_en.ρatke[k] = max(prog_en.ρatke[k], 0.0) @@ -276,7 +271,7 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su ρ_f = aux_gm_f.ρ ae_surf::FT = 1 @inbounds for i in 1:N_up - θ_surf = θ_surface_bc(surf, grid, state, edmf, i, param_set) + θ_surf = θ_surface_bc(surf, grid, state, edmf, i) 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 @@ -292,7 +287,7 @@ function set_edmf_surface_bc(edmf::EDMFModel, grid::Grid, state::State, surf::Su ae_surf -= a_surf end - flux1 = surf.shf / cp + flux1 = surf.ρθ_liq_ice_flux flux2 = surf.ρq_tot_flux zLL::FT = grid.zc[kc_surf].z ustar = surf.ustar @@ -348,28 +343,19 @@ end function uₕ_bcs() return CCO.InterpolateC2F(bottom = CCO.Extrapolate(), top = CCO.Extrapolate()) end -function θ_surface_bc( - surf::SurfaceBase{FT}, - grid::Grid, - state::State, - edmf::EDMFModel, - i::Int, - param_set::APS, -)::FT where {FT} - thermo_params = TCP.thermodynamics_params(param_set) +function θ_surface_bc(surf::SurfaceBase{FT}, grid::Grid, state::State, edmf::EDMFModel, i::Int)::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(thermo_params, ts_gm[kc_surf]) UnPack.@unpack ustar, zLL, oblength, ρLL = surface_helper(surf, grid, state) surf.bflux > 0 || return FT(0) set_src_seed = edmf.set_src_seed a_total = edmf.surface_area a_ = area_surface_bc(surf, edmf, i) - ρθ_liq_ice_flux = surf.shf / c_p # assuming no ql,qi flux + ρθ_liq_ice_flux = surf.ρθ_liq_ice_flux # 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, set_src_seed) @@ -622,7 +608,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, param_set::APS) +function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBase) N_up = n_updrafts(edmf) kc_surf = kc_surface(grid) kf_surf = kf_surface(grid) @@ -700,7 +686,7 @@ function filter_updraft_vars(edmf::EDMFModel, grid::Grid, state::State, surf::Su @. 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].ρ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, param_set) + θ_surf = θ_surface_bc(surf, grid, state, edmf, i) 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 diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl index 5c8cbc084..21480e6b5 100644 --- a/src/TurbulenceConvection.jl +++ b/src/TurbulenceConvection.jl @@ -85,7 +85,7 @@ function debug_state(state, code_location::String) ###### vars_positive = [ - vec(prog_gm.ρe_tot), + vec(prog_gm.ρθ_liq_ice), vec(prog_gm_f.w), vec(prog_up[1].ρarea), vec(prog_up[1].ρaθ_liq_ice), diff --git a/src/diagnostics.jl b/src/diagnostics.jl index 50fad45ff..b020d6eac 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_e_tot_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).e_tot_tendency_precip_formation), + "updraft_thetal_precip" => (; dims = ("zc", "t"), group = "profiles", field = state -> center_aux_bulk(state).θ_liq_ice_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 425d8424a..84a5ea2aa 100644 --- a/src/microphysics_coupling.jl +++ b/src/microphysics_coupling.jl @@ -51,7 +51,6 @@ function precipitation_formation( qs::FT, area::FT, ρ::FT, - z::FT, Δt::Real, ts, precip_fraction, @@ -67,7 +66,6 @@ function precipitation_formation( qr_tendency = FT(0) qs_tendency = FT(0) θ_liq_ice_tendency = FT(0) - e_tot_tendency = FT(0) if area > 0 @@ -77,10 +75,8 @@ function precipitation_formation( c_pm = TD.cp_m(thermo_params, ts) L_v0 = TCP.LH_v0(param_set) L_s0 = TCP.LH_s0(param_set) - I_l = TD.internal_energy_liquid(thermo_params, ts) I_i = TD.internal_energy_ice(thermo_params, ts) I = TD.internal_energy(thermo_params, ts) - Φ = geopotential(param_set, z) if precip_model isa Clima0M qsat = TD.q_vap_saturation(thermo_params, ts) @@ -94,7 +90,6 @@ 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 + Φ) * S_qt end if precip_model isa Clima1M @@ -120,7 +115,6 @@ 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 + Φ) + S_qt_snow * (I_i + Φ) # accretion cloud water + rain S_qr = min(q.liq / Δt, CM1.accretion(microphys_params, liq_type, rain_type, q.liq, qr, ρ)) * precip_fraction @@ -128,7 +122,6 @@ 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 + Φ) # accretion cloud ice + snow S_qs = min(q.ice / Δt, CM1.accretion(microphys_params, ice_type, snow_type, q.ice, qs, ρ)) * precip_fraction @@ -136,7 +129,6 @@ 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 + Φ) # sink of cloud water via accretion cloud water + snow S_qt = @@ -146,7 +138,6 @@ 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 + Φ) else # snow melts, both cloud water and snow become rain α::FT = c_vl / Lf * (T - T_fr) qt_tendency += S_qt @@ -154,7 +145,6 @@ 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 + Φ) end # sink of cloud ice via accretion cloud ice - rain @@ -167,8 +157,6 @@ 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 + Φ) - e_tot_tendency -= S_qr * Lf # accretion rain - snow if T < T_fr @@ -185,16 +173,7 @@ 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, - e_tot_tendency, - qt_tendency, - ql_tendency, - qi_tendency, - qr_tendency, - qs_tendency, - ) + return PrecipFormation{FT}(θ_liq_ice_tendency, qt_tendency, ql_tendency, qi_tendency, qr_tendency, qs_tendency) end diff --git a/src/thermodynamics.jl b/src/thermodynamics.jl index 074ec5c4a..e97f47198 100644 --- a/src/thermodynamics.jl +++ b/src/thermodynamics.jl @@ -12,57 +12,12 @@ function thermo_state_pθq(param_set::APS, p::FT, θ_liq_ice::FT, q_tot::FT, q_l return TD.PhaseNonEquil_pθq(thermo_params, 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 = () - thermo_params = TCP.thermodynamics_params(param_set) - return TD.PhaseEquil_peq(thermo_params, 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) - thermo_params = TCP.thermodynamics_params(param_set) - return TD.PhaseNonEquil_peq(thermo_params, p, e_int, q, config...) -end - -function thermo_state_phq(param_set::APS, p::FT, h::FT, q_tot::FT) where {FT} - config = () - thermo_params = TCP.thermodynamics_params(param_set) - return TD.PhaseEquil_phq(thermo_params, p, h, q_tot, config...) -end - -function thermo_state_phq(param_set::APS, p::FT, h::FT, q_tot::FT, q_liq::FT, q_ice::FT) where {FT} - config = () - q = TD.PhasePartition(q_tot, q_liq, q_ice) - thermo_params = TCP.thermodynamics_params(param_set) - return TD.PhaseNonEquil_phq(thermo_params, p, h, q, config...) -end - function geopotential(param_set, z::Real) FT = eltype(param_set) grav = FT(TCP.grav(param_set)) return grav * z end -# TODO: move to Thermodynamics.jl -function total_enthalpy(param_set::APS, e_tot::FT, ts) where {FT} - thermo_params = TCP.thermodynamics_params(param_set) - Rm = TD.gas_constant_air(thermo_params, ts) - T = TD.air_temperature(thermo_params, ts) - return e_tot + Rm * T -end - -function total_enthalpy(param_set::APS, e_tot::FT, p::FT, ρ::FT) where {FT} - return e_tot + p / ρ -end - -function enthalpy(h_tot::FT, e_kin::FT, e_pot::FT) where {FT} - return h_tot - e_kin - e_pot -end - function enthalpy(mse::FT, e_pot::FT) where {FT} return mse - e_pot end diff --git a/src/types.jl b/src/types.jl index 45e71c4b2..c02927feb 100644 --- a/src/types.jl +++ b/src/types.jl @@ -7,7 +7,6 @@ $(DocStringExtensions.FIELDS) """ Base.@kwdef struct PrecipFormation{FT} θ_liq_ice_tendency::FT - e_tot_tendency::FT qt_tendency::FT ql_tendency::FT qi_tendency::FT @@ -397,7 +396,7 @@ Base.@kwdef struct SurfaceBase{FT} ρq_tot_flux::FT = 0 ρq_liq_flux::FT = 0 ρq_ice_flux::FT = 0 - ρe_tot_flux::FT = 0 + ρθ_liq_ice_flux::FT = 0 ρu_flux::FT = 0 ρv_flux::FT = 0 obukhov_length::FT = 0 diff --git a/src/update_aux.jl b/src/update_aux.jl index 0ee4c743c..50f17c299 100644 --- a/src/update_aux.jl +++ b/src/update_aux.jl @@ -44,12 +44,6 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas ##### ##### center variables ##### - C123 = CCG.Covariant123Vector - @. aux_en.e_kin = LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(wvec(aux_en_f.w)))) / 2 - - @inbounds for i in 1:N_up - @. aux_up[i].e_kin = LA.norm_sqr(C123(prog_gm_uₕ) + C123(Ic(wvec(aux_up_f[i].w)))) / 2 - end @inbounds for k in real_center_indices(grid) ##### @@ -65,7 +59,6 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas aux_up[i].θ_liq_ice[k] = aux_gm.θ_liq_ice[k] aux_up[i].q_tot[k] = aux_gm.q_tot[k] aux_up[i].area[k] = 0 - aux_up[i].e_kin[k] = aux_gm.e_kin[k] end thermo_args = () if edmf.moisture_model isa NonEquilibriumMoisture @@ -79,15 +72,12 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas thermo_args = (aux_up[i].q_liq[k], aux_up[i].q_ice[k]) end ts_up_i = thermo_state_pθq(param_set, p_c[k], aux_up[i].θ_liq_ice[k], aux_up[i].q_tot[k], thermo_args...) - aux_up[i].e_tot[k] = TD.total_energy(thermo_params, ts_up_i, aux_up[i].e_kin[k], e_pot) - aux_up[i].h_tot[k] = total_enthalpy(param_set, aux_up[i].e_tot[k], ts_up_i) end ##### ##### compute bulk ##### aux_bulk.q_tot[k] = 0 - aux_bulk.h_tot[k] = 0 aux_bulk.θ_liq_ice[k] = 0 aux_bulk.area[k] = sum(i -> aux_up[i].area[k], 1:N_up) if aux_bulk.area[k] > 0 @@ -96,12 +86,10 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas a_bulk_k = aux_bulk.area[k] aux_bulk.q_tot[k] += a_k * aux_up[i].q_tot[k] / a_bulk_k aux_bulk.θ_liq_ice[k] += a_k * aux_up[i].θ_liq_ice[k] / a_bulk_k - aux_bulk.h_tot[k] += a_k * aux_up[i].h_tot[k] / a_bulk_k end else aux_bulk.q_tot[k] = aux_gm.q_tot[k] aux_bulk.θ_liq_ice[k] = aux_gm.θ_liq_ice[k] - aux_bulk.h_tot[k] = aux_gm.h_tot[k] end if edmf.moisture_model isa NonEquilibriumMoisture aux_bulk.q_liq[k] = 0 @@ -133,7 +121,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas val1 = 1 / (1 - a_bulk_c) val2 = a_bulk_c * val1 aux_en.q_tot[k] = max(val1 * aux_gm.q_tot[k] - val2 * aux_bulk.q_tot[k], 0) #Yair - this is here to prevent negative QT - aux_en.h_tot[k] = val1 * aux_gm.h_tot[k] - val2 * aux_bulk.h_tot[k] + aux_en.θ_liq_ice[k] = val1 * aux_gm.θ_liq_ice[k] - val2 * aux_bulk.θ_liq_ice[k] if edmf.moisture_model isa NonEquilibriumMoisture aux_en.q_liq[k] = max(val1 * prog_gm.q_liq[k] - val2 * aux_bulk.q_liq[k], 0) aux_en.q_ice[k] = max(val1 * prog_gm.q_ice[k] - val2 * aux_bulk.q_ice[k], 0) @@ -150,11 +138,8 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas error("Something went wrong. The moisture_model options are equilibrium or nonequilibrium") end - h_en = enthalpy(aux_en.h_tot[k], e_pot, aux_en.e_kin[k]) - ts_env[k] = thermo_state_phq(param_set, p_c[k], h_en, aux_en.q_tot[k], thermo_args...) + ts_env[k] = thermo_state_pθq(param_set, p_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k], thermo_args...) ts_en = ts_env[k] - aux_en.θ_liq_ice[k] = TD.liquid_ice_pottemp(thermo_params, ts_en) - aux_en.e_tot[k] = TD.total_energy(thermo_params, ts_en, aux_en.e_kin[k], e_pot) aux_en.T[k] = TD.air_temperature(thermo_params, ts_en) aux_en.θ_virt[k] = TD.virtual_pottemp(thermo_params, ts_en) aux_en.θ_dry[k] = TD.dry_pottemp(thermo_params, ts_en) @@ -328,6 +313,7 @@ function update_aux!(edmf::EDMFModel, grid::Grid, state::State, surf::SurfaceBas ∇0_bcs = (; bottom = CCO.Extrapolate(), top = CCO.Extrapolate()) If0 = CCO.InterpolateC2F(; ∇0_bcs...) + C123 = CCG.Covariant123Vector uₕ_gm = grid_mean_uₕ(state) w_en = aux_en_f.w diff --git a/src/variables.jl b/src/variables.jl index d708c8d6d..037a321be 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -38,11 +38,7 @@ cent_aux_vars_up(FT, local_geometry, edmf) = (; area = FT(0), q_tot = FT(0), θ_liq_ice = FT(0), - e_tot = FT(0), - e_kin = FT(0), - h_tot = 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), @@ -83,8 +79,6 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (; bulk = (; area = FT(0), θ_liq_ice = FT(0), - e_tot = FT(0), - h_tot = FT(0), RH = FT(0), buoy = FT(0), q_tot = FT(0), @@ -92,7 +86,7 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (; q_ice = FT(0), T = FT(0), cloud_fraction = FT(0), - e_tot_tendency_precip_formation = FT(0), + θ_liq_ice_tendency_precip_formation = FT(0), qt_tendency_precip_formation = FT(0), cent_aux_vars_edmf_bulk_moisture(FT, edmf.moisture_model)..., ), @@ -105,9 +99,6 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (; q_liq = FT(0), q_ice = FT(0), θ_liq_ice = FT(0), - e_tot = FT(0), - e_kin = FT(0), - h_tot = FT(0), θ_virt = FT(0), θ_dry = FT(0), RH = FT(0), @@ -119,8 +110,7 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (; Hvar = FT(0), QTvar = FT(0), HQTcov = FT(0), - #θ_liq_ice_tendency_precip_formation = FT(0), - e_tot_tendency_precip_formation = FT(0), + θ_liq_ice_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)), @@ -130,7 +120,6 @@ cent_aux_vars_edmf(::Type{FT}, local_geometry, edmf) where {FT} = (; 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),