prism_calendar.F90

Go to the documentation of this file.

!
! --------------------  PRISM module file calendar.f90 ------------------
!
! $Id

Module PRISM_calendar
  !-----------------------------------------------------------------------
  ! Proleptic Gregorian Calendar routines:
  !
  ! Copyright 2007-2010, Max-Planck-Institute for Meteorology, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
  !
  !
  ! PSMILe_* routines:
  !
  ! Copyright 2007-2010, CERFACS, Toulouse, France.
  ! Copyright 2007-2010, SGI Germany, Munich, Germany.
  ! Copyright 2007-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
  !-----------------------------------------------------------------------
  !+
  !
  ! mo_time_base [module]
  !    routines to handle Julian dates and times
  !
  ! Version:    E5/R1.07+  04-February-2003
  !
  ! Authors:
  !   L. Kornblueh, MPI                prepared basic version
  !   I. Kirchner,  MPI July 2000      large revision
  !   L. Kornblueh, MPI August 2000    adapting for right real kind (dp)
  !                                    and other minor corrections, removed
  !                                    30 day support, removed print overload.
  !   L. Kornblueh, MPI February 2003  adaption to PRISM calendar
  !                                    dates before the change to the gregorian
  !                                    calendar are treated as gregorian as 
  !                                    well to have a consitent view on the 
  !                                    seasons in climte sense - introduce
  !                                    preprocessor flag IAU for the original
  !                                    code (IAU - International Astrophysical
  !                                    Union).
  !   I. Kirchner,  FUB February 2003  code review
  !
  !   R. Redler, NEC-CCRLE, 2003       adaptation to OASIS4 software
  !
  !   R. Redler, NEC-CCRLE, Aug 2004   simple calendars with 360 days per year
  !                                    and 365 days per year added. Julian days
  !                                    calulated by theses calendars are not
  !                                    compatible with those calculated by the
  !                                    proleptic gregorian one.
  !                                    The new calendars need to be revised.
  !
  !-
#undef IAU
!!!#define IAU
  !
  ! References:
  !
  ! Montenbruck, Oliver, "Practical Ephmeris Calculations", Ch. 2, pp 33-40.
  ! The 1992 Astronomical Almanac, page B4.
  !
  ! The Julian Date is defined as the the number of days which have elapsed
  ! since the 1st of January of the year 4713 BC 12:00 Universal Time.
  !
  ! Up to 4th October 1582 AD the Julian Calendar was in force with leap
  ! years every four years, but from then on the Gregorian Calendar carried 
  ! on from 15th October 1582 with the leap years defined by the rule:
  !
  ! Leap year is every year whose yearly number is divisible by four, but
  ! not by a hundred, or is divisible by four hundred.
  !
  ! At midday on 4th October 1582 2299160.5 Julian days had elapsed.
  ! The Gregorian Calendar then carried on at this point from 15th October
  ! 1582, so its beginning occured on the Julian date 2299160.5.
  ! 
  ! Note: the astronomical year -4712 corresponds to the year 4713 BC, the
  !       year 0 to the year 1 BC; thereafter the astronomical year match
  !       the year AD, e.g. year 1 = year 1 AD.
  !
  !       This routines work for the years -5877402 BC until 5868098 AD.  
  !       However, there are not covered dates by the current GRIB edition 1, 
  !       which can handle dates between 1 AD and 25599 AD.
  !
  !       The "Modified Julian Date" is the Julian Date - 2400000.5 and so 
  !       has a zero point of 17th November 1858 AD 00:00 Universal Time.
  ! 
  !       Due to the small time span coverable by the GRIB output there is 
  !       no need to use a "Modified Julian Date".
  !
  ! The Julian day number is stored as two doubles to guarantee sufficient
  ! precision on all machines. The complete value is (day+fraction)
  ! although this addition will sometimes lose precision. Note that day
  ! might not be an integer value and fraction might be greater than one.
  !
  !----------------------------------------------------------------------
  USE PRISM_Constants
  USE PSMILe
 
  IMPLICIT NONE
  !
  ! WORKAROUND:
  ! We assume that SELECTED_REAL_KIND(12,307) corresponds
  ! to double precision on all systems.
  !
  INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12,307)

  LOGICAL, PARAMETER :: eq_months = .false.
  LOGICAL, PARAMETER :: eq_years  = .false.

  PRIVATE

  !----------------------------------------------------------------------
  ! Subroutines/Functions
  !
  ! Note: negative years are valid
  !
  PUBLIC  :: PSMILe_date2ju     ! convert PRISM_Time_Struct into Julian day
  PUBLIC  :: PSMILe_ju2date     ! convert Julian day into PRISM_Time_Struct
  PUBLIC  :: PSMILe_addtime     ! adds time
  PUBLIC  :: PSMILe_subtime     ! subtracts time
  !
  PRIVATE :: Set_JulianDay      ! convert calendar date/time into Julian day
  PRIVATE :: Set_JulianCalendar ! convert Julian day into calendar date/time
  PRIVATE :: Get_JulianYearDay  ! select the day number of the year
  PRIVATE :: Get_JulianYearLen  ! get the length of the year
  PRIVATE :: Get_JulianMonLen   ! get the length of the months
  PRIVATE :: Print_JulianDay    ! print out the date information
  PRIVATE :: sec2frac, frac2sec ! convert sec to fraction and vice versa

  !----------------------------------------------------------------------
  ! Structure declarations, variables
  !
  TYPE, PRIVATE :: julian_date
    !+
    !
    ! julian_date [structure]
    !     day      [real_dp]  (day in calendar)
    !     fraction [real_dp]  (fraction of the day)
    !
    !-
    Double Precision :: day       ! day since 1. January 4713 BC, 12UTC
    Double Precision :: fraction  ! fraction of day
  END TYPE julian_date

  !+
  !
  ! idaylen [integer, constant]
  !    length of a day in seconds
  !
  !-
  INTEGER, PRIVATE, PARAMETER :: IDAYLEN = 86400

