Update Utilities/mean_nemo to handle time-varying masks and make other improvements

[dcalve]

Summary of changes

Add support for time-varying masks in the input data (commit 260077d)

mean_nemo currently cannot handle masked data where the mask varies over time. Like many properties of the input files, it is assumed that the mask of the first time record in the first file is the same as that of all other records. This is not true when NEMO uses open ice shelf cavities, which can have a time-varying mask corresponding to the position of the ice shelf. In such cases, mean_nemo will not correctly exclude masked data from the average and fill values will appear in the result.

This will be the case for the UKESM 1.3 CMIP7 simulations, so annual mean diagnostics will not be correctly calculated (and therefore cannot be submitted) until mean_nemo is able to handle time-varying masks. This is the main change introduced by this development.

The new masked averaging algorithm follows that of numpy.ma.mean- masked points are excluded from the average, with the result being masked where there is no valid data.

Changes:

• Replace time counter scalars for the different data types (ntimes_i1, ntimes_sp, etc) with a single integer type (relying instead on automatic type conversion):
  • 3D and 4D arrays (ntimes_3d, ntimes_4d) for masked averaging
  • A scalar (ntimes) for unmasked averaging
• Implement a masked averaging algorithm
  • The numerator and denominator of the average are accumulated only where the data does not match the _FillValue
  • The numerator is divided by the denominator only where the latter is nonzero, otherwise the result is set to the _FillValue

Fix potential bug in thickness-weighted masked averaging (commit f307cfb)

When performing thickness-weighted averaging, mean_nemo assumes that the cell thickness mask is the same as that of the input data. While this will generally be the case, mean_nemo would not return a correct result if, for example, the cell thickness was masked but the input data was not. This set of changes, to account for differing masks in the input data and cell thickness when calculating masked averages, is therefore largely precautionary.

The average of input data $D$ weighted by cell thickness $C$ is $\frac{\sum{C D}}{\sum{C}}$., where $D$ and/or $C$ may be masked. In calculating the masked average, $\sum{C D}$ and $\sum{C}$ must both use the same mask- the union (bitwise .OR.) of both masks- meaning that $\sum{C}$ must be recalculated for each input data variable, rather than just once as is presently the case. This has a significant performance impact that is counteracted by other optimisations (see below).

Changes:

• Account for the cell thickness mask in the masked average algorithm
  • Either or both of the input data (l_inputdata_ismasked) and cell thickness data (l_cellthick_ismasked) may be masked
  • The mask used by the average is equivalent to the union (bitwise .OR.) of the cell thickness and input data masks, but for efficiency the masking array used in the code (l_mask_4d) is implemented as the intersection (bitwise .AND.) of where the data is not masked
  • The sum over time of the cell thickness is now performed for each thickness-weighted input data variable, rather than just once when calculating the cell thickness mean
  • The _FillValue value of the output data (outputdata_fill_value*) is set to that of the input data (inputdata_fill_value*) if it is masked, otherwise it is set to that of the cell thickness data (cellthick_fill_value*)
    • If the latter is used, we must copy the cell thickness _FillValue to the output file while it is in define mode
    • For this reason, the search for the cell thickness variable has been moved out of the main loop over variables (jv_loop) to before the call to nf90_enddef

Ensure that all STOP errors relating to NetCDF API calls print a human-readable message (commit 13dd91e)

For the most part, error handling of the calls to nf90_* routines only returns exit codes. I've added nf90_strerror calls to provide more informative error messages.

Update compilation scripts to include optimisation flags (commit bb89959)

The first two of the above changes are highly detrimental to the performance of mean_nemo, as they increase the number of calculations. Moving the most intensive calculations into WHERE statements also has an impact on performance due to less optimal vectorisation, although this impact is largely offset by the reduction in number of calculations (masked data is excluded). While this degradation in performance is undesirable, it occurs because some rather strict assumptions about the input data have been relaxed.

This performance degradation could easily be countered by implementing parallelism in the code (there is no MPI or OpenMP support). However, this proved unnecessary as, for some reason, the compilation scripts (compiler.xc40 and compiler.xc40_login) build mean_nemo using only the -O0 flag. Using a reasonable set of optimisation flags for the Cray compiler (-O2 -Ovector1 -hfp0 -hflex_mp=intolerant) counteracts the performance loss of the above changes, resulting in very little difference in run time compared to the original code. These compiler flags have been added to both scripts.

Additionally, both scripts have been updated to include several levels of optimisation (MODE = prod, dev, debug) and to load the correct libraries for the new Cray EX machines.

Remove zeroing of arrays used for storing input data (commit a3896c1)

This was a minor optimisation of the following code patterns:

ALLOCATE(inputdata_1d_i2(outdimlens(dimids(1))))
inputdata_1d_i2(:)=0
iostat = nf90_get_var( ncid, jv, inputdata_1d_i2, start, indimlens)

The zeroing is only necessary if the data read in by nf90_get_var is spatially smaller than inputdata_1d_i2, which should never happen as it is assumed that data has the same non-temporal shape in all input files.

Add basic timing calipers (commit 32d8663)

I used these during development and thought I would include them for future use.

Add script for generating a simple test dataset (commit f8164a3)

The standard testing in place for mean_nemo is limited. The rose-stem suite only checks that mean_nemo is able to successfully process a specific (and old) type of NEMO diagnostic file, but does not check whether the contents are correct. Other than this, a developer might perform regression testing using different sets of model data to check that the mean_nemo output is unchanged, but in practice this tests only a limited part of its functionality. For example, most NEMO diagnostic data is output as single- or double-precision floats, meaning that integer data types are not tested.

In the absence of unit tests, I have added a script (simple_test_data.py) that generates a small set of test data for mean_nemo. Two types of file are generated- input data (intended as input to mean_nemo) and expected results (the results expected from running mean_nemo with the input data)- containing several variables with different data precisions, shapes etc to test different parts of the code. The mean_nemo output can be compared with the expected results using e.g. nccmp:

nccmp -F -f -c 1 -d -m <mean_nemo_output_file> <expected_results_file>

Correctness testing using this simple dataset should complement 'production' testing using actual NEMO diagnostic datasets generated by our suites.