libf/IOIPSL/calendar.f90

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	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	USE strlowercase_m, ONLY : strlowercase
26	USE errioipsl, ONLY : histerr
27	!-
28	PRIVATE
29	PUBLIC :: ymds2ju,ju2ymds,isittime,ioconf_calendar, &
30	ioget_calendar,itau2date, ioconf_startdate
31	!-
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	monthp1 = month-freq
337	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