Opened 13 years ago

Closed 6 years ago

Last modified 6 years ago

#281 closed enhancement (fixed)

Leap Seconds

Reported by: Dave Offiler Owned by: Dave Offiler
Priority: normal Milestone: 9.1
Component: ropp_utils Version: 6.0
Keywords: Date and time, timesince, GPS epoch, GPS week number Cc:

Description

The TimeSince() routine in ROPP_UTILS/datetime/datetime.f90 supports converting calendar date/times to elapsed time in seconds since the GPS epoch (00:00:00 6-Jan-1980) or vice-versa. This option is used within the ROPP_PP module.

Stig points out that TimeSince() does not take into account any accumulated leap seconds which have been applied to UTC since the GPS epoch; there is currently a 15s difference between UTC and GPS times (an offset value contained within the GPS broadcast message).

While Stig believes that the ROPP_PP application has no obvious error caused by neglecting leap seconds, the routine ought to correctly account for them.

It would also be useful to have the functionality to convert GTS times which are given in terms of Week Number since epoch (modulo 1024) and seconds into the week, to UTC calendar date/time (and vice-versa).

References:

Attachments (4)

GRAS_1B_M02_20161006005636Z_20161006005939Z_N_O_20161006020703Z_G17_NN.nc (14.4 MB ) - added by Ian Culverwell 8 years ago.
GRAS_1B_M02_20161006005636Z_20161006005939Z_N_O_20161006020703Z_G17_NN.nc
EPS-SG_Example.bg.nc (2.2 MB ) - added by Ian Culverwell 7 years ago.
EPS-SG_Example.bg.nc
leapsec_plot_p3.gif (54.6 KB ) - added by Ian Culverwell 7 years ago.
leapsec_plot_p3.gif
leapsec_plot_p1.gif (37.8 KB ) - added by Ian Culverwell 6 years ago.
leapsec_plot_p1.gif

Change history (39)

comment:1 by Dave Offiler, 13 years ago

Version: 5.16.0

comment:2 by Dave Offiler, 12 years ago

Status: newaccepted

The ROPP-6 (and earlier) datetime module is basically my genlib all-in-one package datetime.f90 but split into individual routines and with a separate formal interface file. The datetime.f90 module has been substantially re-written this year and amongst other things, now takes proper account of leap seconds (up to June 2012; if any are announced for 2013 and beyond, they would need to be manually added to the internal table - leap second additions are unpredictable and only announced 6-months in advance). Since there have been some API changes & new routines, the new packages is now called timedate.f90 (USE timedate). (Almost all the API changes are concerned with (re-)formatting date/time string representations, not used in ROPP.)

This ticket, and partially #245 (code consolidation), can be actioned by replacing the separate routines in ropp_utils/datetime with the new timedate package. The rest of the ROPP code base would then need to be checked to ensure it remains compatible.

comment:3 by Ian Culverwell, 12 years ago

Milestone: 6.17.0

Ain't gonna happen at 6.1 (if ever). Postponing to ROPP7.0

comment:4 by Dave Offiler, 11 years ago

Milestone: 7.07.1

New timedate.f90 module and use thereof throughout ROPP-land implemented in DO V7.0 branch, but we decided not to merge this complex branch (which also includes several similar F90 module changes with knock-on to most ROPP components) at v7.0. Postponed until v7.1

comment:5 by Ian Culverwell, 11 years ago

Milestone: 7.18.0

comment:6 by Dave Offiler, 10 years ago

Milestone: 8.09.0

Needs re-doing from scratch (new branch) incorporating a few small updates to my GENLIB timedate.f90 module since 2012 (including up-to-date list of announced Leap Seconds). Too late for v8.0 now.

comment:7 by Dave Offiler, 8 years ago

Cc: ian.culverwell@… added

Extracted basic leapseconds() function as new source code file leapseconds.f90 (addeding newly announced l/sec on 31-Dec-2016) and used in timesince() only when converting to/from GPS time. Also updated datetimeprogs.f90 and Makefile.am to build the extra source file. Passed to Ian for review as a stand-alone zip package before (or instead of) svn ci.

comment:8 by Ian Culverwell, 8 years ago

ropp_utils/datetime/leapseconds.f90 has been extended to make it easier to add any subsequent leap seconds, and to allow negative (or even multiple) leap seconds, although none has yet occurred.

A note to remind the ROPP Development Manager to check for additional leap seconds has been added to the ROPP development guide.

This tests out OK, as is shown by running ropp_utils/datetime/lstest.f90 on a few input dates/times. There remains a one second slip if one tries to 'do a return trip' from a time within [0, 1) seconds of one of the leap second datetimes. For example:

Today is 1985  6 30  0 23 59 59 999 2446247.49999998835846781731
Reverse = 1985  6 30  0 23 59 59 999  2446247.49999998835846781731        
in - out = 0.00000000000000000000 msec

Today is 1985  7  1  0  0  0  0   0  2446247.50000000000000000000
Reverse = 1985  6 30  0 23 59 59   0  2446247.49998842598870396614      
in - out = 999.99457597732543945312 msec

