Correct 'T' area code in BUFR Abbreviated Routing Header

[Ian Culverwell]

Owen Lewis points out that there's a problem with the BUFR area code 'T'. Apparently it should mean 45W - 180E/W excluding 0E/W, i.e. the Americas, not 45W - 180E/W including 0E/W, i.e. Europe and Asia.

This error came to light because some ground-based measurements from China were incorrectly assigned the 'T' area code.

We therefore need to change

      IF ( LN >= 315 .OR. LN <= 180 ) THEN       ! Europe/Asia

to

      IF ( LN <= 315 .AND. LN >= 180 ) THEN      ! Americas

in

FUNCTION GTShdrGAD ( Lat, Lon ) RESULT(A2)

in ropp_io/bufr/gtshdrs.f90. The comment section also needs updating.

[Ian Culverwell]

We can test this out by running ropp2bufr_ecodes with the -gi option, since this tells the user the commands needed to run the perl script ropp_io/tools/gtsheaders_bufr.pl, which wraps the BUFR message in the GTS headers.

If we nullify the latitude in the ropp2bufr_eccodes test file, so that it contains

prompt:> ncks -H -Q -vlon,lat ../data/ropp_temp.nc
netcdf ropp_temp {
  dimensions:
    dim_unlim = UNLIMITED ; // (1 currently)

  variables:
    float lat(dim_unlim) ;

    float lon(dim_unlim) ;

  data:
    lat = -9.9999e+07 ;

    lon = -3.89 ;

} // group /

and then pass this through ropp2bufr_eccodes, we find (after adding a few print statements)

prompt:> ../tools/ropp2bufr_eccodes ../data/ropp_temp.nc -gi -o bufr_test2_gi.bfr

---------------------------------------------------------------------
                    ROPP to BUFR (ecCodes) Encoder
---------------------------------------------------------------------

 
WARNING (from ropp2bufr_eccodes):  '-gi' requested, but GTS header encoding not currently available in ecCodes.
   For this functionality either use the perl script at $ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl with the instructions below,
   or use the MetDB or ECMWF BUFR packages.

INFO (from ropp2bufr_eccodes):  Reading  ROPP data from ../data/ropp_temp.nc
INFO (from ropp2bufr_eccodes):  Encoding profile    1 : OC_20090817215807_META_G027_DMI_
 in GTShdrGAD, lat, lon =    -99999000          -3
 LN =          357
 in GTShdrGAD, A2 = T

INFO (from ropp2bufr_eccodes):  To wrap BUFR file in GTS ARH + IPH headers, run:
   perl $ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl --iph bufr_test2_gi.bfr bufr_test2_gi.bfr_gts (for example) IUTT14 EKMI 172158

But the area code 'T' is wrong: the longitude of -3.89 deg E means that this profile is in the 45W-180E/W sector that includes 0E/W, i.e. Europe and Asia. We now think 'T' should apply in the 45W-180E/W sector that excludes 0E/W, i.e. the Americas. By making the change to GTShdrGAD that was suggested above, and rerunning, we find

prompt:> ../tools/ropp2bufr_eccodes ../data/ropp_temp.nc -gi -o bufr_test2_gi.bfr

---------------------------------------------------------------------
                    ROPP to BUFR (ecCodes) Encoder
---------------------------------------------------------------------

 
WARNING (from ropp2bufr_eccodes):  '-gi' requested, but GTS header encoding not currently available in ecCodes.
   For this functionality either use the perl script at $ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl with the instructions below,
   or use the MetDB or ECMWF BUFR packages.

INFO (from ropp2bufr_eccodes):  Reading  ROPP data from ../data/ropp_temp.nc
INFO (from ropp2bufr_eccodes):  Encoding profile    1 : OC_20090817215807_META_G027_DMI_
 in GTShdrGAD, lat, lon =    -99999000          -3
 LN =          357
 in GTShdrGAD, A2 = X

INFO (from ropp2bufr_eccodes):  To wrap BUFR file in GTS ARH + IPH headers, run:
   perl $ROPP_SRC/ropp_io/tools/gtsheaders_bufr.pl --iph bufr_test2_gi.bfr bufr_test2_gi.bfr_gts (for example) IUTX14 EKMI 172158

which is correct.

Conversely, by setting the longitude to -100E, the correct area code of 'T' is suggested (not shown).

[Ian Culverwell]

This change has no effect in the core tests, which don't test the -g or -gi options, or the test folder (ditto), so it seems safe to commit at r6504.

[Ian Culverwell]

Testing this has uncovered some problems with the non-ecCodes tool ropp2bufr. Thus, for example, when we run the above test netCDF file through this exec, we get

