Accelerate grib2bgrasc extraction

[Ian Culverwell]

Axel von Engeln (EUMETSAT) has reported that it takes ~10s to extract a profile from a GRIB file using grib2bgrasc.

The code that packs the 1D GRIB messages containing the background fields into 2D or 3D arrays could probably be speeded up by means for the RESHAPE function. Probably the easiest thing to try first.

I don't think it would be difficult to process lots of profiles from the same grib file at the same time. The difficulty would be restructuring the output namelist, since presumably users would want to convert these profiles to a ROPP multifile at some time. But it could probably be done.

[Ian Culverwell]

The whole thing could be documented a bit better, too. Eg it could be explained that the times/dates in the nml and ROPP files are analysis/data times, and that one needs to add the forecast range to get the validity times.

[Ian Culverwell]

Actually, using RESHAPE actually slows things down a bit! A better bet is to use

CALL grib_get(imess, 'values', values)

instead of

CALL grib_get_data(imess, lats, lons, values)

so that one extracts 1 2D field ('values') rather than 3 ({'lats', 'lons', 'values'}) each pass. As might be expected, this makes the code 3 times faster for very large GRIB fields, as well as less memory-hungry. See r4116 - r4114 for details.

But we can do a little better (4.3 times faster on a 1/4 degree dump, 2x as fast on a 1 degree grid, ~same speed on 4 degree grid) by using

CALL grib_find_nearest(igrib, .FALSE., lat_in, lon_in, lats, lons, values, distances, indices) 

CALL grib_get_element(imess, 'values', indices, values)

instead of

CALL grib_get_data(imess, lats, lons, values)

This way we're only extracting the required 4 neighbouring points per pass. See r4119 - r4114 for details.

(Surprising it's not better, but I guess someone - namely the GRIB API - still has to search through the array and extract the desired values.)

Can't think how to speed this up any more right now.

[Ian Culverwell]

Unfortunately grib_find_nearest doesn't work with g95 or xlf95. It seems to return the nearest neighbour index, and 3 bits of junk.

Contact ECMWF about this.

Add some code to emulate the correct behaviour (I hope) when it fails.

[Ian Culverwell]

Daniel Varela Santoalla (JIRA) [jira@…] developed a fix thatworks for g95 (1.12.0) and xlf95 (1.9.9):

Hi Ian

I've had a look and apparently the problem comes from a bug on our side.

g95 seems to add a double underscore to methods when calling C functions. All other compilers must be adding just a single one (surely gfortran does that). The result is that a different method is called and I don't think the one with the double underscore was properly tested by us. I've tried this change to fortran/grib_fortran.c (comment the last three lines and add a call to the single underscore FOUR method) and now it seems to work.

I will submit the change to our developers and hopefully will be released soon.


int grib_f_find_nearest_four_single__(int* gid,int* is_lsm,
    double* inlats,double* inlons,
    double* outlats,double* outlons,
    double* values,double* distances,
    int* indexes) {

    printf("FORTRAN2 grib_f_find_nearest_four_single__\n");


  return grib_f_find_nearest_four_single_(gid,is_lsm,
    inlats,inlons,outlats,outlons,values,
    distances,indexes);

/*  return grib_f_find_nearest_single_(gid,is_lsm,
    inlats,inlons,outlats,outlons,values,
    distances,indexes); */
}

So the fix should be part of the next release.

However, I cannot compile grib_api_1.12.0 on xlf95:

configure:11536: xlc -o conftest -O2 -w -q64 conftest.c >&5 ld: 0711-317 ERROR: Undefined symbol: .yywrap ld: 0711-345 Use the -bloadmap or -bnoquiet option to obtain more information.

This will need another JIRA query.

[Ian Culverwell]

Daniel Varela Santoalla (JIRA) [ jira@…] came up trumps again:

Hi Ian

The problem seems to be the option "-qsuffix=f=f90", that confuses the macro that tries to guess the right "module include" option letter for the fortran compiler. The test program generated by the macro is an F77 program with .f extension that is then interpreted as F90 and it all breaks down.

If you remove that option it all should be OK.

Cheers
Daniel

That fixes it (and doesn't unfix grib_api_1.9.9). Change to grib_configure_xlf95_aix committed at r4173.

So, we await ECMWF's implementation of their bugfix to grib_api_1.12.0. Meanwhile, my Fortran fix will work (at 1.9.9 and 1.12.0).

[Ian Culverwell]

I think we can close this ticket for now. r4137 fixes the 'nearest neighbour indexing problem' anyway. (For later ROPP releases, when we recommend a later release of GRIB_API which doesn't suffer from these problems, we may consider removing it.)