diff --git a/Project.toml b/Project.toml index f3212eba..fcaa1bf9 100644 --- a/Project.toml +++ b/Project.toml @@ -22,9 +22,11 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed" +TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] @@ -47,7 +49,6 @@ Plots = "1.29" ProgressLogging = "0.1" Random = "1" SpecialFunctions = "2.3" -Statistics = "1.10" TerminalLoggers = "0.1" Thermodynamics = "0.12.6" julia = "1.8" diff --git a/src/Common/initial_condition.jl b/src/Common/initial_condition.jl index ccbae713..b4240e3c 100644 --- a/src/Common/initial_condition.jl +++ b/src/Common/initial_condition.jl @@ -329,6 +329,10 @@ function p3_initial_condition( # use approximate temperature values from # sounding: kin1d/soundings/Tsfc2.txt # (path in p3 fortran github repo) + # although I think it could also be grabbed from + # this website: https://weather.uwyo.edu/upperair/sounding.html + # temperature is a constant gradient with a melting layer at 500 m, + # but of course that could easily be changed T::FT = -0.004 * (z - 500) + 273.15 # temperature p::FT = 990 - (0.1 * z) # pressure _q = TD.PhasePartition(q_tot, q_liq, q_ice) diff --git a/src/Common/ode_utils.jl b/src/Common/ode_utils.jl index 0f4b4b0a..54fbf0bc 100644 --- a/src/Common/ode_utils.jl +++ b/src/Common/ode_utils.jl @@ -141,6 +141,8 @@ function initialise_aux( tmp = similar(ip.q_tot), tmp2 = similar(ip.q_tot), tmp3 = similar(ip.q_tot), + tmp4 = similar(ip.q_tot), + tmp5 = similar(ip.q_tot), tmp_surface = similar(CC.Fields.level(ip.q_tot, 1), Tuple{FT, FT}), ) @@ -200,6 +202,18 @@ function initialise_aux( q_vap = ip.q_vap, ρq_vap = ip.ρq_vap, ) + p3_predicted_eltype = @NamedTuple{ + F_rim::FT, + ρ_rim::FT, + F_liq::FT, + } + p3_predicted = @. p3_predicted_eltype( + tuple( + copy(ip.ρq_rim / (ip.ρq_ice - ip.ρq_liqonice)), + copy(ip.ρq_rim / (ip.ρq_ice - ip.ρq_liqonice) / ip.B_rim), + copy(ip.ρq_liqonice / ip.ρq_ice), + ), + ) velocities = (; term_vel_rai = copy(ip.zero), term_vel_ice = copy(ip.zero), @@ -207,12 +221,12 @@ function initialise_aux( term_vel_N_ice = copy(ip.zero), ) precip_sources_eltype = @NamedTuple{ - q_tot::FT, - q_liq::FT, - q_rai::FT, - q_ice::FT, - q_rim::FT, - q_liqonice::FT, + ρq_tot::FT, + ρq_liq::FT, + ρq_rai::FT, + ρq_ice::FT, + ρq_rim::FT, + ρq_liqonice::FT, N_aer::FT, N_liq::FT, N_rai::FT, @@ -297,7 +311,7 @@ function initialise_aux( if precip isa CloudyPrecip aux = merge(aux, (; cloudy_params, cloudy_variables)) elseif precip isa PrecipitationP3 - aux = merge(aux, (; p3_boundary_condition)) + aux = merge(aux, (; p3_boundary_condition, p3_predicted)) end return aux diff --git a/src/Common/tendency.jl b/src/Common/tendency.jl index 9807ab81..6ed1859b 100644 --- a/src/Common/tendency.jl +++ b/src/Common/tendency.jl @@ -187,26 +187,26 @@ end @. B_rim = max(FT(0), Y.B_rim) # compute predicted quantities - ρ_r = aux.scratch.tmp - F_rim = aux.scratch.tmp2 - F_liq = aux.scratch.tmp3 - # prohibit ρ_r < 0 - @. ρ_r = ifelse(B_rim < eps(FT), eps(FT), max(FT(0), ρq_rim / B_rim)) - # keep F_r in range [0, 1) - @. F_rim = - ifelse((ρq_ice - ρq_liqonice) < eps(FT), FT(0), min(max(FT(0), ρq_rim / (ρq_ice - ρq_liqonice)), 1 - eps(FT))) - # prohibit F_liq < 0 - @. F_liq = ifelse(Y.ρq_ice < eps(FT), FT(0), Y.ρq_liqonice / Y.ρq_ice) + + p3 = ps.p3_params + (; F_rim, ρ_rim, F_liq) = aux.p3_predicted + F_r(ρq_ice, ρq_rim, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_rim / (ρq_ice - ρq_liqonice), 1 - eps(FT))), FT(0)) + ρ_r(ρq_rim, B_rim) = ifelse(B_rim > eps(FT), max(FT(0), min(ρq_rim / B_rim, p3.ρ_l)), FT(0)) + F_l(ρq_ice, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_liqonice / ρq_ice, FT(1 - eps(FT)))), FT(0)) + @. F_rim = F_r(ρq_ice, ρq_rim, ρq_liqonice) + @. ρ_rim = ρ_r(ρq_rim, B_rim) + @. F_liq = F_l(ρq_ice, ρq_liqonice) p3 = ps.p3_params Chen2022 = ps.Chen2022 (; term_vel_ice, term_vel_N_ice, term_vel_rai, term_vel_N_rai) = aux.velocities use_aspect_ratio = true + # TODO - should we be using ρ or ρ_dry in the velocity calculation? @. term_vel_N_ice = - getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_r, F_rim, F_liq, ρ_dry, use_aspect_ratio), 1) + getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_rim, F_rim, F_liq, ρ, use_aspect_ratio), 1) @. term_vel_ice = - getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_r, F_rim, F_liq, ρ_dry, use_aspect_ratio), 2) + getindex(CMP3.ice_terminal_velocity(p3, Chen2022, ρq_ice, N_ice, ρ_rim, F_rim, F_liq, ρ, use_aspect_ratio), 2) # adding 2M rain below: @@ -559,8 +559,267 @@ end end @inline function precompute_aux_precip_sources!(ps::PrecipitationP3, aux) - # TODO [P3] - return nothing + (; dt) = aux.TS + (; thermo_params, air_params, common_params) = aux + FT = eltype(thermo_params) + (; ts, ρ, T) = aux.thermo_variables + (; ρq_liq, ρq_ice, ρq_rim, ρq_liqonice, N_ice, B_rim, ρq_rai, N_rai) = aux.microph_variables + (; F_rim, ρ_rim, F_liq) = aux.p3_predicted + + dLdt_p3 = aux.scratch.tmp + dNdt = aux.scratch.tmp2 + dLdt_rim = aux.scratch.tmp3 + dLdt_liq = aux.scratch.tmp4 + ddt_B_rim = aux.scratch.tmp5 + + p3 = ps.p3_params + type = typeof(p3) + vel = ps.Chen2022 + sb2006 = ps.sb2006 + aps = CMP.AirProperties(FT) + + # helper type and wrapper to populate the tuple with sources + precip_sources_eltype = @NamedTuple{ + ρq_tot::FT, + ρq_liq::FT, + ρq_rai::FT, + ρq_ice::FT, + ρq_rim::FT, + ρq_liqonice::FT, + N_aer::FT, + N_liq::FT, + N_rai::FT, + N_ice::FT, + B_rim::FT, + } + to_sources(args...) = @. precip_sources_eltype(tuple(args...)) + + # zero out the aux sources + @. aux.precip_sources = to_sources(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) + + # calculate temporary variables + + F_r(ρq_ice, ρq_rim, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_rim / (ρq_ice - ρq_liqonice), 1 - eps(FT))), FT(0)) + ρ_r(ρq_rim, B_rim) = ifelse(B_rim > 0, max(FT(0), min(ρq_rim / B_rim, p3.ρ_l)), FT(0)) + F_l(ρq_ice, ρq_liqonice) = ifelse(ρq_ice > eps(FT), max(FT(0), min(ρq_liqonice / ρq_ice, FT(1 - eps(FT)))), FT(0)) + @. F_rim = F_r(ρq_ice, ρq_rim, ρq_liqonice) + @. ρ_rim = ρ_r(ρq_rim, B_rim) + @. F_liq = F_l(ρq_ice, ρq_liqonice) + + # melting: turned off for testing shedding + # if Bool(common_params.precip_sources) + melting = false + if melting # melting works with corresponding CM branch! + + @. dLdt_p3 = limit( + ρq_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_p3_tot, + ) + @. dLdt_rim = limit( + ρq_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_rim, + ) + @. dLdt_liq = limit( + ρq_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_liq, + ) + @. dNdt = limit( + N_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).dNdt_ice, + ) + @. ddt_B_rim= limit( + N_ice, + dt, + CMP3.ice_melt( + p3, + vel, + air_params, + thermo_params, + ρq_ice, + N_ice, + T, + ρ, + F_rim, + ρ_rim, + F_liq, + dt, + ).ddtB_rim, + ) + @. aux.precip_sources += to_sources( + 0, + 0, + dLdt_p3, # source of rain = sink of p3 ice + -dLdt_p3, # sink of p3 ice + -dLdt_rim, # sink of rime + dLdt_liq, # source of L_liq + 0, + 0, + dNdt, # source of rain number + -dNdt, # sink of ice number + -ddt_B_rim, # sink of B_rim + ) + end + shedding = false + if shedding # shedding works with corresponding CM branch! + @. dLdt_p3 = limit( + ρq_ice, + dt, + CMP3.ice_shed( + p3, + ρq_ice, + N_ice, + F_rim, + ρ_rim, + F_liq, + dt, + ).dLdt_p3_tot, + ) + @. dNdt = limit( + N_ice, + dt, + CMP3.ice_shed( + p3, + ρq_ice, + N_ice, + F_rim, + ρ_rim, + F_liq, + dt, + ).dNdt_rai, + ) + @. aux.precip_sources += to_sources( + 0, + 0, + dLdt_p3, # source of rain = sink of p3 ice + -dLdt_p3, # sink of p3 ice + 0, # sink of rime + -dLdt_p3, # sink of L_liq + 0, + 0, + dNdt, # source of rain number + 0, # sink of ice number + 0, # sink of B_rim + ) + end + collisions = true + if collisions # with rain - stable for corresponding CM branch (collision kernel) + # but haven't really been able to test much at all + # TODO - use ρ or ρ_dry? + @. dLdt_p3 = limit( + ρq_rai, + dt, + CMP3.ice_collisions( + sb2006.pdf_r, + p3, + vel, + aps, + thermo_params, + ts, + ρq_ice, + N_ice, + ρq_rai, + N_rai, + ρ, + F_rim, + ρ_rim, + F_liq, + T, + dt, + ).dLdt_p3_tot, + ) + @. dNdt = limit( + N_rai, + dt, + CMP3.ice_collisions( + sb2006.pdf_r, + p3, + vel, + aps, + thermo_params, + ts, + ρq_ice, + N_ice, + ρq_rai, + N_rai, + ρ, + F_rim, + ρ_rim, + F_liq, + T, + dt, + ).dNdt_rai, + ) + @. dLdt_p3 = ifelse(isnan(dLdt_p3), FT(0), dLdt_p3) + @. aux.precip_sources += to_sources( + 0, + 0, + -dLdt_p3, # sink of rain = source of p3 ice + dLdt_p3, # source of p3 ice + dLdt_p3, # source of rime + dLdt_p3, # source of L_liq + 0, + 0, + -dNdt, # sink of rain number + 0, # sink of ice number + dLdt_p3 / FT(900), # source of B_rim + ) + end end # helper function for precomputing aux sources for cloudy @@ -713,6 +972,20 @@ end return dY end @inline function precip_sources_tendency!(ms::MoistureP3, ps::PrecipitationP3, dY, Y, aux, t) + precompute_aux_precip_sources!(ps, aux) + + @. dY.ρq_tot += aux.precip_sources.ρq_tot + @. dY.ρq_liq += aux.precip_sources.ρq_liq + @. dY.ρq_rai += aux.precip_sources.ρq_rai + @. dY.ρq_ice += aux.precip_sources.ρq_ice + @. dY.ρq_rim += aux.precip_sources.ρq_rim + @. dY.ρq_liqonice += aux.precip_sources.ρq_liqonice + @. dY.B_rim += aux.precip_sources.B_rim + @. dY.N_ice += aux.precip_sources.N_ice + @. dY.N_rai += aux.precip_sources.N_rai + @. dY.N_aer += aux.precip_sources.N_aer + @. dY.N_liq += aux.precip_sources.N_liq + return dY end @inline function precip_sources_tendency!(ms::CloudyMoisture, ps::CloudyPrecip, dY, Y, aux, t) diff --git a/test/experiments/KiD_driver/KiD_driver.jl b/test/experiments/KiD_driver/KiD_driver.jl index e061dd02..2716448f 100644 --- a/test/experiments/KiD_driver/KiD_driver.jl +++ b/test/experiments/KiD_driver/KiD_driver.jl @@ -9,9 +9,27 @@ # Get the parameter values for the simulation include("parse_commandline.jl") -if !(@isdefined config) - config = parse_commandline() -end +bound = (; + ice_start = false, + _magnitude = Float64(0.5), + _q_flux = Float64(0.65e-3), + _N_flux = Float64(1e5), + _F_rim = Float64(0.8), + _F_liq = Float64(0.5), + _ρ_r_init = Float64(900), +) + +# if !(@isdefined config) +config = parse_commandline() +config["precipitation_choice"] = "PrecipitationP3" +config["moisture_choice"] = "MoistureP3" +config["n_elem"] = 36 +config["z_max"] = 3000 +config["t_end"] = 75 +config["w1"] = 0 +config["dt"] = Float64(0.5) +config["p3_boundary_condition"] = bound +# end ft_choice = config["FLOAT_TYPE"] diff --git a/test/plotting_utils.jl b/test/plotting_utils.jl index 9c078ae2..d1f9a203 100644 --- a/test/plotting_utils.jl +++ b/test/plotting_utils.jl @@ -97,7 +97,7 @@ function plot_final_aux_profiles(z_centers, aux, precip; output = "output") if precip isa CO.PrecipitationP3 p6 = Plots.plot(N_ice_end .* 1e-6, z_centers, xlabel = "N_ice [1/cm3]", ylabel = "z [m]") - p11 = Plots.plot(N_ice_end .* 1e-6, z_centers, xlabel = "N_ice [1/cm3]", ylabel = "z [m]") + p11 = Plots.plot(N_rai_end .* 1e-6, z_centers, xlabel = "N_rai [1/cm3]", ylabel = "z [m]") p12 = Plots.plot(q_liqonice_end .* 1e3, z_centers, xlabel = "q_liqonice [g/kg]", ylabel = "z [m]") p13 = Plots.plot(q_rim_end .* 1e3, z_centers, xlabel = "q_rim [g/kg]", ylabel = "z [m]") p14 = Plots.plot(B_rim_end, z_centers, xlabel = "B_rim [-]", ylabel = "z [m]") @@ -340,7 +340,7 @@ function plot_timeheight_p3(nc_data_file, precip; output = "output") B_rim_plt = Array(ds.group["profiles"]["B_rim"]) #! format: off p1 = Plots.heatmap(t_plt, z_plt, q_tot_plt .* 1e3, title = "q_tot [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) - p2 = Plots.heatmap(t_plt, z_plt, q_liq_plt .* 1e3, title = "q_liq [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) + p2 = Plots.heatmap(t_plt, z_plt, (ρq_ice_plt .- ρq_liqonice_plt) .* 1e3, title = "ρq_ice (solid ice) [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) p3 = Plots.heatmap(t_plt, z_plt, ρq_ice_plt .* 1e3, title = "ρq_ice [g/m3]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) p4 = Plots.heatmap(t_plt, z_plt, q_rai_plt .* 1e3, title = "q_rai [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis) p5 = Plots.heatmap(t_plt, z_plt, q_vap_plt .* 1e3, title = "q_vap [g/kg]", xlabel = "time [s]", ylabel = "z [m]", color = :viridis)