"""
Copyright (c) 2024 MPI-M, Clara Bayley
----- Microphysics Test Cases -----
File: perform_0dparcel_test_case.py
Project: test_case_0dparcel
Created Date: Wednesday 28th February 2024
Author: Clara Bayley (CB)
Additional Contributors:
-----
Last Modified: Wednesday 4th June 2025
Modified By: CB
-----
License: BSD 3-Clause "New" or "Revised" License
https://opensource.org/licenses/BSD-3-Clause
-----
File Description:
interface called by a test to run the 0-D parcel and then plot the results.
"""
from pathlib import Path
import matplotlib.pyplot as plt
from .run_0dparcel import run_0dparcel
from libs.thermo import formulae
from libs.utility_functions import plot_utilities
[docs]
def plot_0dparcel_thermodynamics(out, binpath, run_name):
"""Plot thermodynamic variables for a 0-D parcel model and save the plots.
This function plots the pressure, density, temperature, and potential temperature(s)
of a run of the 0-D parcel model as a function of time and then saves the plots as a PNG image.
Args:
out (OutputThermodynamics):
OutputThermodynamics object containing the thermodynamic data.
binpath (str):
Path to the directory where the plots will be saved.
run_name (str):
Name of the test run to use in naming saved image.
Raises:
AssertionError: If the specified binpath does not exist or if run_name is empty.
Returns:
None
"""
assert Path(binpath).exists()
assert run_name
print("plotting " + run_name + " and saving plots in: " + str(binpath))
fig, axs = plt.subplots(nrows=2, ncols=3, sharex=True, figsize=(12, 5))
figname = run_name + "_thermodynamics.png"
axs = axs.flatten()
time = out.time.values
plot_utilities.plot_thermodynamics_output_timeseries(axs[0], out, "press")
plot_utilities.plot_thermodynamics_output_timeseries(axs[0].twinx(), out, "rho")
plot_utilities.plot_thermodynamics_output_timeseries(axs[1], out, "temp")
plot_thetas_on_axis(axs[2], time, out.temp, out.press, out.qvap)
plot_mse_on_axis(axs[3], time, out.temp, out.qvap)
plot_relh_on_axis(axs[4], time, out.temp, out.press, out.qvap)
axs[5].remove()
for ax in axs:
ax.set_xlabel(out.time.name + " /" + out.time.units)
fig.tight_layout()
plot_utilities.save_figure(fig, binpath, figname)
[docs]
def plot_0dparcel_massmix_ratios(out, binpath, run_name):
"""Plot mass mixing ratios for a 0-D parcel model and save the plots.
This function plots the mass mixing ratios of water vapor, cloud liquid, cloud ice,
rain, snow, and graupel for a run of the 0-D parcel model as a function of time and then
saves the plots as a PNG image.
Args:
out (OutputThermodynamics):
OutputThermodynamics object containing the mass mixing ratio data.
binpath (str):
Path to the directory where the plots will be saved.
run_name (str):
Name of the test run to use in naming saved image.
Raises:
AssertionError: If the specified binpath does not exist or if run_name is empty.
Returns:
None
"""
assert Path(binpath).exists(), "The specified binpath does not exist."
assert run_name, "The run_name cannot be empty."
print("plotting " + run_name + " and saving plots in: " + str(binpath))
fig, axs = plt.subplots(nrows=2, ncols=4, sharex=True, figsize=(12, 5))
figname = run_name + "_massmix_ratios.png"
axs = axs.flatten()
plot_utilities.plot_thermodynamics_output_timeseries(axs[0], out, "qvap")
plot_utilities.plot_thermodynamics_output_timeseries(axs[1], out, "qcond")
plot_utilities.plot_thermodynamics_output_timeseries(axs[2], out, "qrain")
plot_utilities.plot_thermodynamics_output_timeseries(axs[4], out, "qice")
plot_utilities.plot_thermodynamics_output_timeseries(axs[5], out, "qsnow")
plot_utilities.plot_thermodynamics_output_timeseries(axs[6], out, "qgrau")
qtot_warm = out.qvap.values + out.qcond.values + out.qrain.values
axs[3].plot(out.time.values, qtot_warm)
axs[3].set_ylabel("q$_{v}$ + q$_{cond}$ + q$_{rain}$ / kg/kg")
qtot = qtot_warm + out.qice.values + out.qsnow.values + out.qgrau.values
axs[7].plot(out.time.values, qtot)
axs[7].set_ylabel("q$_{tot}$ / kg/kg")
for ax in axs:
ax.set_xlabel(out.time.name + " /" + out.time.units)
fig.tight_layout()
plot_utilities.save_figure(fig, binpath, figname)
[docs]
def plot_mse_on_axis(ax, time, temp, qvap):
"""Plot moist static energy (MSE) on a specified axis.
This function calculates and plots MSE against time on a specified axis.
Args:
ax (matplotlib.axes.Axes): The (x-y) axis on which to plot the MSE.
time (array-like): Time values (x axis).
temp (OutputVariable): Temperature variable.
qvap (OutputVariable): Mass mixing ratio of water vapour variable.
Returns:
None
"""
mse = formulae.moist_static_energy(temp.values, qvap.values)
ax.plot(time, mse)
ax.set_ylabel("moist static energy /kJ/kg")
[docs]
def plot_thetas_on_axis(ax, time, temp, press, qvap):
"""Plot potential temperature(s) on a specified axis.
This function calculates and plots potential temperature(s) against time on a specified axis.
Args:
ax (matplotlib.axes.Axes): The (x-y) axis on which to plot the potential temperature(s).
time (array-like): Time values (x axis).
temp (OutputVariable): Temperature variable.
press (OutputVariable): Pressure variable.
qvap (OutputVariable): Mass mixing ratio of water vapour variable.
Returns:
None
"""
theta_dry = formulae.dry_potential_temperature(temp.values, press.values)
theta_moist = formulae.moist_equiv_potential_temperature(
temp.values, press.values, qvap.values
)
ax.plot(time, theta_dry, label="dry")
ax.plot(time, theta_moist, label="moist equiv.")
ax.set_ylim(theta_dry[0] - 50, theta_dry[0] + 50)
ax.legend()
ax.set_ylabel("potential temperature /" + temp.units)
[docs]
def plot_relh_on_axis(ax, time, temp, press, qvap):
"""Plot relative humidity (relh) on a specified axis.
This function calculates and plots relh against time on a specified axis.
Args:
ax (matplotlib.axes.Axes): The (x-y) axis on which to plot relh.
time (array-like): Time values (x axis).
temp (OutputVariable): Temperature variable.
press (OutputVariable): Pressure variable.
qvap (OutputVariable): Mass mixing ratio of water vapour variable.
Returns:
None
"""
from metpy.units import units
from metpy import calc
relh = calc.relative_humidity_from_mixing_ratio(
press.values * units.Pa, temp.values * units.kelvin, qvap.values
)
ax.plot(time, relh * 100)
ax.hlines(100, time[0], time[-1], linestyles="--", linewidth=0.8, color="grey")
ax.set_ylabel("relative humidity /%")