timesince.f90

[Ian Culverwell]

There were some discussions about this routine some time ago.

From: Culverwell, Ian
Sent: 11 October 2011 12:28
To: 'Stig Syndergaard'
Cc: Hallgeir Wilhelmsen; Kjartan Kinch; Culverwell, Ian; Offiler, Dave
Subject: RE: NOT URGENT FW: GRAS SAF Helpdesk Request w.r.t. ROPP
Software, [Answered]

Hi Stig,

Dave has asked me to look into this.



1) I didn't know about the UTC/GPS difference. As you say, it could be important. 

As far as I can see, timsince is only called in 2 places where the GPS/UTC anomaly might be important:

ropp_pp_grasrs2ropp, which uses start_time_absdate in the gras RS file to calculate dates and times relative to a datum of "Days since 0001-01-01 00:00:00.00"; and

ropp_pp_preprocess_cosmic.f90, which calculates an occulation start time using a datum of "GPSSEC", ie 00:00:00 6/1/1980.

I suppose the 1st is OK. The second is used to calculate lost carrier flags.  I don’t know whether the mismatch between this and true GPS time is important. (But if it is, why haven't we noticed before?) If so, it wouldn't be difficult to incorporate in the datetime routines. 

All the other uses of timesince refer to a datum of JS200 or JM2000, ie they use UTC.



2) The +1msec error is annoying.  It's partly a precision and partly an algorithm issue.  

The root of the problem is that CalToJul calculates msecs as NINT(msecs past last full sec + 0.5).  Because of rounding errors in the calculations, in your example "msecs past last full sec" = 0.002407550811768 (instead of zero), and therefore NINT rounds up to 1. 

If you shave 0.01msec off your Tsince (ie, Tsince=990919156_wp - 0.00001_wp), CalToJul ends up calculating msec=NINT(0.462174415588379) = 0, and all is well. (Note that the change in the argument to NINT is more than 0.01ms because things get multiplied and divided by 24*60*... en route.)

But even if you're silly enough to work in quadruple precision, to reduce this effect (as I was) you still find a problem (with the original Tsince=990919156_wp) because CalToJul then calculates NINT(msecs past last hour + 0.5) = NINT(0.0+0.5) = NINT(0.5), which is still rounded up to 1.

Both problems can be avoided by calculating msec=NINT(msecs past last full sec + 0.45), which equals 0 in both cases. 

It's not pretty, but the alternative is some fairly major rewriting to prevent rounding errors building up. And anyway, the +0.5 offset should probably be reduced a bit (eg to 0.499999) to prevent 0msec appearing as 1.



3) I think the DOY=150 is as intended. There's a difference between dates and times, isn't there?  When we say month=5 we mean we're in May, not that 5 full months have elapsed and that we're in June. But for times it's the other way round: at 23:19:16 we're in the 24th hour of the day, not the 23rd, but we say HH=23 because 23 full hours have elapsed. It seems to me that you could make a case for the DOY to be either 150 or 151 in this case. As you showed, Timesince chooses 150.

Hope this helps.  Is this a subject suitable for discussion at PT9?

Regards,
Ian.


-----Original Message-----
From: Stig Syndergaard [mailto:ssy@dmi.dk]
Sent: 03 October 2011 18:06
To: Offiler, Dave
Cc: Culverwell, Ian; Hallgeir Wilhelmsen; Kjartan Kinch
Subject: Re: NOT URGENT FW: GRAS SAF Helpdesk Request w.r.t. ROPP Software, [Answered]

Hi Dave,

I see now how clever timesince.f90 really is! I think it could be useful to us here. However, a few issues that I realized when trying it out:

1) It doesn't really convert between UTC time and GPS seconds. UTC and GPS time are currently 15 secs apart, e.g., 2011-10-04 00:00:10 GPS time (which corresponds to 1001721610 GPS seconds since Jan 6 1980 00:00:00) would be 2011-10-03 23:59:55 UTC. The code does not seem to account for the 15 secs difference. I don't know if this is an issue anywhere in ROPP (e.g., when aligning external navbits with the measurements; Kjartan may be able to chime in here). CDAAC uses solely GPS time, which means that their dates and times in file names etc. are not in UTC. But GRAS data dates and times are in UTC, and soon we will also process COSMIC data, so we need to be aware of the difference, and a correct conversion could potentially be useful. To do a correct conversion, one would need a 'living' table with the leap seconds, since a new one is added now and then. Leap seconds can be introduced at the end of June or December, and it is decided about half a year in advance by the International Earth Rotation and Reference System Service (IERS).
(see e.g., http://www.timeanddate.com/time/leapseconds.html)

2) I tried the attached little code and got:

        2011           5          31           0          23          19
           16           1         150