CONTAINS
  !
  !----------------------------------------------------------------------
  !
  subroutine psmile_addtime ( date_in, date_inc, date_out )
    !+
    !
    ! PSMILe_addtime [subroutine]
    ! Adds time to a given date/time in julian days and seconds per day
    ! 1st array element keeps the days, the 2nd keeps the seconds
    !    (
    !    date_in      [double]  input (input date/time)
    !    date_inc     [double]  input (increment date/time)
    !    date_out     [double]  output (resulting date/time)
    !    )
    !
    !-
    !
    IMPLICIT NONE
    !
    Double Precision, INTENT(In)  :: date_in(2)
    Double Precision, INTENT(In)  :: date_inc(2)
    Double Precision, INTENT(Out) :: date_out(2)
    !
    INTEGER               :: days
    !
    days = date_inc(2)/86400.0_dp
    date_out(1) = date_in(1)+date_inc(1)+days
    date_out(2) = date_in(2)+date_inc(2)-(days*86400.0_dp)

  END SUBROUTINE PSMILe_addtime

  subroutine psmile_subtime ( date_in, date_inc, date_out )
    !+
    !
    ! PSMILe_subtime [subroutine]
    ! Subtracts time from a given date/time in julian days and seconds per day
    ! 1st array element keeps the days, the 2nd keeps the seconds
    !    (
    !    date_in      [double]  input (input date/time)
    !    date_inc     [double]  input (increment date/time)
    !    date_out     [double]  output (resulting date/time)
    !    )
    !
    !-
    !
    IMPLICIT NONE
    !
    Double Precision, INTENT(In)  :: date_in(2)
    Double Precision, INTENT(In)  :: date_inc(2)
    Double Precision, INTENT(Out) :: date_out(2)
    !
    INTEGER               :: days
    !
    days = date_inc(2)/86400.0_dp
    date_out(1) = date_in(1)-date_inc(1)-days
    date_out(2) = date_in(2)-date_inc(2)+(days*86400.0_dp)

  END SUBROUTINE PSMILe_subtime

  subroutine psmile_date2ju ( date, julian_day, julian_seconds )
    !+
    !
    ! PSMILe_date2ju [subroutine]
    !    convert PRISM time struct into Julian calendar day and seconds
    !    (
    !    date           [type]    input (PRISM date/time)
    !    julian_day     [double]  output (julian day)
    !    julian_seconds [double]  output (seconds of the julian day)
    !    )
    !
    !-
    !
    IMPLICIT NONE
    !
    !     INPUT
    !
    TYPE (PRISM_Time_Struct), INTENT(In) :: date
    !
    !     OUTPUT
    !
    REAL (PSMILe_float_kind), INTENT(Out) :: julian_day
    REAL (PSMILe_float_kind), INTENT(Out) :: julian_seconds
    !
    !     LOCAL
    !
    TYPE (julian_date) :: my_day

    julian_seconds = real( (date%hour*60.0 + date%minute) * 60.0, dp ) + date%second