prompt:> ../tools/ropp2bufr ../data/ropp_temp.nc -gi -o bufr_test2.bfr                   

---------------------------------------------------------------------
                    ROPP to BUFR (ecbufr) Encoder
---------------------------------------------------------------------

INFO (from ropp2bufr):  Reading  ROPP data from ../data/ropp_temp.nc
INFO (from ropp2bufr):  Encoding profile    1 : OC_20090817215807_META_G027_DMI_
                   ECMWF 
 
      BUFR ENCODING SOFTWARE VERSION -  7.2 
            1 April  2007. 
 
 
 
 Your path for bufr tables is :
 /data/users/idculv/ROPP/data/bufr/       
BUFR TABLES TO BE LOADED  B0000000000000012000.TXT,D0000000000000012000.TXT
 in GTShdrGAD, lat, lon =            0          -4
 in GTShdrGAD, A2 = E
INFO (from ropp2bufr):  Total of  14794 bytes written to bufr_test2.bfr
INFO (from ropp2bufr):  Generated 1 GTS bulletin to bufr_test2.bfr

Note that it has rounded the lon of -3.96 E to -4 E, which is fine. But it has reset the missing data latitude of -9.9999e+07 to zero, which means it returns a 'tropical' area code of 'E', and completely bypasses the possibility of returning area codes of 'N', 'S', 'T' or 'X'. So it doesn't look like this code ever worked with a missing data lon (or lat). Therefore promoting this ticket to ROPP-11 rather than closing it.

[Ian Culverwell]

Note that the justification for this interpretation of the area code of 'T' comes from Fred Branski, former president of CBS and general Grand Poobah¹ of matters GTS, who said:

“The sequencing is intended to provide global coverage quadrants for both northern and southern hemisphere and the tropical belt. The convention for direction is always westerly.

The purpose for "T" was to be able to provide products for the entire North American continent, model, satellite, tropical or whatever. It should follow the westerly direction convention as used throughout the table.”

¹(c) Simon Elliott, EUMETSAT

[Ian Culverwell]

ropp_temp.nc

[Ian Culverwell]

ropp_temp2.nc

[Ian Culverwell]

ropp_temp.nc and ropp_temp2.nc test files attached.

[olewis]

I have found that the when using ropp2bufr with ECMWF bufr that there is arithmetic overflow from the missing values. Missing value is initially set to 1.7E38 which is a float. When the Lat and Lon are passed to the gtshdrs.f90 module they are passed through the NINT function which is changing the float to the nearest integer. The nearest integer for this missing value is -2147483648 which is the smallest integer. When this is passed through the GTShdrGAD function the Lat and Lon values are put through the absolute function to see if they are smaller than the GTSHDR_FLGVAL which is set at 999. With -2137483648 this is the smallest value int but the largest int is 2147483647 and so the negative value does not get changed. This then leads to the wrong area code being assigned.

[olewis]

I have put in a solution to solve the problem for the ROPP2BUFR routine with the ECMWF BUFR when there are missing values in the Lat or the Lon.

The current output is as followed.

{{{prompt:> tools/ropp2bufr data/ropp_temp.nc -gi -o bufr_test.bfr

---

ROPP to BUFR (ecbufr) Encoder

---

INFO (from ropp2bufr): Reading ROPP data from data/ropp_temp.nc INFO (from ropp2bufr): Encoding profile 1 : OC_20090817215807_META_G027_DMI_

ECMWF

BUFR ENCODING SOFTWARE VERSION - 7.2

1 April 2007.

Your path for bufr tables is : /data/users/olewis/ROPP_ECMWF_BUILD/data/bufr/

BUFR TABLES TO BE LOADED B0000000000000012000.TXT,D0000000000000012000.TXT

In ropp2bufr_ec Lat= 1.700000000000000E+038 In ropp2bufr_ec NINT of Lat= -2147483648 In GTShdrGAD, Lat, Lon = -2147483648 -4 In GTShdrGAD, A2 =I

INFO (from ropp2bufr): Total of 14794 bytes written to bufr_test.bfr INFO (from ropp2bufr): Generated 1 GTS bulletin to bufr_test.bfr }}}

The test file has a missing value for the latitude which the ECMWF BUFR uses the Flag value of 1.7E38. When the latitude value is used it is changed to the nearest integer and this value is -2147483648. In the gtshdrs code the latitude is absoluted to check it against a flag value. The -2147483648 does not have a positive equivalent in Fortran and so the ABS function does not do anything and this leads to the area code being assigned incorrectly.