But I was expecting:

        2011           5          31           0          23          19
           16           0         151

Thus, 1 msec off and 1 DOY off. I think the two are unrelated, but I couldn't quite figure out the reason for the 1 msec error. On the DOY, it seems that one would have to do int(Tsince) + 1. Is that intended?

Cheers,
-Stig

On 2011-10-03 15:35, Offiler, Dave wrote:
> Stig/Ian,
>
> The timesince() routine is intended to do all these conversions, so 
> individual routines like gpsseconds(), dayofyear() and julian*(), etc, 
> were removed as redundant. (If you look at the old code, they all did 
> much the same thing with a different base date/time, some returning 
> seconds, some days. As they were all written outside ROPP (ie PES), by 
> different people in different coding styles, they were a bit of a 
> mess. timesince() just brings them all under one 'roof'. The Julian 
> Date algorithm now uses a much more efficient implementation.)
>
> E.g. to convert UTC in year/month/day/hr/min/sec/ms form to 
> GPS-Seconds ('zone' = 0):
>
>     USE  DateTime
>     CHARACTER (LEN=40) :: Tbase = 'GPSSECONDS' !(or abbreviated down 
> to just 'GPS')
>     INTEGER            :: CDT(8) = (/<year>,<month>,<day>,<hour>,
> <zone>,<minute>,<seconds>,<msecs>  /)
>     REAL(dp)           :: Tsince
>     INTEGER            :: inv = 1
>     ...
>     CALL TimeSince ( CDT, Tsince, inv, Base=Tbase )
>
> or explicitly:  ... Base='seconds since 6-Jan-1980 00:00:00', being 
> the starting datum for GPS time.
> Or use a different datum if you don't believe in this one!
>
> Similarly, DayOfYear can be obtained with, e.g.:
>
>     WRITE(Tbase,'(A,I4)') 'JD', CDT(1)   ! Tbase = 'JD<year>'
>     CALL TimeSince ( CDT, Tsince, 1, Base=Tbase)
>     iDOY = INT(Tsince)
>
> The 'absolute' or 'modified' Julian Date can similarly be obtained 
> simply by giving the Base string explicitly, or one of the short-cuts 
> like 'JD2000'. See the header comments in 
> ropp_io/datetime/timesince.f90
>
> You can reverse the conversion by setting Tsince and inv = -1; the 
> converted UTC calendar date/time will be returned in the CDT array.
>
> There's also the datetimeoffset() routine for adding or subtracting 
> calendar date/times. This just uses timesince() to get absolute times, 
> does the maths and converts back.
>
> If there's a conversion you need that timesince() can't do in the 
> present version, please let me have a spec. and I can add it.
>
> HTH
> Dave
>
> --
> David Offiler   Head Satellite Active Sensing Group
> Met Office   FitzRoy Road   Exeter   EX1 3PB   United Kingdom
> Tel: +44 (0)1392 88 6298    Fax: +44 (0)1392 88 5681
> E-mail: dave.offiler@metoffice.gov.uk   http://www.metoffice.gov.uk
>
> -----Original Message-----
> From: Culverwell, Ian
> Sent: 26 September 2011 12:21
> To: Offiler, Dave
> Cc: Culverwell, Ian
> Subject: NOT URGENT FW: GRAS SAF Helpdesk Request w.r.t. ROPP 
> Software, [Answered]
>
>
> Hi Dave,
>
> Just forwarding now, so I don't forget.  ***No need for action until 
> you return to work.***
>
> Enjoy your holiday!
>
> Ian.
>
> -----Original Message-----
> From: Culverwell, Ian
> Sent: 26 September 2011 12:13
> To: 'Stig Syndergaard'
> Cc: Culverwell, Ian
> Subject: RE: GRAS SAF Helpdesk Request w.r.t. ROPP Software, 
> [Answered]
>
>
> Hi Stig,
>
> OK - we think Dave did this (he's listed as author of timesince.f90), 
> and he's not back till next week, but I'll ask him then.
>
> Ian.
>
>
> -----Original Message-----
> From: Stig Syndergaard [mailto:ssy@dmi.dk]
> Sent: 26 September 2011 11:41
> To: Culverwell, Ian
> Subject: Re: GRAS SAF Helpdesk Request w.r.t. ROPP Software, 
> [Answered]
>
> Hi Ian,
>
> I was hoping for a routine that could convert between UTC time and GPS 
> time, but it seems none of timesince.f90 or the old GpsSecond.f90 
> actually does that. Another old subroutine that I think could be 
> useful to us would be DayOfYear.f90. Lets hear what the people you've 
> asked says.
>
> Thanks,
> -Stig
>
> On 2011-09-26 10:23, Culverwell, Ian wrote:
>>
>> Hi Stig,
>>
>> Thanks for drawing my attention to this.
>>
>> I have replied to Glen Carl, saying that 
>> ropp_utils/datetime/timesince.f90, which _is_ in the tarball distro, 
>> should do what he (and you) want.
>>
>> As for _why_ these were changed in going from ropp4.0 to 4.1: the 
>> changelog simply says
>>
>> "The datetime utility library has been reviewed and updated to 
>> improve
>
>> efficiency. Where necessary in other modules (ropp_io, ropp_fm, 
>> ropp_pp), calls to the relevant datetime functions have been
> modified."
>>
>> That's as much as I know.  I've asked those who made the change for 
>> their recollections, and willl pass to you when I receive them.
>>
>> Regards,
>> Ian.
>>
>>
>>
>> -----Original Message-----
>> From: Stig Syndergaard [mailto:ssy@dmi.dk]
>> Sent: 23 September 2011 16:54
>> To: Culverwell, Ian
>> Subject: Re: GRAS SAF Helpdesk Request w.r.t. ROPP Software, 
>> [Answered]
>>
>> Hi Ian,
>>
>> Because I saw the helpdesk question and didn't know of the 
>> GpsSecond.f90 routine before (there are many things I don't know in 
>> ROPP, but I'm learning), I decided to try it out because it could be 
>> useful for us when starting to process CDAAC data (which gives times 
>> in GPS seconds, whereas for GRAS the time is in UTC). But it seems 
>> that GpsSecond.f90 doesn't work with 5.0. The interface is not in 
>> datetimeprogs.f90, and it is not in the makefiles either. It also 
>> calls the function IsLeapYear.f90 and the sub DayOfYear.f90 which are 
>> not in 5.0. I see in 4.0 a few of other datetime routines that are no 
>> longer in 5.0. Do you know why have they been removed?
>>
>> Regards,
>> -Stig
>>
>> On 2011-09-23 09:13, helpdesk@grassaf.org wrote:
>>> The GRAS SAF Helpdesk has received a question from Glen CARL
>>>
>>> The helpdesk subject is:
>>>
>>> ropp_util juliansecond change for upgrade to version 5
>>>
>>> The helpdesk description is:
>>>
>>> I am upgrading to ROPP 5, and I noticed juliansecond has changed in
>> time utilities. What has this function been changed to, and where is 
>> it located?
>>>
>>> I will keep searching, however, I am hoping you have a quick 
>>> reference
>> for the time utility changes that were made for the upgrade to 
>> version
>
>> 5.
>>> Thanks,
>>> Glen
>>>
>>>
>>> Answer:
>>>     Dear Glen,
>>>
>>> Thank you for your interest in ROPP.
>>>
>>> Juliansecond disappeared between ropp4.0 and ropp4.1, in fact. There
>> is a new routine, ropp_utils / datetime / GpsSecond.f90, which 
>> appears
>
>> to do the equivalent calculation, although it uses a different datum
>> time: 0:00 on 6 Jan 1980, rather than 0:00 on 1 Jan 2000.  As far as 
>> I
>
>> can see, the difference in the output is just those 7300 days * 86400 
>> secs.
>>>
>>> Best Regards,
>>> Ian Culverwell, on behalf of the ROPP Development Team.
>>>
>>> This message is sent to the following email addresses:
>>> oliveras@ieec.fcr.es, grasadm@dmi.dk, grassaf@metoffice.gov.uk
>>>