@@ -9,6 +9,7 @@ ENV["GKSwstype"] = "nul"
99import ClimaCorePlots, Plots
1010Plots. GRBackend ()
1111using CairoMakie
12+ using Dates
1213
1314function plot_initial_profiles_comparison (KM; sdm_case = " dry" )
1415 sdm_data = load_sdm_data (sdm_case)
@@ -313,6 +314,99 @@ function plot_animation(nc_data_file; output = "output")
313314 Plots. mp4 (anim, joinpath (path, " animation.mp4" ), fps = 10 )
314315end
315316
317+ function plot_profiles_in_time (nc_data_file; output = " output" , n = 10 )
318+
319+ path = joinpath (@__DIR__ , output)
320+ mkpath (path)
321+
322+ ds = NC. NCDataset (joinpath (@__DIR__ , nc_data_file))
323+ profs = ds. group[" profiles" ]
324+
325+ z_plt = Array (profs[" zc" ])
326+ t_plt = Array (profs[" t" ])
327+ nt = length (t_plt)
328+
329+ # Select n evenly distributed time indices
330+ Δi = round (Int, (nt - 2 ) / n)
331+ time_indices = [1 ; Δi: Δi: (nt - Δi); nt] # explicitly include first and last indices
332+ ni = length (time_indices)
333+ selected_times = t_plt[time_indices]
334+
335+ kg_to_g = 1e3
336+ m⁻³_to_cm⁻³ = 1e-6
337+ q_tot = profs[" q_tot" ][:, time_indices] .* kg_to_g
338+ q_liq = profs[" q_liq" ][:, time_indices] .* kg_to_g
339+ q_ice = profs[" q_ice" ][:, time_indices] .* kg_to_g
340+ q_rai = profs[" q_rai" ][:, time_indices] .* kg_to_g
341+ q_sno = profs[" q_sno" ][:, time_indices] .* kg_to_g
342+ N_liq = profs[" N_liq" ][:, time_indices] .* m⁻³_to_cm⁻³
343+ N_rai = profs[" N_rai" ][:, time_indices] .* m⁻³_to_cm⁻³
344+ SN_liq_act = profs[" SN_liq_act" ][:, time_indices] .* m⁻³_to_cm⁻³
345+
346+ close (ds)
347+
348+ # Create colormap for time progression
349+ # colormap = :viridis
350+ colormap = cgrad (:darkrainbow , (0 : ni) / ni, categorical = true )
351+ # colors = get(colorschemes[colormap], range(0, 1, length = n))
352+ ymax = z_plt[end ] + (z_plt[end ] - z_plt[end - 1 ]) # extend y-axis a bit beyond the last point
353+ lims (x) = iszero (x) ? (nothing , (0 , ymax)) : ((0 , 1.1 * maximum (x)), (0 , ymax))
354+
355+ fig = Figure (size = (2000 , 1500 ))
356+
357+ # Create axes
358+ ax_q_tot = Axis (fig[1 , 1 ]; xlabel = " q_tot [g/kg]" , limits = lims (q_tot), ylabel = " z [m]" )
359+ ax_q_liq = Axis (fig[1 , 2 ]; xlabel = " q_liq [g/kg]" , limits = lims (q_liq))
360+ ax_N_liq = Axis (fig[1 , 3 ]; xlabel = " N_liq [1/cm³]" , limits = lims (N_liq))
361+ ax_q_rai = Axis (fig[2 , 1 ]; xlabel = " q_rai [g/kg]" , limits = lims (q_rai), ylabel = " z [m]" )
362+ ax_N_rai = Axis (fig[2 , 2 ]; xlabel = " N_rai [1/cm³]" , limits = lims (N_rai))
363+ ax_q_ice = Axis (fig[2 , 3 ]; xlabel = " q_ice [g/kg]" , limits = lims (q_ice))
364+ ax_q_sno = Axis (fig[3 , 1 ]; xlabel = " q_sno [g/kg]" , limits = lims (q_sno), ylabel = " z [m]" )
365+ ax_SN_liq_act = Axis (fig[3 , 2 ]; xlabel = " SN_liq_act [1/cm³]" , limits = lims (SN_liq_act))
366+
367+ # Link y-axes
368+ axs = [ax_q_tot, ax_q_liq, ax_N_liq, ax_q_rai, ax_N_rai, ax_q_ice, ax_q_sno]
369+ linkyaxes! (axs... )
370+
371+ # Plot lines for each selected time
372+ args = (; colormap, colorrange = (0 , 1 ))
373+ for (i, t) in enumerate (selected_times)
374+ color = (i - 0.5 ) / ni # Normalize color for colormap
375+ lines! (ax_q_tot, q_tot[:, i], z_plt; args... , color)
376+ lines! (ax_q_liq, q_liq[:, i], z_plt; args... , color)
377+ lines! (ax_N_liq, N_liq[:, i], z_plt; args... , color)
378+ lines! (ax_q_rai, q_rai[:, i], z_plt; args... , color)
379+ lines! (ax_N_rai, N_rai[:, i], z_plt; args... , color)
380+ lines! (ax_q_ice, q_ice[:, i], z_plt; args... , color)
381+ lines! (ax_q_sno, q_sno[:, i], z_plt; args... , color)
382+ lines! (ax_SN_liq_act, SN_liq_act[:, i], z_plt; args... , color)
383+ end
384+
385+ # Add colorbar
386+ tloc = ((1 : ni) .- 0.5 ) / ni # locations offset by 0.5 to center label on the color
387+ tlab = begin
388+ # Canonicalize time labels
389+ can_times = @. canonicalize (Second (selected_times))
390+ # shorten time labels
391+ map (can_times) do t
392+ s = string (t)
393+ s == " empty period" && return " 0 s"
394+ s = replace (s, " seconds" => " s" , " minutes" => " m" , " hours" => " h" , " days" => " d" )
395+ s = replace (s, " second" => " s" , " minute" => " m" , " hour" => " h" , " day" => " d" )
396+ s = replace (s, " milli" => " m" , " micro" => " μ" , " nano" => " n" )
397+ end
398+ end
399+ Colorbar (fig[:, 4 ]; colormap, colorrange = (0 , 1 ), width = 20 , ticks = (tloc, tlab))
400+
401+ # Add title
402+ Label (fig[0 , :], " Vertical Profiles at Multiple Time Steps" , fontsize = 24 )
403+
404+ # Save figure
405+ save (joinpath (path, " profiles_multitime.png" ), fig)
406+
407+ return fig
408+ end
409+
316410function plot_timeheight_p3 (nc_data_file, precip; output = " output" )
317411 path = joinpath (@__DIR__ , output)
318412 mkpath (path)
0 commit comments