Important: ROMS clock/date/calendar revisited

[arango]

I completely reworked the time/date/calendar in ROMS. All the associated routines were put it inside module dateclock.F. The routine get_date.F is replaced with dateclock.F module. This module contains several routines to manage date and calendars:

• caldate: Converts current model time (days) to calendar date. All the returned variables require keyword syntax since they are optional. For example, we can have:

       CALL caldate (tdays(ng), yd_r8=yday)
       CALL caldate (tdays(ng), yd_r8=yday, h_r8=hour)
  

  and so on. Check routine output arguments.

!  On Output:                                                          !
!                                                                      !
!     yy_i          Year including century (integer; OPTIONAL)         !
!     yd_i          Day of the year (integer; OPTIONAL)                !
!     mm_i          Month of the year, 1=Jan, ... (integer; OPTIONAL)  !
!     dd_i          Day of the month (integer; OPTIONAL)               !
!     h_i           Hour of the day (integer; OPTIONAL)                !
!     m_i           Minutes of the hour, 1 - 23 (integer; OPTIONAL)    !
!     s_i           Seconds of the minute (integer; OPTIONAL)          !
!                                                                      !
!     yd_r8         Day of the year (real, fraction; OPTIONAL)         !
!     dd_r8         Day of the month (real, fraction; OPTIONAL)        !
!     h_r8          Hour of the day (real, fraction; OPTION)           !
!     m_r8          Minutes of the hour (real, fraction; OPTION)       !
!     s_r8          Seconds of the minute (real, fraction; OPTIONAL)   !

• datenum: Converts requested date (year, month, day, ...) into a serial number. Similar to Matlab's datenum but the value datenum=0 corresponds to Mar 1, 0000.

• datevec: Converts a given date number to a date vector. It is inverse routine to datenum.

• day_code: Given (month, day, year) it returns a numerical code (0 to 6) for the day of the week.

• get_date: Retuns today date string of the form: DayOfWeak - Month day, year - hh:mm:ss ?M:

  Wednesday - April 19, 2017 -  3:33:09 PM
  

• ref_clock: Sets application time clock/reference and loads its to structure Rclock of TYPE T_CLOCK.

• ROMS_clock: Given (year, month, day, hour, minutes, seconds), this routine returns ROMS clock time since initialization from the reference date. It is used when importing fields from coupled models.

• time_string: Encodes current model time to a string. The time is reported to standard output is now of the form YYYY-MM-DD hh:mm:ss.ss

• yearday: Given (year,month,day) this integer function returns the day of the year.

---

The new routine datenum converts requested date (year, month, day, ...) into a serial date number. It is similar to Matlab's function datenum, except for the base date for day 0:

        Matlab:  datenum(0000,00,00)=0       reference date
                 datenum(0000,01,01)=1

          Here:  datenum(0000,03,01)=0       refecence date: Mar 1, 0000
                 datenum(0000,01,01)=-59

Matlab just has a different reference date for day 0. In October 1582, the Gregorian was adapted with a year length of 365.2425 days. This is coded as:

       365 + 1/4 - 1/100 + 1/400   or   365 + 0.25 - 0.01 + 0.0025

which is used to account for leap years. The base of Mar 1, 0000 is taken for simplicity since the length of February is not fixed.

Although this routine and the Matlab function are not equivalent, they yield the same elapsed time between two dates:

         datenum(2017,1,1)-datenum(0000, 3,1) = 736635
         datenum(2014,3,1)-datenum(1582,10,1) = 158302

In order to correct for the roundoff in the determination of the time (hh:mm:ss.ss), a ROUND function was added (round.F). This module is adapted from H.D. Knoble code (Penn State University) and it is based on a Fuzzy (UFLOOR) or Tolerant Floor function (TFLOOR).

   Usage:  CT is a comparison tolerance

        CT = 3 * EPSILON(X)     That is, CT is about 1 bit on either
                                side of X's mantissa bits.

        Y = ROUND(X, CT) = TFLOOR(X+0.5_r8,CT)

---

As before, the calendar in ROMS is determined by the value of standard input parameter time_ref:

• time_ref > 0, Proleptic Gregorian Calendar: time-units since YYYY-MM-DD hh:mm:ss

• time_ref = 0, Proleptic Gregorian Calendar: time-units since 0001-01-01 00:00:00

• time_ref = -1, 360_day Calendar: time-units since 0001-01-01 00:00:00

• time_ref = -2, Gregorian Calendar: time-units since Jan 1, 4713 BC

The Gregorian Calendar is basically the Julian Calendar corrected on 15 October 1582 when the length of the year became 365.2425 days instead of 365.25 days. The origin (day=0) correspond to Jan 1, 4713 BC. The one used in ROMS is truncated by NASA on May 23, 1968 (offset 2440000 days).

The Proleptic Gregorian Calendar extends backward the date preceding 15 October 1582 with a year length of 365.2425 days. The origin (day=0) correspond to March 1, 0000. This is the one managed with ROMS routined datenum and datevec. This is fine since in ROMS the reference date is time-units since 0001-01-01 00:00:00.

