libf/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	MODULE calendar
2
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	!- - 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	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	julian = julian_day + julian_sec / un_jour
73	!---------------------
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	monthp1 = month - freq
336	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