trunk/IOIPSL/calendar.f90

MODULE calendar

  ! From IOIPSL/src/calendar.f90, version 2.0 2004/04/05 14:47:47

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