Age of air diagnostic porting

[jwarner8]

The age of air diagnostic captures how old air is at each grid-point as a function of forecast lead-time. It provides a measure of dilution within a limited area model. It requires 3D fields at each model output on standard pressure levels. The diagnostic will be applied within the BoM ACCESS models to compare their new wider domain setup to the existing multiple LAM setup (JamesW to visit BoM under the UNITE programme in Feb2023). This issue will capture development, implementation of this diagnostic.

[jwarner8]

Progress: have run through concept and code with James F. Main work identified; returning multiple cubes (cubelist) as part of filter operator, and standard regridding ability (probably iris based, also MET an option). Discussion with Sylvia/James about parallelisation (probably parallelise pressure level in CYLC, and run on all lead-times if cheap enough), and use multiprocessing in python to handle iteration over all grid points (which are independent). Next steps: JamesW to tidy up code, dependencies, write some documentation, create some test data

[jwarner8]

To think further: how to compute grid spacing on grid/regridding, as currently assumes 111km = 1 deg longitude. Not applicable outside tropics. Development of test cases: 1. linear (laminar) flow, 2. linear rising flow with background profile of increasing winds with height. 3. some idealised equatorial Rossby wave flow e.g. this

[jwarner8]

Further testing suggests 0.5 degree is suitable as a standard grid for running the AOA code through currently (without optimisations). 1 degree is quite coarse and we lose a lot of important local detail. Smoothing has been implemented using a Gaussian filter using Scipy for vertical velocity, as we are using instantaneous winds if we sample a larger MCS/storm complex, we might get updraft velocities of 1m/s on 0.5 unsmoothed grid, which integrated over 3hours would reach top of troposphere/surface. So smoothing reduces these grid-points where the back trajectory suddently spans entire atmospheric column. Some further testing needed for optimal sigma arguments for Gaussian smoothing. Test cases have started being created (linear flow, SW flow, and two layer flow with a grid-point set at 1/-1m/s vertical velocity. Visuals for these test cases will be produced to qualitatively evaluate them, and then once established correctly working, correct output will be used to assert value.

[jwarner8]

1. Code now uses the haversine formulae to calculate grid spacing rather than 111km=1deg approximation, which breaks down quickly away from equator. More accurate approaches using geopy but avoiding additional package dependencies.
2. Code now using parallisation across longitude dimension (multiprocessing and partial libs). Scales well across 12 cores, depends on domain size and lead-time length (IO major time drain).
3. Supports cyclic channel.
4. Unit tests mostly done, need to do some visualisations of test cases too.

[jwarner8]

the age of air code scales nicely with multiprocessing. I iterate over longitude points as the pooled processes. Scalability depends on how long the forecast is (more back trajectories), and of course shape of domain (most are square-ish). Beyond about 12 cores performance tends to tail off, guessing due to IO or moving data around/slicing etc. Though we probably wouldn't want to request more than that anyway given queuing. Mem usage scales quite well. We will take pressure level as an argument from the CYLC workflow, rather than run the code on all available pressure levels. Currently working on regrid operator and also unit tests. The diagnostic plays nicely with ACCESS-A data at BoM, on the NCI platform.