#ifdef DEBUGX
    print *, trim(ch_id), ': PSMILe_date2ju: ', date%year, date%month, date%day
    print *, trim(ch_id), ': PSMILe_date2ju: ', date%hour, date%minute, date%second
#endif

    CALL Set_JulianDay ( date%year, date%month, date%day, julian_seconds, my_day )
    !
    ! Shift seconds since Julian days start at noon
    ! For reference: 1. January 1998 00 UTC === Julian Day 2450814, Julian Sec 43200
    !
    if ( julian_seconds >= 43200.0_dp ) then
         julian_seconds = julian_seconds - 43200.0_dp
    else
         julian_seconds = julian_seconds + 43200.0_dp
    endif

    julian_day = my_day%day

#ifdef DEBUGX
    print *, trim(ch_id), ': PSMILe_date2ju: ', julian_day, julian_seconds
#endif

  END SUBROUTINE PSMILe_date2ju
  !
  !----------------------------------------------------------------------
  !
  subroutine psmile_ju2date ( date, julian_day, julian_seconds )
    !+
    !
    ! PSMILe_ju2date [subroutine]
    !    convert Julian calendar day and seconds into PRISM time struct
    !    (
    !    date           [type]    output (PRISM date/time)
    !    julian_day     [double]  input (julian day)
    !    julian_seconds [double]  input (seconds of the julian day)
    !    )
    !
    !-
    !
    IMPLICIT NONE
    !
    !     INPUT
    !
    REAL (PSMILe_float_kind), INTENT(In) :: julian_day
    REAL (PSMILe_float_kind), INTENT(In) :: julian_seconds
    !
    !     OUTPUT
    !
    TYPE (PRISM_Time_Struct), INTENT(Out) :: date
    !
    !     LOCAL
    !
    TYPE (julian_date) :: my_day
    Double Precision           :: sec

    my_day%day      = julian_day
    my_day%fraction = sec2frac(julian_seconds)

    CALL Set_JulianCalendar(my_day, date%year, date%month, date%day, sec)
    !
    ! Seconds are returned but are ignored here since they are truncated
    !
    sec = julian_seconds
    if ( sec < 43200.0_dp ) then
         sec = sec + 43200.0_dp
    else
         sec = sec - 43200.0_dp
    endif

    date%hour   = INT(sec/REAL(3600,dp))
    date%minute = INT((sec-date%hour*REAL(3600,dp))/REAL(60,dp) )
    date%second = sec-date%hour*REAL(3600,dp)-date%minute*REAL(60,dp)
#ifdef DEBUGX
    print *, trim(ch_id), ': PSMILe_ju2date: ', my_day%day, my_day%fraction
    print *, trim(ch_id), ': PSMILe_ju2date: ', date%year, date%month, date%day
    print *, trim(ch_id), ': PSMILe_ju2date: ', date%hour, date%minute, date%second
