diff --git a/Artifacts.toml b/Artifacts.toml index 80d4f0814..a8d6a4acd 100644 --- a/Artifacts.toml +++ b/Artifacts.toml @@ -1,7 +1,7 @@ [AIDA_ice_nucleation] -git-tree-sha1 = "4a298d3558691217adb253b4053611931d62c5c2" +git-tree-sha1 = "91a1a51fb6ebb54282e65de1a2aace5e4ef1c45a" lazy = true [[AIDA_ice_nucleation.download]] - sha256 = "b513219566a3d1125acf1515f0d80b5d212f2bb78acead15137dfff47b503f9d" + sha256 = "df37e510212e328c30f052938537ac29dd217206c266c501c2fb797faa3f4c8a" url = "https://caltech.box.com/shared/static/imo0pyr6o0m9txu8racztmt7lz5n9xf1.gz" diff --git a/papers/Project.toml b/papers/Project.toml index 3438a730f..a60375399 100644 --- a/papers/Project.toml +++ b/papers/Project.toml @@ -5,6 +5,7 @@ ClimaUtilities = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513" CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/papers/ice_nucleation_2024/AIDA_calibrations.jl b/papers/ice_nucleation_2024/AIDA_calibrations.jl index 4adbceb2a..16700a0e9 100644 --- a/papers/ice_nucleation_2024/AIDA_calibrations.jl +++ b/papers/ice_nucleation_2024/AIDA_calibrations.jl @@ -14,20 +14,32 @@ include(joinpath(pkgdir(CM), "papers", "ice_nucleation_2024", "calibration.jl")) # Helper functions function unpack_data(data_file_name) + file_path = CM.ArtifactCalling.AIDA_ice_nucleation(data_file_name) return DelimitedFiles.readdlm(file_path, skipstart = 125) end +function grab_data(unpacked_data) + + AIDA_t_profile = unpacked_data[:, 1] + AIDA_T_profile = unpacked_data[:, 3] + AIDA_P_profile = unpacked_data[:, 2] * 1e2 # hPa to Pa + AIDA_ICNC = unpacked_data[:, 6] .* 1e6 # Nᵢ [m^-3] + AIDA_e = unpacked_data[:, 4] + + data = (; AIDA_t_profile, AIDA_T_profile, AIDA_P_profile, AIDA_ICNC, AIDA_e) + return data +end moving_average(data, n) = - [sum(@view data[i:(i + n - 1)]) / n for i in 1:(length(data) - (n - 1))] + [sum(@view data[i:(i + n)]) / n for i in 1:(length(data) - n)] # Defining data names, start/end times, etc. data_file_names = ["in05_17_aida.edf", "in05_18_aida.edf", "in05_19_aida.edf"] plot_names = ["IN0517", "IN0518", "IN0519"] -end_sim = 25 # Loss func looks at last end_sim timesteps only -start_time_list = [195, 180, 80] .+ 100 # AIDA time starts at -100 seconds. Add freezing onset as needed. -end_time_list = [290, 290, 170] .+ 100 # approximate time freezing stops -moving_average_n = 20 # average every n points -updrafts = [FT(1.5), FT(1.4), FT(5)] # updrafts matching AIDA cooling rate +end_sim = 25 # Loss func looks at last end_sim timesteps only +start_time_list = [195, 180, 80, 50, 35] # freezing onset +end_time_list = [290, 290, 170, 800, 800] # approximate time freezing stops +moving_average_n = 20 # average every n points +updrafts = [FT(1.5), FT(1.4), FT(5), FT(1.5), FT(1.5)] # updrafts matching AIDA cooling rate # Additional definitions tps = TD.Parameters.ThermodynamicsParameters(FT) @@ -42,42 +54,35 @@ for (exp_index, data_file_name) in enumerate(data_file_names) plot_name = plot_names[exp_index] w = updrafts[exp_index] start_time = start_time_list[exp_index] + start_time_index = start_time .+ 100 end_time = end_time_list[exp_index] + end_time_index = end_time .+ 100 ## Moving average to smooth data. - moving_average_start = Int32(start_time + (moving_average_n / 2 - 1)) - moving_average_end = Int32(end_time - (moving_average_n / 2)) - t_max = - (end_time - 100 - (moving_average_n / 2)) - # AIDA time starts at -100 seconds (t = 0 at index 101) - (start_time - 100 + (moving_average_n / 2 - 1)) # duration of simulation + moving_average_start_index = Int32(start_time_index + (moving_average_n / 2)) + moving_average_end_index = Int32(end_time_index - (moving_average_n / 2)) + t_max = moving_average_end_index - moving_average_start_index - ### Check for data file in AIDA_data folder. + ### Check for and grab data in AIDA_data folder. chamber_data = unpack_data(data_file_name) - ICNC = chamber_data[start_time:end_time, 6] .* 1e6 # Nᵢ [m^-3] - t_profile = chamber_data[moving_average_start+1:moving_average_end, 1] .- moving_average_start .+ 100 - T_profile = chamber_data[moving_average_start+1:moving_average_end, 3] - P_profile = chamber_data[moving_average_start+1:moving_average_end, 2] * 1e2 + AIDA_data = grab_data(chamber_data) + (; AIDA_t_profile, AIDA_T_profile, AIDA_P_profile, AIDA_ICNC, AIDA_e) = AIDA_data + t_profile = AIDA_t_profile[moving_average_start_index:moving_average_end_index] .- (moving_average_start_index - 101) + T_profile = AIDA_T_profile[moving_average_start_index:moving_average_end_index] + P_profile = AIDA_P_profile[moving_average_start_index:moving_average_end_index] + ICNC_profile = AIDA_ICNC[moving_average_start_index:moving_average_end_index] + e_profile = AIDA_e[moving_average_start_index:moving_average_end_index] + # S_l_profile = e_profile ./ (TD.saturation_vapor_pressure.(tps, T_profile, TD.Liquid())) params = AIDA_IN05_params(FT, w, t_max, t_profile, T_profile, P_profile) IC = AIDA_IN05_IC(FT, data_file_name) - frozen_frac = ICNC ./ (IC[7] + IC[8] + IC[9]) + Nₜ = IC[7] + IC[8] + IC[9] + frozen_frac = AIDA_ICNC[start_time_index:end_time_index] ./ Nₜ frozen_frac_moving_mean = moving_average(frozen_frac, moving_average_n) - ICNC_moving_avg = moving_average(ICNC, moving_average_n) + ICNC_moving_avg = moving_average(AIDA_ICNC[start_time_index:end_time_index], moving_average_n) Γ = 0.1 * LinearAlgebra.I * (maximum(frozen_frac_moving_mean) - minimum(frozen_frac_moving_mean)) # coeff is an estimated of the noise - ## Calculating some variables from AIDA data. - chamber_S_l = zeros(FT, size(chamber_data[start_time:end_time, 4], 1), size(chamber_data[start_time:end_time, 4], 2)) - for (i, e) in enumerate(chamber_data[start_time:end_time, 4]) - temp = chamber_data[start_time:end_time, 3][i] - eₛ = TD.saturation_vapor_pressure(tps, temp, TD.Liquid()) - chamber_S_l[i] = e / eₛ - end - chamber_r_l = zeros(FT, size(chamber_data[start_time:end_time, 8], 1), size(chamber_data[start_time:end_time, 8], 2)) - for (i, r_l) in enumerate(chamber_data[start_time:end_time, 8] .* 1e-6) - chamber_r_l[i] = r_l < 0 ? FT(0) : r_l - end - ### Calibration. EKI_output = calibrate_J_parameters_EKI(FT, "ABHOM", params, IC, frozen_frac_moving_mean, end_sim, Γ) UKI_output = calibrate_J_parameters_UKI(FT, "ABHOM", params, IC, frozen_frac_moving_mean, end_sim, Γ) @@ -98,23 +103,8 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ## Plotting AIDA data. AIDA_data_fig = MK.Figure(size = (800, 600), fontsize = 24) ax1 = MK.Axis(AIDA_data_fig[1, 1], ylabel = "ICNC [m^-3]", xlabel = "time [s]", title = "AIDA data $plot_name") - MK.lines!( - ax1, - chamber_data[start_time:end_time, 1], - ICNC, - label = "Raw AIDA", - color =:blue, - linestyle =:dash, - linewidth = 2, - ) - MK.lines!( - ax1, - chamber_data[moving_average_start:moving_average_end, 1], - ICNC_moving_avg, - label = "AIDA moving mean", - linewidth = 2.5, - color =:blue, - ) + MK.lines!(ax1, AIDA_t_profile, AIDA_ICNC, label = "Raw AIDA", color =:blue, linestyle =:dash, linewidth = 2) + MK.lines!(ax1, t_profile, ICNC_moving_avg, label = "AIDA moving mean", linewidth = 2.5, color =:blue) MK.axislegend(ax1, framevisible = true, labelsize = 12, position = :rc) MK.save("$plot_name"*"_ICNC.svg", AIDA_data_fig) @@ -139,65 +129,43 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ax_parcel_7 = MK.Axis(calibrated_parcel_fig[1, 3], ylabel = "pressure [Pa]", xlabel = "time [s]") ax_parcel_8 = MK.Axis(calibrated_parcel_fig[2, 3], ylabel = "Nₐ [m^-3]", xlabel = "time [s]") - MK.lines!(ax_parcel_1, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[1, :], label = "EKI Calib Liq", color = :orange) # label = "liquid" - MK.lines!(ax_parcel_1, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[1, :], label = "UKI Calib Liq", color = :fuchsia) # label = "liquid" - MK.lines!(ax_parcel_1, parcel_default.t .+ (moving_average_n / 2), parcel_default[1, :], label = "default", color = :darkorange2) - MK.lines!(ax_parcel_1, EKI_parcel.t .+ (moving_average_n / 2), S_i.(tps, EKI_parcel[3, :], EKI_parcel[1, :]), label = "EKI Calib Ice", color = :green) - # MK.lines!( - # ax_parcel_1, - # chamber_data[start_time:end_time, 1] .- (start_time - 100), - # vec(chamber_S_l), - # label = "chamber", - # color = :blue - # ) - - MK.lines!(ax_parcel_2, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[5, :], color = :orange) - MK.lines!(ax_parcel_2, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[5, :], color = :fuchsia) - MK.lines!(ax_parcel_2, parcel_default.t .+ (moving_average_n / 2), parcel_default[5,:], color = :darkorange2) + MK.lines!(ax_parcel_1, EKI_parcel.t, EKI_parcel[1, :], label = "EKI Calib Liq", color = :orange) # label = "liquid" + MK.lines!(ax_parcel_1, UKI_parcel.t, UKI_parcel[1, :], label = "UKI Calib Liq", color = :fuchsia) # label = "liquid" + MK.lines!(ax_parcel_1, parcel_default.t, parcel_default[1, :], label = "default", color = :darkorange2) + MK.lines!(ax_parcel_1, EKI_parcel.t, S_i.(tps, EKI_parcel[3, :], EKI_parcel[1, :]), label = "EKI Calib Ice", color = :green) + # MK.lines!(ax_parcel_1, t_profile, S_l_profile, label = "chamber", color = :blue) - MK.lines!(ax_parcel_3, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[3, :], color = :orange) - MK.lines!(ax_parcel_3, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[3, :], color = :fuchsia) - MK.lines!(ax_parcel_3, parcel_default.t .+ (moving_average_n / 2), parcel_default[3, :], color = :darkorange2) - MK.lines!( - ax_parcel_3, - chamber_data[start_time:end_time, 1] .- (start_time - 100), - chamber_data[start_time:end_time, 3], - color = :blue, - linestyle =:dash, - ) + MK.lines!(ax_parcel_2, EKI_parcel.t, EKI_parcel[5, :], color = :orange) + MK.lines!(ax_parcel_2, UKI_parcel.t, UKI_parcel[5, :], color = :fuchsia) + MK.lines!(ax_parcel_2, parcel_default.t, parcel_default[5,:], color = :darkorange2) - MK.lines!(ax_parcel_4, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[6, :], color = :orange) - MK.lines!(ax_parcel_4, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[6, :], color = :fuchsia) - MK.lines!(ax_parcel_4, parcel_default.t .+ (moving_average_n / 2), parcel_default[6, :], color = :darkorange2) + MK.lines!(ax_parcel_3, EKI_parcel.t, EKI_parcel[3, :], color = :orange) + MK.lines!(ax_parcel_3, UKI_parcel.t, UKI_parcel[3, :], color = :fuchsia) + MK.lines!(ax_parcel_3, parcel_default.t, parcel_default[3, :], color = :darkorange2) + MK.lines!(ax_parcel_3, t_profile, T_profile, color = :blue, linestyle =:dash) - MK.lines!(ax_parcel_5, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[8, :], color = :orange) - MK.lines!(ax_parcel_5, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[8, :], color = :fuchsia) - MK.lines!(ax_parcel_5, parcel_default.t .+ (moving_average_n / 2), parcel_default[8, :], color = :darkorange2) + MK.lines!(ax_parcel_4, EKI_parcel.t, EKI_parcel[6, :], color = :orange) + MK.lines!(ax_parcel_4, UKI_parcel.t, UKI_parcel[6, :], color = :fuchsia) + MK.lines!(ax_parcel_4, parcel_default.t, parcel_default[6, :], color = :darkorange2) - MK.lines!(ax_parcel_6, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[9, :], color = :orange, label = "EKI") - MK.lines!(ax_parcel_6, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[9, :], color = :fuchsia, label = "UKI") - MK.lines!(ax_parcel_6, parcel_default.t .+ (moving_average_n / 2), parcel_default[9, :], color = :darkorange2, label = "Pre-Calib") - MK.lines!(ax_parcel_6, - chamber_data[start_time:end_time, 1] .- (start_time - 100), - ICNC, - color = :blue, - label = "AIDA", - ) - error = fill(sqrt(Γ[1,1]) * 2, length(chamber_data[start_time:end_time, 1] .- (start_time - 100))) - MK.errorbars!(ax_parcel_6, chamber_data[start_time:end_time, 1] .- (start_time - 100), ICNC ./ (IC[7] + IC[8] + IC[9]), error) + MK.lines!(ax_parcel_5, EKI_parcel.t, EKI_parcel[8, :], color = :orange) + MK.lines!(ax_parcel_5, UKI_parcel.t, UKI_parcel[8, :], color = :fuchsia) + MK.lines!(ax_parcel_5, parcel_default.t, parcel_default[8, :], color = :darkorange2) - MK.lines!(ax_parcel_7, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[2, :], color = :orange) - MK.lines!(ax_parcel_7, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[2, :], color = :fuchsia) - MK.lines!( - ax_parcel_7, - chamber_data[start_time:end_time, 1] .- (start_time - 100), - chamber_data[start_time:end_time, 2] .* 1e2, - color = :blue, - linestyle =:dash, - ) + MK.lines!(ax_parcel_6, EKI_parcel.t, EKI_parcel[9, :], color = :orange, label = "EKI") + MK.lines!(ax_parcel_6, UKI_parcel.t, UKI_parcel[9, :], color = :fuchsia, label = "UKI") + MK.lines!(ax_parcel_6, parcel_default.t, parcel_default[9, :], color = :darkorange2, label = "Pre-Calib") + MK.lines!(ax_parcel_6, t_profile, ICNC_profile, color = :blue, label = "AIDA",) + + error = fill(sqrt(Γ[1,1]) * 2, length(AIDA_t_profile[start_time_index:end_time_index] .- start_time)) + MK.errorbars!(ax_parcel_6, AIDA_t_profile[start_time_index:end_time_index] .- start_time, AIDA_ICNC[start_time_index:end_time_index] ./ Nₜ, error) + + MK.lines!(ax_parcel_7, EKI_parcel.t, EKI_parcel[2, :], color = :orange) + MK.lines!(ax_parcel_7, UKI_parcel.t, UKI_parcel[2, :], color = :fuchsia) + MK.lines!(ax_parcel_7, t_profile, P_profile, color = :blue) #, linestyle =:dash - MK.lines!(ax_parcel_8, EKI_parcel.t .+ (moving_average_n / 2), EKI_parcel[7, :], color = :orange) - MK.lines!(ax_parcel_8, UKI_parcel.t .+ (moving_average_n / 2), UKI_parcel[7, :], color = :fuchsia) + MK.lines!(ax_parcel_8, EKI_parcel.t, EKI_parcel[7, :], color = :orange) + MK.lines!(ax_parcel_8, UKI_parcel.t, UKI_parcel[7, :], color = :fuchsia) MK.axislegend(ax_parcel_6, framevisible = false, labelsize = 16, position = :rb) MK.save("$plot_name"*"_calibrated_parcel_fig.svg", calibrated_parcel_fig) @@ -208,31 +176,31 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ax_compare = MK.Axis(ICNC_comparison_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name") # MK.lines!( # ax_compare, - # parcel_default.t .+ (moving_average_n / 2), - # parcel_default[9, :]./ (IC[7] + IC[8] + IC[9]), + # parcel_default.t, + # parcel_default[9, :]./ Nₜ, # label = "CM.jl Parcel (Pre-Calib)", # linewidth = 2.5, # color =:orangered, # ) MK.lines!( ax_compare, - EKI_parcel.t .+ (moving_average_n / 2), - EKI_parcel[9, :]./ (IC[7] + IC[8] + IC[9]), + EKI_parcel.t, + EKI_parcel[9, :]./ Nₜ, label = "CM.jl Parcel (EKI Calibrated)", linewidth = 2.5, color =:orange, ) MK.lines!( ax_compare, - UKI_parcel.t .+ (moving_average_n / 2), - UKI_parcel[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel.t, + UKI_parcel[9, :]./ Nₜ, label = "CM.jl Parcel (UKI Calibrated)", linewidth = 2.5, color =:fuchsia, ) MK.lines!( ax_compare, - chamber_data[start_time:end_time, 1] .- (start_time - 100), + AIDA_t_profile[start_time_index:end_time_index] .- start_time, frozen_frac, label = "Raw AIDA", linewidth = 2, @@ -241,14 +209,14 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ) MK.lines!( ax_compare, - chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100), + t_profile, frozen_frac_moving_mean, label = "AIDA Moving Avg", linewidth = 2.5, color =:blue ) error = fill(sqrt(Γ[1,1]) * 2, length(frozen_frac_moving_mean)) - MK.errorbars!(ax_compare, chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100), frozen_frac_moving_mean, error) + MK.errorbars!(ax_compare, t_profile, frozen_frac_moving_mean, error) MK.axislegend(ax_compare, framevisible = false, labelsize = 20, position = :rb) MK.save("$plot_name"*"_ICNC_comparison_fig.svg", ICNC_comparison_fig) @@ -265,7 +233,7 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ax_spread = MK.Axis(UKI_spread_fig[1, 1], ylabel = "Frozen Fraction [-]", xlabel = "time [s]", title = "$plot_name") MK.lines!( ax_spread, - chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100), + t_profile, frozen_frac_moving_mean, label = "AIDA Moving Avg", linewidth = 2.5, @@ -273,50 +241,50 @@ for (exp_index, data_file_name) in enumerate(data_file_names) ) MK.lines!( ax_spread, - UKI_parcel_1.t .+ (moving_average_n / 2), - UKI_parcel_1[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel_1.t, + UKI_parcel_1[9, :]./ Nₜ, linewidth = 2.5, color =:grey80, ) MK.lines!( ax_spread, - UKI_parcel_2.t .+ (moving_average_n / 2), - UKI_parcel_2[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel_2.t, + UKI_parcel_2[9, :]./ Nₜ, linewidth = 2.5, color =:grey60, ) MK.lines!( ax_spread, - UKI_parcel_3.t .+ (moving_average_n / 2), - UKI_parcel_3[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel_3.t, + UKI_parcel_3[9, :]./ Nₜ, linewidth = 2.5, color =:grey40, ) MK.lines!( ax_spread, - UKI_parcel_4.t .+ (moving_average_n / 2), - UKI_parcel_4[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel_4.t, + UKI_parcel_4[9, :]./ Nₜ, linewidth = 2.5, color =:grey25, ) MK.lines!( ax_spread, - UKI_parcel_5.t .+ (moving_average_n / 2), - UKI_parcel_5[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel_5.t, + UKI_parcel_5[9, :]./ Nₜ, linewidth = 2.5, color =:grey12, ) MK.lines!( ax_spread, - UKI_parcel.t .+ (moving_average_n / 2), - UKI_parcel[9, :]./ (IC[7] + IC[8] + IC[9]), + UKI_parcel.t, + UKI_parcel[9, :]./ Nₜ, label = "CM.jl Parcel (UKI Calibrated)", linewidth = 2.5, color =:fuchsia, linestyle = :dash, ) error = fill(sqrt(Γ[1,1]) * 2, length(frozen_frac_moving_mean)) - MK.errorbars!(ax_spread, chamber_data[moving_average_start:moving_average_end, 1] .- (start_time - 100), frozen_frac_moving_mean, error) + MK.errorbars!(ax_spread, t_profile, frozen_frac_moving_mean, error) MK.save("$plot_name"*"_UKI_spread_fig.svg", UKI_spread_fig) #! format: on diff --git a/papers/ice_nucleation_2024/AIDA_data_artifact/OutputArtifacts.toml b/papers/ice_nucleation_2024/AIDA_data_artifact/OutputArtifacts.toml index 80d4f0814..3386b8ade 100644 --- a/papers/ice_nucleation_2024/AIDA_data_artifact/OutputArtifacts.toml +++ b/papers/ice_nucleation_2024/AIDA_data_artifact/OutputArtifacts.toml @@ -1,7 +1,6 @@ [AIDA_ice_nucleation] -git-tree-sha1 = "4a298d3558691217adb253b4053611931d62c5c2" -lazy = true +git-tree-sha1 = "91a1a51fb6ebb54282e65de1a2aace5e4ef1c45a" [[AIDA_ice_nucleation.download]] - sha256 = "b513219566a3d1125acf1515f0d80b5d212f2bb78acead15137dfff47b503f9d" + sha256 = "df37e510212e328c30f052938537ac29dd217206c266c501c2fb797faa3f4c8a" url = "https://caltech.box.com/shared/static/imo0pyr6o0m9txu8racztmt7lz5n9xf1.gz" diff --git a/papers/ice_nucleation_2024/AIDA_data_artifact/create_AIDA_artifacts.jl b/papers/ice_nucleation_2024/AIDA_data_artifact/create_AIDA_artifacts.jl index eeb4de62b..1215519e9 100644 --- a/papers/ice_nucleation_2024/AIDA_data_artifact/create_AIDA_artifacts.jl +++ b/papers/ice_nucleation_2024/AIDA_data_artifact/create_AIDA_artifacts.jl @@ -3,10 +3,11 @@ # of a function in CloudMicrophysics that is called by ClimaAtmos). For more info, please # visit the ClimaArtifacts.jl repo on GitHub! -# To use in VS Code, first activate the AIDA_data_artifact environment. +# In terminal, make sure your pwd is in AIDA_data_artifact. Then, +# activate Julia at --project=papers/ice_nucleation_2024/AIDA_data_artifact. +# You are in the right place if you enter the shell (tap ; on keyboard) and +# see that you are in the AIDA_data_artifact folder (NOT CloudMicrophysics.jl) # Then, copy and paste the following script to the Julia REPL. -# To use in terminal, activate Julia at --project=papers/ice_nucleation_2024/AIDA_data_artifact. -# Then, copy and paste the following script in your now opened Julia REPL. # If you do not already have ClimaArtifactsHelper, go into the Julia REPL and type: # Pkg.develop(url="https://github.com/CliMA/ClimaArtifacts.git", subdir="ClimaArtifactsHelper.jl") @@ -25,7 +26,7 @@ create_artifact_guided(data_folder; artifact_name = data_set_name) # After getting an updated "OutputArtifacts.toml" file in this folder, # Add it to the list of artifacts in CloudMicrophysics.jl's Artifact.toml file. -# Add `lazy = true` under the `git-tree-sha1` line in Artifact.toml i fyou would like +# Add `lazy = true` under the `git-tree-sha1` line in Artifact.toml if you would like # the artifact to be downloaded lazily. # To use data, make sure your script has "using LazyArtifacts" at the top. diff --git a/papers/ice_nucleation_2024/calibration_setup.jl b/papers/ice_nucleation_2024/calibration_setup.jl index 15f980689..83b2f3683 100644 --- a/papers/ice_nucleation_2024/calibration_setup.jl +++ b/papers/ice_nucleation_2024/calibration_setup.jl @@ -213,17 +213,17 @@ function AIDA_IN05_IC(FT, data_file) eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) Sₗ = FT(e / eₛ) elseif data_file == "in05_18_aida.edf" - Nₗ = FT(209.46 * 1e6) - Nᵢ = FT(1.53 * 1e6) + Nₗ = FT(210.28 * 1e6) + Nᵢ = FT(1.73 * 1e6) Nₐ = FT(275 * 1e6) - Nₗ - Nᵢ - r₀ = FT(7.03467 * 1e-6) + r₀ = FT(6.99998 * 1e-6) p₀ = FT(873.322 * 1e2) T₀ = FT(237.175) qₗ = FT(Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) # 1.2 should be ρₐ qᵢ = FT(Nᵢ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)) m_l = Nₗ * ρₗ * 4 * π / 3 * r₀^3 m_i = Nᵢ * ρᵢ * 4 * π / 3 * r₀^3 - e = FT(28.1324) + e = FT(28.086) qᵥ = (e / R_v / T₀) / ((p₀ - e) / (R_d * T₀) + e / R_v / T₀ + m_l + m_i) q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) Rₐ = TD.gas_constant_air(tps, q) @@ -232,7 +232,7 @@ function AIDA_IN05_IC(FT, data_file) elseif data_file == "in05_19_aida.edf" Nₐ = FT(0) Nₗ = FT(180 * 1e6) - Nᵢ = FT(0.49 * 1e6) + Nᵢ = FT(0.882 * 1e6) r₀ = FT(6.5e-6) # !!missing in dataset!! p₀ = FT(724.545 * 1e2) T₀ = FT(237.639) diff --git a/parcel/ParcelCommon.jl b/parcel/ParcelCommon.jl index 3c70661f7..b0eb2d0f7 100644 --- a/parcel/ParcelCommon.jl +++ b/parcel/ParcelCommon.jl @@ -13,11 +13,10 @@ S_i(tps, T, S_liq) = ξ(tps, T) * S_liq # Interpolating for cooling/expansion rate function AIDA_rate(model_t, data_t, data) - index = Int32(model_t) + 1 data_at_t = Intp.linear_interpolation(data_t, data) - if index < length(data) - 1 - return (data_at_t(model_t + 2) - data_at_t(model_t + 1)) + if model_t < data_t[end] + return (data_at_t(model_t + 1) - data_at_t(model_t)) else return 0 end diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 15cd58295..2aa3b67c3 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -208,7 +208,7 @@ function benchmark_test(FT) bench_press( CMN.terminal_velocity, (ice, ch2022.small_ice, ρ_air, q_ice), - 350, + 400, ) bench_press(CM1.terminal_velocity, (rain, ch2022.rain, ρ_air, q_rai), 750) bench_press(