In routine caldate, we have:

!
!  The model clock is the elapsed time since reference time of the form
!  'time-units since YYYY-MM-DD hh:mm:ss'.  It is called the Gregorian
!  Proleptic Calendar.
!
      IF (INT(time_ref).gt.0) THEN
        DateNumber=Rclock.Dnumber+CurrentTime         ! fractional days
        DayFraction=DateNumber-AINT(DateNumber)
        HourFraction=DayFraction*24.0_r8
        MinuteFraction=(HourFraction-AINT(HourFraction))*60.0_r8
!
        IsDayUnits=.TRUE.
        CALL datevec (DateNumber, IsDayUnits, MyYear, MyMonth, MyDay,   &
     &                MyHour, MyMinutes, MySeconds)
        MyYday=yearday(MyYear, MyMonth, MyDay)
!
!  The model clock is the elapsed time since reference time of the form
!  'time-units since 0001-01-01 00:00:00'.  It is used in analytical
!  test cases. It has a year length of 365.2425 days (that is, Gregorian
!  Calendar adapted in 15 October 1582). It is called the Gregorian
!  Proleptic Calendar.
!
      ELSE IF (INT(time_ref).eq.0) THEN
        DateNumber=Rclock.Dnumber+CurrentTime         ! fractional days
        DayFraction=DateNumber-AINT(DateNumber)
        HourFraction=DayFraction*24.0_r8
        MinuteFraction=(HourFraction-AINT(HourFraction))*60.0_r8
!
        IsDayUnits=.TRUE.
        CALL datevec (DateNumber, IsDayUnits, MyYear, MyMonth, MyDay,   &
     &                MyHour, MyMinutes, MySeconds)
        MyYday=yearday(MyYear, MyMonth, MyDay)
!
!  The model clock is the elapsed time since reference time of the form
!  'time-units since 0001-01-01 00:00:00'.  It can be used for
!  climatological solutions. It has a year length of 360 days and
!  every month has 30 days.  It is called the 360_day calendar by
!  numerical modelers.
!
      ELSE IF (INT(time_ref).eq.-1) THEN
        DateNumber=Rclock.Dnumber+CurrentTime
        DayFraction=DateNumber-AINT(DateNumber)
        HourFraction=DayFraction*24.0_r8
        MinuteFraction=(HourFraction-AINT(HourFraction))*60.0_r8
!
        MyYear=INT(DateNumber/360.0_r8)
        MyYday=INT(MOD(DateNumber,360.0_r8))
        IF (MyYday.eq.0) THEN
          MyYday=360            ! last day of current year, MOD(x,360)=0
          MyYear=MyYear-1       ! susbract one, new year tomorrow
        END IF
        IF (MOD(MyYday,30).eq.0) THEN
          MyMonth=MyYday/30
          MyDay=30              ! last day of the month, MOD(x,30)=0
        ELSE
          MyMonth=(MyYday/30)+1
          MyDay=MOD(MyYday,30)
        END IF
!
        MySeconds=DayFraction*86400.0_r8
        CT=3.0_r8*EPSILON(MySeconds)           ! comparison tolerance
        MySeconds=ROUND(MySeconds, CT)         ! tolerant round function
        MyHour=INT(MySeconds/3600.0_r8)
        MySeconds=MySeconds-REAL(MyHour*3600,r8)
        MyMinutes=INT(MySeconds/60.0_r8)
        MySeconds=MySeconds-REAL(MyMinutes*60,r8)
!
!  The model clock is the elapsed time since reference time of the form
!  'time-units since 1968-05-23 00:00:00 GMT'. It is a Truncated Julian
!  day introduced by NASA and primarilily used by Astronomers. It has
!  a year length of 365.25 days. It is less used nowadays since the length
!  of the year is 648 seconds less (365.2425) resulting in too many leap
!  years.  So it is correct after 15 October 1582 and it is now called
!  the Gregorian Calendar.
!
      ELSE IF (INT(time_ref).eq.-2) THEN
        DayFraction=CurrentTime-AINT(CurrentTime)
        HourFraction=DayFraction*24.0_r8
        MinuteFraction=(HourFraction-AINT(HourFraction))*60.0_r8
!
        IF (CurrentTime.ge.Rclock.Dnumber) THEN
          jday=INT(CurrentTime)                 ! Origin: Jan 1, 4713 BC
        ELSE
          jday=INT(Rclock.Dnumber+CurrentTime)  ! Truncated Julian Day
        END IF                                  ! add 2440000 offset
        IF (jday.ge.gregorian) THEN
          jalpha=INT(((jday-1867216)-0.25_r8)/36524.25_r8)  ! Gregorian
          ja=jday+1+jalpha-INT(0.25_r8*REAL(jalpha,r8))     ! correction
        ELSE
          ja=jday
        END IF
        jb=ja+1524
        jc=INT(6680.0_r8+(REAL(jb-2439870,r8)-122.1_r8)/365.25_r8)
        jd=365*jc+INT(0.25_r8*REAL(jc,r8))
        je=INT(REAL(jb-jd,r8)/30.6001_r8)
        MyDay=jb-jd-INT(30.6001_r8*REAL(je,r8))
        MyMonth=je-1
        IF (MyMonth.gt.12) MyMonth=MyMonth-12
        MyYear=jc-4715
        IF (MyMonth.gt.2) MyYear=MyYear-1
        IF (MyYear .le.0) MyYear=MyYear-1
        MyYday=yearday(MyYear, MyMonth, MyDay)
