IOIPSL/Calendar/calendar.f

MODULE calendar
  !$Header: /home/ioipsl/CVSROOT/IOIPSL/src/calendar.f90,v 2.0 2004/04/05 14:47:47 adm Exp $
  !-
  !---------------------------------------------------------------------
  !- This is the calendar which going to be used to do all
  !- calculations on time. Three types of calendars are possible :
  !-  - gregorian : The normal calendar. The time origin for the
  !-                julian day in this case is 24 Nov -4713
  !-  - nolap : A 365 day year without leap years.
  !-            The origin for the julian days is in this case 1 Jan 0
  !-  - xxxd : Year of xxx days with month of equal length.
  !-           The origin for the julian days is then also 1 Jan 0
  !- As one can see it is difficult to go from one calendar to the other.
  !- All operations involving julian days will be wrong.
  !- This calendar will lock as soon as possible
  !- the length of the year and  forbid any further modification.
  !-
  !- For the non leap-year calendar the method is still brute force.
  !- We need to find an Integer series which takes care of the length
  !- of the various month. (Jan)
  !-
  !-   un_jour : one day in seconds
  !-   un_an   : one year in days
  !---------------------------------------------------------------------
  USE strlowercase_m, ONLY : strlowercase
  USE errioipsl, ONLY : histerr
  !-
  PRIVATE
  PUBLIC :: ymds2ju,ju2ymds,isittime,ioconf_calendar, &
       ioget_calendar,itau2date, ioconf_startdate
  !-
  INTERFACE ioget_calendar
     MODULE PROCEDURE &
          &    ioget_calendar_real1,ioget_calendar_real2,ioget_calendar_str
  END INTERFACE
  !-
  REAL,PARAMETER :: un_jour = 86400.0
  LOGICAL,SAVE :: lock_startdate = .FALSE.
  !-
  CHARACTER(LEN=30),SAVE :: time_stamp='XXXXXXXXXXXXXXXX'
  !-
  !- Description of calendar
  !-
  CHARACTER(LEN=20),SAVE :: calendar_used="gregorian"
  LOGICAL,SAVE :: lock_unan = .FALSE.
  REAL,SAVE :: un_an = 365.2425
  INTEGER,SAVE :: mon_len(12)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
  !-
  !-
  !-
  CHARACTER(LEN=3),PARAMETER :: &
       &  cal(12) = (/'JAN','FEB','MAR','APR','MAY','JUN', &
       &              'JUL','AUG','SEP','OCT','NOV','DEC'/)
  !-
  REAL,SAVE :: start_day,start_sec

