scale_calendar.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_calendar
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: calendar_setup
  public :: calendar_getdayofyear
  public :: calendar_date2daysec
  public :: calendar_daysec2date
  public :: calendar_ymd2absday
  public :: calendar_hms2abssec
  public :: calendar_adjust_daysec
  public :: calendar_combine_daysec
  public :: calendar_unit2sec
  public :: calendar_sec2unit
  public :: calendar_cfunits2sec
  public :: calendar_date2char
  public :: calendar_get_name

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public, parameter :: i_year  = 1
  integer, public, parameter :: i_month = 2
  integer, public, parameter :: i_day   = 3
  integer, public, parameter :: i_hour  = 4
  integer, public, parameter :: i_min   = 5
  integer, public, parameter :: i_sec   = 6

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: calendar_absday2ymd
  private :: calendar_abssec2hms
  private :: calendar_ymdhms2nd
  private :: calendar_ymdhms2mjd
  private :: checkleap

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical,  private :: calendar_360days = .false.
  logical,  private :: calendar_365days = .false.

  real(DP), private :: calendar_doi  = 365.0_dp 
  real(DP), private :: calendar_hour = 24.0_dp  
  real(DP), private :: calendar_min  = 60.0_dp  
  real(DP), private :: calendar_sec  = 60.0_dp  

  integer,  private, parameter :: i_nonleapyear = 1
  integer,  private, parameter :: i_leapyear    = 2
  integer,  private, parameter :: i_360days     = 3
  integer,  private            :: dayofmonth(12,3)
  data dayofmonth / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, & ! non-leap year
                    31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, & ! leap year
                    30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30  / ! 360 days

  logical,  private :: debug = .false.

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine calendar_setup
    use scale_prc, only: &
       prc_abort
    implicit none

    namelist / param_calendar / &
       calendar_360days, &
       calendar_365days, &
       debug

    integer :: ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("CALENDAR_setup",*) 'Setup'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_calendar,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("CALENDAR_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("CALENDAR_setup",*) 'Not appropriate names in namelist PARAM_CALENDAR. Check!'
       call prc_abort
    endif
    log_nml(param_calendar)

    if    ( calendar_360days ) then
       calendar_doi = 360.0_dp
    elseif( calendar_365days ) then
       calendar_doi = 365.0_dp
    endif

    log_newline
    log_info("CALENDAR_setup",*) 'Calendar settings'
    if    ( calendar_360days ) then
       log_info_cont(*) 'DayOfYear = 360 : ideal setting'
    elseif( calendar_365days ) then
       log_info_cont(*) 'DayOfYear = 365 : ideal setting'
    else
       log_info_cont(*) 'DayOfYear = 365 or 366 : Gregorian calendar'
    endif

    return
  end subroutine calendar_setup

  !-----------------------------------------------------------------------------
  subroutine calendar_getdayofyear( &
       DayOfYear, &
       iyear      )
    implicit none

    real(DP), intent(out) :: DayOfYear ! # of day of year
    integer,  intent(in)  :: iyear     ! current year
    !---------------------------------------------------------------------------

    dayofyear = calendar_doi
    if( checkleap(iyear) ) dayofyear = calendar_doi + 1.0_dp

    return
  end subroutine calendar_getdayofyear

  !-----------------------------------------------------------------------------
  subroutine calendar_date2daysec( &
       absday,     &
       abssec,     &
       ymdhms,     &
       subsec,     &
       offset_year )
    implicit none

    integer,  intent(out) :: absday
    real(DP), intent(out) :: abssec
    integer,  intent(in)  :: ymdhms(6)
    real(DP), intent(in)  :: subsec
    integer,  intent(in)  :: offset_year
    !---------------------------------------------------------------------------

    call calendar_ymd2absday( absday,          & ! [OUT]
                              ymdhms(i_year),  & ! [IN]
                              ymdhms(i_month), & ! [IN]
                              ymdhms(i_day),   & ! [IN]
                              offset_year      ) ! [IN]

    call calendar_hms2abssec( abssec,          & ! [OUT]
                              ymdhms(i_hour),  & ! [IN]
                              ymdhms(i_min),   & ! [IN]
                              ymdhms(i_sec),   & ! [IN]
                              subsec           ) ! [IN]

    return
  end subroutine calendar_date2daysec

  !-----------------------------------------------------------------------------
  subroutine calendar_daysec2date( &
       ymdhms,     &
       subsec,     &
       absday,     &
       abssec,     &
       offset_year )
    implicit none

    integer,  intent(out) :: ymdhms(6)
    real(DP), intent(out) :: subsec
    integer,  intent(in)  :: absday
    real(DP), intent(in)  :: abssec
    integer,  intent(in)  :: offset_year
    !---------------------------------------------------------------------------

    call calendar_absday2ymd( ymdhms(i_year),  & ! [OUT]
                              ymdhms(i_month), & ! [OUT]
                              ymdhms(i_day),   & ! [OUT]
                              absday,          & ! [IN]
                              offset_year      ) ! [IN]

    call calendar_abssec2hms( ymdhms(i_hour),  & ! [OUT]
                              ymdhms(i_min),   & ! [OUT]
                              ymdhms(i_sec),   & ! [OUT]
                              subsec,          & ! [OUT]
                              abssec           ) ! [IN]

    return
  end subroutine calendar_daysec2date

  !-----------------------------------------------------------------------------
  subroutine calendar_ymd2absday( &
       absday, &
       gyear,  &
       gmonth, &
       gday,   &
       oyear   )
    implicit none

    integer, intent(out) :: absday
    integer, intent(in)  :: gyear
    integer, intent(in)  :: gmonth
    integer, intent(in)  :: gday
    integer, intent(in)  :: oyear

    integer :: gyear_mod, gmonth_mod
    integer :: yearday, monthday
    integer :: m, ileap
    !---------------------------------------------------------------------------

    gmonth_mod = mod( gmonth-1, 12 ) + 1
    gyear_mod  = gyear + ( gmonth-gmonth_mod ) / 12

    if ( calendar_360days .OR. calendar_365days ) then
       yearday = int( calendar_doi * ( gyear_mod - oyear ) )
    else
       yearday = int( calendar_doi * ( gyear_mod - oyear ) ) &
               + int( real(gyear_mod-1,kind=DP) /   4.0_DP ) &
               - int( real(gyear_mod-1,kind=DP) / 100.0_dp ) &
               + int( real(gyear_mod-1,kind=DP) / 400.0_dp ) &
               - int( real(oyear    -1,kind=DP) /   4.0_dp ) &
               + int( real(oyear    -1,kind=DP) / 100.0_dp ) &
               - int( real(oyear    -1,kind=DP) / 400.0_dp )
    endif

    ileap = i_nonleapyear
    if( checkleap(gyear_mod) ) ileap = i_leapyear
    if( calendar_360days     ) ileap = i_360days

    monthday = 0
    do m = 1, gmonth_mod-1
       monthday = monthday + dayofmonth(m,ileap)
    enddo

    absday = yearday + monthday + gday - 1

    return
  end subroutine calendar_ymd2absday

  !-----------------------------------------------------------------------------
  subroutine calendar_absday2ymd( &
       gyear,  &
       gmonth, &
       gday,   &
       absday, &
       oyear   )
    implicit none

    integer, intent(out) :: gyear
    integer, intent(out) :: gmonth
    integer, intent(out) :: gday
    integer, intent(in)  :: absday
    integer, intent(in)  :: oyear

    integer :: checkday
    integer :: i, ileap
    !---------------------------------------------------------------------------

    if ( calendar_360days ) then
       gyear = int( real(absday,kind=DP) / 361.0_DP ) + oyear ! first guess
    else
       gyear = int( real(absday,kind=DP) / 366.0_DP ) + oyear ! first guess
    endif

    do i = 1, 1000
       call calendar_ymd2absday( checkday, gyear+1, 1, 1, oyear )
       if( absday < checkday ) exit
       gyear = gyear + 1
    enddo

    ileap = i_nonleapyear
    if( checkleap(gyear) ) ileap = i_leapyear
    if( calendar_360days ) ileap = i_360days

    gmonth = 1
    do i = 1, 1000
       call calendar_ymd2absday( checkday, gyear, gmonth, dayofmonth(gmonth,ileap), oyear )
       if( absday <= checkday ) exit
       gmonth = gmonth + 1
    enddo

    call calendar_ymd2absday( checkday, gyear, gmonth, 1, oyear )
    gday = absday - checkday + 1

    return
  end subroutine calendar_absday2ymd

  !-----------------------------------------------------------------------------
  subroutine calendar_hms2abssec( &
       abssec, &
       hour,   &
       minute, &
       second, &
       subsec  )
    implicit none

    real(DP), intent(out) :: abssec
    integer,  intent(in)  :: hour
    integer,  intent(in)  :: minute
    integer,  intent(in)  :: second
    real(DP), intent(in)  :: subsec
    !---------------------------------------------------------------------------

    abssec = real(hour,  kind=DP) * CALENDAR_MIN * CALENDAR_SEC &
           + real(minute,kind=DP) * CALENDAR_SEC                &
           + real(second,kind=DP)                               &
           + subsec

    return
  end subroutine calendar_hms2abssec

  !-----------------------------------------------------------------------------
  subroutine calendar_abssec2hms( &
       hour,   &
       minute, &
       second, &
       subsec, &
       abssec  )
    implicit none

    integer,  intent(out) :: hour
    integer,  intent(out) :: minute
    integer,  intent(out) :: second
    real(DP), intent(out) :: subsec
    real(DP), intent(in)  :: abssec

    real(DP) :: nsec, nmin, nhour, temp
    !---------------------------------------------------------------------------

    nsec   = real(int(abssec),kind=dp)
    subsec = abssec - nsec

    temp   = mod( nsec, calendar_sec )
    second = int( temp )
    nmin   = ( nsec-temp ) / calendar_sec

    temp   = mod( nmin, calendar_min )
    minute = int( temp )
    nhour  = ( nmin-temp ) / calendar_min

    temp   = mod( nhour, calendar_hour )
    hour   = int( temp )

    return
  end subroutine calendar_abssec2hms

  !-----------------------------------------------------------------------------
  subroutine calendar_adjust_daysec( &
       absday, &
       abssec  )
    implicit none

    integer,  intent(inout) :: absday
    real(DP), intent(inout) :: abssec

    integer :: addday
    !---------------------------------------------------------------------------

    addday = int( abssec / ( calendar_hour * calendar_min * calendar_sec ) )

    absday = absday + addday

    abssec = abssec - real(addday,kind=DP) * CALENDAR_HOUR * CALENDAR_MIN * CALENDAR_SEC

    if ( abssec < 0.0_dp ) then
       absday = absday - 1
       abssec = abssec + calendar_hour * calendar_min * calendar_sec
    endif

    return
  end subroutine calendar_adjust_daysec

  !-----------------------------------------------------------------------------
  function calendar_combine_daysec( absday, abssec ) result(daysec)
    implicit none

    integer,  intent(in) :: absday
    real(DP), intent(in) :: abssec
    real(DP)             :: daysec
    !---------------------------------------------------------------------------

    daysec = real(absday,kind=DP) * CALENDAR_SEC * CALENDAR_MIN * CALENDAR_HOUR &
           + abssec

    return
  end function calendar_combine_daysec

  !-----------------------------------------------------------------------------
  subroutine calendar_unit2sec( &
       second, &
       value,  &
       unit    )
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP),         intent(out) :: second
    real(DP),         intent(in)  :: value
    character(len=*), intent(in)  :: unit
    !---------------------------------------------------------------------------

    select case(unit)
    case('MSEC')
       second = value * 1.e-3_dp
    case('SEC', 'seconds')
       second = value
    case('MIN')
       second = value * calendar_sec
    case('HOUR')
       second = value * calendar_sec * calendar_min
    case('DAY')
       second = value * calendar_sec * calendar_min * calendar_hour
    case default
       log_error("CALENDAR_unit2sec",*) 'Unsupported UNIT: ', trim(unit), ', ', value
       call prc_abort
    endselect

    return
  end subroutine calendar_unit2sec

  !-----------------------------------------------------------------------------
  subroutine calendar_sec2unit( &
     value,  &
     second, &
     unit    )
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP),         intent(out) :: value
    real(DP),         intent( in) :: second
    character(len=*), intent( in) :: unit
    !---------------------------------------------------------------------------

    select case(trim(unit))
    case('MSEC', 'msec')
       value = second / 1.0e-3_dp
    case('SEC', 'seconds', 'sec', 's')
       value = second
    case('MIN', 'mins', 'min')
       value = second / calendar_sec
    case('HOUR', 'hours', 'hour', 'h')
       value = second / (calendar_sec * calendar_min)
    case('DAY', 'days', 'day')
       value = second / (calendar_sec * calendar_min * calendar_hour)
    case default
       log_error("CALENDAR_sec2unit",*) 'Unsupported UNIT: ', trim(unit), ', ', value
       call prc_abort
    endselect

  end subroutine calendar_sec2unit

  !-----------------------------------------------------------------------------
  function calendar_cfunits2sec( cftime, cfunits, offset_year, startdaysec ) result( sec )
    use scale_prc, only: &
       prc_abort
    implicit none

    real(DP),         intent(in) :: cftime
    character(len=*), intent(in) :: cfunits
    integer,          intent(in) :: offset_year
    real(DP),         intent(in), optional :: startdaysec
    real(DP)                     :: sec

    character(len=H_MID) :: tunit
    character(len=H_MID) :: buf

    integer  :: date(6)
    integer  :: day
    real(DP) :: sec0

    integer :: l

    intrinsic index
    !---------------------------------------------------------------------------

    l = index( cfunits, " since " )
    if ( l > 1 ) then ! untis is under the CF convension
       tunit = cfunits(1:l-1)
       buf   = cfunits(l+7:)

       l = index(buf,"-")
       if ( l /= 5 ) then
          log_error("CALENDAR_CFunits2sec",*) 'units for time is invalid (year)'
          log_error_cont(*) trim(cfunits)
          log_error_cont(*) trim(buf)
          call prc_abort
       end if
       read(buf(1:4),*) date(1) ! year
       buf = buf(6:)

       l = index(buf,"-")
       if ( l /= 3 ) then
          log_error("CALENDAR_CFunits2sec",*) 'units for time is invalid (month)'
          log_error_cont(*) trim(cfunits)
          log_error_cont(*) trim(buf)
          call prc_abort
       end if
       read(buf(1:2),*) date(2) ! month
       buf = buf(4:)

       l = index(buf," ")
       if ( l /= 3 ) then
          log_error("CALENDAR_CFunits2sec",*) 'units for time is invalid (day)'
          log_error_cont(*) trim(cfunits)
          log_error_cont(*) trim(buf)
          call prc_abort
       end if
       read(buf(1:2),*) date(3) ! day
       buf = buf(4:)

       l = index(buf,":")
       if ( l /= 3 ) then
          log_error("CALENDAR_CFunits2sec",*) 'units for time is invalid (hour)'
          log_error_cont(*) trim(cfunits)
          log_error_cont(*) trim(buf)
          call prc_abort
       end if
       read(buf(1:2),*) date(4) ! hour
       buf = buf(4:)

       l = index(buf,":")
       if ( l /= 3 ) then
          log_error("CALENDAR_CFunits2sec",*) 'units for time is invalid (min)'
          log_error_cont(*) trim(cfunits)
          log_error_cont(*) trim(buf)
          call prc_abort
       end if
       read(buf(1:2),*) date(5) ! min
       buf = buf(4:)

       if ( len_trim(buf) /= 2 ) then
          log_error("CALENDAR_CFunits2sec",*) 'units for time is invalid (sec)'
          log_error_cont(*) trim(cfunits)
          log_error_cont(*) trim(buf)
          log_error_cont(*) len_trim(buf)
          call prc_abort
       end if
       read(buf(1:2),*) date(6) ! sec

       call calendar_date2daysec( day,        & ! (out)
                                  sec0,       & ! (out)
                                  date(:),    & ! (in)
                                  0.0_dp,     & ! (in)
                                  offset_year ) ! (in)

       sec0 = calendar_combine_daysec( day, sec0 )
    else
       tunit = cfunits
       if ( present(startdaysec) ) then
          sec0 = startdaysec
       else
          sec0 = 0.0_dp
       end if
    end if

    call calendar_unit2sec( sec, cftime, tunit )

    sec = sec0 + sec

    return
  end function calendar_cfunits2sec

  !-----------------------------------------------------------------------------
  subroutine calendar_date2char( &
       chardate, &
       ymdhms,   &
       subsec    )
    implicit none

    character(len=27), intent(out) :: chardate
    integer,           intent(in)  :: ymdhms(6)
    real(DP),          intent(in)  :: subsec
    !---------------------------------------------------------------------------

    write(chardate,'(I4.4,A1,I2.2,A1,I2.2,A1,I2.2,A1,I2.2,A1,I2.2,A2,F6.3)')   &
                       ymdhms(1),'/',ymdhms(2),'/',ymdhms(3),' ',  &
                       ymdhms(4),':',ymdhms(5),':',ymdhms(6),' +', &
                       subsec

    return
  end subroutine calendar_date2char

  !-----------------------------------------------------------------------------
  function checkleap( iyear )
    implicit none

    integer, intent(in) :: iyear
    logical             :: checkleap

    integer :: check4, check100, check400
    !---------------------------------------------------------------------------

    check4   = mod(iyear,4  )
    check100 = mod(iyear,100)
    check400 = mod(iyear,400)

    checkleap = .false.
    if( check4   == 0 ) checkleap = .true.
    if( check100 == 0 ) checkleap = .false.
    if( check400 == 0 ) checkleap = .true.

    if( calendar_360days ) checkleap = .false.
    if( calendar_365days ) checkleap = .false.

  end function checkleap

  !-----------------------------------------------------------------------------
  subroutine calendar_ymdhms2nd( &
       nd,     &
       ymdhms, &
       oyear   )
    implicit none

    real(DP), intent(out) :: nd
    integer,  intent(in)  :: ymdhms(6)
    integer,  intent(in)  :: oyear

    integer :: absday, absday_jan1
    !---------------------------------------------------------------------------

    call calendar_ymd2absday( absday,          & ! [OUT]
                              ymdhms(i_year),  & ! [IN]
                              ymdhms(i_month), & ! [IN]
                              ymdhms(i_day),   & ! [IN]
                              oyear            ) ! [IN]

    call calendar_ymd2absday( absday_jan1,     & ! [OUT]
                              ymdhms(i_year),  & ! [IN]
                              1,               & ! [IN]
                              1,               & ! [IN]
                              oyear            ) ! [IN]

    nd = absday - absday_jan1 + 1.0_dp

    return
  end subroutine calendar_ymdhms2nd

  !-----------------------------------------------------------------------------
  subroutine calendar_ymdhms2mjd( &
       mjd,    &
       ymdhms, &
       oyear   )
    implicit none

    real(DP), intent(out) :: mjd
    integer,  intent(in)  :: ymdhms(6)
    integer,  intent(in)  :: oyear

    integer :: y, m, mjd_day
    !---------------------------------------------------------------------------

    if (      ymdhms(i_month) == 1 &
         .OR. ymdhms(i_month) == 2 ) then
       y = ymdhms(i_year)  - 1
       m = ymdhms(i_month) + 2
    else
       y = ymdhms(i_year)
       m = ymdhms(i_month)
    endif

    mjd_day = int( 365.25_dp * y )                  & ! year
            + int( y/400.0_dp ) - int( y/100.0_dp ) & ! leap year
            + int( 30.59_dp * m-2 )                 & ! month
            + ymdhms(i_day)                         & ! day
            + 678912                                ! constant

    mjd     = real(mjd_day,kind=DP)       & ! day
            + ymdhms(i_hour) /    24.0_DP & ! hour
            + ymdhms(i_min)  /  1440.0_DP & ! min
            + ymdhms(i_sec)  / 86400.0_DP   ! sec

    return
  end subroutine calendar_ymdhms2mjd

  subroutine calendar_get_name(name)
    character(len=*), intent(out) :: name

    if    ( calendar_360days ) then
       name = "360_day"
    elseif( calendar_365days ) then
       name = "365_day"
    else
       name = "gregorian"
    endif

    return
  end subroutine calendar_get_name

end module scale_calendar