Hi there, Will.
Nothing much I can say about point from 2 to 5. However, I have already made something similar to your point 1, when it came to calculating AOD (Aerosol Optical Depth).
For that, you need individual layer thicknesses. To compute that, I used (staggered) geopotential layer interfaaces (N+1, staggered) to compute thickess (top-bottom). With that I cross multiplied it with extinction coefficients (both N vertical levels). Then summing up and obtaining AOD. So, in some way, this might be useful to what you are trying to achieve.
Xarray is partcicularly good at this task (dimensions, .diff() and cross products along given dimensions) and also allows to equip newly calculated products with metadata to keep track of each step.
Hope this helps.
Happy coding.
Marco