CONTAINS

  SUBROUTINE ymds2ju (year,month,day,sec,julian)

    IMPLICIT NONE

    INTEGER,INTENT(IN) :: year,month,day
    REAL,INTENT(IN)    :: sec

    REAL,INTENT(OUT) :: julian

    INTEGER :: julian_day
    REAL    :: julian_sec
    !---------------------------------------------------------------------
    CALL ymds2ju_internal (year,month,day,sec,julian_day,julian_sec)

    julian = julian_day+julian_sec / un_jour
    !---------------------
  END SUBROUTINE ymds2ju

  !===

  SUBROUTINE ymds2ju_internal (year,month,day,sec,julian_day,julian_sec)
    !---------------------------------------------------------------------
    !- Converts year, month, day and seconds into a julian day

    !- In 1968 in a letter to the editor of Communications of the ACM
    !- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
    !- and Thomas C. Van Flandern presented such an algorithm.

    !- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm

    !- In the case of the Gregorian calendar we have chosen to use
    !- the Lilian day numbers. This is the day counter which starts
    !- on the 15th October 1582.
    !- This is the day at which Pope Gregory XIII introduced the
    !- Gregorian calendar.
    !- Compared to the true Julian calendar, which starts some
    !- 7980 years ago, the Lilian days are smaler and are dealt with
    !- easily on 32 bit machines. With the true Julian days you can only
    !- the fraction of the day in the real part to a precision of
    !- a 1/4 of a day with 32 bits.
    !---------------------------------------------------------------------
    IMPLICIT NONE

    INTEGER,INTENT(IN) :: year,month,day
    REAL,INTENT(IN)    :: sec

    INTEGER,INTENT(OUT) :: julian_day
    REAL,INTENT(OUT)    :: julian_sec

    INTEGER :: jd,m,y,d,ml
    !---------------------------------------------------------------------
    lock_unan = .TRUE.

    m = month
    y = year
    d = day

    !- We deduce the calendar from the length of the year as it
    !- is faster than an INDEX on the calendar variable.

    !- Gregorian
    IF ( (un_an > 365.0).AND.(un_an < 366.0) ) THEN
       jd = (1461*(y+4800+INT(( m-14 )/12)))/4 &
            &      +(367*(m-2-12*(INT(( m-14 )/12))))/12 &
            &      -(3*((y+4900+INT((m-14)/12))/100))/4 &
            &      +d-32075
       jd = jd-2299160
       !- No leap or All leap
    ELSE IF (ABS(un_an-365.0) <= EPSILON(un_an) .OR. &
         &   ABS(un_an-366.0) <= EPSILON(un_an)) THEN
       ml = SUM(mon_len(1:m-1))
       jd = y*INT(un_an)+ml+(d-1)
       !- Calendar with regular month
    ELSE
       ml = INT(un_an)/12
       jd = y*INT(un_an)+(m-1)*ml+(d-1)
    ENDIF

    julian_day = jd
    julian_sec = sec
    !------------------------------
  END SUBROUTINE ymds2ju_internal
  !-
  !===
  !-
  SUBROUTINE ju2ymds (julian,year,month,day,sec)
    !---------------------------------------------------------------------
    IMPLICIT NONE

    REAL,INTENT(IN) :: julian

    INTEGER,INTENT(OUT) :: year,month,day
    REAL,INTENT(OUT)    :: sec

    INTEGER :: julian_day
    REAL    :: julian_sec
    !---------------------------------------------------------------------
    julian_day = INT(julian)
    julian_sec = (julian-julian_day)*un_jour

    CALL ju2ymds_internal(julian_day,julian_sec,year,month,day,sec)
    !---------------------
  END SUBROUTINE ju2ymds
  !-
  !===
  !-
  SUBROUTINE ju2ymds_internal (julian_day,julian_sec,year,month,day,sec)
    !---------------------------------------------------------------------
    !- This subroutine computes from the julian day the year,
    !- month, day and seconds

    !- In 1968 in a letter to the editor of Communications of the ACM
    !- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
    !- and Thomas C. Van Flandern presented such an algorithm.

    !- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm

    !- In the case of the Gregorian calendar we have chosen to use
    !- the Lilian day numbers. This is the day counter which starts
    !- on the 15th October 1582. This is the day at which Pope
    !- Gregory XIII introduced the Gregorian calendar.
    !- Compared to the true Julian calendar, which starts some 7980
    !- years ago, the Lilian days are smaler and are dealt with easily
    !- on 32 bit machines. With the true Julian days you can only the
    !- fraction of the day in the real part to a precision of a 1/4 of
    !- a day with 32 bits.
    !---------------------------------------------------------------------
    IMPLICIT NONE

    INTEGER,INTENT(IN) :: julian_day
    REAL,INTENT(IN)    :: julian_sec

    INTEGER,INTENT(OUT) :: year,month,day
    REAL,INTENT(OUT)    :: sec

    INTEGER :: l,n,i,jd,j,d,m,y,ml
    INTEGER :: add_day
    !---------------------------------------------------------------------
    lock_unan = .TRUE.

    jd = julian_day
    sec = julian_sec
    IF (sec > un_jour) THEN
       add_day = INT(sec/un_jour)
       sec = sec-add_day*un_jour
       jd = jd+add_day
    ENDIF

    !- Gregorian
    IF ( (un_an > 365.0).AND.(un_an < 366.0) ) THEN
       jd = jd+2299160
   
       l = jd+68569
       n = (4*l)/146097
       l = l-(146097*n+3)/4
       i = (4000*(l+1))/1461001
       l = l-(1461*i)/4+31
       j = (80*l)/2447
       d = l-(2447*j)/80
       l = j/11
       m = j+2-(12*l)
       y = 100*(n-49)+i+l
       !- No leap or All leap
    ELSE IF (ABS(un_an-365.0) <= EPSILON(un_an) .OR. &
         &   ABS(un_an-366.0) <= EPSILON(un_an) ) THEN
       y = jd/INT(un_an)
       l = jd-y*INT(un_an)
       m = 1
       ml = 0
       DO WHILE (ml+mon_len(m) <= l)
          ml = ml+mon_len(m)
          m = m+1
       ENDDO
       d = l-ml+1
       !- others
    ELSE
       ml = INT(un_an)/12
       y = jd/INT(un_an)
       l = jd-y*INT(un_an)
       m = (l/ml)+1
       d = l-(m-1)*ml+1
    ENDIF

    day = d
    month = m
    year = y
    !------------------------------
  END SUBROUTINE ju2ymds_internal
  !-
  !===
  !-
  REAL FUNCTION itau2date (itau,date0,deltat)
    !---------------------------------------------------------------------
    !- This function transforms itau into a date. The date whith which
    !- the time axis is going to be labeled

    !- INPUT
    !-   itau   : current time step
    !-   date0  : Date at which itau was equal to 0
    !-   deltat : time step between itau s

    !- OUTPUT
    !-   itau2date : Date for the given itau
    !---------------------------------------------------------------------
    IMPLICIT NONE

    INTEGER  :: itau
    REAL     :: date0,deltat
    !---------------------------------------------------------------------
    itau2date = REAL(itau)*deltat/un_jour+date0
    !---------------------
  END FUNCTION itau2date
  !-
  !===
  !-
  SUBROUTINE isittime &
       &  (itau,date0,dt,freq,last_action,last_check,do_action)
    !---------------------------------------------------------------------
    !- This subroutine checks the time has come for a given action.
    !- This is computed from the current time-step(itau).
    !- Thus we need to have the time delta (dt), the frequency
    !- of the action (freq) and the last time it was done
    !- (last_action in units of itau).
    !- In order to extrapolate when will be the next check we need
    !- the time step of the last call (last_check).

    !- The test is done on the following condition :
    !- the distance from the current time to the time for the next
    !- action is smaller than the one from the next expected
    !- check to the next action.
    !- When the test is done on the time steps simplifactions make
    !- it more difficult to read in the code.
    !- For the real time case it is easier to understand !
    !---------------------------------------------------------------------
    IMPLICIT NONE

    INTEGER,INTENT(IN) :: itau
    REAL,INTENT(IN)    :: dt,freq
    INTEGER,INTENT(IN) :: last_action,last_check
    REAL,INTENT(IN)    :: date0

    LOGICAL,INTENT(OUT)  :: do_action

    REAL :: dt_action,dt_check
    REAL :: date_last_act,date_next_check,date_next_act, &
         &        date_now,date_mp1,date_mpf
    INTEGER :: year,month,monthp1,day,next_check_itau,next_act_itau
    INTEGER :: yearp,dayp
    REAL :: sec,secp
    LOGICAL :: check = .FALSE.
    !---------------------------------------------------------------------
    IF (check) THEN
       WRITE(*,*) &
            &    "isittime 1.0 ",itau,date0,dt,freq,last_action,last_check
    ENDIF

    IF (last_check >= 0) THEN
       dt_action = (itau-last_action)*dt
       dt_check = (itau-last_check)*dt
       next_check_itau = itau+(itau-last_check)
   
       !-- We are dealing with frequencies in seconds and thus operation
       !-- can be done on the time steps.
   
       IF (freq > 0) THEN
          IF (ABS(dt_action-freq) <= ABS(dt_action+dt_check-freq)) THEN
             do_action = .TRUE.
          ELSE
             do_action = .FALSE.
          ENDIF
      
          !---- Here we deal with frequencies in month and work on julian days.
      
       ELSE
          date_now = itau2date (itau,date0,dt)
          date_last_act = itau2date (last_action,date0,dt)
          CALL ju2ymds (date_last_act,year,month,day,sec)
          monthp1 = month - freq
          yearp = year
      
          !---- Here we compute what logically should be the next month
      
          IF (month >= 13) THEN
             yearp = year+1
             monthp1 = monthp1-12
          ENDIF
          CALL ymds2ju (year,monthp1,day,sec,date_mpf)
      
          !---- But it could be that because of a shorter month or a bad
          !---- starting date that we end up further than we should be.
          !---- Thus we compute the first day of the next month.
          !---- We can not be beyond this date and if we are close
          !---- then we will take it as it is better.
      
          monthp1 = month+ABS(freq)
          yearp=year
          IF (monthp1 >= 13) THEN
             yearp = year+1
             monthp1 = monthp1 -12
          ENDIF
          dayp = 1
          secp = 0.0
          CALL ymds2ju (yearp,monthp1,dayp,secp,date_mp1)
      
          !---- If date_mp1 is smaller than date_mpf or only less than 4 days
          !---- larger then we take it. This needed to ensure that short month
          !---- like February do not mess up the thing !
      
          IF (date_mp1-date_mpf < 4.) THEN
             date_next_act = date_mp1
          ELSE
             date_next_act = date_mpf
          ENDIF
          date_next_check = itau2date (next_check_itau,date0,dt)
      
          !---- Transform the dates into time-steps for the needed precisions.
      
          next_act_itau = &
               &      last_action+INT((date_next_act-date_last_act)*(un_jour/dt))
          !-----
          IF (   ABS(itau-next_act_itau) &
               &        <= ABS( next_check_itau-next_act_itau)) THEN
             do_action = .TRUE.
             IF (check) THEN
                WRITE(*,*) &
                     &         'ACT-TIME : itau, next_act_itau, next_check_itau : ', &
                     &         itau,next_act_itau,next_check_itau
                CALL ju2ymds (date_now,year,month,day,sec)
                WRITE(*,*) 'ACT-TIME : y, m, d, s : ',year,month,day,sec
                WRITE(*,*) &
                     &         'ACT-TIME : date_mp1, date_mpf : ',date_mp1,date_mpf
             ENDIF
          ELSE
             do_action = .FALSE.
          ENDIF
       ENDIF
   
       IF (check) THEN
          WRITE(*,*) "isittime 2.0 ", &
               &     date_next_check,date_next_act,ABS(dt_action-freq), &
               &     ABS(dt_action+dt_check-freq),dt_action,dt_check, &
               &     next_check_itau,do_action
       ENDIF
    ELSE
       do_action=.FALSE.
    ENDIF
    !----------------------
  END SUBROUTINE isittime
  !-
  !===
  !-
  SUBROUTINE ioconf_calendar (str)
    !---------------------------------------------------------------------
    !- This routine allows to configure the calendar to be used.
    !- This operation is only allowed once and the first call to
    !- ymds2ju or ju2ymsd will lock the current configuration.
    !- the argument to ioconf_calendar can be any of the following :
    !-  - gregorian : This is the gregorian calendar (default here)
    !-  - noleap    : A calendar without leap years = 365 days
    !-  - xxxd      : A calendar of xxx days (has to be a modulo of 12)
    !-                with 12 month of equal length
    !---------------------------------------------------------------------
    IMPLICIT NONE

    CHARACTER(LEN=*),INTENT(IN) :: str

    INTEGER :: leng,ipos
    CHARACTER(LEN=10) :: str10
    !---------------------------------------------------------------------

    ! 1.0 Clean up the sring !

    CALL strlowercase (str)

    IF (.NOT.lock_unan) THEN
       !---
       lock_unan=.TRUE.
       !---
       SELECT CASE(str)
       CASE('gregorian')
          calendar_used = 'gregorian'
          un_an = 365.2425
          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
       CASE('standard')
          calendar_used = 'gregorian'
          un_an = 365.2425
          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
       CASE('proleptic_gregorian')
          calendar_used = 'gregorian'
          un_an = 365.2425
          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
       CASE('noleap')
          calendar_used = 'noleap'
          un_an = 365.0
          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
       CASE('365_day')
          calendar_used = 'noleap'
          un_an = 365.0
          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
       CASE('365d')
          calendar_used = 'noleap'
          un_an = 365.0
          mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
       CASE('all_leap')
          calendar_used = 'all_leap'
          un_an = 366.0
          mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
       CASE('366_day')
          calendar_used = 'all_leap'
          un_an = 366.0
          mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
       CASE('366d')
          calendar_used = 'all_leap'
          un_an = 366.0
          mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
       CASE DEFAULT
          ipos = INDEX(str,'d')
          IF (ipos == 4) THEN
             READ(str(1:3),'(I3)') leng
             IF ( (MOD(leng,12) == 0).AND.(leng > 1) ) THEN
                calendar_used = str
                un_an = leng
                mon_len(:) = leng
             ELSE
                CALL histerr (3,'ioconf_calendar', &
                     &         'The length of the year as to be a modulo of 12', &
                     &         'so that it can be divided into 12 month of equal length', &
                     &         str)
             ENDIF
          ELSE
             CALL histerr (3,'ioconf_calendar', &
                  &       'Unrecognized input, please ceck the man pages.',str,' ')
          ENDIF
       END SELECT
    ELSE
       WRITE(str10,'(f10.4)') un_an
       CALL histerr (2,'ioconf_calendar', &
            &   'The calendar was already used or configured. You are not', &
            &   'allowed to change it again. '// &
            &   'The following length of year is used :',str10)
    ENDIF
    !-----------------------------
  END SUBROUTINE ioconf_calendar
  !-
  !===
  !-
  SUBROUTINE ioget_calendar_str (str)
    !---------------------------------------------------------------------
    !- This subroutine returns the name of the calendar used here.
    !- Three options exist :
    !-  - gregorian : This is the gregorian calendar (default here)
    !-  - noleap    : A calendar without leap years = 365 days
    !-  - xxxd      : A calendar of xxx days (has to be a modulo of 12)
    !-                with 12 month of equal length

    !- This routine will lock the calendar.
    !- You do not want it to change after your inquiry.
    !---------------------------------------------------------------------
    IMPLICIT NONE

    CHARACTER(LEN=*),INTENT(OUT) :: str
    !---------------------------------------------------------------------
    lock_unan = .TRUE.

    str = calendar_used
    !--------------------------------
  END SUBROUTINE ioget_calendar_str
  !-
  !===
  !-
  SUBROUTINE ioget_calendar_real1 (long_an)
    !---------------------------------------------------------------------
    !- This subroutine returns the name of the calendar used here.
    !- Three options exist :
    !-  - gregorian : This is the gregorian calendar (default here)
    !-  - noleap    : A calendar without leap years = 365 days
    !-  - xxxd      : A calendar of xxx days (has to be a modulo of 12)
    !-                with 12 month of equal length

    !- This routine will lock the calendar.
    !- You do not want it to change after your inquiry.
    !---------------------------------------------------------------------
    IMPLICIT NONE

    REAL,INTENT(OUT) :: long_an
    !---------------------------------------------------------------------
    lock_unan = .TRUE.

    long_an = un_an
    !----------------------------------
  END SUBROUTINE ioget_calendar_real1
  !-
  !===
  !-
  SUBROUTINE ioget_calendar_real2 (long_an,long_jour)
    !---------------------------------------------------------------------
    !- This subroutine returns the name of the calendar used here.
    !- Three options exist :
    !-  - gregorian : This is the gregorian calendar (default here)
    !-  - noleap    : A calendar without leap years = 365 days
    !-  - xxxd      : A calendar of xxx days (has to be a modulo of 12)
    !-                with 12 month of equal length

    !- This routine will lock the calendar.
    !- You do not want it to change after your inquiry.
    !---------------------------------------------------------------------
    IMPLICIT NONE

    REAL,INTENT(OUT) :: long_an,long_jour
    !---------------------------------------------------------------------
    lock_unan = .TRUE.

    long_an = un_an
    long_jour = un_jour
    !----------------------------------
  END SUBROUTINE ioget_calendar_real2

