libaccessom2 doesn't deal with netcdf unpacking properly (we think...) in ERA-5 runs

[rmholmes]

As described on the ERA-5 forcing issue I think libaccessom2 may have an issue dealing with netcdf unpacking across file boundaries. I'll summarize the problem here.

The problem occurs when transitioning between two months (the ERA-5 forcing is stored in monthly files), best demonstrated by plotting daily minimum wind stress at 92W, 0N from a 1deg_era5_iaf run spanning 1980-01-01 to 1980-05-01:

There is a large burst of negative wind stress in the first day of April in the "raw" run (this causes all sorts of crazy stuff...). The add_offset netcdf packing value in the ERA-5 10m zonal winds file is particularly anomalous for March of this year (listed below per month of the files in /g/data/rt52/era5/single-levels/reanalysis/10u/1980/

                u10:add_offset = -3.54318567240813 ;
                u10:add_offset = 0.856332909292688 ;
                u10:add_offset = -32.1081480318141 ;
                u10:add_offset = -0.761652898754254 ;
                u10:add_offset = -0.10650583633675 ;
                u10:add_offset = -2.55211599669929 ;
...

If I change the netcdf packing values in the single March 1980 10m winds file (using the below python) and rerun, then I remove the burst of wind stress ("Altered packing" run above). This confirms to me that it is a packing issue.

file_in = '/g/data/rt52/era5/single-levels/reanalysis/10u/1980/10u_era5_oper_sfc_19800301-19800331.nc'
file_out = '/g/data/e14/rmh561/access-om2/input/ERA-5/IAF/10u/1980/10u_era5_oper_sfc_19800301-19800331.nc'
DS = xr.open_dataset(file_in)
encoding = {}
scale = 0.000966930321007164 # Apr 1980 value
offset = -0.761652898754254 # Apr 1980 value
encoding['u10'] = {'scale_factor': scale, 'add_offset': offset, 'dtype': 'int16'}
DS.to_netcdf(file_out,encoding=encoding)

Yes, the packing in the ERA-5 files is weird. But in any case, libaccessom2 should be able to deal with the variable packing. Xarray in python can, as shown by this plot of the time series of 10m zonal wind at the same point from the original file:

I've had a quick look through the code and am none the wiser. As @aekiss said, the netcdf unpacking seems to be handled by the netcdf library, so I don't understand how there can be a problem. Clearly it only affects the times between months when an interpolation has to be done. The rest of the month is fine.

[rmholmes]

Note: It doesn't look like this error will have affected the JRA-55 forced runs because these don't seem to use netcdf packing (no add_offset or scale_factor in the forcing files).

[access-hive-bot]

This issue has been mentioned on ACCESS Hive Community Forum. There might be relevant details there:

https://forum.access-hive.org.au/t/help-from-access-nri-on-libaccessom2-netcdf-packing-issue-for-cosima-era-5-runs/379/1

[aidanheerdegen]

Where is your experiment directory located @rmholmes?

[rmholmes]

/home/561/rmh561/access-om2/1deg_era5_iaf for the run with the issue and /home/561/rmh561/access-om2/1deg_era5_iaf_perror for the altered packing run.

Note this comment about potentially dropping this project, although I'd still love to know what the issue is!

[aidanheerdegen]

Note COSIMA/access-om2#242 (comment), although I'd still love to know what the issue is!

Even when this is working for MOM6 that is likely through CMEPS or some other nupoc framework, so this problem will still need to be fixed.

For future reference, the code for this in a separate branch of libaccessom2. The specific lines where the scale_factor and offset are applied are here:

https://github.com/COSIMA/libaccessom2/blob/242-era5-support/libutil/src/util.F90#L96-L103

I'd start with adding some logging to make sure those values are as expected.

It is strange that the first time it is read the values are bad, then they're ok.

The other place to look is this ncvar_read_data routine, which has the caching Nic added

https://github.com/COSIMA/libaccessom2/blob/242-era5-support/libforcing/src/ncvar.F90#L303-L343

An off-by-one index error is a candidate, but I can't see it having quite the catastrophic effects seen in the model.

[aekiss]

I've had a close look at subroutine read_data
https://github.com/COSIMA/libaccessom2/blob/70e2c46800fd9/libutil/src/util.F90#L46-L107
and ncvar_read_data
https://github.com/COSIMA/libaccessom2/blob/242-era5-support/libforcing/src/ncvar.F90#L303-L343
and didn't see anything amiss. But maybe somebody with sharper eyes would do better.

The order of operations in unpacking the data is also correct - see https://docs.unidata.ucar.edu/netcdf-c/current/attribute_conventions.html

[rmholmes]

I started adding some logging the other day but didn't get far. Hopefully I can find some time this week.

[russfiedler]

It looked good to me too. I think there may need to be some testing with artificial fields to see what's happening.

[rmholmes]

Please go ahead and do some testing if you can. I'm sure you can both make progress more quickly than I can.

[aekiss]

silly question - will these

libaccessom2/libutil/src/util.F90

Line 96 in 70e2c46

status = nf90_get_att(ncid, varid, "scale_factor", scale_factor)

libaccessom2/libutil/src/util.F90

Line 100 in 70e2c46

status = nf90_get_att(ncid, varid, "add_offset", offset)

return double precision? In C and f77 there are nc_get_att_double and nf_get_att_double but AFAICS nf90_get_att_double does not exist, so I guess types are handled automatically by nf90_get_att ...?

[aekiss]

so it would seem https://docs.unidata.ucar.edu/netcdf-fortran/current/f90-attributes.html#f90-get-attributes-values-nf90_get_att

[russfiedler]

My worry is that there is a chance that the compiler could perform the operations out of order since operations are on the inside of the if/endif constructs and it might not understand. The operation should be performed as dataout=scale_factor * dataout + offset anyway and not separately. It's also more efficient to do it in one go.
I also detest the fact that scale_factor and offset are never initialised to 1 and 0 respectively.

scale_factor = 1
offset = 0
status1=nf90_get_att(,,,scale_factor)
status2=nf90_get_att(,,,offset)
if(status1==nf90_noerr .or. status2 ==noerr) dataout=scale_factor*dataout + offset

[aekiss]

Ah yes that looks a lot safer and more efficient

[rmholmes]

I'll give it a go.

[rmholmes]

No luck unfortunately:

I'm putting in a PR anyway.