Check that grib2bgrasc works OK with MIR-interpolated GRIB files

[Ian Culverwell]

Joe Nielsen at DMI draws attention to a change in the Product Generation Software to be used by ECMWF from the end of February 2019, as described at ​https://confluence.ecmwf.int/display/UDOC/New+Product+Generation+software. In particular, there will be a change in the longitude encoding of GRIB1 files. ERA-interim files are in GRIB1 format, and DMI use these, so we need to check that grib2bgrasc works OK on them.

[Ian Culverwell]

ei_25_2018020112_000.grib

[Ian Culverwell]

od_25_2018020112_000.grib

[Ian Culverwell]

ei_25_2018020112_000_MIR.grib

[Ian Culverwell]

ei_25_2018020112_000_lon_eq_0.log

[Ian Culverwell]

ei_25_2018020112_000_MIR_lon_eq_0.log

[Ian Culverwell]

ei_25_2018020112_000_lon_eq_180.log

[Ian Culverwell]

ei_25_2018020112_000_MIR_lon_eq_180.log

[Ian Culverwell]

Hans Gleisner has provided two existing GRIB files, for comparison. ei_25_2018020112_000.grib, a 2.5x2.5 deg ERA-I file, has the following parameters:

ei_25_2018020112_000.grib
  editionNumber = 1;
  longitudeOfFirstGridPointInDegrees = -180;
  longitudeOfLastGridPointInDegrees = 177.5;
  iDirectionIncrementInDegrees = 2.5;
  iScansNegatively = 0;
  latitudeOfFirstGridPointInDegrees = 90;
  latitudeOfLastGridPointInDegrees = -90;
  jDirectionIncrementInDegrees = 2.5;
  jScansPositively = 0;
Latitude, Longitude, Value, shortName
   90.000 -180.000 1.1549090385e+01 lnsp
   90.000 -177.500 1.1549090385e+01 lnsp
   90.000 -175.000 1.1549090385e+01 lnsp
   90.000 -172.500 1.1549090385e+01 lnsp
...
  -90.000  170.000 1.1121310234e+01 lnsp
  -90.000  172.500 1.1121310234e+01 lnsp
  -90.000  175.000 1.1121310234e+01 lnsp
  -90.000  177.500 1.1121310234e+01 lnsp

I downloaded the equivalent file, ei_25_2018020112_000_MIR.grib, from the ECMWF web interface. This appears to be produced by the new product generation software, because its metadata say

ei_25_2018020112_000_MIR.grib
  editionNumber = 1;
  longitudeOfFirstGridPointInDegrees = 0;
  longitudeOfLastGridPointInDegrees = 357.5;
  iDirectionIncrementInDegrees = 2.5;
  iScansNegatively = 0;
  latitudeOfFirstGridPointInDegrees = 90;
  latitudeOfLastGridPointInDegrees = -90;
  jDirectionIncrementInDegrees = 2.5;
  jScansPositively = 0;
Latitude, Longitude, Value, shortName
   90.000    0.000 1.1549090385e+01 lnsp
   90.000    2.500 1.1549090385e+01 lnsp
   90.000    5.000 1.1549090385e+01 lnsp
   90.000    7.500 1.1549090385e+01 lnsp
...
  -90.000  350.000 1.1121310234e+01 lnsp
  -90.000  352.500 1.1121310234e+01 lnsp
  -90.000  355.000 1.1121310234e+01 lnsp
  -90.000  357.500 1.1121310234e+01 lnsp

The important point is that grib2bgrasc gives the same results, within rounding error, with both files. See the attached log files from the two runs, with input lons of 0 and 180. These also agree with what you get by examining the grib file 'manually'.

[Ian Culverwell]

We can check these longitudes (which are the results of running grib_get_data $IFILE -p shortName -w shortName='lnsp' through the files) by plotting the surface geopotential.

ei_25_2018020112_000.grib:

ei_25_2018020112_000_MIR.grib:

[Ian Culverwell]

ei_25_2018020112_000.png

[Ian Culverwell]

ei_25_2018020112_000_MIR.png

[Ian Culverwell]

