From 2b64ce499957db85713b7474a11ba818ffbb84a0 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 23 Mar 2023 18:39:32 -0700 Subject: [PATCH 01/16] Integration of number density with the 2m parametrization of seifert and beheng 2006 --- .buildkite/pipeline.yml | 8 ++ Project.toml | 1 - src/calibrateKiD/KiDUtils.jl | 10 ++ src/calibrateKiD/ReferenceModels.jl | 4 + src/driveKiD/EquationTypes.jl | 4 + src/driveKiD/KiD_model.jl | 33 ++++- src/driveKiD/initial_condition.jl | 8 ++ src/driveKiD/tendency.jl | 128 ++++++++++++++++++ test/Project.toml | 2 +- test/experiments/KiD_driver/KiD_driver.jl | 8 +- .../KiD_driver/parse_commandline.jl | 8 +- test/experiments/calibrations/config.jl | 2 +- test/plotting_utils.jl | 6 +- test/unit_tests/unit_test.jl | 8 ++ 14 files changed, 217 insertions(+), 13 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 6ef6bad0..0f6f1fd7 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -92,6 +92,14 @@ steps: command: "julia --color=yes --project=test test/experiments/KiD_driver/KiD_driver.jl --moisture_choice=EquilibriumMoisture --prognostic_vars=RhodTQ --precipitation_choice=Precipitation1M --rain_formation_scheme_choice LD2004 --prescribed_Nd 1e8" artifact_paths: "test/experiments/KiD_driver/Output_EquilibriumMoisture_RhodTQ_Precipitation1M_LD2004/figures/*" + - label: ":crystal_ball: Experiments: equil (fixed rhod and T) + SB2006 " + command: "julia --color=yes --project=test test/experiments/KiD_driver/KiD_driver.jl --moisture_choice=EquilibriumMoisture --prognostic_vars=RhodTQ --precipitation_choice=Precipitation2M --rain_formation_scheme_choice SB2006 --prescribed_Nd 1e8" + artifact_paths: "test/experiments/KiD_driver/Output_EquilibriumMoisture_RhodTQ_Precipitation2M_SB2006/figures/*" + + - label: ":crystal_ball: Experiments: non equil (fixed rhod and T) + SB2006 " + command: "julia --color=yes --project=test test/experiments/KiD_driver/KiD_driver.jl --moisture_choice=NonEquilibriumMoisture --prognostic_vars=RhodTQ --precipitation_choice=Precipitation2M --rain_formation_scheme_choice SB2006 --prescribed_Nd 1e8" + artifact_paths: "test/experiments/KiD_driver/Output_NonEquilibriumMoisture_RhodTQ_Precipitation2M_SB2006/figures/*" + - label: ":crystal_ball: Experiments: perfect-model calibration + EKP " command: "julia --color=yes --project=test test/experiments/calibrations/run_calibration.jl" artifact_paths: "output/*" diff --git a/Project.toml b/Project.toml index 6fd7aac2..7712422e 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,6 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] ClimaCore = "0.10" ClimaCorePlots = "0.2" -CloudMicrophysics = "0.9" DiffEqBase = "6.75" EnsembleKalmanProcesses = "0.13" NCDatasets = "0.12" diff --git a/src/calibrateKiD/KiDUtils.jl b/src/calibrateKiD/KiDUtils.jl index 5281f5fc..4ea3abd2 100644 --- a/src/calibrateKiD/KiDUtils.jl +++ b/src/calibrateKiD/KiDUtils.jl @@ -155,6 +155,10 @@ function get_variable_data_from_ODE(u, ρ::Vector{Float64}, var::String) elseif var == "rlr" _qtot = parent(u.ρq_tot) ./ ρ output = (parent(u.ρq_liq) .+ parent(u.ρq_rai)) ./ ρ ./ (1 .- _qtot) + elseif var == "Nl" + output = parent(u.N_liq) + elseif var == "Nr" + output = parent(u.N_rai) elseif var == "rho" output = ρ elseif var == "rain averaged terminal velocity" @@ -240,6 +244,12 @@ function equation_types(moisture_choice::String, precipitation_choice::String, r else error("Invalid rain formation choice: $rain_formation_choice") end + elseif precipitation_choice == "Precipitation2M" + if rain_formation_choice == "SB2006" + precip = Precipitation2M(CMT.SB2006Type()) + else + error("Invalid rain formation choice: $rain_formation_choice") + end else error("Invalid precipitation choice: $precipitation_choice") end diff --git a/src/calibrateKiD/ReferenceModels.jl b/src/calibrateKiD/ReferenceModels.jl index eabdc6d3..ef1b002d 100644 --- a/src/calibrateKiD/ReferenceModels.jl +++ b/src/calibrateKiD/ReferenceModels.jl @@ -143,6 +143,10 @@ function get_single_obs_field( _data = (_data_pysdm["qc"] .+ _data_pysdm["qr"]) .* 1e-3 elseif var == "rtr" _data = _data_pysdm["qv"] .+ (_data_pysdm["qc"] .+ _data_pysdm["qr"]) .* 1e-3 + elseif var == "Nl" + _data = _data_pysdm["na"] + _data_pysdm["nc"] + elseif var == "Nr" + _data = _data_pysdm["nr"] else _data = _data_pysdm[var] end diff --git a/src/driveKiD/EquationTypes.jl b/src/driveKiD/EquationTypes.jl index 91507493..f79e88a9 100644 --- a/src/driveKiD/EquationTypes.jl +++ b/src/driveKiD/EquationTypes.jl @@ -19,6 +19,7 @@ export NonEquilibriumMoisture_ρdTq export NoPrecipitation export Precipitation0M export Precipitation1M +export Precipitation2M export OneMomentRainFormation @@ -42,5 +43,8 @@ struct Precipitation0M <: AbstractPrecipitationStyle end struct Precipitation1M{PT} <: AbstractPrecipitationStyle rain_formation::PT end +struct Precipitation2M{PT} <: AbstractPrecipitationStyle + rain_formation::PT +end struct OneMomentRainFormation <: AbstractRainFormationStyle end diff --git a/src/driveKiD/KiD_model.jl b/src/driveKiD/KiD_model.jl index 62ff3070..4b7a2a36 100644 --- a/src/driveKiD/KiD_model.jl +++ b/src/driveKiD/KiD_model.jl @@ -79,6 +79,24 @@ function initialise_state(::NonEquilibriumMoisture, ::Precipitation1M, initial_p ρq_sno = initial_profiles.ρq_sno, ) end +function initialise_state(::EquilibriumMoisture, ::Precipitation2M, initial_profiles) + return CC.Fields.FieldVector(; + ρq_tot = initial_profiles.ρq_tot, + ρq_rai = initial_profiles.ρq_rai, + N_liq = initial_profiles.N_liq, + N_rai = initial_profiles.N_rai, + ) +end +function initialise_state(::NonEquilibriumMoisture, ::Precipitation2M, initial_profiles) + return CC.Fields.FieldVector(; + ρq_tot = initial_profiles.ρq_tot, + ρq_liq = initial_profiles.ρq_liq, + ρq_ice = initial_profiles.ρq_ice, + ρq_rai = initial_profiles.ρq_rai, + N_liq = initial_profiles.N_liq, + N_rai = initial_profiles.N_rai, + ) +end """ Interface to ODE solver. It initializes the auxiliary state. @@ -119,7 +137,12 @@ function initialise_aux(FT, ip, params, TS, Stats, face_space, moisture) q_ice = ip.q_ice, ts = ts, ), - precip_variables = CC.Fields.FieldVector(; q_rai = ip.q_rai, q_sno = ip.q_sno), + precip_variables = CC.Fields.FieldVector(; + q_rai = ip.q_rai, + q_sno = ip.q_sno, + N_liq = ip.N_liq, + N_rai = ip.N_rai, + ), moisture_sources = CC.Fields.FieldVector(; S_q_liq = ip.S_ql_moisture, S_q_ice = ip.S_qi_moisture), precip_sources = CC.Fields.FieldVector(; S_q_tot = ip.S_qt_precip, @@ -127,8 +150,14 @@ function initialise_aux(FT, ip, params, TS, Stats, face_space, moisture) S_q_ice = ip.S_qi_precip, S_q_rai = ip.S_qr_precip, S_q_sno = ip.S_qs_precip, + S_N_liq = ip.S_Nl_precip, + S_N_rai = ip.S_Nr_precip, + ), + precip_velocities = CC.Fields.FieldVector(; + term_vel_rai = term_vel_rai, + term_vel_sno = term_vel_sno, + term_vel_N_rai = term_vel_rai, ), - precip_velocities = CC.Fields.FieldVector(; term_vel_rai = term_vel_rai, term_vel_sno = term_vel_sno), prescribed_velocity = CC.Fields.FieldVector(; ρw = ρw, ρw0 = ρw0), params = params, q_surf = q_surf, diff --git a/src/driveKiD/initial_condition.jl b/src/driveKiD/initial_condition.jl index f0e62ec8..8cfac58a 100644 --- a/src/driveKiD/initial_condition.jl +++ b/src/driveKiD/initial_condition.jl @@ -132,6 +132,8 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F q_sno::FT = FT(0.0) ρq_rai::FT = q_rai * ρ ρq_sno::FT = q_sno * ρ + N_liq::FT = Parameters.prescribed_Nd(params) + N_rai::FT = FT(0) S_ql_moisture::FT = FT(0.0) S_qi_moisture::FT = FT(0.0) @@ -141,6 +143,8 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F S_qi_precip::FT = FT(0.0) S_qr_precip::FT = FT(0.0) S_qs_precip::FT = FT(0.0) + S_Nl_precip::FT = FT(0) + S_Nr_precip::FT = FT(0) return (; ρ, @@ -159,6 +163,8 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F q_ice, q_rai, q_sno, + N_liq, + N_rai, S_ql_moisture, S_qi_moisture, S_qt_precip, @@ -166,5 +172,7 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F S_qi_precip, S_qr_precip, S_qs_precip, + S_Nl_precip, + S_Nr_precip, ) end diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index e1aacc9d..4c61ccfe 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -28,6 +28,12 @@ end @. dY.ρq_rai = FT(0) @. dY.ρq_sno = FT(0) end +@inline function zero_tendencies!(::Precipitation2M, dY, Y, aux, t) + FT = eltype(Y.ρq_tot) + @. dY.ρq_rai = FT(0) + @. dY.N_liq = FT(0) + @. dY.N_rai = FT(0) +end """ Helper functions for broadcasting and populating the auxiliary state @@ -294,6 +300,66 @@ end return (; S_q_tot, S_q_liq, S_q_ice, S_q_rai, S_q_sno, S_q_vap, term_vel_rai, term_vel_sno) end +@inline function precip_helper_sources_2M!(params, q_tot, q_liq, q_rai, N_liq, N_rai, T, ρ, dt, ps) + + microphys_params = KP.microphysics_params(params) + + rf = ps.rain_formation + + FT = eltype(q_tot) + + S_q_rai::FT = FT(0) + S_q_tot::FT = FT(0) + S_q_liq::FT = FT(0) + S_q_vap::FT = FT(0) # added for testing + S_N_liq::FT = FT(0) + S_N_rai::FT = FT(0) + + q = TD.PhasePartition(q_tot, q_liq, FT(0)) + + # autoconversion liquid to rain + tmp = CM2.autoconversion(microphys_params, rf, q.liq, q_rai, ρ, N_liq) + S_qr = min(max(0, q.liq / dt), tmp.dq_rai_dt) + S_q_rai += S_qr + S_q_tot -= S_qr + S_q_liq -= S_qr + S_N_liq += -min(max(0, N_liq / dt), -tmp.dN_liq_dt) + S_N_rai += min(max(0, N_liq / dt), tmp.dN_rai_dt) + + # self_collection + tmp_l = CM2.liquid_self_collection(microphys_params, rf, q.liq, ρ, -S_qr) + tmp_r = CM2.rain_self_collection(microphys_params, rf, q_rai, ρ, N_rai) + S_N_liq += -min(max(0, N_liq / dt), -tmp_l) + S_N_rai += -min(max(0, N_rai / dt), -tmp_r) + + # rain breakup + tmp = CM2.rain_breakup(microphys_params, rf, q_rai, ρ, N_rai, tmp_r) + S_N_rai += min(max(0, N_rai / dt), tmp_r) + + # accretion cloud water + rain + tmp = CM2.accretion(microphys_params, rf, q.liq, q_rai, ρ, N_liq) + S_qr = min(max(0, q.liq / dt), tmp.dq_rai_dt) + S_q_rai += S_qr + S_q_tot -= S_qr + S_q_liq -= S_qr + S_N_liq += -min(max(0, N_liq / dt), -tmp.dN_liq_dt) + S_N_rai += min(max(0, N_rai / dt), tmp.dN_rai_dt) + + # evaporation + tmp = CM2.evaporation(microphys_params, rf, q, q_rai, ρ, N_rai, T) + S_Nr = -min(max(0, N_rai / dt), -tmp[1]) + S_qr = -min(max(0, q_rai / dt), -tmp[2]) + S_q_rai += S_qr + S_q_tot -= S_qr + S_q_vap -= S_qr + S_N_rai += S_Nr + + vt = CM2.terminal_velocity(microphys_params, CMT.SB2006Type(), q_rai, ρ, N_rai) + term_vel_N_rai = vt[1] + term_vel_rai = vt[2] + + return (; S_q_tot, S_q_liq, S_q_rai, S_N_liq, S_N_rai, S_q_vap, term_vel_rai, term_vel_N_rai) +end """ Precompute the auxiliary values @@ -475,6 +541,35 @@ end @. aux.precip_velocities.term_vel_rai = CC.Geometry.WVector(If(tmp.term_vel_rai) * FT(-1)) @. aux.precip_velocities.term_vel_sno = CC.Geometry.WVector(If(tmp.term_vel_sno) * FT(-1)) end +@inline function precompute_aux_precip!(ps::Precipitation2M, dY, Y, aux, t) + FT = eltype(Y.ρq_rai) + + aux.precip_variables.q_rai = Y.ρq_rai ./ aux.moisture_variables.ρ + aux.precip_variables.N_liq = Y.N_liq + aux.precip_variables.N_rai = Y.N_rai + + tmp = @. precip_helper_sources_2M!( + aux.params, + aux.moisture_variables.q_tot, + aux.moisture_variables.q_liq, + aux.precip_variables.q_rai, + aux.precip_variables.N_liq, + aux.precip_variables.N_rai, + aux.moisture_variables.T, + aux.moisture_variables.ρ, + aux.TS.dt, + ps, + ) + aux.precip_sources.S_q_rai = tmp.S_q_rai + aux.precip_sources.S_q_tot = tmp.S_q_tot + aux.precip_sources.S_q_liq = tmp.S_q_liq + aux.precip_sources.S_N_liq = tmp.S_N_liq + aux.precip_sources.S_N_rai = tmp.S_N_rai + + If = CC.Operators.InterpolateC2F(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) + @. aux.precip_velocities.term_vel_N_rai = CC.Geometry.WVector(If(tmp.term_vel_N_rai) * FT(-1)) + @. aux.precip_velocities.term_vel_rai = CC.Geometry.WVector(If(tmp.term_vel_rai) * FT(-1)) +end """ Advection Equation: ∂ϕ/dt = -∂(vΦ) @@ -561,6 +656,31 @@ end return dY end +@inline function advection_tendency!(::Precipitation2M, dY, Y, aux, t) + FT = eltype(Y.ρq_tot) + + If = CC.Operators.InterpolateC2F() + ∂ = CC.Operators.DivergenceF2C(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) + + @. dY.N_rai += + -∂( + (aux.prescribed_velocity.ρw / If(aux.moisture_variables.ρ) + aux.precip_velocities.term_vel_N_rai) * + If(Y.N_rai), + ) + @. dY.ρq_rai += + -∂( + (aux.prescribed_velocity.ρw / If(aux.moisture_variables.ρ) + aux.precip_velocities.term_vel_rai) * + If(Y.ρq_rai), + ) + + fcc = CC.Operators.FluxCorrectionC2C(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) + @. dY.N_rai += + fcc((aux.prescribed_velocity.ρw / If(aux.moisture_variables.ρ) + aux.precip_velocities.term_vel_N_rai), Y.N_rai) + @. dY.ρq_rai += + fcc((aux.prescribed_velocity.ρw / If(aux.moisture_variables.ρ) + aux.precip_velocities.term_vel_rai), Y.ρq_rai) + + return dY +end """ Additional source terms @@ -596,3 +716,11 @@ end return dY end +@inline function sources_tendency!(::Precipitation2M, dY, Y, aux, t) + + @. dY.ρq_rai += aux.moisture_variables.ρ * aux.precip_sources.S_q_rai + @. dY.N_liq += aux.precip_sources.S_N_liq + @. dY.N_rai += aux.precip_sources.S_N_rai + + return dY +end diff --git a/test/Project.toml b/test/Project.toml index bd7acc0f..1298520f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -26,7 +26,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" CLIMAParameters = "0.6" ClimaCore = "0.10" ClimaCorePlots = "0.2" -CloudMicrophysics = "0.9" +CloudMicrophysics = "0.10" DiffEqBase = "6.75" NCDatasets = "0.12" Plots = "1.29" diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index ff83aa8e..b4732f76 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -69,6 +69,12 @@ elseif precipitation_choice == "Precipitation1M" else error("Invalid rain formation choice: $rain_formation_choice") end +elseif precipitation_choice == "Precipitation2M" + if rain_formation_choice == "SB2006" + precip = KID.Precipitation2M(CMT.SB2006Type()) + else + error("Invalid rain formation choice: $rain_formation_choice") + end else error("Invalid precipitation choice: $precipitation_choice") end @@ -83,7 +89,7 @@ face_coord = CC.Fields.coordinate_field(face_space) # Initialize the netcdf output Stats struct output_folder = string("Output_", moisture_choice, "_", prognostics_choice, "_", precipitation_choice) -if precipitation_choice == "Precipitation1M" +if precipitation_choice in ["Precipitation1M", "Precipitation2M"] output_folder = output_folder * "_" * rain_formation_choice end if opts["qtot_flux_correction"] diff --git a/test/experiments/KiD_driver/parse_commandline.jl b/test/experiments/KiD_driver/parse_commandline.jl index b66daf3c..9668f284 100644 --- a/test/experiments/KiD_driver/parse_commandline.jl +++ b/test/experiments/KiD_driver/parse_commandline.jl @@ -15,15 +15,15 @@ function parse_commandline() arg_type = String default = "RhoThetaQ" "--precipitation_choice" - help = "Precipitation model choice: NoPrecipitation, Precipitation0M, Precipitation1M" + help = "Precipitation model choice: NoPrecipitation, Precipitation0M, Precipitation1M, Precipitation2M" arg_type = String default = "Precipitation1M" "--rain_formation_scheme_choice" - help = "Rain formation scheme choice: CliMA_1M, KK2000, B1994, TC1980, LD2004" + help = "Rain formation scheme choice: CliMA_1M, KK2000, B1994, TC1980, LD2004 for Precipitation1M; and SB2006 for Precipitation2M" arg_type = String default = "CliMA_1M" "--prescribed_Nd" - help = "Prescribed number of cloud droplets (used in KK2000, B1994, TC1980 and LD2004 rain formation schemes" + help = "Prescribed number of cloud droplets (used in KK2000, B1994, TC1980, LD2004 and SB2006 rain formation schemes)" arg_type = Float64 default = Float64(1e8) "--plotting_flag" @@ -82,7 +82,7 @@ function parse_commandline() "--p0" help = "Pressure at the surface [pa]" arg_type = Real - default = 100700.0 + default = 100000.0 end return AP.parse_args(s) diff --git a/test/experiments/calibrations/config.jl b/test/experiments/calibrations/config.jl index 2401d367..33cb8001 100644 --- a/test/experiments/calibrations/config.jl +++ b/test/experiments/calibrations/config.jl @@ -100,7 +100,7 @@ function get_model_config(params_calib_names::Array{String}) config["model"] = "KiD" config["moisture_choice"] = "NonEquilibriumMoisture" config["precipitation_choice"] = "Precipitation1M" - # Define rain formation choice: "CliMA_1M", "KK2000", "B1994", "TC1980", "LD2004" + # Define rain formation choice: "CliMA_1M", "KK2000", "B1994", "TC1980", "LD2004", "SB2006" config["rain_formation_choice"] = "CliMA_1M" config["z_min"] = 0.0 config["z_max"] = 3000.0 diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index d9d9e697..bd74c6b6 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -127,9 +127,9 @@ function plot_timeheight(nc_data_file; output = "output") q_tot_plt = collect(ds.group["profiles"]["q_tot"]) q_liq_plt = collect(ds.group["profiles"]["q_liq"]) q_rai_plt = collect(ds.group["profiles"]["q_rai"]) - p1 = Plots.heatmap(t_plt, z_plt, q_tot_plt, title = "q_tot [g/kg]", xlabel = "time [s]", ylabel = "z [m]") - p2 = Plots.heatmap(t_plt, z_plt, q_liq_plt, title = "q_liq [g/kg]", xlabel = "time [s]", ylabel = "z [m]") - p3 = Plots.heatmap(t_plt, z_plt, q_rai_plt, title = "q_rai [g/kg]", xlabel = "time [s]", ylabel = "z [m]") + p1 = Plots.heatmap(t_plt, z_plt, q_tot_plt .* 1e3, title = "q_tot [g/kg]", xlabel = "time [s]", ylabel = "z [m]") + p2 = Plots.heatmap(t_plt, z_plt, q_liq_plt .* 1e3, title = "q_liq [g/kg]", xlabel = "time [s]", ylabel = "z [m]") + p3 = Plots.heatmap(t_plt, z_plt, q_rai_plt .* 1e3, title = "q_rai [g/kg]", xlabel = "time [s]", ylabel = "z [m]") p = Plots.plot( p1, p2, diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index e4bed311..58da4e79 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -169,6 +169,8 @@ end θ_dry = [440.0, 450], q_rai = [0.0, 0.0], q_sno = [0.0, 0.0], + N_liq = [0.0, 0.0], + N_rai = [0.0, 0.0], S_ql_moisture = [0.0, 0.0], S_qi_moisture = [0.0, 0.0], S_qt_precip = [0.0, 0.0], @@ -176,6 +178,8 @@ end S_qi_precip = [0.0, 0.0], S_qr_precip = [0.0, 0.0], S_qs_precip = [0.0, 0.0], + S_Nl_precip = [0.0, 0.0], + S_Nr_precip = [0.0, 0.0], ) space, face_space = KID.make_function_space(Float64, 0, 100, 10) @@ -252,6 +256,8 @@ end θ_dry = [440.0, 450], q_rai = [0.0, 0.0], q_sno = [0.0, 0.0], + N_liq = [0.0, 0.0], + N_rai = [0.0, 0.0], S_ql_moisture = [0.0, 0.0], S_qi_moisture = [0.0, 0.0], S_qt_precip = [0.0, 0.0], @@ -259,6 +265,8 @@ end S_qi_precip = [0.0, 0.0], S_qr_precip = [0.0, 0.0], S_qs_precip = [0.0, 0.0], + S_Nl_precip = [0.0, 0.0], + S_Nr_precip = [0.0, 0.0], ) space, face_space = KID.make_function_space(Float64, 0, 100, 10) aux = KID.initialise_aux(Float64, ip, params, 0.0, 0.0, face_space, KID.EquilibriumMoisture_ρθq()) From 82a33e895afba388d1f33ae538cce0e66ed31b3b Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 6 Apr 2023 14:48:16 -0700 Subject: [PATCH 02/16] add compat entries to make registrator happy --- Project.toml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7712422e..08d37e27 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Kinematic1D" uuid = "a8a3bb05-6cc5-4342-abf6-5cc636bd2b35" authors = ["Climate Modeling Alliance"] -version = "0.3.2" +version = "0.3.0" [deps] ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884" @@ -28,8 +28,13 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" ClimaCore = "0.10" ClimaCorePlots = "0.2" DiffEqBase = "6.75" +Distributions = "0.25" EnsembleKalmanProcesses = "0.13" +Interpolations = "0.14" +JLD2 = "0.4" NCDatasets = "0.12" +Optim = "1.7" +OrdinaryDiffEq = "6.29" Plots = "1.29" ProgressLogging = "0.1" TerminalLoggers = "0.1" From 4fa80c82b742e6d183315ea2df0adfb81f1dfd23 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 6 Apr 2023 14:51:45 -0700 Subject: [PATCH 03/16] fix number density initialization inconsistency and a minor bug --- src/driveKiD/initial_condition.jl | 3 ++- src/driveKiD/tendency.jl | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/driveKiD/initial_condition.jl b/src/driveKiD/initial_condition.jl index 8cfac58a..db506d76 100644 --- a/src/driveKiD/initial_condition.jl +++ b/src/driveKiD/initial_condition.jl @@ -111,6 +111,7 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F ρ::FT = ρ_profile(z) ρ_dry::FT = SDM_ρ_dry_of_ρ(ρ, q_vap) + ρ_SDP = 1.225 # assuming no cloud condensate in the initial profile θ_liq_ice::FT = θ_std # TODO - compute this based on TS @@ -132,7 +133,7 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F q_sno::FT = FT(0.0) ρq_rai::FT = q_rai * ρ ρq_sno::FT = q_sno * ρ - N_liq::FT = Parameters.prescribed_Nd(params) + N_liq::FT = Parameters.prescribed_Nd(params) * ρ_dry / ρ_SDP N_rai::FT = FT(0) S_ql_moisture::FT = FT(0.0) diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index 4c61ccfe..caaa62a4 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -346,7 +346,7 @@ end S_N_rai += min(max(0, N_rai / dt), tmp.dN_rai_dt) # evaporation - tmp = CM2.evaporation(microphys_params, rf, q, q_rai, ρ, N_rai, T) + tmp = CM2.rain_evaporation(microphys_params, rf, q, q_rai, ρ, N_rai, T) S_Nr = -min(max(0, N_rai / dt), -tmp[1]) S_qr = -min(max(0, q_rai / dt), -tmp[2]) S_q_rai += S_qr @@ -354,7 +354,7 @@ end S_q_vap -= S_qr S_N_rai += S_Nr - vt = CM2.terminal_velocity(microphys_params, CMT.SB2006Type(), q_rai, ρ, N_rai) + vt = CM2.rain_terminal_velocity(microphys_params, CMT.SB2006Type(), q_rai, ρ, N_rai) term_vel_N_rai = vt[1] term_vel_rai = vt[2] From be27addbc1bf2b6895e2f938cd7378d02a8d4954 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 6 Apr 2023 14:52:59 -0700 Subject: [PATCH 04/16] minor --- test/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Project.toml b/test/Project.toml index 1298520f..1e254a91 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,7 +23,7 @@ Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] -CLIMAParameters = "0.6" +CLIMAParameters = "0.7" ClimaCore = "0.10" ClimaCorePlots = "0.2" CloudMicrophysics = "0.10" From 60b036fba9c3d624689633baba8ac95b7e310fd0 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Tue, 11 Apr 2023 09:44:30 -0700 Subject: [PATCH 05/16] aerosol activation, ARG model, added --- src/Kinematic1D.jl | 2 + src/calibrateKiD/KiDUtils.jl | 5 +++ src/calibrateKiD/ReferenceModels.jl | 4 +- src/driveKiD/KiD_model.jl | 3 ++ src/driveKiD/Parameters.jl | 7 ++++ src/driveKiD/helper_functions.jl | 37 +++++++++++++++++++ src/driveKiD/initial_condition.jl | 8 +++- src/driveKiD/tendency.jl | 26 +++++++++++-- test/create_parameters.jl | 17 ++++++++- test/experiments/KiD_driver/KiD_driver.jl | 3 ++ .../KiD_driver/parse_commandline.jl | 12 ++++++ test/experiments/calibrations/config.jl | 3 ++ test/unit_tests/unit_test.jl | 1 + 13 files changed, 121 insertions(+), 7 deletions(-) diff --git a/src/Kinematic1D.jl b/src/Kinematic1D.jl index 5d73a2d4..55c5f310 100644 --- a/src/Kinematic1D.jl +++ b/src/Kinematic1D.jl @@ -21,6 +21,8 @@ const CMNe = CM.MicrophysicsNonEq const CM0 = CM.Microphysics0M const CM1 = CM.Microphysics1M const CM2 = CM.Microphysics2M +const CMAM = CM.AerosolModel +const CMAA = CM.AerosolActivation # modules used only for calibrateKiD using Distributions diff --git a/src/calibrateKiD/KiDUtils.jl b/src/calibrateKiD/KiDUtils.jl index 4ea3abd2..06a5d413 100644 --- a/src/calibrateKiD/KiDUtils.jl +++ b/src/calibrateKiD/KiDUtils.jl @@ -159,6 +159,8 @@ function get_variable_data_from_ODE(u, ρ::Vector{Float64}, var::String) output = parent(u.N_liq) elseif var == "Nr" output = parent(u.N_rai) + elseif var == "Na" + output = parent(u.N_aer) elseif var == "rho" output = ρ elseif var == "rain averaged terminal velocity" @@ -285,6 +287,9 @@ function create_parameter_set(FT, model_settings::Dict, params_cal::Dict) precip_sinks = 1, prescribed_Nd = model_settings["Nd"], qtot_flux_correction = Int(model_settings["qtot_flux_correction"]), + r_dry = model_settings["r_dry"], + std_dry = model_settings["std_dry"], + κ = model_settings["κ"], microphys_params, ) return param_set diff --git a/src/calibrateKiD/ReferenceModels.jl b/src/calibrateKiD/ReferenceModels.jl index ef1b002d..ce5b67d7 100644 --- a/src/calibrateKiD/ReferenceModels.jl +++ b/src/calibrateKiD/ReferenceModels.jl @@ -144,9 +144,11 @@ function get_single_obs_field( elseif var == "rtr" _data = _data_pysdm["qv"] .+ (_data_pysdm["qc"] .+ _data_pysdm["qr"]) .* 1e-3 elseif var == "Nl" - _data = _data_pysdm["na"] + _data_pysdm["nc"] + _data = _data_pysdm["nc"] elseif var == "Nr" _data = _data_pysdm["nr"] + elseif var == "Na" + _data = _data_pysdm["na"] else _data = _data_pysdm[var] end diff --git a/src/driveKiD/KiD_model.jl b/src/driveKiD/KiD_model.jl index 4b7a2a36..4304c245 100644 --- a/src/driveKiD/KiD_model.jl +++ b/src/driveKiD/KiD_model.jl @@ -85,6 +85,7 @@ function initialise_state(::EquilibriumMoisture, ::Precipitation2M, initial_prof ρq_rai = initial_profiles.ρq_rai, N_liq = initial_profiles.N_liq, N_rai = initial_profiles.N_rai, + N_aer = initial_profiles.N_aer, ) end function initialise_state(::NonEquilibriumMoisture, ::Precipitation2M, initial_profiles) @@ -95,6 +96,7 @@ function initialise_state(::NonEquilibriumMoisture, ::Precipitation2M, initial_p ρq_rai = initial_profiles.ρq_rai, N_liq = initial_profiles.N_liq, N_rai = initial_profiles.N_rai, + N_aer = initial_profiles.N_aer, ) end @@ -143,6 +145,7 @@ function initialise_aux(FT, ip, params, TS, Stats, face_space, moisture) N_liq = ip.N_liq, N_rai = ip.N_rai, ), + aerosol_variables = CC.Fields.FieldVector(; N_aer = ip.N_aer, N_aer_0 = ip.N_aer_0, S_N_aer = ip.S_Na), moisture_sources = CC.Fields.FieldVector(; S_q_liq = ip.S_ql_moisture, S_q_ice = ip.S_qi_moisture), precip_sources = CC.Fields.FieldVector(; S_q_tot = ip.S_qt_precip, diff --git a/src/driveKiD/Parameters.jl b/src/driveKiD/Parameters.jl index 730b68ea..6026c7fb 100644 --- a/src/driveKiD/Parameters.jl +++ b/src/driveKiD/Parameters.jl @@ -22,6 +22,9 @@ Base.@kwdef struct KinematicParameters{FT, MP} <: AKP precip_sinks::Int prescribed_Nd::FT qtot_flux_correction::Int + r_dry::FT + std_dry::FT + κ::FT microphys_params::MP end thermodynamics_params(ps::AKP) = CM.Parameters.thermodynamics_params(ps.microphys_params) @@ -33,6 +36,10 @@ t1(ps::AKP) = ps.t1 precip_sources(ps::AKP) = ps.precip_sources precip_sinks(ps::AKP) = ps.precip_sinks prescribed_Nd(ps::AKP) = ps.prescribed_Nd +#aerosol parameters for 2M scheme +r_dry(ps::AKP) = ps.r_dry +std_dry(ps::AKP) = ps.std_dry +κ(ps::AKP) = ps.κ # Forward parameters to Thermodynamics const TDPS = TD.Parameters.ThermodynamicsParameters diff --git a/src/driveKiD/helper_functions.jl b/src/driveKiD/helper_functions.jl index ec166824..7ee9533e 100644 --- a/src/driveKiD/helper_functions.jl +++ b/src/driveKiD/helper_functions.jl @@ -4,3 +4,40 @@ @inline function ρw_helper(t, w1, t1) return t < t1 ? w1 * sin(pi * t / t1) : 0.0 end + +""" + Returns the number of new activated aerosol particles and updates aerosol number density +""" +@inline function aerosol_activation_helper!(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw, dt) + + microphys_params = KP.microphysics_params(params) + thermo_params = CM.Parameters.thermodynamics_params(microphys_params) + + FT = eltype(q_tot) + S_Nl::FT = FT(0) + S_Na::FT = FT(0) + + q = TD.PhasePartition(q_tot, q_liq, FT(0)) + S::FT = TD.supersaturation(thermo_params, q, ρ, T, TD.Liquid()) + + if (S < FT(0)) + return (; S_Nl, S_Na) + end + + r_dry = KP.r_dry(params) + std_dry = KP.std_dry(params) + κ = KP.κ(params) + w = ρw / ρ + + aerosol_distribution = CMAM.AerosolDistribution((CMAM.Mode_κ(r_dry, std_dry, N_aer_0, FT(1), FT(1), FT(0), κ, 1),)) + N_act = CMAA.N_activated_per_mode(microphys_params, aerosol_distribution, T, p, w, q)[1] + + if isnan(N_act) + return (; S_Nl, S_Na) + end + + S_Nl = max(0, N_act - (N_aer_0 - N_aer)) + S_Na = -S_Nl + + return (; S_Nl, S_Na) +end \ No newline at end of file diff --git a/src/driveKiD/initial_condition.jl b/src/driveKiD/initial_condition.jl index db506d76..9d7fa969 100644 --- a/src/driveKiD/initial_condition.jl +++ b/src/driveKiD/initial_condition.jl @@ -133,8 +133,10 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F q_sno::FT = FT(0.0) ρq_rai::FT = q_rai * ρ ρq_sno::FT = q_sno * ρ - N_liq::FT = Parameters.prescribed_Nd(params) * ρ_dry / ρ_SDP + N_liq::FT = FT(0) N_rai::FT = FT(0) + N_aer_0::FT = Parameters.prescribed_Nd(params) * ρ_dry / ρ_SDP + N_aer::FT = N_aer_0 S_ql_moisture::FT = FT(0.0) S_qi_moisture::FT = FT(0.0) @@ -146,6 +148,7 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F S_qs_precip::FT = FT(0.0) S_Nl_precip::FT = FT(0) S_Nr_precip::FT = FT(0) + S_Na::FT = FT(0) return (; ρ, @@ -166,6 +169,8 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F q_sno, N_liq, N_rai, + N_aer, + N_aer_0, S_ql_moisture, S_qi_moisture, S_qt_precip, @@ -175,5 +180,6 @@ function init_1d_column(::Type{FT}, params, ρ_profile, z; dry = false) where {F S_qs_precip, S_Nl_precip, S_Nr_precip, + S_Na, ) end diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index caaa62a4..8424ac88 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -33,6 +33,7 @@ end @. dY.ρq_rai = FT(0) @. dY.N_liq = FT(0) @. dY.N_rai = FT(0) + @. dY.N_aer = FT(0) end """ @@ -323,8 +324,9 @@ end S_q_rai += S_qr S_q_tot -= S_qr S_q_liq -= S_qr - S_N_liq += -min(max(0, N_liq / dt), -tmp.dN_liq_dt) - S_N_rai += min(max(0, N_liq / dt), tmp.dN_rai_dt) + S_Nr = min(max(0, N_liq / 2 / dt), tmp.dN_rai_dt) + S_N_rai += S_Nr + S_N_liq -= 2*S_Nr # self_collection tmp_l = CM2.liquid_self_collection(microphys_params, rf, q.liq, ρ, -S_qr) @@ -343,7 +345,7 @@ end S_q_tot -= S_qr S_q_liq -= S_qr S_N_liq += -min(max(0, N_liq / dt), -tmp.dN_liq_dt) - S_N_rai += min(max(0, N_rai / dt), tmp.dN_rai_dt) + S_N_rai += FT(0) # evaporation tmp = CM2.rain_evaporation(microphys_params, rf, q, q_rai, ρ, N_rai, T) @@ -544,10 +546,25 @@ end @inline function precompute_aux_precip!(ps::Precipitation2M, dY, Y, aux, t) FT = eltype(Y.ρq_rai) + aux.aerosol_variables.N_aer = Y.N_aer + tmp = @. aerosol_activation_helper!( + aux.params, + aux.moisture_variables.q_tot, + aux.moisture_variables.q_liq, + aux.aerosol_variables.N_aer, + aux.aerosol_variables.N_aer_0, + aux.moisture_variables.T, + aux.moisture_variables.p, + aux.moisture_variables.ρ, + ClimaCore.Operators.InterpolateF2C().(aux.prescribed_velocity.ρw.components.data.:1), + aux.TS.dt, + ) + @. Y.N_liq += tmp.S_Nl + @. Y.N_aer += tmp.S_Na + aux.precip_variables.q_rai = Y.ρq_rai ./ aux.moisture_variables.ρ aux.precip_variables.N_liq = Y.N_liq aux.precip_variables.N_rai = Y.N_rai - tmp = @. precip_helper_sources_2M!( aux.params, aux.moisture_variables.q_tot, @@ -721,6 +738,7 @@ end @. dY.ρq_rai += aux.moisture_variables.ρ * aux.precip_sources.S_q_rai @. dY.N_liq += aux.precip_sources.S_N_liq @. dY.N_rai += aux.precip_sources.S_N_rai + @. dY.N_aer += aux.aerosol_variables.S_N_aer return dY end diff --git a/test/create_parameters.jl b/test/create_parameters.jl index 85b55b2e..b49fed9c 100644 --- a/test/create_parameters.jl +++ b/test/create_parameters.jl @@ -15,6 +15,9 @@ function create_parameter_set( precip_sinks = 1, qtot_flux_correction = 0, prescribed_Nd = 100 * 1e6, + r_dry = 0.04 * 1e-6, + std_dry = 1.4, + κ = 0.9, ) FT = CP.float_type(toml_dict) override_file = joinpath(out_dir, "override_dict.toml") @@ -79,6 +82,18 @@ function create_parameter_set( println(io, "alias = \"prescribed_Nd\"") println(io, "value = " * string(prescribed_Nd)) println(io, "type = \"float\"") + println(io, "[r_dry]") + println(io, "alias = \"r_dry\"") + println(io, "value = " * string(r_dry)) + println(io, "type = \"float\"") + println(io, "[std_dry]") + println(io, "alias = \"std_dry\"") + println(io, "value = " * string(std_dry)) + println(io, "type = \"float\"") + println(io, "[κ]") + println(io, "alias = \"κ\"") + println(io, "value = " * string(κ)) + println(io, "type = \"float\"") end toml_dict = CP.create_toml_dict(FT; override_file, dict_type="alias") isfile(override_file) && rm(override_file; force=true) @@ -97,7 +112,7 @@ function create_parameter_set( ) MP = typeof(microphys_params) - aliases = ["w1", "t1", "p0", "precip_sources", "precip_sinks", "qtot_flux_correction", "prescribed_Nd"] + aliases = ["w1", "t1", "p0", "precip_sources", "precip_sinks", "qtot_flux_correction", "prescribed_Nd", "r_dry", "std_dry", "κ"] pairs = CP.get_parameter_values!(toml_dict, aliases, "Kinematic1D") param_set = KID.Parameters.KinematicParameters{FTD, MP}(; pairs..., microphys_params) diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index b4732f76..06ed4ff4 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -113,6 +113,9 @@ params = create_parameter_set( Int(opts["precip_sinks"]), Int(opts["qtot_flux_correction"]), FT(opts["prescribed_Nd"]), + FT(opts["r_dry"]), + FT(opts["std_dry"]), + FT(opts["kappa"]), ) # Solve the initial value problem for density profile diff --git a/test/experiments/KiD_driver/parse_commandline.jl b/test/experiments/KiD_driver/parse_commandline.jl index 9668f284..5acdcfcc 100644 --- a/test/experiments/KiD_driver/parse_commandline.jl +++ b/test/experiments/KiD_driver/parse_commandline.jl @@ -83,6 +83,18 @@ function parse_commandline() help = "Pressure at the surface [pa]" arg_type = Real default = 100000.0 + "--r_dry" + help = "aerosol distribution mean radius for aerosol activation calculations in 2M schemes [m]" + arg_type = Real + default = 0.04 * 1e-6 + "--std_dry" + help = "aerosol distribution standard deviation for aerosol activation calucaulations in 2M schemes" + arg_type = Real + default = 1.4 + "--kappa" + help = "hygroscopicity of aerosols for aerosol activation calucaulations in 2M schemes" + arg_type = Real + default = 0.9 end return AP.parse_args(s) diff --git a/test/experiments/calibrations/config.jl b/test/experiments/calibrations/config.jl index 33cb8001..aa69a895 100644 --- a/test/experiments/calibrations/config.jl +++ b/test/experiments/calibrations/config.jl @@ -115,6 +115,9 @@ function get_model_config(params_calib_names::Array{String}) config["p0"] = 100000.0 config["Nd"] = 100 * 1e6 config["qtot_flux_correction"] = false + config["r_dry"] = 0.04 * 1e-6 + config["std_dry"] = 1.4 + config["κ"] = 0.9 config["filter"] = KD.make_filter_props( config["n_elem"], config["t_calib"]; diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 58da4e79..959a11c1 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -38,6 +38,7 @@ thermo_params = KP.thermodynamics_params(params) @test KID.NoPrecipitation <: KID.AbstractPrecipitationStyle @test KID.Precipitation0M <: KID.AbstractPrecipitationStyle @test KID.Precipitation1M <: KID.AbstractPrecipitationStyle + @test KID.Precipitation2M <: KID.AbstractPrecipitationStyle end From 678a49d1088b4194a382b26f427a456f7afe7729 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Tue, 11 Apr 2023 16:08:57 -0700 Subject: [PATCH 06/16] style format --- Project.toml | 1 + src/driveKiD/helper_functions.jl | 4 ++-- src/driveKiD/tendency.jl | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 08d37e27..60368e38 100644 --- a/Project.toml +++ b/Project.toml @@ -27,6 +27,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] ClimaCore = "0.10" ClimaCorePlots = "0.2" +CloudMicrophysics = "0.10" DiffEqBase = "6.75" Distributions = "0.25" EnsembleKalmanProcesses = "0.13" diff --git a/src/driveKiD/helper_functions.jl b/src/driveKiD/helper_functions.jl index 7ee9533e..16d68359 100644 --- a/src/driveKiD/helper_functions.jl +++ b/src/driveKiD/helper_functions.jl @@ -29,7 +29,7 @@ end κ = KP.κ(params) w = ρw / ρ - aerosol_distribution = CMAM.AerosolDistribution((CMAM.Mode_κ(r_dry, std_dry, N_aer_0, FT(1), FT(1), FT(0), κ, 1),)) + aerosol_distribution = CMAM.AerosolDistribution((CMAM.Mode_κ(r_dry, std_dry, N_aer_0, FT(1), FT(1), FT(0), κ, 1),)) N_act = CMAA.N_activated_per_mode(microphys_params, aerosol_distribution, T, p, w, q)[1] if isnan(N_act) @@ -40,4 +40,4 @@ end S_Na = -S_Nl return (; S_Nl, S_Na) -end \ No newline at end of file +end diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index 8424ac88..0d90747a 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -326,7 +326,7 @@ end S_q_liq -= S_qr S_Nr = min(max(0, N_liq / 2 / dt), tmp.dN_rai_dt) S_N_rai += S_Nr - S_N_liq -= 2*S_Nr + S_N_liq -= 2 * S_Nr # self_collection tmp_l = CM2.liquid_self_collection(microphys_params, rf, q.liq, ρ, -S_qr) From 624fca40651b5edbca4cd495354fc923dd862dd4 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 12:39:34 -0700 Subject: [PATCH 07/16] minor --- test/create_parameters.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/create_parameters.jl b/test/create_parameters.jl index b49fed9c..95138ae4 100644 --- a/test/create_parameters.jl +++ b/test/create_parameters.jl @@ -90,8 +90,8 @@ function create_parameter_set( println(io, "alias = \"std_dry\"") println(io, "value = " * string(std_dry)) println(io, "type = \"float\"") - println(io, "[κ]") - println(io, "alias = \"κ\"") + println(io, "[kappa]") + println(io, "alias = \"kappa\"") println(io, "value = " * string(κ)) println(io, "type = \"float\"") end From a4210b7f8aa2ca8334058a3874a1a10e19a24402 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 13:17:41 -0700 Subject: [PATCH 08/16] minor --- test/create_parameters.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/create_parameters.jl b/test/create_parameters.jl index 95138ae4..2f549244 100644 --- a/test/create_parameters.jl +++ b/test/create_parameters.jl @@ -91,7 +91,7 @@ function create_parameter_set( println(io, "value = " * string(std_dry)) println(io, "type = \"float\"") println(io, "[kappa]") - println(io, "alias = \"kappa\"") + println(io, "alias = \"κ\"") println(io, "value = " * string(κ)) println(io, "type = \"float\"") end From 97967f3bd1e289361f8b44078a6830922d512f5a Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 13:53:44 -0700 Subject: [PATCH 09/16] minor --- test/unit_tests/calibration_pipeline/config.jl | 3 +++ test/unit_tests/unit_test.jl | 6 ++++++ 2 files changed, 9 insertions(+) diff --git a/test/unit_tests/calibration_pipeline/config.jl b/test/unit_tests/calibration_pipeline/config.jl index fa71bcb5..1d2fde6e 100644 --- a/test/unit_tests/calibration_pipeline/config.jl +++ b/test/unit_tests/calibration_pipeline/config.jl @@ -82,6 +82,9 @@ function get_model_config(params_calib_names::Array{String}) config["p0"] = 99000.0 config["Nd"] = 50 * 1e6 config["qtot_flux_correction"] = false + config["r_dry"] = 0.04 * 1e-6 + config["std_dry"] = 1.4 + config["κ"] = 0.9 config["filter"] = KD.make_filter_props(config["n_elem"], config["t_calib"]; apply = false) fixed_parameters = create_fixed_parameter_set(params_calib_names) config["thermo_params"] = fixed_parameters.thermo_params diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 959a11c1..123ab9a0 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -172,6 +172,8 @@ end q_sno = [0.0, 0.0], N_liq = [0.0, 0.0], N_rai = [0.0, 0.0], + N_aer_0 = [0.0, 0.0], + N_aer = [0.0, 0.0], S_ql_moisture = [0.0, 0.0], S_qi_moisture = [0.0, 0.0], S_qt_precip = [0.0, 0.0], @@ -181,6 +183,7 @@ end S_qs_precip = [0.0, 0.0], S_Nl_precip = [0.0, 0.0], S_Nr_precip = [0.0, 0.0], + S_Na = [0.0, 0.0], ) space, face_space = KID.make_function_space(Float64, 0, 100, 10) @@ -259,6 +262,8 @@ end q_sno = [0.0, 0.0], N_liq = [0.0, 0.0], N_rai = [0.0, 0.0], + N_aer_0 = [0.0, 0.0], + N_aer = [0.0, 0.0], S_ql_moisture = [0.0, 0.0], S_qi_moisture = [0.0, 0.0], S_qt_precip = [0.0, 0.0], @@ -268,6 +273,7 @@ end S_qs_precip = [0.0, 0.0], S_Nl_precip = [0.0, 0.0], S_Nr_precip = [0.0, 0.0], + S_Na = [0.0, 0.0], ) space, face_space = KID.make_function_space(Float64, 0, 100, 10) aux = KID.initialise_aux(Float64, ip, params, 0.0, 0.0, face_space, KID.EquilibriumMoisture_ρθq()) From 05d225efcbd280e4c1595f1de0a401e889c09bff Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 16:05:43 -0700 Subject: [PATCH 10/16] tests --- test/unit_tests/unit_test.jl | 71 +++++++++++++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 2 deletions(-) diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 123ab9a0..078f1f82 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -10,6 +10,7 @@ import CLIMAParameters import ClimaCore import Thermodynamics import Kinematic1D +import CloudMicrophysics const kid_dir = pkgdir(Kinematic1D) include(joinpath(kid_dir, "test", "create_parameters.jl")) @@ -20,6 +21,7 @@ const CC = ClimaCore const TD = Thermodynamics const KID = Kinematic1D const KP = KID.Parameters +const CMT = CloudMicrophysics.CommonTypes # Create all the params boxes and overwrite the defaults to match PySDM toml_dict = CP.create_toml_dict(Float64; dict_type = "alias") @@ -68,7 +70,7 @@ end @test_throws Exception KID.initialise_state(KID.AbstractMoistureStyle(), KID.AbstractPrecipitationStyle(), 0) - initial_profiles = (; ρq_tot = 0, ρq_liq = 0, ρq_ice = 0, ρq_rai = 0, ρq_sno = 0) + initial_profiles = (; ρq_tot = 0, ρq_liq = 0, ρq_ice = 0, ρq_rai = 0, ρq_sno = 0, N_liq = 0, N_rai = 0, N_aer = 0) state = KID.initialise_state(KID.EquilibriumMoisture_ρθq(), KID.NoPrecipitation(), initial_profiles) @test state isa CC.Fields.FieldVector @@ -88,6 +90,14 @@ end @test LA.norm(state.ρq_rai) == 0 @test LA.norm(state.ρq_sno) == 0 + state = KID.initialise_state(KID.EquilibriumMoisture_ρθq(), KID.Precipitation2M(CMT.SB2006Type()), initial_profiles) + @test state isa CC.Fields.FieldVector + @test LA.norm(state.ρq_tot) == 0 + @test LA.norm(state.ρq_rai) == 0 + @test LA.norm(state.N_liq) == 0 + @test LA.norm(state.N_rai) == 0 + @test LA.norm(state.N_aer) == 0 + state = KID.initialise_state(KID.EquilibriumMoisture_ρdTq(), KID.NoPrecipitation(), initial_profiles) @test state isa CC.Fields.FieldVector @test LA.norm(state.ρq_tot) == 0 @@ -106,6 +116,15 @@ end @test LA.norm(state.ρq_rai) == 0 @test LA.norm(state.ρq_sno) == 0 + state = + KID.initialise_state(KID.EquilibriumMoisture_ρdTq(), KID.Precipitation2M(CMT.SB2006Type()), initial_profiles) + @test state isa CC.Fields.FieldVector + @test LA.norm(state.ρq_tot) == 0 + @test LA.norm(state.ρq_rai) == 0 + @test LA.norm(state.N_liq) == 0 + @test LA.norm(state.N_rai) == 0 + @test LA.norm(state.N_aer) == 0 + state = KID.initialise_state(KID.NonEquilibriumMoisture_ρθq(), KID.NoPrecipitation(), initial_profiles) @test state isa CC.Fields.FieldVector @test LA.norm(state.ρq_tot) == 0 @@ -130,6 +149,17 @@ end @test LA.norm(state.ρq_rai) == 0 @test LA.norm(state.ρq_sno) == 0 + state = + KID.initialise_state(KID.NonEquilibriumMoisture_ρθq(), KID.Precipitation2M(CMT.SB2006Type()), initial_profiles) + @test state isa CC.Fields.FieldVector + @test LA.norm(state.ρq_tot) == 0 + @test LA.norm(state.ρq_liq) == 0 + @test LA.norm(state.ρq_ice) == 0 + @test LA.norm(state.ρq_rai) == 0 + @test LA.norm(state.N_liq) == 0 + @test LA.norm(state.N_rai) == 0 + @test LA.norm(state.N_aer) == 0 + state = KID.initialise_state(KID.NonEquilibriumMoisture_ρdTq(), KID.NoPrecipitation(), initial_profiles) @test state isa CC.Fields.FieldVector @test LA.norm(state.ρq_tot) == 0 @@ -154,6 +184,17 @@ end @test LA.norm(state.ρq_rai) == 0 @test LA.norm(state.ρq_sno) == 0 + state = + KID.initialise_state(KID.NonEquilibriumMoisture_ρdTq(), KID.Precipitation2M(CMT.SB2006Type()), initial_profiles) + @test state isa CC.Fields.FieldVector + @test LA.norm(state.ρq_tot) == 0 + @test LA.norm(state.ρq_liq) == 0 + @test LA.norm(state.ρq_ice) == 0 + @test LA.norm(state.ρq_rai) == 0 + @test LA.norm(state.N_liq) == 0 + @test LA.norm(state.N_rai) == 0 + @test LA.norm(state.N_aer) == 0 + end @testset "Initialise aux" begin @@ -196,6 +237,7 @@ end @test aux.moisture_variables isa CC.Fields.FieldVector @test aux.precip_variables isa CC.Fields.FieldVector @test aux.moisture_sources isa CC.Fields.FieldVector + @test aux.aerosol_variables isa CC.Fields.FieldVector @test aux.precip_sources isa CC.Fields.FieldVector @test aux.precip_velocities isa CC.Fields.FieldVector @test aux.prescribed_velocity isa CC.Fields.FieldVector @@ -208,6 +250,7 @@ end @test aux.moisture_variables isa CC.Fields.FieldVector @test aux.precip_variables isa CC.Fields.FieldVector @test aux.moisture_sources isa CC.Fields.FieldVector + @test aux.aerosol_variables isa CC.Fields.FieldVector @test aux.precip_sources isa CC.Fields.FieldVector @test aux.precip_velocities isa CC.Fields.FieldVector @test aux.prescribed_velocity isa CC.Fields.FieldVector @@ -220,6 +263,7 @@ end @test aux.moisture_variables isa CC.Fields.FieldVector @test aux.precip_variables isa CC.Fields.FieldVector @test aux.moisture_sources isa CC.Fields.FieldVector + @test aux.aerosol_variables isa CC.Fields.FieldVector @test aux.precip_sources isa CC.Fields.FieldVector @test aux.precip_velocities isa CC.Fields.FieldVector @test aux.prescribed_velocity isa CC.Fields.FieldVector @@ -232,6 +276,7 @@ end @test aux.moisture_variables isa CC.Fields.FieldVector @test aux.precip_variables isa CC.Fields.FieldVector @test aux.moisture_sources isa CC.Fields.FieldVector + @test aux.aerosol_variables isa CC.Fields.FieldVector @test aux.precip_sources isa CC.Fields.FieldVector @test aux.precip_velocities isa CC.Fields.FieldVector @test aux.prescribed_velocity isa CC.Fields.FieldVector @@ -300,6 +345,19 @@ end KID.zero_tendencies!(KID.Precipitation1M(KID.OneMomentRainFormation()), dY, Y, aux, 1.0) @test LA.norm(dY) == 0 + dY = (; + ρq_tot = [10.0, 13.0], + ρq_liq = [5.0, 10.0], + ρq_ice = [44.0, -42.0], + ρq_rai = [1.0, 2.0], + N_liq = [1e8, 1e7], + N_rai = [1e4, 1e4], + N_aer = [1e9, 1e8], + ) + KID.zero_tendencies!(KID.NonEquilibriumMoisture_ρθq(), dY, Y, aux, 1.0) + KID.zero_tendencies!(KID.Precipitation2M(CMT.SB2006Type()), dY, Y, aux, 1.0) + @test LA.norm(dY) == 0 + end @testset "Tendency helper functions" begin @@ -317,6 +375,8 @@ end q_ice = ρq_ice / ρ q_rai = ρq_rai / ρ q_sno = ρq_sno / ρ + N_liq = 1e8 + N_rai = 1e4 T = 280.0 q = TD.PhasePartition(q_tot, q_liq, q_ice) @@ -370,13 +430,20 @@ end @test tmp.S_q_tot == tmp.S_q_liq + tmp.S_q_ice @test tmp.S_q_tot == -(tmp.S_q_rai + tmp.S_q_sno) - ps = KID.Precipitation1M(KID.OneMomentRainFormation()) tmp = KID.precip_helper_sources_1M!(params, ts, q_tot, q_liq, q_ice, q_rai, q_sno, T, ρ, dt, ps) @test !isnan(tmp.S_q_tot .+ tmp.S_q_liq .+ tmp.S_q_ice .+ tmp.S_q_rai .+ tmp.S_q_sno) @test tmp.S_q_tot ≈ tmp.S_q_liq + tmp.S_q_ice + tmp.S_q_vap @test tmp.S_q_tot ≈ -(tmp.S_q_rai + tmp.S_q_sno) + ps = KID.Precipitation2M(CMT.SB2006Type()) + tmp = KID.precip_helper_sources_2M!(params, q_tot, q_liq, q_rai, N_liq, N_rai, T, ρ, dt, ps) + @test !isnan(tmp.S_q_tot .+ tmp.S_q_liq .+ tmp.S_q_rai) + @test !isnan(tmp.S_N_liq .+ tmp.S_N_rai) + @test !isnan(tmp.term_vel_rai .+ tmp.term_vel_N_rai) + @test tmp.S_q_tot ≈ tmp.S_q_liq + tmp.S_q_vap + @test tmp.S_q_tot ≈ -tmp.S_q_rai + end @testset "precompute_aux, advection_tendency and sources_tendency" begin From 1dbe6c8409e5221c1e4216ebf3b2ce403c3efd19 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 16:08:15 -0700 Subject: [PATCH 11/16] minor --- src/driveKiD/helper_functions.jl | 2 +- src/driveKiD/tendency.jl | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/driveKiD/helper_functions.jl b/src/driveKiD/helper_functions.jl index 16d68359..981ba38b 100644 --- a/src/driveKiD/helper_functions.jl +++ b/src/driveKiD/helper_functions.jl @@ -8,7 +8,7 @@ end """ Returns the number of new activated aerosol particles and updates aerosol number density """ -@inline function aerosol_activation_helper!(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw, dt) +@inline function aerosol_activation_helper!(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw) microphys_params = KP.microphysics_params(params) thermo_params = CM.Parameters.thermodynamics_params(microphys_params) diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index 0d90747a..c6ea6210 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -557,7 +557,6 @@ end aux.moisture_variables.p, aux.moisture_variables.ρ, ClimaCore.Operators.InterpolateF2C().(aux.prescribed_velocity.ρw.components.data.:1), - aux.TS.dt, ) @. Y.N_liq += tmp.S_Nl @. Y.N_aer += tmp.S_Na From 5930ee0ceba3bd74af1a674949f62b83cf3bcabb Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 16:20:12 -0700 Subject: [PATCH 12/16] tests --- src/driveKiD/helper_functions.jl | 2 +- src/driveKiD/tendency.jl | 2 +- test/unit_tests/unit_test.jl | 20 ++++++++++++++++++++ 3 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/driveKiD/helper_functions.jl b/src/driveKiD/helper_functions.jl index 981ba38b..4ded338d 100644 --- a/src/driveKiD/helper_functions.jl +++ b/src/driveKiD/helper_functions.jl @@ -8,7 +8,7 @@ end """ Returns the number of new activated aerosol particles and updates aerosol number density """ -@inline function aerosol_activation_helper!(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw) +@inline function aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw) microphys_params = KP.microphysics_params(params) thermo_params = CM.Parameters.thermodynamics_params(microphys_params) diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index c6ea6210..f6019b86 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -547,7 +547,7 @@ end FT = eltype(Y.ρq_rai) aux.aerosol_variables.N_aer = Y.N_aer - tmp = @. aerosol_activation_helper!( + tmp = @. aerosol_activation_helper( aux.params, aux.moisture_variables.q_tot, aux.moisture_variables.q_liq, diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 078f1f82..1ab2465e 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -508,5 +508,25 @@ end end +@testset "aerosol activation for 2M schemes" begin + #setup + q_tot = 1e-2 + q_liq = 1e-3 + N_aer = 5e7 + N_aer_0 = 1e8 + T = 280.0 + p = 1e5 + ρ = 1.0 + ρw = 2.0 + + #action + tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw) + + #test + @test tmp isa NamedTuple + @test tmp.S_Nl ≈ -tmp.S_Na + @test 0 < tmp.S_Nl < N_aer_0 +end + # calibration pipeline unit tests include("./calibration_pipeline/run_calibration_pipeline_unit_tests.jl") From fdaa387d3c84f21e33baebbcc5330b609f52c148 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 13 Apr 2023 18:09:22 -0700 Subject: [PATCH 13/16] tests --- test/unit_tests/unit_test.jl | 101 +++++++++++++++++++++++++++++++---- 1 file changed, 90 insertions(+), 11 deletions(-) diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index 1ab2465e..db52a909 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -446,20 +446,22 @@ end end -@testset "precompute_aux, advection_tendency and sources_tendency" begin +@testset "advection_tendency and sources_tendency" begin space, face_space = KID.make_function_space(Float64, 0, 100, 5) coord = CC.Fields.coordinate_field(space) ρ_profile = KID.ρ_ivp(Float64, params) init = map(coord -> KID.init_1d_column(Float64, params, ρ_profile, coord.z), coord) + t = 13.0 + # eq aux = KID.initialise_aux(Float64, init, params, 0.0, 0.0, face_space, KID.EquilibriumMoisture_ρθq()) + ρw = 0.0 + @. aux.prescribed_velocity.ρw = CC.Geometry.WVector.(ρw) + aux.prescribed_velocity.ρw0 = ρw + Y = KID.initialise_state(KID.EquilibriumMoisture_ρθq(), KID.NoPrecipitation(), init) dY = Y / 10 - t = 13.0 - - @test_throws Exception KID.precompute_aux_thermo!(KID.AbstractMoistureStyle(), dY, Y, aux, t) - @test_throws Exception KID.precompute_aux_precip!(KID.AbstractPrecipitationStyle(), dY, Y, aux, t) @test_throws Exception advection_tendency!(KID.AbstractMoistureStyle(), dY, Y, aux, t) @test_throws Exception advection_tendency!(KID.AbstractPrecipitationStyle(), dY, Y, aux, t) @@ -467,8 +469,50 @@ end @test_throws Exception KID.sources_tendency!(KID.AbstractMoistureStyle(), dY, Y, aux, t) @test_throws Exception KID.sources_tendency!(KID.AbstractPrecipitationStyle(), dY, Y, aux, t) - KID.advection_tendency!(KID.EquilibriumMoisture_ρθq(), dY, Y, aux, 1.0) - @test dY == Y / 10 + KID.advection_tendency!(KID.EquilibriumMoisture_ρθq(), dY, Y, aux, t) + @test dY ≈ Y / 10 atol = eps(Float64) * 10 + + Y = KID.initialise_state(KID.EquilibriumMoisture_ρθq(), KID.Precipitation2M(CMT.SB2006Type()), init) + dY = Y / 10 + KID.advection_tendency!(KID.Precipitation2M(CMT.SB2006Type()), dY, Y, aux, t) + @test dY ≈ Y / 10 atol = eps(Float64) * 10 + + # Non-eq + aux = KID.initialise_aux(Float64, init, params, 0.0, 0.0, face_space, KID.NonEquilibriumMoisture_ρθq()) + ρw = 0.0 + @. aux.prescribed_velocity.ρw = CC.Geometry.WVector.(ρw) + aux.prescribed_velocity.ρw0 = ρw + + Y = KID.initialise_state(KID.NonEquilibriumMoisture_ρθq(), KID.NoPrecipitation(), init) + dY = Y / 10 + + KID.advection_tendency!(KID.NonEquilibriumMoisture_ρθq(), dY, Y, aux, t) + @test dY ≈ Y / 10 atol = eps(Float64) * 10 + + + Y = KID.initialise_state(KID.NonEquilibriumMoisture_ρθq(), KID.Precipitation2M(CMT.SB2006Type()), init) + dY = Y / 10 + KID.advection_tendency!(KID.Precipitation2M(CMT.SB2006Type()), dY, Y, aux, t) + @test dY ≈ Y / 10 atol = eps(Float64) * 10 + +end + +@testset "precompute aux" begin + + space, face_space = KID.make_function_space(Float64, 0, 100, 5) + coord = CC.Fields.coordinate_field(space) + ρ_profile = KID.ρ_ivp(Float64, params) + init = map(coord -> KID.init_1d_column(Float64, params, ρ_profile, coord.z), coord) + t = 13.0 + + # eq + aux = KID.initialise_aux(Float64, init, params, 0.0, 0.0, face_space, KID.EquilibriumMoisture_ρθq()) + + Y = KID.initialise_state(KID.EquilibriumMoisture_ρθq(), KID.NoPrecipitation(), init) + dY = Y / 10 + + @test_throws Exception KID.precompute_aux_thermo!(KID.AbstractMoistureStyle(), dY, Y, aux, t) + @test_throws Exception KID.precompute_aux_precip!(KID.AbstractPrecipitationStyle(), dY, Y, aux, t) KID.precompute_aux_thermo!(KID.EquilibriumMoisture_ρθq(), dY, Y, aux, t) @test aux.moisture_variables.ρ_dry == aux.moisture_variables.ρ .- Y.ρq_tot @@ -484,13 +528,10 @@ end @test aux.moisture_variables.θ_dry == TD.dry_pottemp.(thermo_params, aux.moisture_variables.T, aux.moisture_variables.ρ_dry) + # Non-eq aux = KID.initialise_aux(Float64, init, params, 0.0, 0.0, face_space, KID.NonEquilibriumMoisture_ρθq()) Y = KID.initialise_state(KID.NonEquilibriumMoisture_ρθq(), KID.NoPrecipitation(), init) dY = Y / 10 - t = 13.0 - - KID.advection_tendency!(KID.NonEquilibriumMoisture_ρθq(), dY, Y, aux, 1.0) - @test dY == Y / 10 KID.precompute_aux_thermo!(KID.NonEquilibriumMoisture_ρθq(), dY, Y, aux, t) @test aux.moisture_variables.ρ_dry == aux.moisture_variables.ρ .- Y.ρq_tot @@ -506,6 +547,29 @@ end @test aux.moisture_variables.θ_dry == TD.dry_pottemp.(thermo_params, aux.moisture_variables.T, aux.moisture_variables.ρ_dry) + aux = KID.initialise_aux(Float64, init, params, (; dt = 2.0), 0.0, face_space, KID.NonEquilibriumMoisture_ρθq()) + Y = KID.initialise_state(KID.NonEquilibriumMoisture_ρθq(), KID.Precipitation2M(CMT.SB2006Type()), init) + dY = Y / 10 + + KID.precompute_aux_precip!(KID.Precipitation2M(CMT.SB2006Type()), dY, Y, aux, t) + tmp = @. KID.precip_helper_sources_2M!( + aux.params, + aux.moisture_variables.q_tot, + aux.moisture_variables.q_liq, + aux.precip_variables.q_rai, + aux.precip_variables.N_liq, + aux.precip_variables.N_rai, + aux.moisture_variables.T, + aux.moisture_variables.ρ, + aux.TS.dt, + KID.Precipitation2M(CMT.SB2006Type()), + ) + @test aux.precip_sources.S_q_rai ≈ tmp.S_q_rai atol = eps(Float64) * 10 + @test aux.precip_sources.S_q_tot ≈ tmp.S_q_tot atol = eps(Float64) * 10 + @test aux.precip_sources.S_q_liq ≈ tmp.S_q_liq atol = eps(Float64) * 10 + @test aux.precip_sources.S_N_liq ≈ tmp.S_N_liq atol = eps(Float64) * 10 + @test aux.precip_sources.S_N_rai ≈ tmp.S_N_rai atol = eps(Float64) * 10 + end @testset "aerosol activation for 2M schemes" begin @@ -526,6 +590,21 @@ end @test tmp isa NamedTuple @test tmp.S_Nl ≈ -tmp.S_Na @test 0 < tmp.S_Nl < N_aer_0 + + #action + tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, 290.0, p, ρ, ρw) + + #test + @test tmp.S_Nl ≈ 0.0 atol = eps(Float64) + @test tmp.S_Na ≈ 0.0 atol = eps(Float64) + + #action + tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, 0.0) + + #test + @test tmp.S_Nl ≈ 0.0 atol = eps(Float64) + @test tmp.S_Na ≈ 0.0 atol = eps(Float64) + end # calibration pipeline unit tests From 969bad94651b0822f117cad75d25308116a96501 Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Fri, 14 Apr 2023 18:44:59 -0700 Subject: [PATCH 14/16] minor --- test/experiments/calibrations/grid_search/grid_search.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/experiments/calibrations/grid_search/grid_search.jl b/test/experiments/calibrations/grid_search/grid_search.jl index ffd06258..3e150244 100644 --- a/test/experiments/calibrations/grid_search/grid_search.jl +++ b/test/experiments/calibrations/grid_search/grid_search.jl @@ -26,7 +26,8 @@ model_settings = get_model_config(param_names) config["model"] = model_settings param_names = collect(keys(parameters)) truth = KID.get_obs!(config) -ref_stats = KID.ReferenceStatistics(truth, config["statistics"]) +ref_stats_list = KID.make_ref_stats_list(obs, config["statistics"], KID.get_numbers_from_config(config)...) +ref_stats = KID.combine_ref_stats(ref_stats_list) NC.Dataset(output_file_name, "c") do ds for (combination_num, (param_name_1, param_name_2)) in enumerate(combinations(param_names, 2)) From bdcabe4c8eee59a3dcf5c406c1b98f5efa648a8b Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Fri, 14 Apr 2023 18:58:16 -0700 Subject: [PATCH 15/16] minor --- test/experiments/calibrations/grid_search/grid_search.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/experiments/calibrations/grid_search/grid_search.jl b/test/experiments/calibrations/grid_search/grid_search.jl index 3e150244..aaf12674 100644 --- a/test/experiments/calibrations/grid_search/grid_search.jl +++ b/test/experiments/calibrations/grid_search/grid_search.jl @@ -26,7 +26,7 @@ model_settings = get_model_config(param_names) config["model"] = model_settings param_names = collect(keys(parameters)) truth = KID.get_obs!(config) -ref_stats_list = KID.make_ref_stats_list(obs, config["statistics"], KID.get_numbers_from_config(config)...) +ref_stats_list = KID.make_ref_stats_list(truth, config["statistics"], KID.get_numbers_from_config(config)...) ref_stats = KID.combine_ref_stats(ref_stats_list) NC.Dataset(output_file_name, "c") do ds From 5bdee22f8450da2294cda3472ba66518a00b78fb Mon Sep 17 00:00:00 2001 From: Sajjad Azimi Date: Thu, 20 Apr 2023 18:35:39 -0700 Subject: [PATCH 16/16] changing aerosol sources; using an approximate rate instead of updating the state Y --- src/driveKiD/helper_functions.jl | 4 ++-- src/driveKiD/tendency.jl | 7 ++++--- test/unit_tests/unit_test.jl | 7 ++++--- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/driveKiD/helper_functions.jl b/src/driveKiD/helper_functions.jl index 4ded338d..055935bf 100644 --- a/src/driveKiD/helper_functions.jl +++ b/src/driveKiD/helper_functions.jl @@ -8,7 +8,7 @@ end """ Returns the number of new activated aerosol particles and updates aerosol number density """ -@inline function aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw) +@inline function aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw, dt) microphys_params = KP.microphysics_params(params) thermo_params = CM.Parameters.thermodynamics_params(microphys_params) @@ -36,7 +36,7 @@ end return (; S_Nl, S_Na) end - S_Nl = max(0, N_act - (N_aer_0 - N_aer)) + S_Nl = max(0, N_act - (N_aer_0 - N_aer)) / dt S_Na = -S_Nl return (; S_Nl, S_Na) diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl index f6019b86..e83f0679 100644 --- a/src/driveKiD/tendency.jl +++ b/src/driveKiD/tendency.jl @@ -557,9 +557,10 @@ end aux.moisture_variables.p, aux.moisture_variables.ρ, ClimaCore.Operators.InterpolateF2C().(aux.prescribed_velocity.ρw.components.data.:1), + aux.TS.dt, ) - @. Y.N_liq += tmp.S_Nl - @. Y.N_aer += tmp.S_Na + aux.aerosol_variables.S_N_aer = tmp.S_Na + aux.precip_sources.S_N_liq = tmp.S_Nl aux.precip_variables.q_rai = Y.ρq_rai ./ aux.moisture_variables.ρ aux.precip_variables.N_liq = Y.N_liq @@ -579,7 +580,7 @@ end aux.precip_sources.S_q_rai = tmp.S_q_rai aux.precip_sources.S_q_tot = tmp.S_q_tot aux.precip_sources.S_q_liq = tmp.S_q_liq - aux.precip_sources.S_N_liq = tmp.S_N_liq + @. aux.precip_sources.S_N_liq += tmp.S_N_liq aux.precip_sources.S_N_rai = tmp.S_N_rai If = CC.Operators.InterpolateC2F(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate()) diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl index db52a909..5c3ea07d 100644 --- a/test/unit_tests/unit_test.jl +++ b/test/unit_tests/unit_test.jl @@ -582,9 +582,10 @@ end p = 1e5 ρ = 1.0 ρw = 2.0 + dt = 1.0 #action - tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw) + tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw, dt) #test @test tmp isa NamedTuple @@ -592,14 +593,14 @@ end @test 0 < tmp.S_Nl < N_aer_0 #action - tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, 290.0, p, ρ, ρw) + tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, 290.0, p, ρ, ρw, dt) #test @test tmp.S_Nl ≈ 0.0 atol = eps(Float64) @test tmp.S_Na ≈ 0.0 atol = eps(Float64) #action - tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, 0.0) + tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, 0.0, dt) #test @test tmp.S_Nl ≈ 0.0 atol = eps(Float64)