Opened 8 years ago
Closed 8 years ago
#724 closed upgrade (Done)
Important: ROMS clock/date/calendar revisited
Reported by: | arango | Owned by: | arango |
---|---|---|---|
Priority: | major | Milestone: | Release ROMS/TOMS 3.7 |
Component: | Nonlinear | Version: | 3.7 |
Keywords: | Cc: |
Description
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