END MODULE calendar
1	guez	30	MODULE calendar
2			!$Header: /home/ioipsl/CVSROOT/IOIPSL/src/calendar.f90,v 2.0 2004/04/05 14:47:47 adm Exp $
3			!-
4			!---------------------------------------------------------------------
5			!- This is the calendar which going to be used to do all
6			!- calculations on time. Three types of calendars are possible :
7			!- - gregorian : The normal calendar. The time origin for the
8			!- julian day in this case is 24 Nov -4713
9			!- - nolap : A 365 day year without leap years.
10			!- The origin for the julian days is in this case 1 Jan 0
11			!- - xxxd : Year of xxx days with month of equal length.
12			!- The origin for the julian days is then also 1 Jan 0
13			!- As one can see it is difficult to go from one calendar to the other.
14			!- All operations involving julian days will be wrong.
15			!- This calendar will lock as soon as possible
16			!- the length of the year and forbid any further modification.
17			!-
18			!- For the non leap-year calendar the method is still brute force.
19			!- We need to find an Integer series which takes care of the length
20			!- of the various month. (Jan)
21			!-
22			!- un_jour : one day in seconds
23			!- un_an : one year in days
24			!---------------------------------------------------------------------
25	guez	32	USE strlowercase_m, ONLY : strlowercase
26	guez	30	USE errioipsl, ONLY : histerr
27			!-
28			PRIVATE
29	guez	32	PUBLIC :: ymds2ju,ju2ymds,isittime,ioconf_calendar, &
30			ioget_calendar,itau2date, ioconf_startdate
31	guez	30	!-
32			INTERFACE ioget_calendar
33			MODULE PROCEDURE &
34			& ioget_calendar_real1,ioget_calendar_real2,ioget_calendar_str
35			END INTERFACE
36			!-
37			REAL,PARAMETER :: un_jour = 86400.0
38			LOGICAL,SAVE :: lock_startdate = .FALSE.
39			!-
40			CHARACTER(LEN=30),SAVE :: time_stamp='XXXXXXXXXXXXXXXX'
41			!-
42			!- Description of calendar
43			!-
44			CHARACTER(LEN=20),SAVE :: calendar_used="gregorian"
45			LOGICAL,SAVE :: lock_unan = .FALSE.
46			REAL,SAVE :: un_an = 365.2425
47			INTEGER,SAVE :: mon_len(12)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
48			!-
49			!-
50			!-
51			CHARACTER(LEN=3),PARAMETER :: &
52			& cal(12) = (/'JAN','FEB','MAR','APR','MAY','JUN', &
53			& 'JUL','AUG','SEP','OCT','NOV','DEC'/)
54			!-
55			REAL,SAVE :: start_day,start_sec
56
57			CONTAINS
58
59			SUBROUTINE ymds2ju (year,month,day,sec,julian)
60
61			IMPLICIT NONE
62
63			INTEGER,INTENT(IN) :: year,month,day
64			REAL,INTENT(IN) :: sec
65
66			REAL,INTENT(OUT) :: julian
67
68			INTEGER :: julian_day
69			REAL :: julian_sec
70			!---------------------------------------------------------------------
71			CALL ymds2ju_internal (year,month,day,sec,julian_day,julian_sec)
72
73			julian = julian_day+julian_sec / un_jour
74			!---------------------
75			END SUBROUTINE ymds2ju
76
77			!===
78
79			SUBROUTINE ymds2ju_internal (year,month,day,sec,julian_day,julian_sec)
80			!---------------------------------------------------------------------
81			!- Converts year, month, day and seconds into a julian day
82
83			!- In 1968 in a letter to the editor of Communications of the ACM
84			!- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
85			!- and Thomas C. Van Flandern presented such an algorithm.
86
87			!- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm
88
89			!- In the case of the Gregorian calendar we have chosen to use
90			!- the Lilian day numbers. This is the day counter which starts
91			!- on the 15th October 1582.
92			!- This is the day at which Pope Gregory XIII introduced the
93			!- Gregorian calendar.
94			!- Compared to the true Julian calendar, which starts some
95			!- 7980 years ago, the Lilian days are smaler and are dealt with
96			!- easily on 32 bit machines. With the true Julian days you can only
97			!- the fraction of the day in the real part to a precision of
98			!- a 1/4 of a day with 32 bits.
99			!---------------------------------------------------------------------
100			IMPLICIT NONE
101
102			INTEGER,INTENT(IN) :: year,month,day
103			REAL,INTENT(IN) :: sec
104
105			INTEGER,INTENT(OUT) :: julian_day
106			REAL,INTENT(OUT) :: julian_sec
107
108			INTEGER :: jd,m,y,d,ml
109			!---------------------------------------------------------------------
110			lock_unan = .TRUE.
111
112			m = month
113			y = year
114			d = day
115
116			!- We deduce the calendar from the length of the year as it
117			!- is faster than an INDEX on the calendar variable.
118
119			!- Gregorian
120			IF ( (un_an > 365.0).AND.(un_an < 366.0) ) THEN
121			jd = (1461*(y+4800+INT(( m-14 )/12)))/4 &
122			& +(367(m-2-12(INT(( m-14 )/12))))/12 &
123			& -(3*((y+4900+INT((m-14)/12))/100))/4 &
124			& +d-32075
125			jd = jd-2299160
126			!- No leap or All leap
127			ELSE IF (ABS(un_an-365.0) <= EPSILON(un_an) .OR. &
128			& ABS(un_an-366.0) <= EPSILON(un_an)) THEN
129			ml = SUM(mon_len(1:m-1))
130			jd = y*INT(un_an)+ml+(d-1)
131			!- Calendar with regular month
132			ELSE
133			ml = INT(un_an)/12
134			jd = yINT(un_an)+(m-1)ml+(d-1)
135			ENDIF
136
137			julian_day = jd
138			julian_sec = sec
139			!------------------------------
140			END SUBROUTINE ymds2ju_internal
141			!-
142			!===
143			!-
144			SUBROUTINE ju2ymds (julian,year,month,day,sec)
145			!---------------------------------------------------------------------
146			IMPLICIT NONE
147
148			REAL,INTENT(IN) :: julian
149
150			INTEGER,INTENT(OUT) :: year,month,day
151			REAL,INTENT(OUT) :: sec
152
153			INTEGER :: julian_day
154			REAL :: julian_sec
155			!---------------------------------------------------------------------
156			julian_day = INT(julian)
157			julian_sec = (julian-julian_day)*un_jour
158
159			CALL ju2ymds_internal(julian_day,julian_sec,year,month,day,sec)
160			!---------------------
161			END SUBROUTINE ju2ymds
162			!-
163			!===
164			!-
165			SUBROUTINE ju2ymds_internal (julian_day,julian_sec,year,month,day,sec)
166			!---------------------------------------------------------------------
167			!- This subroutine computes from the julian day the year,
168			!- month, day and seconds
169
170			!- In 1968 in a letter to the editor of Communications of the ACM
171			!- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
172			!- and Thomas C. Van Flandern presented such an algorithm.
173
174			!- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm
175
176			!- In the case of the Gregorian calendar we have chosen to use
177			!- the Lilian day numbers. This is the day counter which starts
178			!- on the 15th October 1582. This is the day at which Pope
179			!- Gregory XIII introduced the Gregorian calendar.
180			!- Compared to the true Julian calendar, which starts some 7980
181			!- years ago, the Lilian days are smaler and are dealt with easily
182			!- on 32 bit machines. With the true Julian days you can only the
183			!- fraction of the day in the real part to a precision of a 1/4 of
184			!- a day with 32 bits.
185			!---------------------------------------------------------------------
186			IMPLICIT NONE
187
188			INTEGER,INTENT(IN) :: julian_day
189			REAL,INTENT(IN) :: julian_sec
190
191			INTEGER,INTENT(OUT) :: year,month,day
192			REAL,INTENT(OUT) :: sec
193
194			INTEGER :: l,n,i,jd,j,d,m,y,ml
195			INTEGER :: add_day
196			!---------------------------------------------------------------------
197			lock_unan = .TRUE.
198
199			jd = julian_day
200			sec = julian_sec
201			IF (sec > un_jour) THEN
202			add_day = INT(sec/un_jour)
203			sec = sec-add_day*un_jour
204			jd = jd+add_day
205			ENDIF
206
207			!- Gregorian
208			IF ( (un_an > 365.0).AND.(un_an < 366.0) ) THEN
209			jd = jd+2299160
210
211			l = jd+68569
212			n = (4*l)/146097
213			l = l-(146097*n+3)/4
214			i = (4000*(l+1))/1461001
215			l = l-(1461*i)/4+31
216			j = (80*l)/2447
217			d = l-(2447*j)/80
218			l = j/11
219			m = j+2-(12*l)
220			y = 100*(n-49)+i+l
221			!- No leap or All leap
222			ELSE IF (ABS(un_an-365.0) <= EPSILON(un_an) .OR. &
223			& ABS(un_an-366.0) <= EPSILON(un_an) ) THEN
224			y = jd/INT(un_an)
225			l = jd-y*INT(un_an)
226			m = 1
227			ml = 0
228			DO WHILE (ml+mon_len(m) <= l)
229			ml = ml+mon_len(m)
230			m = m+1
231			ENDDO
232			d = l-ml+1
233			!- others
234			ELSE
235			ml = INT(un_an)/12
236			y = jd/INT(un_an)
237			l = jd-y*INT(un_an)
238			m = (l/ml)+1
239			d = l-(m-1)*ml+1
240			ENDIF
241
242			day = d
243			month = m
244			year = y
245			!------------------------------
246			END SUBROUTINE ju2ymds_internal
247			!-
248			!===
249			!-
250			REAL FUNCTION itau2date (itau,date0,deltat)
251			!---------------------------------------------------------------------
252			!- This function transforms itau into a date. The date whith which
253			!- the time axis is going to be labeled
254
255			!- INPUT
256			!- itau : current time step
257			!- date0 : Date at which itau was equal to 0
258			!- deltat : time step between itau s
259
260			!- OUTPUT
261			!- itau2date : Date for the given itau
262			!---------------------------------------------------------------------
263			IMPLICIT NONE
264
265			INTEGER :: itau
266			REAL :: date0,deltat
267			!---------------------------------------------------------------------
268			itau2date = REAL(itau)*deltat/un_jour+date0
269			!---------------------
270			END FUNCTION itau2date
271			!-
272			!===
273			!-
274			SUBROUTINE isittime &
275			& (itau,date0,dt,freq,last_action,last_check,do_action)
276			!---------------------------------------------------------------------
277			!- This subroutine checks the time has come for a given action.
278			!- This is computed from the current time-step(itau).
279			!- Thus we need to have the time delta (dt), the frequency
280			!- of the action (freq) and the last time it was done
281			!- (last_action in units of itau).
282			!- In order to extrapolate when will be the next check we need
283			!- the time step of the last call (last_check).
284
285			!- The test is done on the following condition :
286			!- the distance from the current time to the time for the next
287			!- action is smaller than the one from the next expected
288			!- check to the next action.
289			!- When the test is done on the time steps simplifactions make
290			!- it more difficult to read in the code.
291			!- For the real time case it is easier to understand !
292			!---------------------------------------------------------------------
293			IMPLICIT NONE
294
295			INTEGER,INTENT(IN) :: itau
296			REAL,INTENT(IN) :: dt,freq
297			INTEGER,INTENT(IN) :: last_action,last_check
298			REAL,INTENT(IN) :: date0
299
300			LOGICAL,INTENT(OUT) :: do_action
301
302			REAL :: dt_action,dt_check
303			REAL :: date_last_act,date_next_check,date_next_act, &
304			& date_now,date_mp1,date_mpf
305			INTEGER :: year,month,monthp1,day,next_check_itau,next_act_itau
306			INTEGER :: yearp,dayp
307			REAL :: sec,secp
308			LOGICAL :: check = .FALSE.
309			!---------------------------------------------------------------------
310			IF (check) THEN
311			WRITE(,) &
312			& "isittime 1.0 ",itau,date0,dt,freq,last_action,last_check
313			ENDIF
314
315			IF (last_check >= 0) THEN
316			dt_action = (itau-last_action)*dt
317			dt_check = (itau-last_check)*dt
318			next_check_itau = itau+(itau-last_check)
319
320			!-- We are dealing with frequencies in seconds and thus operation
321			!-- can be done on the time steps.
322
323			IF (freq > 0) THEN
324			IF (ABS(dt_action-freq) <= ABS(dt_action+dt_check-freq)) THEN
325			do_action = .TRUE.
326			ELSE
327			do_action = .FALSE.
328			ENDIF
329
330			!---- Here we deal with frequencies in month and work on julian days.
331
332			ELSE
333			date_now = itau2date (itau,date0,dt)
334			date_last_act = itau2date (last_action,date0,dt)
335			CALL ju2ymds (date_last_act,year,month,day,sec)
336	guez	35	monthp1 = month - freq
337	guez	30	yearp = year
338
339			!---- Here we compute what logically should be the next month
340
341			IF (month >= 13) THEN
342			yearp = year+1
343			monthp1 = monthp1-12
344			ENDIF
345			CALL ymds2ju (year,monthp1,day,sec,date_mpf)
346
347			!---- But it could be that because of a shorter month or a bad
348			!---- starting date that we end up further than we should be.
349			!---- Thus we compute the first day of the next month.
350			!---- We can not be beyond this date and if we are close
351			!---- then we will take it as it is better.
352
353			monthp1 = month+ABS(freq)
354			yearp=year
355			IF (monthp1 >= 13) THEN
356			yearp = year+1
357			monthp1 = monthp1 -12
358			ENDIF
359			dayp = 1
360			secp = 0.0
361			CALL ymds2ju (yearp,monthp1,dayp,secp,date_mp1)
362
363			!---- If date_mp1 is smaller than date_mpf or only less than 4 days
364			!---- larger then we take it. This needed to ensure that short month
365			!---- like February do not mess up the thing !
366
367			IF (date_mp1-date_mpf < 4.) THEN
368			date_next_act = date_mp1
369			ELSE
370			date_next_act = date_mpf
371			ENDIF
372			date_next_check = itau2date (next_check_itau,date0,dt)
373
374			!---- Transform the dates into time-steps for the needed precisions.
375
376			next_act_itau = &
377			& last_action+INT((date_next_act-date_last_act)*(un_jour/dt))
378			!-----
379			IF ( ABS(itau-next_act_itau) &
380			& <= ABS( next_check_itau-next_act_itau)) THEN
381			do_action = .TRUE.
382			IF (check) THEN
383			WRITE(,) &
384			& 'ACT-TIME : itau, next_act_itau, next_check_itau : ', &
385			& itau,next_act_itau,next_check_itau
386			CALL ju2ymds (date_now,year,month,day,sec)
387			WRITE(,) 'ACT-TIME : y, m, d, s : ',year,month,day,sec
388			WRITE(,) &
389			& 'ACT-TIME : date_mp1, date_mpf : ',date_mp1,date_mpf
390			ENDIF
391			ELSE
392			do_action = .FALSE.
393			ENDIF
394			ENDIF
395
396			IF (check) THEN
397			WRITE(,) "isittime 2.0 ", &
398			& date_next_check,date_next_act,ABS(dt_action-freq), &
399			& ABS(dt_action+dt_check-freq),dt_action,dt_check, &
400			& next_check_itau,do_action
401			ENDIF
402			ELSE
403			do_action=.FALSE.
404			ENDIF
405			!----------------------
406			END SUBROUTINE isittime
407			!-
408			!===
409			!-
410			SUBROUTINE ioconf_calendar (str)
411			!---------------------------------------------------------------------
412			!- This routine allows to configure the calendar to be used.
413			!- This operation is only allowed once and the first call to
414			!- ymds2ju or ju2ymsd will lock the current configuration.
415			!- the argument to ioconf_calendar can be any of the following :
416			!- - gregorian : This is the gregorian calendar (default here)
417			!- - noleap : A calendar without leap years = 365 days
418			!- - xxxd : A calendar of xxx days (has to be a modulo of 12)
419			!- with 12 month of equal length
420			!---------------------------------------------------------------------
421			IMPLICIT NONE
422
423			CHARACTER(LEN=*),INTENT(IN) :: str
424
425			INTEGER :: leng,ipos
426			CHARACTER(LEN=10) :: str10
427			!---------------------------------------------------------------------
428
429			! 1.0 Clean up the sring !
430
431			CALL strlowercase (str)
432
433			IF (.NOT.lock_unan) THEN
434			!---
435			lock_unan=.TRUE.
436			!---
437			SELECT CASE(str)
438			CASE('gregorian')
439			calendar_used = 'gregorian'
440			un_an = 365.2425
441			mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
442			CASE('standard')
443			calendar_used = 'gregorian'
444			un_an = 365.2425
445			mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
446			CASE('proleptic_gregorian')
447			calendar_used = 'gregorian'
448			un_an = 365.2425
449			mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
450			CASE('noleap')
451			calendar_used = 'noleap'
452			un_an = 365.0
453			mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
454			CASE('365_day')
455			calendar_used = 'noleap'
456			un_an = 365.0
457			mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
458			CASE('365d')
459			calendar_used = 'noleap'
460			un_an = 365.0
461			mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
462			CASE('all_leap')
463			calendar_used = 'all_leap'
464			un_an = 366.0
465			mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
466			CASE('366_day')
467			calendar_used = 'all_leap'
468			un_an = 366.0
469			mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
470			CASE('366d')
471			calendar_used = 'all_leap'
472			un_an = 366.0
473			mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
474			CASE DEFAULT
475			ipos = INDEX(str,'d')
476			IF (ipos == 4) THEN
477			READ(str(1:3),'(I3)') leng
478			IF ( (MOD(leng,12) == 0).AND.(leng > 1) ) THEN
479			calendar_used = str
480			un_an = leng
481			mon_len(:) = leng
482			ELSE
483			CALL histerr (3,'ioconf_calendar', &
484			& 'The length of the year as to be a modulo of 12', &
485			& 'so that it can be divided into 12 month of equal length', &
486			& str)
487			ENDIF
488			ELSE
489			CALL histerr (3,'ioconf_calendar', &
490			& 'Unrecognized input, please ceck the man pages.',str,' ')
491			ENDIF
492			END SELECT
493			ELSE
494			WRITE(str10,'(f10.4)') un_an
495			CALL histerr (2,'ioconf_calendar', &
496			& 'The calendar was already used or configured. You are not', &
497			& 'allowed to change it again. '// &
498			& 'The following length of year is used :',str10)
499			ENDIF
500			!-----------------------------
501			END SUBROUTINE ioconf_calendar
502			!-
503			!===
504			!-
505			SUBROUTINE ioget_calendar_str (str)
506			!---------------------------------------------------------------------
507			!- This subroutine returns the name of the calendar used here.
508			!- Three options exist :
509			!- - gregorian : This is the gregorian calendar (default here)
510			!- - noleap : A calendar without leap years = 365 days
511			!- - xxxd : A calendar of xxx days (has to be a modulo of 12)
512			!- with 12 month of equal length
513
514			!- This routine will lock the calendar.
515			!- You do not want it to change after your inquiry.
516			!---------------------------------------------------------------------
517			IMPLICIT NONE
518
519			CHARACTER(LEN=*),INTENT(OUT) :: str
520			!---------------------------------------------------------------------
521			lock_unan = .TRUE.
522
523			str = calendar_used
524			!--------------------------------
525			END SUBROUTINE ioget_calendar_str
526			!-
527			!===
528			!-
529			SUBROUTINE ioget_calendar_real1 (long_an)
530			!---------------------------------------------------------------------
531			!- This subroutine returns the name of the calendar used here.
532			!- Three options exist :
533			!- - gregorian : This is the gregorian calendar (default here)
534			!- - noleap : A calendar without leap years = 365 days
535			!- - xxxd : A calendar of xxx days (has to be a modulo of 12)
536			!- with 12 month of equal length
537
538			!- This routine will lock the calendar.
539			!- You do not want it to change after your inquiry.
540			!---------------------------------------------------------------------
541			IMPLICIT NONE
542
543			REAL,INTENT(OUT) :: long_an
544			!---------------------------------------------------------------------
545			lock_unan = .TRUE.
546
547			long_an = un_an
548			!----------------------------------
549			END SUBROUTINE ioget_calendar_real1
550			!-
551			!===
552			!-
553			SUBROUTINE ioget_calendar_real2 (long_an,long_jour)
554			!---------------------------------------------------------------------
555			!- This subroutine returns the name of the calendar used here.
556			!- Three options exist :
557			!- - gregorian : This is the gregorian calendar (default here)
558			!- - noleap : A calendar without leap years = 365 days
559			!- - xxxd : A calendar of xxx days (has to be a modulo of 12)
560			!- with 12 month of equal length
561
562			!- This routine will lock the calendar.
563			!- You do not want it to change after your inquiry.
564			!---------------------------------------------------------------------
565			IMPLICIT NONE
566
567			REAL,INTENT(OUT) :: long_an,long_jour
568			!---------------------------------------------------------------------
569			lock_unan = .TRUE.
570
571			long_an = un_an
572			long_jour = un_jour
573			!----------------------------------
574			END SUBROUTINE ioget_calendar_real2
575
576			END MODULE calendar