Bug in DATUM_HMSL()

[Dave Offiler]

While reviewing the EGM96 calculations in ROPP, I found a small bug in DATUM_HMSL() (part of earth.f90 in ROPP_UTILS). This is documented in RSR-16 (in draft).

The effect is to introduce a sinusoidal error in undulation values (compared with the 'official' NMA F477.F program) within a few degrees of the equator, having a short wavelength of 1 degree and plus/minus 3mm peak (and a longer on/off wavelength). (See attached image) This is of no consequence to real-world usage in ROPP, but is clearly a systematic error.

Investigation pinned this down to the initialisation in the DATUM_HMSL() code which loads the coefficients and corrections values from the two EGM96 input files. E.g. for the former:

    DO WHILE ( ierr == 0 )
      READ ( UNIT=PUNIT, &
              FMT=*,     &
           IOSTAT=ierr ) N, M, C, S
      IF ( ierr /= 0 .OR.     &
           ( N  == NMAX .AND. &                ! (1.1)
             M  == NDEG ) ) EXIT               ! (1.1)

      N = ( N * ( N + 1 ) ) / 2 + M + 1
      HC(N) = C
      HS(N) = S
    END DO

The purpose of the IF statement is to terminate the loop on any I/O error, on EoF (when ierr is negative) or when the maximum degree & order are found, so as to avoid going 'off the end' of the arrays in memory.

The bug is the fact that when the latter test is TRUE, the loop terminates, but without having stored the last pair of values in the arrays (these having been previously initialised to zero).

Removing this test corrects the undulation errors, such that values completely agree with the NMA calculation to <<1mm (see RSR-16). Note that simply changing the loop to:

    DO
      READ ( UNIT=PUNIT, &
              FMT=*,     &
           IOSTAT=ierr ) N, M, C, S
      IF ( ierr1 /= 0 ) EXIT
      N = ( N * ( N + 1 ) ) / 2 + M + 1
      HC(N) = C
      HS(N) = S
    END DO

is not the whole story, as ierr is later tested for zero before the loading is flagged as OK. This test needs to also accept negative (EoF) values. It also now makes the assumption that the file has no records following the M=360 and N=360 values (presumably an issue at the time this test was introduced).

The same fix is needed for the corrections file loading. NB Fix only tested with ifort v12 and v14.

Given that this is a minor bug of no practical import, and (if the recommendations in RSR-16 are accepted) DATUM_HMSL() will anyway be replaced by a completely new routine, this fix is not being checked into SVN for now. Should the proposed new routine not be accepted as a replacement, then the fixed code needs to to tested on the usual range of compilers and will then be checked in for ROPP-8 (or perhaps v7.1).

[Dave Offiler]

Image showing systematic 3mm undulation errors

[Ian Culverwell]

I've taken over the implementation of this.

[Ian Culverwell]

Compiles OK with ifort, and it's clear that the N=M=360 harmonic was missing before:

ROPP7.1

... (from grib2bgrasc):  lon_in =   -0.50000
... (from grib2bgrasc):  lat_in =    0.00000
... (from grib2bgrasc):  file_out = cntl.nml
INFO (from grib2bgrasc):  Reading file /net/data/nwp1/idculv/ROPP/ropp_src/branches/dev/Share/ROPP80_prototype/ropp_io/data/an_20121001000000_T+0.grib
 ropp7.1 reading coeffs
 NUM, HC(NUM), HS(NUM) =        65341  0.000000000000000E+000  0.000000000000000E+000
 ropp7.1 reading corrns
 NUM, CC(NUM), CS(NUM) =        65341  0.000000000000000E+000  0.000000000000000E+000

ropp7.1 undulation =     0.17156920889114871898E+02

ROPP8.0

... (from grib2bgrasc):  lon_in =     -0.50000
... (from grib2bgrasc):  lat_in =      0.00000

INFO (from grib2bgrasc):  Reading file /net/data/nwp1/idculv/ROPP/ropp_src/branches/dev/Share/ROPP80_prototype/ropp_io/data/an_20121001000000_T+0.grib
... (from extract_prof_from_GRIB):  GRIB Edition 2
... (from extract_prof_from_GRIB):  ECMWF generatingProcessIdentifier =  142
 ropp8.0 reading coeffs
 NUM, HC(NUM), HS(NUM) =        65341 -4.475163896780000E-025 -8.302249455249999E-011
 ropp8.0 reading corrns
 NUM, CC(NUM), CS(NUM) =        65341  2.054457448700123E-018  7.252469661987566E-003

ropp8.0 undulation =     0.17156920889114875450E+02

This generates +/- 3mm differences, varying with 1 deg wavelength, as might be expected by omitting an M=360 spherical harmonic (all these at lat=0):

At lon = 0.0
ropp7.1 undulation =     0.17161578498270731785E+02
ropp8.0 undulation =     0.17161578498270731785E+02
ie difference in undulation[0]=0.

At lon = 0.25
ropp7.1 undulation =     0.17082670300284675591E+02
ropp8.0 undulation =     0.17079671915399114113E+02
ie difference in undulation[0]=-0.00299835 = -3mm, as DO found.  

At lon = 0.5
ropp7.1 undulation =     0.17012674571632835807E+02
ropp8.0 undulation =     0.17012674571632835807E+02
ie difference in undulation[0]=0.

At lon = 0.75
ropp7.1 undulation =     0.17012642029719213355E+02
ropp8.0 undulation =     0.17015640414604771280E+02
ie difference in undulation[0]=+0.00299835 = +3mm, as DO found.  

At lon = 1.0
ropp7.1 undulation =     0.16997073488848911893E+02
ropp8.0 undulation =     0.16997073488848915446E+02
ie difference in undulation[0]=0.

So it's on a 1 deg wavelength, as you might expect by omitting N=M=360.

For -ve lons:

lon=-0.25 ==>  undulation[0]=0.00299835 
lon=-0.50 ==>  undulation[0]=0
lon=-0.75 ==>  undulation[0]=-0.00299835 
lon=-1.00 ==>  undulation[0]=0

I don't, however, see any evidence of the ~10 deg amplitude modulation: I get the same 3 mm oscillation at 1 deg freq all along the equator. I assume something else was responsible for this, as M=360 _is_ just sin(2pi lon/1 deg).

[Ian Culverwell]

Committed change to ROPP8.0 at r4163.

[Ian Culverwell]

Closing ticket.

[Ian Culverwell]

Testing with nagfor revealed a nasty bug in the new logic: it might be trying to read from a file that it couldn't open (eg because it doesn't exist).

This has been fixed at r4192.

We also take the opportunity to issue warning messages if it can't open COEF_file or CORR_file, and to print out the undulation, to draw attention to the fact that it might have failed to find one.

[Ian Culverwell]

Reclosing ticket.

[Ian Culverwell]

No, really.