Here are the surface geopotentials for Hans's quartet of files: Note the 'ringing' in the EI-25 data - probably Gibbs' oscillations resulting from interpolation from the underlying spectral model onto a lat-lon grid that is so coarse that the geopotential field has sharp discontinuities. This ringing has been almost completely silenced in the other files. (The OD files come from the operational data stream from ECMWF to DMI, whereas the EI data are obtained from MARS retrievals. And the effect should be much smaller at 1 deg anyway. So that might explain it.)

The surface geopotentials in my retrievals of the analagous files are The big difference is the OD-25 run, which rings as noisily as EI-25 in my case. This makes it difficult to compare grib2bgrasc on these files (although the temperatures and specific humidities aren't too bad.)

It turns out that the geops in the two EI-10 files match exactly as far as one can tell (because of the shift of the longitudes you can't difference the fields), and grib2bgrasc therefore gives very similar results (not shown).

There is some difference in the geopotentials and the other fields at OD-10, but they are close enough to test grib2bgrasc. The results follow.

[Ian Culverwell]

all_xx_2018020112_000_geop.png

[Ian Culverwell]

all_xx_2018020112_000_MIR_geop.png

[Ian Culverwell]

(Incidentally, I tried transplanting the geop field from Hans's od_10_28020112_000.grib into my od_10_28020112_000_MIR.grib, to get a closer comparison. (grib_copy is the trick.) This turns out to be a bad idea, because Hans's geopotentials are on a different grid (-180 < lons < 180) to the other fields in my file (0 < lons < 360), and grib2bgrasc assumes that all fields are on the same grid. We therefore end up comparing geops at the same point and temps at different points, 180 deg apart, or vice-versa. Either way it's pants.)

[Ian Culverwell]

od_10_2018020112_000_lon_eq_0.log

[Ian Culverwell]

od_10_2018020112_000_MIR_lon_eq_0.log

[Ian Culverwell]

od_10_2018020112_000_lon_eq_180.log

[Ian Culverwell]

od_10_2018020112_000_MIR_lon_eq_180.log

[Ian Culverwell]

The metadata of Hans's 1 deg operational data file, od_10_2018020112_000.grib, say

od_10_2018020112_000.grib
  editionNumber = 2;
  longitudeOfFirstGridPointInDegrees = 180;
  longitudeOfLastGridPointInDegrees = 179;
  iDirectionIncrementInDegrees = 1;
  iScansNegatively = 0;
  latitudeOfFirstGridPointInDegrees = 90;
  latitudeOfLastGridPointInDegrees = -90;
  jDirectionIncrementInDegrees = 1;
  jScansPositively = 0;
Latitude, Longitude, Value, shortName
   90.000 -180.000 1.1545543671e+01 lnsp
   90.000 -179.000 1.1545543671e+01 lnsp
   90.000 -178.000 1.1545543671e+01 lnsp
   90.000 -177.000 1.1545543671e+01 lnsp
...
  -90.000  176.000 1.1146862030e+01 lnsp
  -90.000  177.000 1.1146862030e+01 lnsp
  -90.000  178.000 1.1146862030e+01 lnsp
  -90.000  179.000 1.1146862030e+01 lnsp

and this is borne out by the geop pictures above, which start at lon=(+/-)180 deg.

The metadata in the equivalent file downloaded from the ECMWF web api, od_10_2018020112_000_MIR.grib, say

od_10_2018020112_000_MIR.grib
  editionNumber = 2;
  longitudeOfFirstGridPointInDegrees = 0;
  longitudeOfLastGridPointInDegrees = 359;
  iDirectionIncrementInDegrees = 1;
  iScansNegatively = 0;
  latitudeOfFirstGridPointInDegrees = 90;
  latitudeOfLastGridPointInDegrees = -90;
  jDirectionIncrementInDegrees = 1;
  jScansPositively = 0;
Latitude, Longitude, Value, shortName
   90.000    0.000 1.1545681000e+01 lnsp
   90.000    1.000 1.1545681000e+01 lnsp
   90.000    2.000 1.1545681000e+01 lnsp
   90.000    3.000 1.1545681000e+01 lnsp
...
  -90.000  356.000 1.1142818451e+01 lnsp
  -90.000  357.000 1.1142818451e+01 lnsp
  -90.000  358.000 1.1142818451e+01 lnsp
  -90.000  359.000 1.1142818451e+01 lnsp

and the geopotentials plotted above do indeed go from 0 to 360.

[Ian Culverwell]

The temps and shums in the two files are close but ot identical, perhaps because in Hans's file these fields undergo the same processing that silences the ringing in the geop field. But they are close enough to test that grib2bgrasc works the same on both files. See, for example, the od_10_*.log log files, which print the extracted profile data as well as some diagnostics. The results are certainly a lot closer than they are in the inadvertent sensitivity study I made above by using a shifted geop field (also attached as od_10_2018020112_000_MIR_transplanted_geop_ lon_eq_0.txt).

[Ian Culverwell]

od_10_2018020112_000_MIR_transplanted_geop_lon_eq_0.log

[Ian Culverwell]

Summary

• grib2bgrasc has worked OK on all the files I've thrown at it.

• A shift in longitudes from [-180, 180) to [0, 360), as expected with the new product generation software, makes no difference to its behaviour.

• We therefore foresee no difficulties when MIR becomes the default MARS client at ECMWF on 29 January 2019.

• There is evidently some unrecorded/unexpected processing in DMI's stream of operational data from ECMWF, because the geopotentials are much less oscillatory than their counterparts from the ERA-I analyses, as retrieved from MARS.

• grib2bgrasc assumes all fields in a grib file are on the same grid. This should be borne in mind.

[Johannes K. Nielsen]

Memo on the ripple issue: Marc Schwaerz noted that MARS allows retrieval of geopotential either as a surface parameter or as a model level (ML) parameter. In the retrieval of ERA-I data we use "LEVT = ML" which is the correct choice, so the ripples issue not connected to that.

[Johannes K. Nielsen]

The source of ringing has now been identified. ​https://confluence.ecmwf.int/display/UDOC/Post-processing+keywords#Post-processingkeywords-grid

it appears that the default MARS behaviour is to truncate the spectral representation at sufficient high resolution for the requested grid before interpolation. I.e. all the higher order modes are left out of the interpolation.

What we need to do is to insert the keyword line: "resol=av" in the MARS call ("AV" means archived value in all keywords), like.

RETRIEVE,

TARGET='myFile', USE=INFREQUENT, STREAM=OPER, CLASS=EI, EXPVER=1, GRID=2.5/2.5, RESOL = AV, AREA=90.0/-180.0/-90.0/177.5, DATE=$DATE, TIME=$TIME, TYPE=AN, STEP=0, LEVT=ML, LEVE=1, PARA=Z

I tried to do that, and the ringing is gone. Completely.

I suggest that we do not change anything in the ERA-I data set at DMI and that we let ICDR1 run to the end of ERA-I with ringing. Then for the next episode, ERA5, we put RESOL = AV, and we should be safe.

jkn

[Johannes K. Nielsen]

As of the 20th March 2019 we switched to use MIR operationally at DMI, with no observed anomalies, possible bias jumps will be monitored in operational reports, so you can close this ticket seen from DMI point of view.

[Ian Culverwell]

I cannot reproduce this. I access the MARS archive through the web interface, not by submitting a retrieval script, and the web interface doesn't allow me to specify RESOL = AV as well as the grid spacing. It does give me the option of retrieving something called Default (as archived), but that just gives me 256 x 256 complex numbers: spherical harmonic components, presumably.

Since DMI are happy with the solution above, I'll close the ticket.

NB: There is a python-API to the archive, which I could try in future if necessary. This involves downloading and building some libraries, but it does allow users to write their own retrieval scripts, like DMI use.

[Ian Culverwell]

ei_25_2018020112_000_on_archived_grid.grib

[Ian Culverwell]

ei_25_2018020112_000_on_archived_grid.nc

[Ian Culverwell]

ei_25_2018020112_000_on_2.5deg_grid.grib

[Ian Culverwell]

(Added ei_25_2018020112_000_on_archived_grid.grib = ERA-I geopotential on 'archived grid', as retrieved from MARS using the web interface. This contains 256x256 complex numbers. Also added ei_25_2018020112_000_on_archived_grid.nc = the result of converting this file to netCDF format using xconv. (You can't ask for data in netCDF format in the archived format directly from MARS.) It has just concatenated the real and imaginary parts of the other file - not much use. I also attach ei_25_2018020112_000_on_2.5deg_grid.grib = what you get when you ask for geop on 2.5 deg grid, because - hilariously - this is different to what I got back in January. Need to stop pulling at these threads.)