#endif
  END SUBROUTINE PSMILe_ju2date
  !
  !----------------------------------------------------------------------
  !
  SUBROUTINE Set_JulianDay(ky, km, kd, sec, my_day)
    !+
    !
    ! Set_JulianDay  [subroutine]
    !    convert year, month, day, seconds into Julian calendar day
    !    (
    !    year   [integer] input (calendar year)
    !    month  [integer] input (month of the year)
    !    day    [integer] input (day of the month)
    !    second [dble]    input (seconds of the day)
    !    date   [julian_date] output (Julian day)
    !    )
    !
    !-
    !
    INTEGER, INTENT(IN)  :: ky
    INTEGER, INTENT(IN)  :: km
    INTEGER, INTENT(IN)  :: kd
    Double Precision, INTENT(IN) :: sec
    !
    ! for reference: 1. January 1998 00 UTC === Julian Day 2450814.5
    !
    TYPE (julian_date), INTENT(out) :: my_day

    INTEGER :: ib, iy, im, idmax
    Double Precision :: zd, zsec

    Double Precision :: yr_len
    Double Precision :: mon_len

    Double Precision :: mon_lens(12)

    DATA mon_lens / 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 /

    IF ( sec < 0 .OR. sec > 86400.0 ) THEN
      WRITE(*,*) 'mo_time_base: Set_JulianDay', ' invalid number of seconds', sec
      call psmile_abort
    ENDIF
    zsec = sec2frac(sec)

    IF ( eq_years .and. eq_months ) THEN
          mon_len = 30.0000_dp
          yr_len  = 360.0000_dp
          iy = ky
          im = km
    ELSE IF ( eq_years .and. .not. eq_months ) THEN
          yr_len  = 365.00_dp
          iy = ky
          im = km
    ELSE
       IF (km <= 2) THEN
          iy = ky-1
          im = km+12
       ELSE 
          iy = ky
          im = km
       ENDIF
       mon_len = 30.6001_dp
       yr_len  = 365.25_dp
    ENDIF

    ib = INT(iy/400)-INT(iy/100)

#ifdef IAU
    IF (ky > 1582 .OR. (ky == 1582 .AND. km > 10 &
         .OR. (km == 10 .AND. kd >= 15))) THEN
      ! 15th October 1582 AD or later no changes

    ELSE IF (ky == 1582 .AND. km == 10 .AND. (4 < kd .AND. kd < 15) ) THEN
      ! a small pice from the history:
      !     Papst Gregor XIII corrected the calendar,
      !     he skipped 10 days between the 4th and the 15th of October 1582
      !
      WRITE (*, '(a29,a5,i2.2,a1,i2.2,a1,i4.4,a)') &
           'mpi-calendar: Set_JulianDay: '         &
           'date ', km, '/', kd, '/', ky, ' not valid'
      call psmile_abort
    ELSE 

      ! 4th October 1582 AD or earlier
      ib = -2

    ENDIF
#endif

    ! check the length of the month
    idmax = Get_JulianMonLen (ky, km)

    IF (kd < 1 .OR. idmax < kd) then
      WRITE(*,'(a,a,i10,i10)') 'mo_time_base: Set_JulianDay' ,' day in months invalid', kd, idmax
      call psmile_abort
    ENDIF

    IF ( eq_years .and. eq_months ) THEN

       zd = FLOOR(yr_len*iy)+INT(mon_len*(im-1)) &
            +REAL(kd,dp)-1.0_dp+zsec

    ELSE IF ( eq_years .and. .not. eq_months ) THEN

       zd = FLOOR(yr_len*iy)+INT(mon_lens(im)) &
            +REAL(kd,dp)+zsec

    ELSE

       zd = FLOOR(yr_len*iy)+INT(mon_len*(im+1)) &
            +REAL(ib,dp)+1720996.5_dp+REAL(kd,dp)+zsec

    ENDIF

    my_day%day      = AINT(zd,dp)
    my_day%fraction = zd-AINT(zd,dp)

  END SUBROUTINE Set_JulianDay

  SUBROUTINE Set_JulianCalendar(my_day, ky, km, kd, sec)
    !+
    !
    ! Set_JulianCalendar [subroutine]
    !    convert Julian date into year, months, day, seconds of day
    !    (
    !    date   [julian_date]  input (Julian day)
    !    year   [integer] output (calendar year)
    !    month  [integer] output (month of the year)
    !    day    [integer] output (day of the month)
    !    second [dble]    output (seconds of the day)
    !    )
    !
    !-
    TYPE (julian_date), INTENT(IN) :: my_day
    INTEGER, INTENT(OUT)           :: ky
    INTEGER, INTENT(OUT)           :: km
    INTEGER, INTENT(OUT)           :: kd
    Double Precision, INTENT(OUT)          :: sec
    
    Double Precision :: za, zb, zc, zd, ze, zf, zday, zfrac, zsec

    Double Precision :: yr_len
    Double Precision :: mon_len

    zday  = my_day%day
    zfrac = my_day%fraction

    IF ( eq_years .and. eq_months ) THEN

       mon_len = 30.0000_dp
       yr_len  = 360.0000_dp
       if ( zfrac >= 0.5 ) zfrac = zfrac - 0.5_dp
       za = FLOOR(zday+zfrac+0.5_dp)
       ky = FLOOR(za/yr_len)
       zb = za - (ky*yr_len)
       km = FLOOR(zb/mon_len) + 1
       zc = zb - ((km-1)*mon_len)
       kd = zc + 1
       RETURN

    ELSE IF ( eq_years .and. .not. eq_months ) THEN

       yr_len  = 365.0000_dp
       if ( zfrac >= 0.5 ) zfrac = zfrac - 0.5_dp
       za = FLOOR(zday+zfrac+0.5_dp)
       ky = FLOOR(za/yr_len)
       zb = za - (ky*yr_len)
       CALL Get_MonDay ( zb, km, kd )
       RETURN

    ELSE

       mon_len = 30.6001_dp
       yr_len  = 365.25_dp

    ENDIF

    za = FLOOR(zday+zfrac+0.5_dp)
    zb = FLOOR((za-1867216.25_dp)/36524.25_dp)
    zc = za+zb-FLOOR(zb/4)+1525