Today is 1985  7  1  0  0  0  0 999  2446247.50001156236976385117
Reverse = 1985  6 30  0 23 59 59 999  2446247.49999998835846781731      
in - out = 999.99457597732543945312 msec

Today is 1985  7  1  0  0  0  1   0  2446247.50001157401129603386
Reverse = 1985  7  1  0  0  0  1   0  2446247.50001157401129603386        
in - out = 0.00000000000000000000 msec

Fixing this would probably require making ropp_utils/datetime/timesince.f90 RECURSIVE, and the number of times when it matters (27 seconds in 44.5 years ~ 20 times in a billion), and the relative harmlessness of making an error on such occasions means that we can probably safely ignore this contingency.

comment:9 by Ian Culverwell, 8 years ago

In fact, it wouldn't be possible to fix this, because there's a fundamental ambiguity for times near a leapsecond day. For example:

cdt1 = (/ 1985,  6, 30,  0, 23, 59, 59, 500  /)

and

cdt1 = (/ 1985,  7,  1,  0,  0,  0,  0, 500  /)

both return

 gps1 =  173059197.500

from a call to

CALL timesince ( cdt1, gps1, 1, base='gps')

This is because at the first time, 12 leap seconds have elapsed, which must therefore be deducted from the GPS time. The second time, however, one second later, occurs after the introduction of a leap second on 0Z/01-07-1985, and therefore requires 13 leap seconds to be deducted. Hence both instants have the same GPS time.

There is therefore no 'correct' CDT corresponding to this GPS time (or indeed any GPS time within 1 sec of a leapsecond time). The chosen algorithm chooses the smaller of the two possible CDTs which correspond to such a given 'bordeline' GPS time, which is why we see the above 'round trip' differences for CDTs of 00Z+[0,1) sec 01/07/1985 that were reported above. Round trip differences for 00Z-(0,1] sec are within +/- 1 msec, as they are for general times.

comment:10 by Ian Culverwell, 8 years ago

This checks out OK, so commit at r5072.

comment:11 by Ian Culverwell, 8 years ago

Resolution: fixed
Status: acceptedclosed

Ticket owner says

Ian - I'm OK, that's all that can be done on this ticket. /D

so closing ticket.

comment:12 by Ian Culverwell, 8 years ago

Resolution: fixed
Status: closedreopened

Reopening this ticket as I'm not happy with several aspects:

  • Only 9 leap seconds, I think, were added before 00Z/01-01-1980, not 10, as used in ropp_utils/datetime/timesince.f90. (Obviously easily fixed.)
  • A base time of 'GPS' is used in ropp_pp/preprocess/ropp_pp_preprocess_cosmic.f90:
      CALL TimeSince ( DT8, start_time, 1, Base="GPSSEC" )
    

The leapsecond change therefore affects a lot of the ropp_pp tests, and we get a lot of warning messages like

WARNING (from ropp_io_read_ncdf_get):  'start_time' and yr/mo/dy/hr/mn/sc/ms timestamps differ by 15.000 seconds - using yr/../ms timestamp

All these need checking, in principle.

  • Axel has identified an issue with the start_time in eum2ropp, which this change signally fails to fix:
    Hi Ian,
    
    I am not sure to what extent it is a problem at other 
    places of ROPP, but when running e.g. the stuff I 
    added, the eum2ropp, with the attached input / output 
    files, I get from the input level_1b group the 
    following fields for the start time:
    
         utc_start_absdate = 6123 ;
         utc_start_abstime = 3396.67 ;
    
    When I convert that to Julian seconds, without taking 
    into account any leap seconds, I get: 
    
    IDL> print, 6123.D * 86400.D + 3396.67D, format='(F20.10)'
    529030596.6700000167
    
    Which is also the number I find in the ROPP output file for eum2ropp:
    
    start_time = 529030596.670009 ;
    
    But this somehow cannot be correct, since my Julian 
    second conversion is without leap seconds, thus the 
    real number should be 4s more. So I obviously didn't 
    take that into account back then, but assume that ROPP 
    is affected in general.
    
    Cheers,
    
                Axel
    

Axel is therefore expecting start_time in Julian secs to be in advance of the given start_time, which is calculated without accounting for leapseconds, and also to be relative to the quoted base time (01/01/2000 in this case). Subtracting leapseconds from GPS time, always relative to the start of the GPS epoch on 01/01/1980, doesn't do the job.

For background information, I checked that in the attached file (dated 06/10/2016):

/data/level_1a/utc_start_absdate = 6123
/data/level_1a/utc_start_abstime = 3396.67 

/data/level_1a/gps_start_absdate = 6123
/data/level_1a/gps_start_abstime = 3413.67

gps_start_abstime is 17 secs (= the accumulated leap secs from the start of the GPS era on 1/1/1980) larger than utc_start_abstime, which makes sense to me. The leapsecond change reduces GPS time compared to UTC. It also forces the base time to be 00Z/01-01-1980, so if, for the record, I replace

    CALL TimeSince ( DT8, time, 1, Base="JS2000" )