{{{prompt:> tools/ropp2bufr data/ropp_temp.nc -gi -o bufr_test.bfr

---

ROPP to BUFR (ecbufr) Encoder

---

INFO (from ropp2bufr): Reading ROPP data from data/ropp_temp.nc INFO (from ropp2bufr): Encoding profile 1 : OC_20090817215807_META_G027_DMI_

ECMWF

BUFR ENCODING SOFTWARE VERSION - 7.2

1 April 2007.

Your path for bufr tables is : /data/users/olewis/ROPP_ECMWF_BUILD/data/bufr/

BUFR TABLES TO BE LOADED B0000000000000012000.TXT,D0000000000000012000.TXT

In ropp2bufr_ec Lat= 1.700000000000000E+038 In ropp2bufr_ec NINT of Lat= -2147483648 In GTShdrGAD, Lat, Lon = -2147483647 -4 In GTShdrGAD, A2 =X

INFO (from ropp2bufr): Total of 14794 bytes written to bufr_test.bfr INFO (from ropp2bufr): Generated 1 GTS bulletin to bufr_test.bfr }}}

With the change which checks if the values of longitude and latitude in ropp2bufr_ec are equalt to -2147483648 and if so adds 1 to put the absolute value of this into the acceptable range. This leads to the correct area code being produced of X for global when the Longitude or Latitude is missing.

[olewis]

build_ropp log file

[olewis]

The log file from the build of ROPP with this change has been attached to this ticket. The change to the code is in ropp2bufr_ec.f90 and is adding {{{IF (Lat == -2147483648) THEN

Lat = Lat+1

END IF IF (Lon == -2147483648) THEN

Lon = Lon+1

END IF

}}} before GTShdrARH is called.

[Ian Culverwell]

When I had a look at this, by adding

print*,'ABS(-2147483648) = ', ABS(-2147483648)
print*,'ABS(-2147483647) = ', ABS(-2147483647)

to a fortran file, gfortran complained that

gfortran -I../build -O2 -I/data/users/idculv/ROPP/gfortran
/include -c -o ropp_1dvar_dbangle.o  ropp_1dvar_dbangle.f90
ropp_1dvar_dbangle.f90:432.45:

print*,'ABS(-2147483648) = ', ABS(-2147483648)
                                             1
Error: Integer too big for its kind at (1). This check can be
disabled with the option -fno-range-check"

So I don't think -2147483648 is legal fortran. Perhaps the fix should read

IF (Lon < -2147483647) Lon = -2147483647

etc.

(When I compile with -fno-range-check I find

 ABS(-2147483648) =  -2147483648
 ABS(-2147483647) =   2147483647

as expected.)

[olewis]

Replying to idculv:

When I had a look at this, by adding

print*,'ABS(-2147483648) = ', ABS(-2147483648)
print*,'ABS(-2147483647) = ', ABS(-2147483647)

to a fortran file, gfortran complained that

gfortran -I../build -O2 -I/data/users/idculv/ROPP/gfortran
/include -c -o ropp_1dvar_dbangle.o  ropp_1dvar_dbangle.f90
ropp_1dvar_dbangle.f90:432.45:

print*,'ABS(-2147483648) = ', ABS(-2147483648)
                                             1
Error: Integer too big for its kind at (1). This check can be
disabled with the option -fno-range-check"

So I don't think -2147483648 is legal fortran. Perhaps the fix should read

IF (Lon < -2147483647) Lon = -2147483647

etc.

(When I compile with -fno-range-check I find

 ABS(-2147483648) =  -2147483648
 ABS(-2147483647) =   2147483647

as expected.)

I have now tried this option with gfortran and this compiler does not treat the NINT overflow in the same way as the other compilers. When using NINT on 1.7E38 it does not roll over to -2147483648. Instead the value comes out as 0 which is incorrect!

I think another solution may be required to solve this. I think we should add a variable to ropp2bufr_ec to define the ecmwf missing value flag and then have IF statements to act accordingly to this.

[olewis]

     REAL(dp), PARAMETER :: MISSING = 1.7E38_dp    ! Real    missing data flag value (ECMWF)
     IF (ABS(Values(30) - MISSING) < 1.0e-5_dp) THEN 
          Lat = -2147483647
      ELSE
          Lat = NINT(Values(30))
      END IF
      
      IF (ABS(Values(31) - MISSING) < 1.0e-5_dp) THEN
          Lon = -2147483647
      ELSE
          Lon = NINT(Values(31))
      END IF

