diff --git a/config/model_configs/diagnostic_edmfx_dycoms_rf02_2M_box.yml b/config/model_configs/diagnostic_edmfx_dycoms_rf02_2M_box.yml new file mode 100644 index 0000000000..d3592f4926 --- /dev/null +++ b/config/model_configs/diagnostic_edmfx_dycoms_rf02_2M_box.yml @@ -0,0 +1,41 @@ +initial_condition: DYCOMS_RF02 +subsidence: DYCOMS +scm_coriolis: DYCOMS_RF02 +rad: DYCOMS +surface_setup: DYCOMS_RF02 +turbconv: diagnostic_edmfx +implicit_diffusion: true +approximate_linear_solve_iters: 2 +prognostic_tke: true +edmfx_upwinding: first_order +edmfx_entr_model: "Generalized" +edmfx_detr_model: "Generalized" +edmfx_nh_pressure: true +edmfx_sgs_mass_flux: true +edmfx_sgs_diffusive_flux: true +moist: nonequil +cloud_model: "quadrature_sgs" +precip_model: "2M" +call_cloud_diagnostics_per_stage: true +config: box +x_max: 1e8 +y_max: 1e8 +x_elem: 2 +y_elem: 2 +z_elem: 30 +z_max: 1500 +z_stretch: false +dt: 10secs +t_end: 6hours +dt_save_state_to_disk: 10mins +toml: [toml/diagnostic_edmfx_2M.toml] +netcdf_interpolation_num_points: [8, 8, 30] +diagnostics: + - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr, cdnc, ncra] + period: 10mins + - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, husraup, hussnup, tke, ncraup, cdncup] + period: 10mins + - short_name: [entr, detr, lmix, bgrad, strain, edt, evu] + period: 10mins + - short_name: [husra, hussn] + period: 10mins diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 8a25bd5f12..78146331c9 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -1403,6 +1403,7 @@ EDMFBoxPlotsWithPrecip = Union{ DiagEDMFBoxPlotsWithPrecip = Union{ Val{:diagnostic_edmfx_dycoms_rf02_box}, + Val{:diagnostic_edmfx_dycoms_rf02_2M_box}, Val{:diagnostic_edmfx_rico_box}, Val{:diagnostic_edmfx_trmm_box}, } @@ -1511,7 +1512,20 @@ function make_plots( ("husra", "hussn", "husraup", "hussnup", "husraen", "hussnen") end elseif sim_type isa DiagEDMFBoxPlotsWithPrecip - precip_names = ("husra", "hussn", "husraup", "hussnup") + if sim_type isa Val{:diagnostic_edmfx_dycoms_rf02_2M_box} + precip_names = ( + "husra", + "hussn", + "husraup", + "hussnup", + "cdnc", + "ncra", + "cdncup", + "ncraup", + ) + else + precip_names = ("husra", "hussn", "husraup", "hussnup") + end else precip_names = () end diff --git a/src/cache/diagnostic_edmf_precomputed_quantities.jl b/src/cache/diagnostic_edmf_precomputed_quantities.jl index d7003b6d96..fa18c697db 100644 --- a/src/cache/diagnostic_edmf_precomputed_quantities.jl +++ b/src/cache/diagnostic_edmf_precomputed_quantities.jl @@ -204,7 +204,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( obukhov_length_sfc_halflevel = Fields.field_values(obukhov_length) if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment (; ᶜq_liqʲs, ᶜq_iceʲs, ᶜq_raiʲs, ᶜq_snoʲs) = p.precomputed ᶜq_liq = @. lazy(specific(Y.c.ρq_liq, Y.c.ρ)) @@ -223,6 +224,21 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( ρ_flux_q_ice_sfc_halflevel = FT(0) ρ_flux_q_rai_sfc_halflevel = FT(0) ρ_flux_q_sno_sfc_halflevel = FT(0) + + if microphysics_model isa Microphysics2Moment + (; ᶜn_liqʲs, ᶜn_raiʲs) = p.precomputed + + ᶜn_liq = @. lazy(specific(Y.c.ρn_liq, Y.c.ρ)) + ᶜn_rai = @. lazy(specific(Y.c.ρn_rai, Y.c.ρ)) + + n_liq_int_level = Fields.field_values(Fields.level(ᶜn_liq, 1)) + n_rai_int_level = Fields.field_values(Fields.level(ᶜn_rai, 1)) + + # TODO consider adding them to p.precomputed.sfc_conditions + # Though they will always be zero. + ρ_flux_n_liq_sfc_halflevel = FT(0) + ρ_flux_n_rai_sfc_halflevel = FT(0) + end end # boundary condition @@ -278,7 +294,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( ) if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment ᶜq_liqʲ = ᶜq_liqʲs.:($j) ᶜq_iceʲ = ᶜq_iceʲs.:($j) ᶜq_raiʲ = ᶜq_raiʲs.:($j) @@ -333,10 +350,41 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_bottom_bc!( obukhov_length_sfc_halflevel, local_geometry_int_halflevel, ) + if microphysics_model isa Microphysics2Moment + ᶜn_liqʲ = ᶜn_liqʲs.:($j) + ᶜn_raiʲ = ᶜn_raiʲs.:($j) + + n_liqʲ_int_level = Fields.field_values(Fields.level(ᶜn_liqʲ, 1)) + n_raiʲ_int_level = Fields.field_values(Fields.level(ᶜn_raiʲ, 1)) + + @. n_liqʲ_int_level = sgs_scalar_first_interior_bc( + z_int_level - z_sfc_halflevel, + ρ_int_level, + FT(turbconv_params.surface_area), + n_liq_int_level, + buoyancy_flux_sfc_halflevel, + ρ_flux_n_liq_sfc_halflevel, + ustar_sfc_halflevel, + obukhov_length_sfc_halflevel, + local_geometry_int_halflevel, + ) + @. n_raiʲ_int_level = sgs_scalar_first_interior_bc( + z_int_level - z_sfc_halflevel, + ρ_int_level, + FT(turbconv_params.surface_area), + n_rai_int_level, + buoyancy_flux_sfc_halflevel, + ρ_flux_n_rai_sfc_halflevel, + ustar_sfc_halflevel, + obukhov_length_sfc_halflevel, + local_geometry_int_halflevel, + ) + end end if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment set_diagnostic_edmfx_draft_quantities_level!( thermo_params, tsʲ_int_level, @@ -469,20 +517,29 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( (; ᶠu³⁰, ᶜK⁰) = p.precomputed if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment ᶜq_liq = @. lazy(specific(Y.c.ρq_liq, Y.c.ρ)) ᶜq_ice = @. lazy(specific(Y.c.ρq_ice, Y.c.ρ)) ᶜq_rai = @. lazy(specific(Y.c.ρq_rai, Y.c.ρ)) ᶜq_sno = @. lazy(specific(Y.c.ρq_sno, Y.c.ρ)) (; ᶜq_liqʲs, ᶜq_iceʲs, ᶜq_raiʲs, ᶜq_snoʲs) = p.precomputed + if microphysics_model isa Microphysics2Moment + ᶜn_liq = @. lazy(specific(Y.c.ρn_liq, Y.c.ρ)) + ᶜn_rai = @. lazy(specific(Y.c.ρn_rai, Y.c.ρ)) + + (; ᶜn_liqʲs, ᶜn_raiʲs) = p.precomputed + end end thermo_params = CAP.thermodynamics_params(params) microphys_0m_params = CAP.microphysics_0m_params(params) microphys_1m_params = CAP.microphysics_1m_params(params) + microphysics_2m_params = CAP.microphysics_2m_params(params) cloud_params = CAP.microphysics_cloud_params(params) turbconv_params = CAP.turbconv_params(params) + aerosol_params = cloud_params.aml ᶠΦ = p.scratch.ᶠtemp_scalar @. ᶠΦ = CAP.grav(params) * ᶠz @@ -545,7 +602,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( dz_prev_level = Fields.field_values(Fields.level(ᶜdz, i - 1)) if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment q_liq_level = Fields.field_values(Fields.level(ᶜq_liq, i)) q_liq_prev_level = Fields.field_values(Fields.level(ᶜq_liq, i - 1)) @@ -557,6 +615,16 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( q_sno_level = Fields.field_values(Fields.level(ᶜq_sno, i)) q_sno_prev_level = Fields.field_values(Fields.level(ᶜq_sno, i - 1)) + + if microphysics_model isa Microphysics2Moment + n_liq_level = Fields.field_values(Fields.level(ᶜn_liq, i)) + n_liq_prev_level = + Fields.field_values(Fields.level(ᶜn_liq, i - 1)) + + n_rai_level = Fields.field_values(Fields.level(ᶜn_rai, i)) + n_rai_prev_level = + Fields.field_values(Fields.level(ᶜn_rai, i - 1)) + end end local_geometry_prev_level = Fields.field_values( @@ -583,7 +651,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( ᶜS_q_totʲ = p.precomputed.ᶜSqₜᵖʲs.:($j) end if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment ᶜq_liqʲ = ᶜq_liqʲs.:($j) ᶜq_iceʲ = ᶜq_iceʲs.:($j) ᶜq_raiʲ = ᶜq_raiʲs.:($j) @@ -595,7 +664,14 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( ᶜS_q_snoʲ = p.precomputed.ᶜSqₛᵖʲs.:($j) ᶜSᵖ = p.scratch.ᶜtemp_scalar - ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2 + ᶜSᵖ_2 = p.scratch.ᶜtemp_scalar_2 + if microphysics_model isa Microphysics2Moment + ᶜn_liqʲ = ᶜn_liqʲs.:($j) + ᶜn_raiʲ = ᶜn_raiʲs.:($j) + + ᶜS_n_liqʲ = p.precomputed.ᶜSnₗᵖʲs.:($j) + ᶜS_n_raiʲ = p.precomputed.ᶜSnᵣᵖʲs.:($j) + end end ρaʲ_level = Fields.field_values(Fields.level(ᶜρaʲ, i)) @@ -636,7 +712,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( Ref(nothing) end if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment q_liqʲ_level = Fields.field_values(Fields.level(ᶜq_liqʲ, i)) q_iceʲ_level = Fields.field_values(Fields.level(ᶜq_iceʲ, i)) @@ -662,8 +739,23 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( Fields.field_values(Fields.level(ᶜS_q_snoʲ, i - 1)) Sᵖ_prev_level = Fields.field_values(Fields.level(ᶜSᵖ, i - 1)) - Sᵖ_snow_prev_level = - Fields.field_values(Fields.level(ᶜSᵖ_snow, i - 1)) + Sᵖ_2_prev_level = + Fields.field_values(Fields.level(ᶜSᵖ_2, i - 1)) + + if microphysics_model isa Microphysics2Moment + n_liqʲ_level = Fields.field_values(Fields.level(ᶜn_liqʲ, i)) + n_raiʲ_level = Fields.field_values(Fields.level(ᶜn_raiʲ, i)) + + n_liqʲ_prev_level = + Fields.field_values(Fields.level(ᶜn_liqʲ, i - 1)) + n_raiʲ_prev_level = + Fields.field_values(Fields.level(ᶜn_raiʲ, i - 1)) + + S_n_liqʲ_prev_level = + Fields.field_values(Fields.level(ᶜS_n_liqʲ, i - 1)) + S_n_raiʲ_prev_level = + Fields.field_values(Fields.level(ᶜS_n_raiʲ, i - 1)) + end end tke_prev_level = Fields.field_values(Fields.level(ᶜtke⁰, i - 1)) @@ -777,7 +869,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( # Rain formation from the updrafts compute_precipitation_sources!( Sᵖ_prev_level, - Sᵖ_snow_prev_level, + Sᵖ_2_prev_level, S_q_liqʲ_prev_level, S_q_iceʲ_prev_level, S_q_raiʲ_prev_level, @@ -834,7 +926,85 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( TD.air_temperature(thermo_params, tsʲ_prev_level), dt, ) + elseif moisture_model isa NonEquilMoistModel && + microphysics_model isa Microphysics2Moment + + # 2M warm rain formation + compute_warm_precipitation_sources_2M!( + Sᵖ_prev_level, + Sᵖ_2_prev_level, + S_n_liqʲ_prev_level, + S_n_raiʲ_prev_level, + S_q_liqʲ_prev_level, + S_q_raiʲ_prev_level, + ρʲ_prev_level, + n_liqʲ_prev_level, + n_raiʲ_prev_level, + q_totʲ_prev_level, + q_liqʲ_prev_level, + q_iceʲ_prev_level, + q_raiʲ_prev_level, + q_snoʲ_prev_level, + tsʲ_prev_level, + dt, + microphysics_2m_params, + thermo_params, + ) + # Snow not implemented + @. S_q_iceʲ_prev_level = FT(0) + @. S_q_snoʲ_prev_level = FT(0) + # Cloud sources from the updrafts + @. S_q_liqʲ_prev_level += cloud_sources( + cloud_params.liquid, + thermo_params, + q_totʲ_prev_level, + q_liqʲ_prev_level, + q_iceʲ_prev_level, + q_raiʲ_prev_level, + q_snoʲ_prev_level, + ρʲ_prev_level, + TD.air_temperature(thermo_params, tsʲ_prev_level), + dt, + ) + @. S_q_iceʲ_prev_level += cloud_sources( + cloud_params.ice, + thermo_params, + q_totʲ_prev_level, + q_liqʲ_prev_level, + q_iceʲ_prev_level, + q_raiʲ_prev_level, + q_snoʲ_prev_level, + ρʲ_prev_level, + TD.air_temperature(thermo_params, tsʲ_prev_level), + dt, + ) + # Aerosol activation + # TODO - right now I'm hijacking the + seasalt_num_prev_level = aerosol_params.c₀_seasalt + seasalt_mean_radius_prev_level = aerosol_params.α_seasalt + sulfate_num_prev_level = FT(0) + + @. S_n_liqʲ_prev_level += aerosol_activation_sources( + seasalt_num_prev_level, + seasalt_mean_radius_prev_level, + sulfate_num_prev_level, + q_totʲ_prev_level, + q_liqʲ_prev_level + q_raiʲ_prev_level, + q_iceʲ_prev_level + q_snoʲ_prev_level, + n_liqʲ_prev_level + n_raiʲ_prev_level, + ρʲ_prev_level, + max( + 0, + w_component.( + Geometry.WVector.(u³ʲ_data_prev_halflevel) + ), + ), + (microphysics_2m_params,), + thermo_params, + tsʲ_prev_level, + dt, + ) end u³ʲ_datau³ʲ_data = p.scratch.temp_data_level @@ -1047,7 +1217,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( ### ### 1-moment microphysics tracers ### - if microphysics_model isa Microphysics1Moment && + if microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment && moisture_model isa NonEquilMoistModel # TODO - loop over tracres ρaʲu³ʲ_dataq_ = p.scratch.temp_data_level_3 @@ -1170,6 +1341,72 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( q_sno_level, ρaʲu³ʲ_dataq_ / ρaʲu³ʲ_data, ) + + ### + ### Additional 2-moment microphysics tracers + ### + if microphysics_model isa Microphysics2Moment + # Cloud droplets + @. ρaʲu³ʲ_dataq_ = + diag_edmf_advection( + local_geometry_halflevel.J, + local_geometry_prev_halflevel.J, + ρaʲ_prev_level, + u³ʲ_data_prev_halflevel, + n_liqʲ_prev_level, + ) + + entr_detr( + local_geometry_halflevel.J, + local_geometry_prev_level.J, + ρaʲ_prev_level, + entrʲ_prev_level, + detrʲ_prev_level, + turb_entrʲ_prev_level, + n_liq_prev_level, + n_liqʲ_prev_level, + ) + + microphysics_sources( + local_geometry_halflevel.J, + local_geometry_prev_level.J, + ρaʲ_prev_level, + S_n_liqʲ_prev_level, + ) + @. n_liqʲ_level = ifelse( + kill_updraft, + n_liq_level, + ρaʲu³ʲ_dataq_ / ρaʲu³ʲ_data, + ) + # Rain drops + @. ρaʲu³ʲ_dataq_ = + diag_edmf_advection( + local_geometry_halflevel.J, + local_geometry_prev_halflevel.J, + ρaʲ_prev_level, + u³ʲ_data_prev_halflevel, + n_raiʲ_prev_level, + ) + + entr_detr( + local_geometry_halflevel.J, + local_geometry_prev_level.J, + ρaʲ_prev_level, + entrʲ_prev_level, + detrʲ_prev_level, + turb_entrʲ_prev_level, + n_rai_prev_level, + n_raiʲ_prev_level, + ) + + microphysics_sources( + local_geometry_halflevel.J, + local_geometry_prev_level.J, + ρaʲ_prev_level, + S_n_raiʲ_prev_level, + ) + @. n_raiʲ_level = ifelse( + kill_updraft, + n_rai_level, + ρaʲu³ʲ_dataq_ / ρaʲu³ʲ_data, + ) + end end # set updraft to grid-mean if vertical velocity is too small @@ -1186,7 +1423,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( ifelse(kill_updraft_2, h_tot_level - K_level, mseʲ_level) @. q_totʲ_level = ifelse(kill_updraft_2, q_tot_level, q_totʲ_level) - if microphysics_model isa Microphysics1Moment && + if microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment && moisture_model isa NonEquilMoistModel @. q_liqʲ_level = ifelse(kill_updraft_2, q_liq_level, q_liqʲ_level) @@ -1196,6 +1434,12 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( ifelse(kill_updraft_2, q_rai_level, q_raiʲ_level) @. q_snoʲ_level = ifelse(kill_updraft_2, q_sno_level, q_snoʲ_level) + if microphysics_model isa Microphysics2Moment + @. n_liqʲ_level = + ifelse(kill_updraft_2, n_liq_level, n_liqʲ_level) + @. n_raiʲ_level = + ifelse(kill_updraft_2, n_rai_level, n_raiʲ_level) + end end end @@ -1206,7 +1450,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!( local_geometry_halflevel, ) if moisture_model isa NonEquilMoistModel && - microphysics_model isa Microphysics1Moment + microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment set_diagnostic_edmfx_draft_quantities_level!( thermo_params, tsʲ_level, @@ -1310,7 +1555,8 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!( S_q_totʲ_level = Fields.field_values(Fields.level(ᶜS_q_totʲ, i_top)) @. S_q_totʲ_level = 0 end - if microphysics_model isa Microphysics1Moment && + if microphysics_model isa Microphysics1Moment || + microphysics_model isa Microphysics2Moment && moisture_model isa NonEquilMoistModel ᶜS_q_liqʲ = p.precomputed.ᶜSqₗᵖʲs.:($j) ᶜS_q_iceʲ = p.precomputed.ᶜSqᵢᵖʲs.:($j) @@ -1324,6 +1570,16 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_top_bc!( @. S_q_iceʲ_level = 0 @. S_q_raiʲ_level = 0 @. S_q_snoʲ_level = 0 + if microphysics_model isa Microphysics2Moment + ᶜS_n_liqʲ = p.precomputed.ᶜSnₗᵖʲs.:($j) + ᶜS_n_raiʲ = p.precomputed.ᶜSnᵣᵖʲs.:($j) + S_n_liqʲ_level = + Fields.field_values(Fields.level(ᶜS_n_liqʲ, i_top)) + S_n_raiʲ_level = + Fields.field_values(Fields.level(ᶜS_n_raiʲ, i_top)) + @. S_n_liqʲ_level = 0 + @. S_n_raiʲ_level = 0 + end end end return nothing @@ -1467,7 +1723,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita microphys_1m_params, thermo_params, ) - # Rain sinks from the updrafts + # Rain sinks from the environment compute_precipitation_sinks!( ᶜSᵖ, ᶜSqᵣᵖ⁰, @@ -1483,7 +1739,7 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita microphys_1m_params, thermo_params, ) - # Cloud formation from the updrafts + # Cloud formation from the environment @. ᶜSqₗᵖ⁰ += cloud_sources( cloud_params.liquid, thermo_params, @@ -1510,3 +1766,103 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita ) return nothing end +NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipitation!( + Y, + p, + t, + microphysics_model::Microphysics2Moment, +) + thermo_params = CAP.thermodynamics_params(p.params) + microphysics_2m_params = CAP.microphysics_2m_params(p.params) + cloud_params = CAP.microphysics_cloud_params(p.params) + turbconv_params = CAP.turbconv_params(p.params) + aerosol_params = cloud_params.aml + FT = eltype(thermo_params) + + (; dt) = p + + (; ᶜts, ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜSnₗᵖ⁰, ᶜSnᵣᵖ⁰) = p.precomputed + ᶜSᵖ = p.scratch.ᶜtemp_scalar + ᶜSᵖ_2 = p.scratch.ᶜtemp_scalar_2 + + # TODO - check if this is the correct velocity + if p.atmos.turbconv_model isa DiagnosticEDMFX + (; ᶜu⁰) = p.precomputed + elseif p.atmos.turbconv_model isa EDOnlyEDMFX + ᶜu⁰ = ᶜu + end + + # 2M warm rain formation + compute_warm_precipitation_sources_2M!( + ᶜSᵖ, + ᶜSᵖ_2, + ᶜSnₗᵖ⁰, + ᶜSnᵣᵖ⁰, + ᶜSqₗᵖ⁰, + ᶜSqᵣᵖ⁰, + Y.c.ρ, + specific.(Y.c.ρn_liq, Y.c.ρ), + specific.(Y.c.ρn_rai, Y.c.ρ), + specific.(Y.c.ρq_tot, Y.c.ρ), + specific.(Y.c.ρq_liq, Y.c.ρ), + specific.(Y.c.ρq_ice, Y.c.ρ), + specific.(Y.c.ρq_rai, Y.c.ρ), + specific.(Y.c.ρq_sno, Y.c.ρ), + ᶜts, + dt, + microphysics_2m_params, + thermo_params, + ) + + # Snow not implemented + @. ᶜSqᵢᵖ⁰ = FT(0) + @. ᶜSqₛᵖ⁰ = FT(0) + + # Cloud formation from the environment + @. ᶜSqₗᵖ⁰ += cloud_sources( + cloud_params.liquid, + thermo_params, + specific(Y.c.ρq_tot, Y.c.ρ), + specific(Y.c.ρq_liq, Y.c.ρ), + specific(Y.c.ρq_ice, Y.c.ρ), + specific(Y.c.ρq_rai, Y.c.ρ), + specific(Y.c.ρq_sno, Y.c.ρ), + Y.c.ρ, + TD.air_temperature(thermo_params, ᶜts), + dt, + ) + @. ᶜSqᵢᵖ⁰ += cloud_sources( + cloud_params.ice, + thermo_params, + specific(Y.c.ρq_tot, Y.c.ρ), + specific(Y.c.ρq_liq, Y.c.ρ), + specific(Y.c.ρq_ice, Y.c.ρ), + specific(Y.c.ρq_rai, Y.c.ρ), + specific(Y.c.ρq_sno, Y.c.ρ), + Y.c.ρ, + TD.air_temperature(thermo_params, ᶜts), + dt, + ) + + # Aerosol activation + # TODO - make it a user specified free parameter + seasalt_num_prev_level = aerosol_params.c₀_seasalt + seasalt_mean_radius_prev_level = aerosol_params.α_seasalt + sulfate_num_prev_level = FT(0) + + @. ᶜSnₗᵖ⁰ += aerosol_activation_sources( + seasalt_num_prev_level, + seasalt_mean_radius_prev_level, + sulfate_num_prev_level, + specific(Y.c.ρq_tot, Y.c.ρ), + specific(Y.c.ρq_liq, Y.c.ρ) + specific(Y.c.ρq_rai, Y.c.ρ), + specific(Y.c.ρq_ice, Y.c.ρ) + specific(Y.c.ρq_sno, Y.c.ρ), + specific(Y.c.ρn_liq, Y.c.ρ) + specific(Y.c.ρn_rai, Y.c.ρ), + Y.c.ρ, + max(0, w_component.(Geometry.WVector.(ᶜu⁰))), + (microphysics_2m_params,), + thermo_params, + ᶜts, + dt, + ) +end diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl index 3dd66e8bcb..15c69436a8 100644 --- a/src/cache/precipitation_precomputed_quantities.jl +++ b/src/cache/precipitation_precomputed_quantities.jl @@ -411,7 +411,8 @@ function set_precipitation_cache!( ::Microphysics2Moment, ::DiagnosticEDMFX, ) - error("Not implemented yet") + # Nothing needs to be done on the grid mean. The Sources are computed + # in edmf sub-domains. return nothing end function set_precipitation_cache!( diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index 05e2972609..096f8b6046 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -138,12 +138,12 @@ function precomputed_quantities(Y, atmos) ᶜSqᵢᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSqᵣᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSqₛᵖʲs = similar(Y.c, NTuple{n, FT}), - ᶜwₗʲs = similar(Y.c, NTuple{n, FT}), - ᶜwᵢʲs = similar(Y.c, NTuple{n, FT}), - ᶜwᵣʲs = similar(Y.c, NTuple{n, FT}), - ᶜwₛʲs = similar(Y.c, NTuple{n, FT}), - ᶜwₜʲs = similar(Y.c, NTuple{n, FT}), - ᶜwₕʲs = similar(Y.c, NTuple{n, FT}), + ᶜwₗʲs = similar(Y.c, NTuple{n, FT}), # TODO I don't think I need those + ᶜwᵢʲs = similar(Y.c, NTuple{n, FT}), # + ᶜwᵣʲs = similar(Y.c, NTuple{n, FT}), # + ᶜwₛʲs = similar(Y.c, NTuple{n, FT}), # + ᶜwₜʲs = similar(Y.c, NTuple{n, FT}), # + ᶜwₕʲs = similar(Y.c, NTuple{n, FT}), # ᶜSqₗᵖ⁰ = similar(Y.c, FT), ᶜSqᵢᵖ⁰ = similar(Y.c, FT), ᶜSqᵣᵖ⁰ = similar(Y.c, FT), @@ -157,8 +157,8 @@ function precomputed_quantities(Y, atmos) ᶜSqₛᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSnₗᵖʲs = similar(Y.c, NTuple{n, FT}), ᶜSnᵣᵖʲs = similar(Y.c, NTuple{n, FT}), - ᶜwₗʲs = similar(Y.c, NTuple{n, FT}), - ᶜwᵢʲs = similar(Y.c, NTuple{n, FT}), + ᶜwₗʲs = similar(Y.c, NTuple{n, FT}), # TODO - I don't think I need those + ᶜwᵢʲs = similar(Y.c, NTuple{n, FT}), # for diagnostci EDMF ᶜwᵣʲs = similar(Y.c, NTuple{n, FT}), ᶜwₛʲs = similar(Y.c, NTuple{n, FT}), ᶜwₜʲs = similar(Y.c, NTuple{n, FT}), @@ -206,8 +206,16 @@ function precomputed_quantities(Y, atmos) ᶜq_iceʲs = similar(Y.c, NTuple{n, FT}), ᶜq_raiʲs = similar(Y.c, NTuple{n, FT}), ᶜq_snoʲs = similar(Y.c, NTuple{n, FT}), + ) : + atmos.microphysics_model isa Microphysics2Moment ? + (; + ᶜq_liqʲs = similar(Y.c, NTuple{n, FT}), + ᶜq_iceʲs = similar(Y.c, NTuple{n, FT}), + ᶜq_raiʲs = similar(Y.c, NTuple{n, FT}), + ᶜq_snoʲs = similar(Y.c, NTuple{n, FT}), + ᶜn_liqʲs = similar(Y.c, NTuple{n, FT}), + ᶜn_raiʲs = similar(Y.c, NTuple{n, FT}), ) : (;) - diagnostic_sgs_quantities = atmos.turbconv_model isa DiagnosticEDMFX ? (; diff --git a/src/diagnostics/core_diagnostics.jl b/src/diagnostics/core_diagnostics.jl index 63eb0742cc..04b637c511 100644 --- a/src/diagnostics/core_diagnostics.jl +++ b/src/diagnostics/core_diagnostics.jl @@ -915,7 +915,7 @@ add_diagnostic_variable!( standard_name = "number_concentration_of_cloud_liquid_water_particles_in_air", units = "m^-3", comments = """ - This is calculated as the number of cloud liquid water droplets in the grid + This is calculated as the number of cloud liquid water droplets in the grid cell divided by the cell volume. """, compute! = compute_cdnc!, @@ -1479,7 +1479,7 @@ function compute_cape!(out, state, cache, time) ᶜbuoyancy = cache.scratch.ᶜtemp_scalar ᶜbuoyancy .= g .* (parcel_Tv .- env_Tv) ./ env_Tv - # restrict to tropospheric buoyancy (generously below 20km) TODO: integrate from LFC to LNB + # restrict to tropospheric buoyancy (generously below 20km) TODO: integrate from LFC to LNB FT = Spaces.undertype(axes(ᶜbuoyancy)) ᶜbuoyancy .= ᶜbuoyancy .* @@ -1553,7 +1553,7 @@ function compute_rwp!( state, cache, time, - moisture_model::T, + microphysics_model::T, ) where {T <: Union{Microphysics1Moment, Microphysics2Moment}} if isnothing(out) out = zeros(axes(Fields.level(state.f, half))) diff --git a/src/diagnostics/edmfx_diagnostics.jl b/src/diagnostics/edmfx_diagnostics.jl index ba072e4add..a4eceeeb9c 100644 --- a/src/diagnostics/edmfx_diagnostics.jl +++ b/src/diagnostics/edmfx_diagnostics.jl @@ -411,7 +411,7 @@ function compute_cdncup!( state, cache, time, - microphysics_model_model::Microphysics2Moment, + microphysics_model::Microphysics2Moment, turbconv_model::PrognosticEDMFX, ) if isnothing(out) @@ -420,6 +420,20 @@ function compute_cdncup!( out .= (state.c.sgsʲs.:1).n_liq end end +function compute_cdncup!( + out, + state, + cache, + time, + microphysics_model::Microphysics2Moment, + turbconv_model::DiagnosticEDMFX, +) + if isnothing(out) + return (cache.precomputed.ᶜn_liqʲs.:1) + else + out .= (cache.precomputed.ᶜn_liqʲs.:1) + end +end add_diagnostic_variable!( short_name = "cdncup", @@ -554,7 +568,7 @@ function compute_husraup!( state, cache, time, - moisture_model::Microphysics1Moment, + moisture_model::Union{Microphysics1Moment, Microphysics2Moment}, turbconv_model::DiagnosticEDMFX, ) if isnothing(out) @@ -608,6 +622,20 @@ function compute_ncraup!( out .= (state.c.sgsʲs.:1).n_rai end end +function compute_ncraup!( + out, + state, + cache, + time, + moisture_model::Microphysics2Moment, + turbconv_model::DiagnosticEDMFX, +) + if isnothing(out) + return (cache.precomputed.ᶜn_raiʲs.:1) + else + out .= (cache.precomputed.ᶜn_raiʲs.:1) + end +end add_diagnostic_variable!( short_name = "ncraup", @@ -661,7 +689,7 @@ function compute_hussnup!( state, cache, time, - moisture_model::Microphysics1Moment, + moisture_model::Union{Microphysics1Moment, Microphysics2Moment}, turbconv_model::DiagnosticEDMFX, ) if isnothing(out) diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index d8c439882e..679588b9d6 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -222,7 +222,30 @@ function precipitation_tendency!( microphysics_model::Microphysics2Moment, turbconv_model::DiagnosticEDMFX, ) - error("Not implemented yet") + # Source terms from EDMFX environment + (; ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜSnₗᵖ⁰, ᶜSnᵣᵖ⁰) = p.precomputed + # Source terms from EDMFX updrafts + (; ᶜSqₗᵖʲs, ᶜSqᵢᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs, ᶜSnₗᵖʲs, ᶜSnᵣᵖʲs) = p.precomputed + (; ᶜρaʲs) = p.precomputed + + # Update from environment precipitation sources + @. Yₜ.c.ρq_liq += Y.c.ρ * ᶜSqₗᵖ⁰ + @. Yₜ.c.ρq_ice += Y.c.ρ * ᶜSqᵢᵖ⁰ + @. Yₜ.c.ρq_rai += Y.c.ρ * ᶜSqᵣᵖ⁰ + @. Yₜ.c.ρq_sno += Y.c.ρ * ᶜSqₛᵖ⁰ + @. Yₜ.c.ρn_liq += Y.c.ρ * ᶜSnₗᵖ⁰ + @. Yₜ.c.ρn_rai += Y.c.ρ * ᶜSnᵣᵖ⁰ + + # Update from the updraft precipitation sources + n = n_mass_flux_subdomains(p.atmos.turbconv_model) + for j in 1:n + @. Yₜ.c.ρq_liq += ᶜρaʲs.:($$j) * ᶜSqₗᵖʲs.:($$j) + @. Yₜ.c.ρq_ice += ᶜρaʲs.:($$j) * ᶜSqᵢᵖʲs.:($$j) + @. Yₜ.c.ρq_rai += ᶜρaʲs.:($$j) * ᶜSqᵣᵖʲs.:($$j) + @. Yₜ.c.ρq_sno += ᶜρaʲs.:($$j) * ᶜSqₛᵖʲs.:($$j) + @. Yₜ.c.ρn_liq += ᶜρaʲs.:($$j) * ᶜSnₗᵖʲs.:($$j) + @. Yₜ.c.ρn_rai += ᶜρaʲs.:($$j) * ᶜSnᵣᵖʲs.:($$j) + end end function precipitation_tendency!( Yₜ, diff --git a/toml/diagnostic_edmfx_2M.toml b/toml/diagnostic_edmfx_2M.toml new file mode 100644 index 0000000000..751a602ab5 --- /dev/null +++ b/toml/diagnostic_edmfx_2M.toml @@ -0,0 +1,32 @@ +[condensation_evaporation_timescale] +value = 100 + +[sublimation_deposition_timescale] +value = 100 + +[entr_inv_tau] +value = 0.002 + +[entr_coeff] +value = 0 + +[detr_inv_tau] +value = 0 + +[detr_vertdiv_coeff] +value = 0.6 + +[detr_buoy_coeff] +value = 0.12 + +[min_area_limiter_scale] +value = 0 + +[max_area_limiter_scale] +value = 0 + +[seasalt_calibration_coefficient] +value = 1e-8 + +[reference_seasalt_aerosol_mass_concentration] +value = 1e10