Misprocessing of DTocc information in UCAR files

[Ian Culverwell]

Axel von Engeln (EUM) reports:

Some more trouble with new UCAR files, an example for the trouble 
is attached.

It has second of 59.506... which is converted in 
ropp_io_read_ncdf_get.f90 with an NINT to seconds, leading to 60 
seconds here.

CALL ncdf_getatt('second', sec)
DTocc%Second = NINT(sec)
DTocc%Msec = 0

This in turn seems to lead to a start time of 0 seconds in the ROPP 
output file when running through ucar2ropp (or I believe this is 
the reason, I didn't go fully into the details of it just saw that 
the DTocc has 60 seconds and my ROPP output file has a start time 
of 0s).

Axel's file is attached. Clearly we need a better breakdown of second into dtocc%second and dtocc%msec.

[Ian Culverwell]

atmPrf_MTPA.2013.001.01.55.G30_2016.0120_nc

[Ian Culverwell]

By replacing

  DTocc%Second = NINT(sec)
  DTocc%Msec = 0

by

  DTocc%Second = INT(sec)
  DTocc%Msec = NINT(1000.0*(sec - DTocc%Second))

in SUBROUTINE gettime in ropp_io_read_ncdf_get.f90, the results of passing the attached file through ucar2ropp change from

ncks -H -Q -vyear,month,day,hour,minute,second,msec,start_time atmPrf_MTPA.2013.001.01.55.G30_2016.0120_test.nc
day[0]=1 

hour[0]=1 

minute[0]=55 

month[0]=1 

msec[0]=0 

second[0]=99 

start_time[0]=0 

year[0]=2013 

to

ncks -H -Q -vyear,month,day,hour,minute,second,msec,start_time atmPrf_MTPA.2013.001.01.55.G30_2016.0120_test.nc
day[0]=1 

hour[0]=1 

minute[0]=55 

month[0]=1 

msec[0]=507 

second[0]=59 

start_time[0]=410320559.507 

year[0]=2013 

which looks OK, given that second = 59.5067520141602 in the input file. (Check: 01/01/2013 ~ 410 000 000 secs after 01/01/2000, and the ROPP start_time:units = "seconds since 2000-01-01 00:00:00".)

[Ian Culverwell]

The root of the problem, as is often the case, is range-checking. In this case, NINT(59.5067520141602) = 60, which is outside the valid range of seconds ([0, 59]). Range-checking therefore sets ROdata%DTocc%Second to 99. The code to define the start_time (Sec 1.5.1 of ropp_io_write_ncdf_put.f90) says

  IF (isroppinrange(data%dtocc)) THEN
    DT8 = (/data%dtocc%year,data%dtocc%month, data%dtocc%day,0,  &
            data%dtocc%hour,data%dtocc%minute,data%dtocc%second, &
            data%dtocc%msec/)
    CALL TimeSince ( DT8, time, 1, Base="JS2000" ) 
  ELSE
    time = 0.0_wp
  ENDIF
  CALL ncdf_putvar('start_time', time, rec = irec)

isroppinrange(data%dtocc) checks that every element of data%dtocc, including data%dtocc%second, is in the required range, and this obviously fails in this case. Hence time, and start_time, are set to zero.

[Ian Culverwell]

This change has been tested on the dataset above, and with earlier COSMIC files whose seconds were artificially set to 59.9. (Note that for these files second was always a whole number in real form, like 25., so the problem would not have arisen. This also, presumably, explains the original design decision to set DTocc%Msec to zero, and to round DTocc%Second to the nearest integer, not the one below.)

[Ian Culverwell]

The change was committed at r4989.

[Ian Culverwell]

Stig Syndergaard (DMI) suggests that a better fix would be to replace

DTocc%Msec = NINT(1000.0*(sec - DTocc%Second))

by

DTocc%Msec = INT(1000.0_wp*(sec - DTocc%Second))

I agree!

[Ian Culverwell]

Checked that the above correction works if sec=59.999. But sec is defined as single precision in the calling routine; the above trick therefore fails if (for example) second = 59.99999999, since this is rounded to sec = 60.0 before we even start mucking around with DTocc. Defining

  REAL(wp)                           :: sec

in SUBROUTINE gettime(DTocc) allows second = 59.99999999 and even second = 59.99999999999999 to work, although it fails (ie, DTocc%Second gets set to 99, DTocc%Msec gets set to 0, and start_time gets set to 0) if second = 59.999999999999999 - a contingency with which we can probably live.

[Ian Culverwell]

There is one slight infelicity, which is that 59.9_dp is stored in the computer as 59.8999999999999985789145285, which results in msec being set to 899 according to the rounding down rule employed by INT function. Perhaps this was why went for NINT originally?

[Ian Culverwell]

We can get around this last objection by saying

  IF ( sec - DTocc%Second > 0.999_wp ) THEN
    DTocc%Msec =  INT(1000.0_wp*(sec - DTocc%Second))
  ELSE
    DTocc%Msec = NINT(1000.0_wp*(sec - DTocc%Second))
  END IF

Then second=59.9 ==> second=59 and msec=900.

[Ian Culverwell]

This change has been committed at r5051. Closing ticket as 'fixed'.