by

    CALL TimeSince ( DT8, time, 1, Base="GPS" )

in ropp_io/ropp/ropp_io_write_ncdf_put.f90, I get:

start_time[0]=1159750580.67 ( ~ no. of secs since 01-01-1980)

instead of

start_time = 529030600.67 ( ~ no. of secs since 01-01-2000)

That's not what is wanted.

by Ian Culverwell, 8 years ago

GRAS_1B_M02_20161006005636Z_20161006005939Z_N_O_20161006020703Z_G17_NN.nc

comment:13 by Ian Culverwell, 8 years ago

Axel's file attached.

comment:14 by Ian Culverwell, 8 years ago

Cc: ian.culverwell@… removed

Actions:

1) Liaise with Axel and Stig to check that r5072 does what they want.

2) Check that the effects on the ropp_pp test files makes sense.

3) Revert, temporarily, r5072 at r5081.

comment:15 by Dave Offiler, 8 years ago

While there have indeed been 9 leap seconds added between 1972 and the start of GPS Time in 1980, according to Wikipedia[ref.1] (my source for the original LeapSeconds() code), there was already a 10 second difference before 1972:

In 1972, the leap-second system was introduced so that the broadcast UTC seconds could be made exactly equal to the standard SI second, while still maintaining the UTC time of day and changes of UTC date synchronized with those of UT1 (the solar time standard that superseded GMT).[10] By then, the UTC clock was already 10 seconds behind TAI, which had been synchronized with UT1 in 1958, but had been counting true SI seconds since then. After 1972, both clocks have been ticking in SI seconds, so the difference between their readouts at any time is 10 seconds plus the total number of leap seconds that have been applied to UTC (36 seconds in July 2015).

So I believe the the 10 secs fixed offset is correct. and not to be confused with the 9 seconds added leap secs between 1972 and 1980.

Ref1: https://en.wikipedia.org/wiki/Leap_second

comment:16 by Ian Culverwell, 8 years ago

Milestone: 9.010.0

What a mess. Defer to ROPP10.0 when I should have more strength to face it.

comment:17 by Ian Culverwell, 7 years ago

In connection (perhaps) with this, Axel von Engeln writes

I was just sending this wonderful background data that  is attached through ROPP
(using 9.0 version). And I  noted ROPP was doing terrible things with it. First,
it  complains about the "'start_time' and yr/mo/dy/hr/mn/sc/ms timestamps differ
by 1.000 seconds" and then corrects the start_time (from actually being correct
with leap second to being incorrect...), but the ropp_fm also seems to set the
start_time and the time variables to the same value (independent on the leap
second being there). Is that a bug or should it be like that?

The file to which he refers, EPS-SG_Example.bg.nc, is attached.

by Ian Culverwell, 7 years ago

Attachment: EPS-SG_Example.bg.nc added

EPS-SG_Example.bg.nc

comment:18 by Ian Culverwell, 7 years ago

Axel adds the following.

With some examples from our Python Julian seconds to UTC conversion:

>>> dt.js2utc(242903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2007-09-12 09:07:05.73' at f634f8>
>>> dt.js2utc(243903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2007-09-23 22:53:45.73' at f66150>
>>> dt.js2utc(242903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2007-09-12 09:07:05.73' at f66198>
>>> dt.js2utc(252903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-01-06 02:53:45.73' at f66108>
>>> dt.js2utc(262903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-04-30 20:40:25.73' at f660c0>
>>> dt.js2utc(272903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-08-24 14:27:05.73' at f66078>
>>> dt.js2utc(282903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-12-18 08:13:45.73' at f634f8>
>>> dt.js2utc(292903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2009-04-13 02:00:24.73' at f66198>
>>> dt.js2utc(302903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2009-08-06 19:47:04.73' at f66108>
>>> dt.js2utc(312903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2009-11-30 13:33:44.73' at f660c0>
>>> dt.js2utc(322903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2010-03-26 07:20:24.73' at f66078>
>>> dt.js2utc(332903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2010-07-20 01:07:04.73' at f634f8>
>>> dt.js2utc(342903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2010-11-12 18:53:44.73' at f66198>
>>> dt.js2utc(352903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2011-03-08 12:40:24.73' at f66108>
>>> dt.js2utc(402903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2012-10-07 05:33:43.73' at f66150>
>>> dt.js2utc(502903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2015-12-08 15:20:22.73' at f66078>
>>> dt.js2utc(552903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2017-07-09 08:13:41.73' at f660c0>

comment:19 by Ian Culverwell, 7 years ago

Milestone: 10.09.1

comment:20 by Ian Culverwell, 7 years ago

The situation is summarised here: leapsec_plot_p3.gif

by Ian Culverwell, 7 years ago

Attachment: leapsec_plot_p3.gif added

leapsec_plot_p3.gif

comment:21 by Ian Culverwell, 7 years ago

Before the first leap second was added at 0Z 01/07/1972, the three timescales (TAI2, GPS and UTC) were in agreement. (TAI was already adrift of UTC by 10 seconds; TAI2 = TAI - 10 secs takes account of this.) As soon as the first leap second was added, UTC started lagging behind TAI2. GPS time wasn't really defined, but if it had been, it would have matched UTC. By the time the GPS clock started, at 0Z 06/01/1980, TAI2 time was 9 leap seconds ahead of UTC/GPS time. Like TAI2, GPS takes no account of leap seconds, and therefore remains 9 seconds behind TAI2 (i.e. 19 seconds behind TAI), while UTC continues to fall behind both of them as further leap seconds accumulate.

Three tools to facilitate the conversion between GPS and UTC times have been written and committed at r5471.

  • ropp_utils/datetime/leapseconds.f90 calculates the total number of leapseconds that have accrued up to a given UTC. In other words, it provides xU-xT as a function of yU. Thus, for example, leapseconds(0Z 01/01/2000) = 9 + 13 = 22.
  • ropp_utils/datetime/gpstoutc.f90 converts between the two date formats, i.e. it maps yU to yG and yG to yU. For example, gpstoutc returns a GPS datetime of 00:00:13 01/12/2005 when given an input UTC datetime of 00:00:00 01/12/2005, showing that the 22-9 = 13 leapseconds difference between the two clocks at this date have been correctly accounted for. When this GPS time is passed through gpstoutc in the other direction, it correctly returns the original UTC of 00:00:00 01/12/2005.

There remains a slight infelicity when dealing with datetimes close to the leapsecond insertion times. For example, 23:59:60 UTC 31/12/2005 returns 00:00:14 GPS 01/01/2006, rather than 00:00:13 GPS 01/01/2006. (This is because 23:59:60 UTC 31/12/2005 is equivalent to 00:00:00 UTC 01/01/2006 as far as CalToJul is concerned, and by that time 14 leapseconds have been added.) Similarly, 00:00:13 GPS 01/01/2006 maps to 23:59:59 UTC 31/12/2005, rather than 23:59:60 UTC 31/12/2005 (which would fail range-checking anyway). The computational cost of encoding checks and fixes for this behaviour outweighs the marginal benefit of solving a problem that will be encountered only very rarely - historically, for about 17 seconds in 37 years, or about 14 times per billion. This decision could be revisited according to users' wishes.

  • ropp_utils/datetime/augmentutc.f90 adds a specified number of days to a given UTC and returns the resulting UTC. In other words, it maps from yU(xU) to yU(xU+jday). It therefore emulates Axel's Python tool dt.js2utc. Thus, for example, augmentutc(0Z 01/01/2000, 242903226.727/86400.d0) = 09:07:05.727 12/09/20107, which matches Axel's first example.

comment:22 by Ian Culverwell, 7 years ago

Here is some evidence that these tools are doing the right job.

leapseconds.f90 See the above picture, as well as:

        UTC                     leapseconds
----------------------------------------------
 1980 01  1   0  0   0  0   0        9
 2000 01  1   0  0   0  0   0       22
 2016 10  6   0  0  56 36 670       26

comment:23 by Ian Culverwell, 7 years ago

gpstoutc.f90

Around the turn of a non-leapsecond year:

  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:12.000 GPS 01/01/2001 **** 
  ------------------------------------------ 
GPS =   2001     1     1     0     0     0    12     0
 Converting this to UTC
UTC =   2000    12    31     0    23    59    59     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    12     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:13.000 GPS 01/01/2001 **** 
  ------------------------------------------ 
GPS =   2001     1     1     0     0     0    13     0
 Converting this to UTC
UTC =   2001     1     1     0     0     0     0     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    13     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:14.000 GPS 01/01/2001 **** 
  ------------------------------------------ 
GPS =   2001     1     1     0     0     0    14     0
 Converting this to UTC
UTC =   2001     1     1     0     0     0     1     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    14     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 23:59:59.000 UTC 31/12/2000 **** 
  ------------------------------------------ 
UTC =   2000    12    31     0    23    59    59     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    12     0
 Converting this to UTC
UTC =   2000    12    31     0    23    59    59     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 23:59:60.000 UTC 31/12/2000 **** 
  ------------------------------------------ 
UTC =   2000    12    31     0    23    59    60     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    13     0
 Converting this to UTC
UTC =   2001     1     1     0     0     0     0     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:00.000 UTC 01/01/2001 **** 
  ------------------------------------------ 
UTC =   2001     1     1     0     0     0     0     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    13     0
 Converting this to UTC
UTC =   2001     1     1     0     0     0     0     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:01.000 UTC 01/01/2001 **** 
  ------------------------------------------ 
UTC =   2001     1     1     0     0     0     1     0
 Converting this to GPS
GPS =   2001     1     1     0     0     0    14     0
 Converting this to UTC
UTC =   2001     1     1     0     0     0     1     0

We see that the mappings work in both directions. ( UTC = 2000 12 31 0 23 59 60 is equivalent to UTC = 2001 1 1 0 0 0 0 as far as CalToJul is concerned.)

Around the turn of a leapsecond year:

  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:12.000 GPS 01/01/2006 **** 
  ------------------------------------------ 
GPS =   2006     1     1     0     0     0    12     0
 Converting this to UTC
UTC =   2005    12    31     0    23    59    59     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    12     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:13.000 GPS 01/01/2006 **** 
  ------------------------------------------ 
GPS =   2006     1     1     0     0     0    13     0
 Converting this to UTC
UTC =   2005    12    31     0    23    59    59     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    12     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:14.000 GPS 01/01/2006 **** 
  ------------------------------------------ 
GPS =   2006     1     1     0     0     0    14     0
 Converting this to UTC
UTC =   2006     1     1     0     0     0     0     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    14     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 23:59:59.000 UTC 31/12/2005 **** 
  ------------------------------------------ 
UTC =   2005    12    31     0    23    59    59     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    12     0
 Converting this to UTC
UTC =   2005    12    31     0    23    59    59     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 23:59:60.000 UTC 31/12/2005 **** 
  ------------------------------------------ 
UTC =   2005    12    31     0    23    59    60     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    14     0
 Converting this to UTC
UTC =   2006     1     1     0     0     0     0     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:00.000 UTC 01/01/2006 **** 
  ------------------------------------------ 
UTC =   2006     1     1     0     0     0     0     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    14     0
 Converting this to UTC
UTC =   2006     1     1     0     0     0     0     0
  ------------------------------------------ 
  **** Testing GPStoUTC on 00:00:01.000 UTC 01/01/2006 **** 
  ------------------------------------------ 
UTC =   2006     1     1     0     0     0     1     0
 Converting this to GPS
GPS =   2006     1     1     0     0     0    15     0
 Converting this to UTC
UTC =   2006     1     1     0     0     0     1     0

We see the two minor problems mentioned above, but otherwise all is OK.

For Axel's test file GRAS_1B_M02_20161006005636Z_20161006005939Z_N_O_20161006020703Z_G17_NN.nc, we have

 ------------------------------------------ 
  **** Testing GPStoUTC on 00:56:36.670 UTC 06/10/2016 **** 
  ------------------------------------------ 
UTC =   2016    10     6     0     0    56    36   670
 Converting this to GPS
GPS =   2016    10     6     0     0    56    53   670
 Converting this to UTC
UTC =   2016    10     6     0     0    56    36   670

The 17 seconds difference is as expected.

We can also confirm that TimeSince takes no account of leapseconds, because (5*366 + 12*365 + 279)*86400 + (56*60 + 53.670) = 529030613.67 clock seconds have elapsed between 0Z/06/01/1980 and 00:56:53.670 GPS 06/10/2016 = 00:56:36.670 UTC 06/10/2016.

comment:24 by Ian Culverwell, 7 years ago

augmentutc.f90

As well as the example in the introduction, we have:

jsec, UTC = 242903226.7269999981  2007     9    12     0     9     7     5   727
jsec, UTC = 243903226.7269999981  2007     9    23     0    22    53    45   727
jsec, UTC = 242903226.7269999981  2007     9    12     0     9     7     5   727
jsec, UTC = 252903226.7269999981  2008     1     6     0     2    53    45   727
jsec, UTC = 262903226.7269999981  2008     4    30     0    20    40    25   727
jsec, UTC = 272903226.7269999981  2008     8    24     0    14    27     5   728
jsec, UTC = 282903226.7269999981  2008    12    18     0     8    13    45   728
jsec, UTC = 292903226.7269999981  2009     4    13     0     2     0    24   727
jsec, UTC = 302903226.7269999981  2009     8     6     0    19    47     4   727
jsec, UTC = 312903226.7269999981  2009    11    30     0    13    33    44   727
jsec, UTC = 322903226.7269999981  2010     3    26     0     7    20    24   727
jsec, UTC = 332903226.7269999981  2010     7    20     0     1     7     4   727
jsec, UTC = 342903226.7269999981  2010    11    12     0    18    53    44   728
jsec, UTC = 352903226.7269999981  2011     3     8     0    12    40    24   728
jsec, UTC = 402903226.7269999981  2012    10     7     0     5    33    43   727
jsec, UTC = 502903226.7269999981  2015    12     8     0    15    20    22   727
jsec, UTC = 552903226.7269999981  2017     7     9     0     8    13    41   727

These match the values in Axel's table, reproduced here for convenience:

>>> dt.js2utc(242903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2007-09-12 09:07:05.73' at f634f8>
>>> dt.js2utc(243903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2007-09-23 22:53:45.73' at f66150>
>>> dt.js2utc(242903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2007-09-12 09:07:05.73' at f66198>
>>> dt.js2utc(252903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-01-06 02:53:45.73' at f66108>
>>> dt.js2utc(262903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-04-30 20:40:25.73' at f660c0>
>>> dt.js2utc(272903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-08-24 14:27:05.73' at f66078>
>>> dt.js2utc(282903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2008-12-18 08:13:45.73' at f634f8>
>>> dt.js2utc(292903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2009-04-13 02:00:24.73' at f66198>
>>> dt.js2utc(302903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2009-08-06 19:47:04.73' at f66108>
>>> dt.js2utc(312903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2009-11-30 13:33:44.73' at f660c0>
>>> dt.js2utc(322903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2010-03-26 07:20:24.73' at f66078>
>>> dt.js2utc(332903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2010-07-20 01:07:04.73' at f634f8>
>>> dt.js2utc(342903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2010-11-12 18:53:44.73' at f66198>
>>> dt.js2utc(352903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2011-03-08 12:40:24.73' at f66108>
>>> dt.js2utc(402903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2012-10-07 05:33:43.73' at f66150>
>>> dt.js2utc(502903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2015-12-08 15:20:22.73' at f66078>
>>> dt.js2utc(552903226.727, ref = dt.DateTime(2000, 1, 1, 0, 0, 0.)) <mx.DateTime.DateTime object for '2017-07-09 08:13:41.73' at f660c0>

(There is evidently an occasional +/- 1 msec wobble, due to rounding error, which we can live with.)

comment:25 by Ian Culverwell, 7 years ago

Brief email discussion with Axel confirms that he wants the start_time increased to account for the extra leapseconds:

On your last question, yes I would be happy if ROPP tells 
me it is 529030600.67 seconds since 2000.

He also confirms that the gps_start_absdate/time are apparently in error (or inconsistent, at least), and that Chris M will look at it in the PFS.

comment:26 by Ian Culverwell, 7 years ago

Changes to ropp_io_write_ncdf_put.f90check out OK on the above example:

ncks -H -Q -vstart_time,time_offset data/GRAS_1B_M02_20161006005636Z_20161006005939Z_N_O_20161006020703Z_G17_NN_test1.nc
start_time[0]=529030600.67

cf ROPP9.0 control:

ncks -H -Q -vstart_time,time_offset data/GRAS_1B_M02_20161006005636Z_20161006005939Z_N_O_20161006020703Z_G17_NN_cntl1.nc
start_time[0]=529030596.67

Note that we have to modify ropp_io_read_ncdf_get.f90 too, because it compares start_time and DTocc, and they are no longer consistent unless leapseconds are taken into account. This means that when earlier files, which have incorrect start_times, are read, they will now provoke a warning that the start time and the UTC date stamp are inconsistent, because they differ by the number of leapseconds that have accrued since 0Z 01/01/2000. In such cases ROPP will use the UTC stamp, which is still correct, so the warning may be safely ignored.

comment:27 by Ian Culverwell, 7 years ago

Code changes submitted at r5477.

ROPP IO UG changes submitted at r5478.

Change Log updated at r5479.

comment:28 by Ian Culverwell, 7 years ago

Axel was also complaining that the ROPP9.0 forward model was screwing up his file EPS-SG_Example.bg.nc. It certainly is (just looking at the first profile):

ncks -H -Q -vstart_time,time,time_offset data/EPS-SG_Example.bg_prof1.nc
start_time[0]=242901664.901 

time[0]=242901762.227 

time_offset[0]=97.326 


-----------------------------------------------------------------------
                     ROPP Forward Model
-----------------------------------------------------------------------


INFO (from ropp_fm_bg2ro_1d):  Processing profile    1 of      1
 
WARNING (from ropp_io_read_ncdf_get):  'start_time' and yr/mo/dy/hr/mn/sc/ms timestamps differ by 1.000 seconds - using yr/../ms timestamp
INFO (from ropp_fm_bg2ro_1d):  (BG_20070912084103_MSA1_G002_EUME) 

... (from ropp_io_ascend):  Level2b (geop,temp,shum,pres) data found to be in reverse order (descending height order) 
  Reversing data direction.
... (from ropp_io_ascend):  Note 1dVar users must ensure that the error description is appropriate for profiles in increasing height order. 
 See User Guide for more details 

ncks -H -Q -vstart_time,time,time_offset data/EPS-SG_Example.bg_prof1_cntl1.nc
start_time[0]=242901663.901 

time[0]=242901663.901 

time_offset[0]=97.326 

/data/users/idculv/ROPP/ropp_src/branches/dev/Share/ROPP90_prototype/ropp_fm/tools/../../ropp_io/tools/ropp2ropp data/EPS-SG_Example.bg_prof1_cntl1.nc -d -o data/EPS-SG_Example.bg_prof1_cntl2.nc

---------------------------------------------------------------------
                ROPP-to-ROPP generic netCDF tool
---------------------------------------------------------------------

INFO (from ropp2ropp):  Reading data/EPS-SG_Example.bg_prof1_cntl1.nc
... (from ropp_io_thin):  Thinning method:  [v3.1]
... (from ropp2ropp):  Ensuring all profiles are in ascending height order
INFO (from ropp2ropp):  Profile    1 : BG_20070912084103_MSA1_G002_EUME
... (from ropp2ropp):  FmtVer:  ROPP I/O V1.1
... (from ropp2ropp):  OccID :  BG_20070912084103_MSA1_G002_EUME
... (from ropp2ropp):  LEOID :  MSA1
... (from ropp2ropp):  GNSID :  G002
... (from ropp2ropp):  STNID :
... (from ropp2ropp):  ProCen:  EUMETSAT
... (from ropp2ropp):  ProS/W:  UNKNOWN
... (from ropp2ropp):  DTocc :  08:41:03.901 12-Sep-2007
... (from ropp2ropp):  DTpro :  14:11:31.000 28-Aug-2017
... (from ropp2ropp):  GEOref:     97.326    28.5    10.8
... (from ropp2ropp):           6352868.0 26017.7  4949.5 -6011.9     0.0    28.3
... (from ropp2ropp):             -9999.9   -9999.9   -9999.9   -9999.9   -9999.9   -9999.9
... (from ropp2ropp):             -9999.9   -9999.9   -9999.9   -9999.9   -9999.9   -9999.9
... (from ropp2ropp):  PCD   :  0100000000000110
... (from ropp2ropp):  Qual  :    100.0
... (from ropp2ropp):  PODmth:  Not applicable
... (from ropp2ropp):  LVLmth:  ROPP forward ROPP forward Occultation Simulator and IDL based ECMW
... (from ropp2ropp):  THNmth:   [v3.1]
... (from ropp2ropp):  SW-ID :  UNKNOWN
... (from ropp2ropp):  B/G   :  ECMWF ANA-SG, CAT: 3  06:00 12-Sep-2007
... (from ropp2ropp):  ModLev:  HYBRID ECMWF
... (from ropp2ropp):  Levels:    L1a  L1b  L2a  L2b  L2c  L2d
... (from ropp2ropp):  Npts  :      0  923  300  137    1  138
... (from ropp2ropp):  Missng:      T    F    F    F    F    F
INFO (from ropp2ropp):  Writing data/EPS-SG_Example.bg_prof1_cntl2.nc

ncks -H -Q -vstart_time,time,time_offset data/EPS-SG_Example.bg_prof1_cntl2.nc
start_time[0]=242901663.901 

time[0]=242901663.901 

time_offset[0]=97.326 

One leapsec was introduced at the end of 31/12/2005, and this means that the given UTC of 08:41:03.901 12-Sep-2007 is (2*366 + 5*365 + 31+28+31+30+31+30+31+31+11)*86400.0 + 8*3600 + 41*60 + 3.901 + 1 = 242901664.9009999931 secs after 0Z 01/01/2000.

This is the start_time in the original file, which ROPP9.0 misinterprets as being 1 sec out of synch with DTocc.

This file also has time = start_time + time_offset, which makes sense. ROPP keeps the time_offset of 97.326 sec, but sets start_time and time to be the same. Why?

From the ROPP input files:

	double start_time(dim_unlim) ;
		start_time:long_name = "Starting time for the occultation" ;
		start_time:units = "seconds since 2000-01-01 00:00:00" ;

	double time(dim_unlim) ;
		time:long_name = "Reference time for the occultation" ;
		time:units = "seconds since 2000-01-01 00:00:00" ;

	float time_offset(dim_unlim) ;
		time_offset:long_name = "Time offset for georeferencing (since start of occ.)" ;
		time_offset:units = "seconds" ;
		time_offset:valid_range = -100., 239.999;

So time = start_time + time_offset looks correct to me.

As far as I can see, time is read but never used in ROPP, so we could safely set it to a more sensible value.

Putting in a little fix to this effect at r5480, we get this result when we pass the same file through ROPP9.1:

ncks -H -Q -vstart_time,time,time_offset data/EPS-SG_Example.bg_prof1.nc
start_time[0]=242901664.901 

time[0]=242901762.227 

time_offset[0]=97.326 


-----------------------------------------------------------------------
                     ROPP Forward Model
-----------------------------------------------------------------------


INFO (from ropp_fm_bg2ro_1d):  Processing profile    1 of      1


INFO (from ropp_fm_bg2ro_1d):  (BG_20070912084103_MSA1_G002_EUME) 

... (from ropp_io_ascend):  Level2b (geop,temp,shum,pres) data found to be in reverse order (descending height order) 
  Reversing data direction.
... (from ropp_io_ascend):  Note 1dVar users must ensure that the error description is appropriate for profiles in increasing height order. 
 See User Guide for more details 

ncks -H -Q -vstart_time,time,time_offset data/EPS-SG_Example.bg_prof1_test1.nc
start_time[0]=242901664.901 

time[0]=242901762.227 

time_offset[0]=97.326 

/data/users/idculv/ROPP/ropp_src/branches/dev/Share/ic_utcgps/ropp_fm/tools/../../ropp_io/tools/ropp2ropp data/EPS-SG_Example.bg_prof1_test1.nc -d -o data/EPS-SG_Example.bg_prof1_test2.nc

---------------------------------------------------------------------
                ROPP-to-ROPP generic netCDF tool
---------------------------------------------------------------------

INFO (from ropp2ropp):  Reading data/EPS-SG_Example.bg_prof1_test1.nc
... (from ropp_io_thin):  Thinning method:  [v3.1]
... (from ropp2ropp):  Ensuring all profiles are in ascending height order
INFO (from ropp2ropp):  Profile    1 : BG_20070912084103_MSA1_G002_EUME
... (from ropp2ropp):  FmtVer:  ROPP I/O V1.1
... (from ropp2ropp):  OccID :  BG_20070912084103_MSA1_G002_EUME
... (from ropp2ropp):  LEOID :  MSA1
... (from ropp2ropp):  GNSID :  G002
... (from ropp2ropp):  STNID :
... (from ropp2ropp):  ProCen:  EUMETSAT
... (from ropp2ropp):  ProS/W:  UNKNOWN
... (from ropp2ropp):  DTocc :  08:41:03.901 12-Sep-2007
... (from ropp2ropp):  DTpro :  14:11:31.000 28-Aug-2017
... (from ropp2ropp):  GEOref:     97.326    28.5    10.8
... (from ropp2ropp):           6352868.0 26017.7  4949.5 -6011.9     0.0    28.3
... (from ropp2ropp):             -9999.9   -9999.9   -9999.9   -9999.9   -9999.9   -9999.9
... (from ropp2ropp):             -9999.9   -9999.9   -9999.9   -9999.9   -9999.9   -9999.9
... (from ropp2ropp):  PCD   :  0100000000000110
... (from ropp2ropp):  Qual  :    100.0
... (from ropp2ropp):  PODmth:  Not applicable
... (from ropp2ropp):  LVLmth:  ROPP forward ROPP forward Occultation Simulator and IDL based ECMW
... (from ropp2ropp):  THNmth:   [v3.1]
... (from ropp2ropp):  SW-ID :  UNKNOWN
... (from ropp2ropp):  B/G   :  ECMWF ANA-SG, CAT: 3  06:00 12-Sep-2007
... (from ropp2ropp):  ModLev:  HYBRID ECMWF
... (from ropp2ropp):  Levels:    L1a  L1b  L2a  L2b  L2c  L2d
... (from ropp2ropp):  Npts  :      0  923  300  137    1  138
... (from ropp2ropp):  Missng:      T    F    F    F    F    F
INFO (from ropp2ropp):  Writing data/EPS-SG_Example.bg_prof1_test2.nc

ncks -H -Q -vstart_time,time,time_offset data/EPS-SG_Example.bg_prof1_test2.nc
start_time[0]=242901664.901 

time[0]=242901762.227 

time_offset[0]=97.326

We see that the start_time, time and time_offset of the original dataset are preserved. And there are no warnings about mismatched start_time and Dtocc.

comment:29 by Ian Culverwell, 7 years ago

I/O UG changes at r5481, and updated Change Log at r5482.

comment:30 by Ian Culverwell, 7 years ago

Updates to the reference files, to avoid those annoying

WARNING (from ropp_io_read_ncdf_get):  'start_time' and yr/mo/dy/hr/mn/sc/ms timestamps differ by 2.000 seconds - using yr/../ms timestamp

messages, were committed at r5483.

comment:31 by Ian Culverwell, 6 years ago

Resolution: fixed
Status: reopenedclosed

Made some further updates to warning messages, so that they are only triggered if the time difference is > 30 secs, and only then if msg_MODE = VerboseMode:

  IF (isroppinrange(data%DTocc)) THEN
    IF ( isinrange(start_time, (/ 0.001_dp, now_time /)) ) 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, dtocc_time, 1, Base="JS2000")
      dtocc_time = dtocc_time + (LeapSeconds(DT8) - 22)*1.D0 ! 22 leapseconds had elapsed by 0Z 01/01/2000
      time = ABS(start_time - dtocc_time)
      IF ( time > 0.0005_dp ) THEN
        WRITE ( cval, FMT="(F15.3)" ) time
        IF ( time < 30.0_dp ) THEN ! Probably a result of not accounting for leap-seconds, so reduce severity of message to DIAG
          CALL message(msg_diag, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // &
                       "timestamps differ by " // &
                       TRIM(ADJUSTL(cval))//" seconds (probably a leap-second issue) - using yr/../ms timestamp")
        ELSE
          CALL message(msg_warn, "'start_time' and yr/mo/dy/hr/mn/sc/ms " // &
                       "timestamps differ by " // &
                       TRIM(ADJUSTL(cval))//" seconds - using yr/../ms timestamp")
        END IF
      END IF
    END IF
...
  END IF

comment:32 by Ian Culverwell, 6 years ago

Resolution: fixed
Status: closedreopened

Oops, prematurely closed.

comment:33 by Ian Culverwell, 6 years ago

Resolution: fixed
Status: reopenedclosed

This now gets through the test folder without too much bellyaching, so close the ticket.

comment:34 by Ian Culverwell, 6 years ago

(For the record, the ic_utcgps mods were merged into ROPP91_prototype at r5624.)

comment:35 by Ian Culverwell, 6 years ago

A rather better figure, constructed for the ROPP9.1 Change Log (in the absence of a ROPP_UTILS User Guide to describe the tools) is this: leapsec_plot_p1.gif

by Ian Culverwell, 6 years ago

Attachment: leapsec_plot_p1.gif added

leapsec_plot_p1.gif

Note: See TracTickets for help on using tickets.