Moisture budget decomposition

[nora-fahrenbach]

Hello everyone,

I wanted to ask if anyone here has used MetPy to perform a moisture budget decomposition before including an evaluation of the thermodynamic $\delta TH$, mean circulation dynamic $\delta MCD$ and transient eddy $\delta TE$ term (e.g. something similar in papers from Seager et al. 2010 and 2011.)

The equation is as follows:
$\rho_w g \delta (P-E) = \delta TH + \delta MCD + \delta TE - \delta S$
where
$\delta TH = - \int_{0}^{p_s} \nabla \cdot (\bar{u} [\delta \bar{q}]) dp$
$\delta MCD = - \int_{0}^{p_s} \nabla \cdot ([\delta \bar{u}] \bar{q}) dp$
$\delta TE = - \int_{0}^{p_s} \nabla \cdot \delta (\overline{u' q'}) dp$

with overbars indicating monthly means and primes indicating departures from the monthly mean.

I have CMIP6 data for precipitation, evaporation, 3D wind components u and v and 3D specific humidity q, with units time/plev/lat/lon or time/lat/lon respectively. I want to then perform a moisture budget decomposition for a certain timeslice (i.e. time will not be a variable anymore) to get the values for the thermodynamic, mean circulation dynamic and transient eddy term.

How would I e.g. code up the nabla of the terms and then the integral over pressure levels?
Any ideas how to approach this would be very much appreciated :-)