#ifdef IAU
    IF (za < 2299161.0_dp) zc = za+1524.0_dp
#endif

    zd = FLOOR((zc-122.1_dp)/365.25_dp)
    ze = FLOOR(365.25_dp*zd)
    zf = FLOOR((zc-ze)/30.6001_dp)

    kd = INT(zc-ze- FLOOR(30.6001_dp*zf))
    km = INT(zf-1-12*FLOOR(zf/14))
    ky = INT(zd-4715-((7+km)/10))

    ! convert fraction in seconds of the day
    IF (zfrac < -0.5_dp) THEN
       zsec = 1.5_dp + zfrac
    ELSE IF (zfrac < 0.5_dp) THEN
       zsec = 0.5_dp + zfrac
    ELSE
       zsec = zfrac - 0.5_dp
    END IF
    sec = frac2sec(zsec)

  END SUBROUTINE Set_JulianCalendar

  SUBROUTINE Get_JulianYearDay(my_day, kjy, kd, sec)
    !+
    !
    ! Get_JulianYearDay [subroutine]
    !    get Julian year, day of year and seconds of day from Julian date
    !    (
    !    date   [julian_date]  input (Julian day)
    !    year   [integer] output (Julian calendar year)
    !    day    [integer] output (day of the year)
    !    second [dble]    output (seconds of the day)
    !    )
    !
    !-
    ! Day of Astronomical Year (approximated)
    !
    ! it is used for the correct annual cycle
    ! remember the mismatch between calendar and
    ! astronomical year accumulated until 1582
    !
    ! this version gives adjustment until 1583
    !
    TYPE (julian_date), INTENT(IN) :: my_day
    INTEGER, INTENT(OUT)           :: kjy
    INTEGER, INTENT(OUT)           :: kd
    Double Precision, INTENT(OUT)          :: sec

    TYPE (julian_date) :: first_jan, present_day
    INTEGER :: iy, id, im

    ! AD 1 is year 4713 in Julian calender, 31st of December 1997 the Julian 
    ! day 2450814.0 is elapsed at 12 UTC and the Julian year is 6711, so:
    ! year_len = 2450814. / 6710. = 365.25

    ! get the date of the present julian day
    CALL Set_JulianCalendar(my_day, iy, im, id, sec)
    
    ! find the first of January 00UTC
    CALL Set_JulianDay(iy,  1,  1, 0.0_dp, first_jan)

    ! find the present Julian Day at 00UTC
    CALL Set_JulianDay(iy, im, id, 0.0_dp, present_day)

    ! the day in the year is given
    kd = present_day%day-first_jan%day+1

    ! get the Julian year
    kjy = iy+4712

  END SUBROUTINE Get_JulianYearDay

  INTEGER FUNCTION Get_JulianYearLen(ky)
    !+
    !
    ! Get_JulianYearLen  [function, integer]
    !    get the length of a Julian year in days
    !    (
    !    year [integer] input (Calendar year)
    !    )
    !
    !-
    INTEGER, INTENT(in) :: ky

    INTEGER :: ylen

    IF ( eq_months .and. eq_years ) THEN
        Get_JulianYearLen = 360
        RETURN
    ENDIF
    IF ( eq_years ) THEN
        Get_JulianYearLen = 365
        RETURN
    ENDIF

    IF (ky == 1582) THEN
       ylen = 355
    ELSE IF ( (MOD(ky,4)==0 .AND. MOD(ky,100)/=0) .OR. MOD(ky,400)==0 ) THEN
       ylen = 366
    ELSE
       ylen = 365
    END IF
    Get_JulianYearLen = ylen

  END FUNCTION Get_JulianYearLen

  INTEGER FUNCTION Get_JulianMonLen(ky, km)
    !+
    !
    ! Get_JulianMonLen [function, integer]
    !    get the length of a months in a Julian year
    !    (
    !    year  [integer] input (Calendar year)
    !    month [integer] input (month of the year)
    !    )
    !
    !-
    INTEGER, INTENT(in) :: km, ky

    INTEGER :: idmax

    IF ( eq_months ) THEN
        Get_JulianMonLen = 30
        RETURN
    ENDIF

    SELECT CASE(km)
    CASE(1,3,5,7,8,10,12);  idmax = 31
    CASE(4,6,9,11);         idmax = 30
    CASE(2)
      IF ( (MOD(ky,4)==0 .AND. MOD(ky,100)/=0) .OR. MOD(ky,400)==0 ) THEN
        ! leap year found
        idmax = 29
      ELSE
        idmax = 28
      END IF

    CASE default
      WRITE(*,'(a,a)') 'mo_time_base: Get_JulianMonLen', ' month invalid', km
      call psmile_abort
    END SELECT
    Get_JulianMonLen = idmax

  END FUNCTION Get_JulianMonLen

  SUBROUTINE Print_JulianDay(my_day)
    !+
    !
    ! Print_JulianDay [subroutine] interface [print]
    !   print Julian/model day on standard output
    !   (
    !   date [julian_date] input (model/Julian calendar day)
    !   )
    !
    !-
    TYPE (julian_date), INTENT(in) :: my_day

    WRITE(*,*) 'Print_JulianDay: ', my_day%day+my_day%fraction

  END SUBROUTINE Print_JulianDay

  FUNCTION sec2frac(sec) RESULT (zfrac)
    !+
    !
    ! sec2frac [function, real_dp]
    !    convert seconds of day into fraction
    !    (
    !    second [dble] input (seconds of the day)
    !    )
    !
    !-
    Double Precision, INTENT(in) :: sec
    Double Precision             :: zfrac

    zfrac = sec/REAL(IDAYLEN,dp)

  END FUNCTION Sec2frac

  FUNCTION frac2sec(zfrac) RESULT (sec)
    !+
    !
    ! frac2sec [function, dble]
    !    convert fraction of day into seconds
    !    (
    !    fraction [real_dp] input (fraction of the day)
    !    )
    !
    !-
    Double Precision, INTENT(in) :: zfrac
    Double Precision             :: sec

    sec = INT(zfrac*REAL(IDAYLEN,dp))

  END FUNCTION frac2sec

  SUBROUTINE Get_MonDay ( days, km, kd )
    Double Precision, INTENT(IN) :: days
    INTEGER, INTENT(OUT) :: km 
    INTEGER, INTENT(OUT) :: kd

    INTEGER              :: m

    Double Precision :: mon_lens(12)

    DATA mon_lens / 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 /

    DO m = 1, 11
       IF ( days <= mon_lens(m+1) ) EXIT
    END DO

    kd = FLOOR(days-mon_lens(m))
    km = m

  END SUBROUTINE Get_MonDay

END MODULE PRISM_calendar