This is the solution that has been implemented. All the compilers build these changes successfully and the change has succeeded in fixing the issue of missing values being reassigned to 0 and therefore the wrong area code being assigned to the BUFR.

[olewis]

The results from the different compilers for the area results are shown below. These were run using the ropp_temp.nc file. The file was run through the ropp2bufr routine built with the ECMWF BUFR and then the od -c tool was used on the outputted BUFR file to get the header info from the binary files.

ifort control I U T I 1 4 E K M I

ifort I U T X 1 4 E K M I

ifort16 control I U T I 1 4 E K M I

ifort16 I U T X 1 4 E K M I

gfortran control I U T E 1 4 E K M I

gfortran I U T X 1 4 E K M I

g95 control I U T X 1 4 E K M I

g95 I U T X 1 4 E K M I

pgf95 control I U T I 1 4 E K M I

pgf95 I U T X 1 4 E K M I

Nagfor61 is already known to not work with ECMWF BUFR.

These results demonstrate that the solution fixes the issue with the NINT function changing the missing value flag to acceptable solutions for the latitude/longitude and therefore fixing the observation to an area code which it is not located in. All compilers now produce the global character X which is the value that should be produced when either the latitude o longitude values are missing.

[Ian Culverwell]

ropp_temp3.nc

[Ian Culverwell]

Tried to reproduce this with gfortran.

Control on ropp_temp.nc (lon=-3.89, lat=-9.9999e7): I   U   T   E   1   4
Control on ropp_temp2.nc (lon=-100, lat=-9.9999e7): I   U   T   F   1   4

These would be correct if the profile were in the tropics.

When I try the fix, I get

Test on ropp_temp.nc (lon=-3.89, lat=-9.9999e7): I   U   T   X   1   4
Test on ropp_temp2.nc (lon=-100, lat=-9.9999e7): I   U   T   T   1   4

X for the first is good. Why have we got T in the second one? Because FUNCTION GTShdrGAD ( Lat, Lon ) says

      IF ( LN <= 315 .AND. LN >= 180 ) THEN    ! Americas
        A2 = "T"
      ELSE                                     ! Global
        A2 = "X"
      END IF

But this is only activated if ABS(Lat) >= GTSHDR_FLGVAL, whereas it should be activated if Lat > 0 (N Hem).

(Note that ​https://community.wmo.int/table-c3 says

    T    45°W - 180° northern hemisphere

and this is what it says in the ROM SAF BUFR doc too.)

When I fix this in ropp_io/bufr/gtshdrs.f90, I get

Test on ropp_temp.nc (lon=-3.89, lat=-9.9999e7): I   U   T   X   1   4
Test on ropp_temp2.nc (lon=-100, lat=-9.9999e7): I   U   T   X   1   4

which makes more sense, because the latitude is still undefined.

On a new file (attached), sited over N America, I get

Control on ropp_temp3.nc (lon=-100, lat=45):     I   U   T   B   1   4
Test on ropp_temp3.nc (lon=-100, lat=45):        I   U   T   T   1   4

B in the first case is further evidence that the gtshdr.f90 algorithm is not working; T in the second is correct.

[Ian Culverwell]

Analagous changes need to be made to eum2ropp_ec.f90. These, and the changes to gtshdrs.f90 have been made at r6859.

[Ian Culverwell]

Check that eum2ropp (built on ECBUFR libs) also works OK now.

Cntl on ../data/eum_test.n4 (lon=162.36, lat=-49.36): I   U   T   K   1   4
Cntl on eum_temp.nc (lon=162.36, lat=-9.9999e7)     : I   U   T   G   1   4
Cntl on eum_temp2.nc (lon=-100, lat=-9.9999e7)      : I   U   T   F   1   4
Cntl on eum_temp3.nc (lon=-100, lat=45)             : I   U   T   B   1   4

The first is correct; the other 3 are wrong.

Test on ../data/eum_test.n4 (lon=162.36, lat=-49.36): I   U   T   K   1   4
Test on eum_temp.nc (lon=162.36, lat=-9.9999e7)     : I   U   T   X   1   4
Test on eum_temp2.nc (lon=-100, lat=-9.9999e7)      : I   U   T   X   1   4
Test on eum_temp3.nc (lon=-100, lat=45)             : I   U   T   T   1   4

All 4 are correct.

[Ian Culverwell]

Cosmetic change to eum2bufr_ecodes.f90 made at r6860. (Will make it into the reference manual.)

[Ian Culverwell]

Closing ticket as 'fixed' for now.

[Ian Culverwell]

The two meanings of NH etc have been spelled out in the ROM SAF BUFR doc at r6861.