!
        MySeconds=DayFraction*86400.0_r8
        CT=3.0_r8*EPSILON(MySeconds)           ! comparison tolerance
        MySeconds=ROUND(MySeconds, CT)         ! tolerant round function
        MyHour=INT(MySeconds/3600.0_r8)
        MySeconds=MySeconds-REAL(MyHour*3600,r8)
        MyMinutes=INT(MySeconds/60.0_r8)
        MySeconds=MySeconds-REAL(MyMinutes*60,r8)
      END IF

---

The reporting of ROMS clock is improved in the standard output file. It provides information about the grid, REGRID interpolation, time in days, and date/clock:

 NLM: GET_STATE - Read state initial conditions,                   2014-01-01 00:00:00.00
                   (Grid 01, t = 2922.0000, File: ini_doppio_20140101_atmpress_NAVD88.nc, Rec=0001, Index=1)

...

 NLM: GET_STATE - Read state initial conditions,                   2014-01-01 00:00:00.00
                   (Grid 02, t = 2922.0000, File: ini_pioneer_20140101_atmpress_NAVD88_v3bathy.nc, Rec=0001, Index=1)

...

    GET_2DFLD   - surface u-wind component,                        2014-01-01 00:00:00.00
                   (Grid=01, Rec=0000001, Index=1, File: Uwind_nam_3hourly_MAB_and_GoM_2014.nc)
                   (Tmin=       2922.0000 Tmax=       3287.0000)      t =       2922.0000
                   (Min = -1.66847594E+00 Max =  1.65528660E+01)      regrid = T
    GET_2DFLD   - surface v-wind component,                        2014-01-01 00:00:00.00
                   (Grid=01, Rec=0000001, Index=1, File: Vwind_nam_3hourly_MAB_and_GoM_2014.nc)
                   (Tmin=       2922.0000 Tmax=       3287.0000)      t =       2922.0000
                   (Min = -7.48610586E+00 Max =  5.12102658E+00)      regrid = T

...

    GET_2DFLD   - surface u-wind component,                        2014-01-01 00:00:00.00
                   (Grid=02, Rec=0000001, Index=1, File: Uwind_nam_3hourly_MAB_and_GoM_2014.nc)
                   (Tmin=       2922.0000 Tmax=       3287.0000)      t =       2922.0000
                   (Min =  3.47775363E+00 Max =  1.34733235E+01)      regrid = T
    GET_2DFLD   - surface v-wind component,                        2014-01-01 00:00:00.00
                   (Grid=02, Rec=0000001, Index=1, File: Vwind_nam_3hourly_MAB_and_GoM_2014.nc)
                   (Tmin=       2922.0000 Tmax=       3287.0000)      t =       2922.0000
                   (Min = -4.34088251E+00 Max =  5.16844070E+00)      regrid = T

...

 TIME-STEP YYYY-MM-DD hh:mm:ss.ss  KINETIC_ENRG   POTEN_ENRG    TOTAL_ENRG    NET_VOLUME  Grid
                     C => (i,j,k)       Cu            Cv            Cw         Max Speed

         0 2014-01-01 00:00:00.00  2.479700E-02  1.895958E+04  1.895961E+04  2.057379E+15  01
                     (128,010,40)  6.586887E-02  8.830755E-02  0.000000E+00  2.146392E+00
      DEF_HIS     - creating  history      file, Grid 01: doppio_his_0001.nc
      WRT_HIS     - wrote history     fields (Index=1,1) in record = 0000001               01
      DEF_STATION - creating  stations     file, Grid 01: doppio_sta.nc

...

         0 2014-01-01 00:00:00.00  9.826892E-03  1.355795E+04  1.355796E+04  2.279106E+14  02
                     (067,001,40)  6.153040E-02  3.379189E-02  0.000000E+00  1.549213E+00
      DEF_HIS     - creating  history      file, Grid 02: pioneer_his_0001.nc
      WRT_HIS     - wrote history     fields (Index=1,1) in record = 0000001               02
      DEF_STATION - creating  stations     file, Grid 02: pioneer_sta.nc

The reporting of the time-step has a i10 field format and the reported counter is reset ever 10 billion time-steps in diag.F, as follows:

          IF (Master) THEN    ! restart counter after 10 billion steps
            WRITE (stdout,30) MOD(iic(ng)-1,10000000000), time_code(ng),&
#ifdef NESTING
     &                        avgke, avgpe, avgkp, volume, ng
 30         FORMAT (i10,1x,a,4(1pe14.6),2x,i2.2)
#else
     &                        avgke, avgpe, avgkp, volume
 30         FORMAT (i10,1x,a,4(1pe14.6))
#endif
            ...
          END IF