id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
345	Bug in DATUM_HMSL()	Dave Offiler	Ian Culverwell	"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).

"	defect	closed	minor	8.0	ropp_utils	7.0	fixed	EGM96