[Ian Culverwell]

Sorry: one further tweak: MISSING can be replaced by the identically valued RVIND, which is in the module ropp2bufr, which is already being USEd. Checks out OK, so commit this at r6862.

[Ian Culverwell]

Stig Syndergaard (DMI) raised some complaints about this (nearly two years after the work was done):

We are in the process of getting Spire data ready from ROM SAF. They will be disseminated via the GTS, and we have been working on the routing header recently. Yesterday I found that ROPP 11 sets the ’T’ area code in the routing header when the location is within 45W - 180W in NH, as in Table 13 here: ​https://www.romsaf.org/romsaf_bufr.pdf. I expected only the letters A-L.

I found that it was changed in gtshdrs.f90 here: https://trac.romsaf.org/ropp/changeset/6859. I also read the relevant parts in this ticket: https://trac.romsaf.org/ropp/ticket/691. Somewhere near the end of that ticket you wrote: “But this is only activated if ABS(Lat) >= GTSHDR_FLGVAL, whereas it should be activated if Lat > 0 (N Hem).” And you changed the code accordingly. But is it correct? In the description of FUNCTION GTShdrGAD (line 1041 of gtshdrs.f90) it says:

When both Lat and Lon values are within their normal ranges, the appropriate A-L code is determined; if |Lon|>=999 , N or S is selected from the sign of Lat; if |Lat|>=999 and -180<=Lon<=-45, (sector 180-315 deg E for the Americas) then T is returned, otherwise X; and if both |Lat|>=999 and |Lon|>=999, X is returned.

So, I think the intention was that it should be letters A-L when (Lat, Lon) is within normal ranges, and only letter T when ABS(Lat) >= GTSHDR_FLGVAL. I think it was meant as a way to set the T area code instead of some of the A-L. Would you agree with this?

I said that I thiought T meant NHem, not unspecified hemisphere, and that Dave's comments should not be considered sacrosanct.

Stig countered:

If you are right, then how would a user of this subroutine specify that they want the A-L codes and not the T code? Note, that the T area overlap fully the B area and partly the A, E, and F areas. I believe T is meant as an alternative to these codes, not as the default to overwrite them. I don’t think the note by Fred Branski goes against that.

As I understand the subroutine, the user can get the T area code by inputting Lat >= 999. It is up to the user to find out how to do that although the real latitude is in the Northern hemisphere. It is just a way of getting the T code instead of the A-L codes. I imagine the user would say something like:

IF (I want T code and (Lat, Lon) is in the right range) THEN
  GTShdrGAD(999, Lon)
ELSE
  GTShdrGAD(Lat, Lon)
ENDIF

And:

The new code in ROPP 11 does not make as much sense to me because it does not allow me to get all A-L codes. We have been using those for Metop for many years and I have not heard any wishes to change that. So A-L is what we prefer to continue to give users (at least we have not discussed otherwise).

The T code may be preferred by other providers, and that is alright. They can get it by the example I wrote (using the earlier code), and I think it is obvious from the comments in the code that this was the original intention. I can understand that you think it should be by other more clear input options, and that would be fine with me too, but the code was what it was, and the documentation in the code of intentionally using 999 to control these different outputs was clear to me. However, the intended functionality is broken in ROPP 11, as I see it.

And:

I looked into gtshdrs.f90 a little closer, and found that also the GTShdrARH_arr subroutine needs update regarding the wrong interpretation of the T area being the North Americas and not Europe/Asia. Obviously Dave Offiler who implemented all this thought that 45W-180 meant 45W-180E whereas the correct interpretation is 45W-180W - I agree with you on that. Although we don’t use GTShdrARH_arr in our operational processing, I did the change in our dmi_trunk, where I also reverted GTShdrGAD to what is was earlier (https://trac.romsaf.org/ropp/changeset/7201/ropp_src/branches/dev/Share/dmi_trunk_11.0/ropp_io/bufr/gtshdrs.f90).

Note that in GTShdrARH_arr the Lat is set to GTSHDR_FLGVAL if we are in the 'T' area, and it then calls GTShdrARH_sca which calls GTShdrGAD. The GTShdrARH_arr subroutine takes an array of (lat,lon) values and finds a common header with ’T’ or ’N’ or ’S’ if one of A-L does not apply to all (lat,lon) pairs.

I suggest that the issue be raised for discussion and possible inclusion in ROPP-12.0, using the new proposal/justification/suggested fix/test/review/documentation development model for ROPP that we have discussed before, and which is used pretty much everywhere except the ROM SAF.

Reopening ticket and reassigning it to ROPP-12.0.