time_interp.inc

!***********************************************************************
!*                   GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
 
!> @addtogroup time_interp_mod
!> @{
 
 !> @brief Calculates the fractional time into the current year
 subroutine time_interp_frac_ ( Time, weight )
 
   type(time_type),   intent(in)  :: Time
   real(FMS_TI_KIND_), intent(out) :: weight !< fractional time
 
   integer         :: yr, mo, dy, hour, minute, second
   type(time_type) :: Year_beg, Year_end
 
 
   if ( .not. module_is_initialized ) call time_interp_init
 
!  ---- compute fractional time of year -----
 
     call get_date (time, yr, mo, dy, hour, minute, second)
 
     year_beg = set_date(yr  , 1, 1)
     year_end = set_date(yr+1, 1, 1)
 
     weight = real( (time - year_beg) // (year_end - year_beg) , kind=fms_ti_kind_)
 
 end subroutine time_interp_frac_
 
 
 !> @brief Calculates fractional time between mid points of consecutive years
 subroutine time_interp_year_ ( Time, weight, year1, year2 )
 
   type(time_type),   intent(in)  :: Time
   real(FMS_TI_KIND_), intent(out) :: weight !< fractional time between midpoints of year1 and year2
   integer        ,   intent(out) :: year1, year2
 
   integer :: yr, mo, dy, hour, minute, second
   type (time_type) :: Mid_year, Mid_year1, Mid_year2
 
 
   if ( .not. module_is_initialized ) call time_interp_init()
 
      call get_date (time, yr, mo, dy, hour, minute, second)
 
    ! mid point of current year
      mid_year = year_midpt(yr)
 
      if ( time >= mid_year ) then
    ! current time is after mid point of current year
           year1  = yr
           year2  = yr+1
           mid_year2 = year_midpt(year2)
           weight = real( (time - mid_year) // (mid_year2 - mid_year) , kind=fms_ti_kind_ )
      else
    ! current time is before mid point of current year
           year2  = yr
           year1  = yr-1
           mid_year1 = year_midpt(year1)
           weight = real( (time - mid_year1) // (mid_year - mid_year1), kind=fms_ti_kind_ )
      endif
 
 end subroutine time_interp_year_
 
 !> @brief Calculates fractional time between mid points of consecutive months
 subroutine time_interp_month_ ( Time, weight, year1, year2, month1, month2 )
 
   type(time_type), intent(in)  :: Time
   real(FMS_TI_KIND_)           , intent(out) :: weight
   integer        , intent(out) :: year1, year2, month1, month2
 
   integer :: yr, mo, dy, hour, minute, second,  &
              mid_month, cur_month, mid1, mid2
 
   if ( .not. module_is_initialized ) call time_interp_init()
 
      call get_date (time, yr, mo, dy, hour, minute, second)
 
    ! mid point of current month in seconds
      mid_month = days_in_month(time) * halfday
    ! time into current month in seconds
      cur_month = second + secmin*minute + sechour*hour + secday*(dy-1)
 
      if ( cur_month >= mid_month ) then
    ! current time is after mid point of current month
           year1  = yr;  month1 = mo
           year2  = yr;  month2 = mo+1
           if (month2 > monyear)  then
              year2 = year2+1;  month2 = 1
           endif
           mid1 = mid_month
           mid2 = days_in_month(set_date(year2,month2,2)) * halfday
           weight = real(cur_month - mid1, fms_ti_kind_) / real(mid1+mid2, fms_ti_kind_)
      else
    ! current time is before mid point of current month
           year2  = yr;  month2 = mo
           year1  = yr;  month1 = mo-1
           if (month1 < 1)  then
              year1 = year1-1;  month1 = monyear
           endif
           if (year1>0) then
              mid1 = days_in_month(set_date(year1,month1,2)) * halfday
           else
              ! this can happen if we are at the beginning of year 1. In this case
              ! use December 0001 to calculate the duration of December 0000.
              ! This should work for all calendars
              mid1 = days_in_month(set_date(1,month1,2)) * halfday
           endif
           mid2 = mid_month
           weight = real(cur_month + mid1, fms_ti_kind_) / real(mid1+mid2, fms_ti_kind_)
      endif
 
 end subroutine time_interp_month_
 
 !> @brief Calculates fractional time between mid points of consecutive days
 subroutine time_interp_day_ ( Time, weight, year1, year2, month1, month2, day1, day2 )
 
   type(time_type), intent(in)  :: Time
   real(FMS_TI_KIND_), intent(out) :: weight
   integer        , intent(out) :: year1, year2, month1, month2, day1, day2
 
   integer :: yr, mo, dy, hour, minute, second, sday
 
   if ( .not. module_is_initialized ) call time_interp_init()
 
      call get_date (time, yr, mo, dy, hour, minute, second)
 
    ! time into current day in seconds
      sday = second + secmin*minute + sechour*hour
 
      if ( sday >= halfday ) then
    ! current time is after mid point of day
           year1 = yr;  month1 = mo;  day1 = dy
           year2 = yr;  month2 = mo;  day2 = dy + 1
           weight  = real(sday - halfday, fms_ti_kind_) / real(secday, fms_ti_kind_)
 
           if (day2 > days_in_month(time)) then
               month2 = month2 + 1
               day2 = 1
               if (month2 > monyear) then
                    month2 = 1;  year2 = year2+1
               endif
           endif
      else
    ! current time is before mid point of day
           year2 = yr;  month2 = mo   ;  day2 = dy
           year1 = yr;  month1 = mo;  day1 = dy - 1
           weight  = real(sday + halfday,fms_ti_kind_) / real(secday,fms_ti_kind_)
 
           if (day1 < 1) then
               month1 = month1 - 1
               if (month1 < 1) then
                   month1 = monyear;  year1 = year1-1
               endif
               day1 = days_in_month(set_date(year1,month1,2))
           endif
      endif
 
 end subroutine time_interp_day_
 
 !> Part of the time_interp interface, calculates for cyclical data
 !! Time_beg and Time_end mark a repeating period
 !!
 !! Finds mid points and fractional weight for a time perioid
subroutine time_interp_modulo_(Time, Time_beg, Time_end, Timelist, weight, index1, index2, &
                              correct_leap_year_inconsistency, err_msg)
type(time_type), intent(in)  :: Time !< a specific time value
type(time_type), intent(in)  :: Time_beg !< begining of period to search with
type(time_type), intent(in)  :: Time_end !< end of period to search with
type(time_type), intent(in)  :: Timelist(:) !< ascending time values to search between
real(FMS_TI_KIND_)           , intent(out) :: weight
integer        , intent(out) :: index1, index2 !< indices of bounding time values within Timelist
logical, intent(in), optional :: correct_leap_year_inconsistency!< When true turns on a kluge for an
                                !! inconsistency which may occur in a special case.
                                !! When the modulo time period (i.e. Time_end - Time_beg) is a
                                !! whole number of years and is not a multiple of 4, and the calendar
                                !! in use has leap years, then it is likely that the interpolation
                                !! will involve mapping a common year onto a leap year. In this case
                                !! it is often desirable, but not absolutely necessary, to use data
                                !! for Feb 28 of the leap year when it is mapped onto a common year.
character(len=*), intent(out), optional :: err_msg
 
  type(time_type) :: Period, T
  integer :: is, ie,i1,i2
  integer :: ys,ms,ds,hs,mins,ss ! components of the starting date
  integer :: ye,me,de,he,mine,se ! components of the ending date
  integer :: yt,mt,dt,ht,mint,st ! components of the current date
  integer :: dt1                 ! temporary value for day
  integer :: n                   ! size of Timelist
  integer :: stdoutunit
  logical :: correct_lyr, calendar_has_leap_years, do_the_lyr_correction
  integer, parameter :: kindl = fms_ti_kind_
 
  if ( .not. module_is_initialized ) call time_interp_init
  if( present(err_msg) ) err_msg = ''
 
  stdoutunit = stdout()
  n = size(timelist)
 
  if (time_beg>=time_end) then
     if(fms_error_handler('time_interp_modulo', &
     'end of the specified time loop interval must be later than its beginning',err_msg)) return
  endif
 
  calendar_has_leap_years = (get_calendar_type() == julian .or. get_calendar_type() == gregorian)
 
  period = time_end-time_beg ! period of the time axis
 
  if(present(correct_leap_year_inconsistency)) then
    correct_lyr = correct_leap_year_inconsistency
  else
    correct_lyr = .false.
  endif
 
  ! bring the requested time inside the specified time period
  t = time
 
  do_the_lyr_correction = .false.
 
  ! Determine if the leap year correction needs to be done.
  ! It never needs to be done unless 3 conditions are met:
  ! 1) We are using a calendar with leap years
  ! 2) optional argument correct_leap_year_inconsistency is present and equals .true.
  ! 3) The modulo time period is an integer number of years
  ! If all of these are true then set do_the_lyr_correction to .true.
 
  if(calendar_has_leap_years .and. correct_lyr) then
    call get_date(time_beg,ys,ms,ds,hs,mins,ss)
    call get_date(time_end,ye,me,de,he,mine,se)
    if(ms==me.and.ds==de.and.hs==he.and.mins==mine.and.ss==se) then
      ! whole number of years
      do_the_lyr_correction = .true.
    endif
  endif
 
  if(do_the_lyr_correction) then
     call get_date(t,yt,mt,dt,ht,mint,st)
     yt = ys+modulo(yt-ys,ye-ys)
     dt1 = dt
     ! If it is Feb 29, but we map into a common year, use Feb 28
     if(mt==2.and.dt==29.and..not.leap_year(set_date(yt,1,1))) dt1=28
     t = set_date(yt,mt,dt1,ht,mint,st)
     if (t < time_beg) then
       ! the requested time is within the first year,
       ! but before the starting date. So we shift it to the last year.
       if(mt==2.and.dt==29.and..not.leap_year(set_date(ye,1,1))) dt=28
       t = set_date(ye,mt,dt,ht,mint,st)
     endif
  else
     do while ( t >= time_end )
        t = t-period
     enddo
     do while ( t < time_beg )
        t = t+period
     enddo
  endif
 
  ! find indices of the first and last records in the Timelist that are within
  ! the requested time period.
  if (time_end<=timelist(1).or.time_beg>=timelist(n)) then
     if(get_calendar_type() == no_calendar) then
       call print_time(time_beg,    'Time_beg'    )
       call print_time(time_end,    'Time_end'    )
       call print_time(timelist(1), 'Timelist(1)' )
       call print_time(timelist(n), 'Timelist(n)' )
     else
       call print_date(time_beg,    'Time_beg'    )
       call print_date(time_end,    'Time_end'    )
       call print_date(timelist(1), 'Timelist(1)' )
       call print_date(timelist(n), 'Timelist(n)' )
     endif
     write(stdoutunit,*)'where n = size(Timelist) =',n
     if(fms_error_handler('time_interp_modulo', &
     'the entire time list is outside the specified time loop interval',err_msg)) return
  endif
 
  call bisect(timelist,time_beg,index1=i1,index2=i2)
  if (i1 < 1) then
     is = 1 ! Time_beg before lower boundary
  else if (time_beg == timelist(i1)) then
     is = i1 ! Time_beg right on the lower boundary
  else
     is = i2 ! Time_beg inside the interval or on upper boundary
  endif
  call bisect(timelist,time_end,index1=i1,index2=i2)
  if (time_end > timelist(i1)) then
    ie = i1
  else if (time_end == timelist(i1)) then
    if(time_beg == timelist(is)) then
      ! Timelist includes time levels at both the lower and upper ends of the period.
      ! The endpoints of Timelist specify the same point in the cycle.
      ! This ambiguity is resolved by ignoring the last time level.
      ie = i1-1
    else
      ie = i1
    endif
  else
!   This should never happen because bisect does not return i1 such that Time_end < Timelist(i1)
  endif
  if (is>=ie) then
     if(get_calendar_type() == no_calendar) then
       call print_time(time_beg,    'Time_beg   =')
       call print_time(time_end,    'Time_end   =')
       call print_time(timelist(1), 'Timelist(1)=')
       call print_time(timelist(n), 'Timelist(n)=')
     else
       call print_date(time_beg,    'Time_beg   =')
       call print_date(time_end,    'Time_end   =')
       call print_date(timelist(1), 'Timelist(1)=')
       call print_date(timelist(n), 'Timelist(n)=')
     endif
     write(stdoutunit,*)'where n = size(Timelist) =',n
     write(stdoutunit,*)'is =',is,'ie =',ie
     if(fms_error_handler('time_interp_modulo', &
     'error in calculation of time list bounds within the specified time loop interval',err_msg)) return
  endif
 
  ! handle special cases:
  if( t>=timelist(ie) ) then
     ! time is after the end of the portion of the time list within the requested period
     index1 = ie;   index2 = is
     weight = real((t-timelist(ie))//(period-(timelist(ie)-timelist(is))), fms_ti_kind_ )
  else if (t<timelist(is)) then
     ! time is before the beginning of the portion of the time list within the requested period
     index1 = ie;   index2 = is
     weight = 1.0_kindl - real(((timelist(is)-t)//(period-(timelist(ie)-timelist(is)))), fms_ti_kind_ )
  else
     call bisect(timelist,t,index1,index2)
     weight = real((t-timelist(index1)) // (timelist(index2)-timelist(index1)), fms_ti_kind_ )
  endif
 
end subroutine time_interp_modulo_
 
subroutine time_interp_list_ ( Time, Timelist, weight, index1, index2, modtime, err_msg )
type(time_type)  , intent(in)  :: Time, Timelist(:)
real(FMS_TI_KIND_)             , intent(out) :: weight
integer          , intent(out) :: index1, index2
integer, optional, intent(in)  :: modtime
character(len=*), intent(out), optional :: err_msg
 
integer :: n, hr, mn, se, mtime
type(time_type) :: T, Ts, Te, Td, Period, Time_mod
character(len=:),allocatable :: terr, tserr, teerr
integer, parameter :: kindl = fms_ti_kind_
 
  if ( .not. module_is_initialized ) call time_interp_init
 
  if( present(err_msg) ) err_msg = ''
 
  weight = 0.0_kindl; index1 = 0; index2 = 0
  n = size(timelist(:))
 
! setup modular time axis?
  mtime = none
  if (present(modtime)) then
     mtime = modtime
     time_mod = (timelist(1)+timelist(n))/2
     call get_date (time_mod, yrmod, momod, dymod, hr, mn, se)
     mod_leapyear = leap_year(time_mod)
  endif
 
! set period for modulo axis
  select case (mtime)
     case (none)
       ! do nothing
     case (year)
         period = set_time(0,days_in_year(time_mod))
     case (month)
       ! month length must be equal
         if (days_in_month(time_mod) /= days_in_month(time)) then
            if(fms_error_handler('time_interp_list','modulo months must have same length',err_msg)) return
         endif
         period = set_time(0,days_in_month(time_mod))
     case (day)
         period = set_time(0,1)
     case default
         if(fms_error_handler('time_interp_list','invalid value for argument modtime',err_msg)) return
  end select
 
! If modulo time is in effect and Timelist spans a time interval exactly equal to
! the modulo period, then the endpoints of Timelist specify the same point in the cycle.
! This ambiguity is resolved by ignoring the last time level.
  if (mtime /= none .and. timelist(size(timelist))-timelist(1) == period) then
     n = size(timelist) - 1
  else
     n = size(timelist)
  endif
 
! starting and ending times from list
  ts = timelist(1)
  te = timelist(n)
  td = te-ts
  t  = set_modtime(time,mtime)
 
! Check that Timelist does not span a time interval greater than the modulo period
  if (mtime /= none) then
     if (td > period) then
        if(fms_error_handler('time_interp_list','period of list exceeds modulo period',err_msg)) return
     endif
  endif
 
! time falls on start or between start and end list values
  if ( t >= ts .and. t < te ) then
     call bisect(timelist(1:n),t,index1,index2)
     weight = real( (t-timelist(index1)) // (timelist(index2)-timelist(index1)), fms_ti_kind_)
 
! time falls before starting list value
  else if ( t < ts ) then
     if (mtime == none) then
        call time_list_error(t,terr)
        call time_list_error(ts,tserr)
        call time_list_error(te,teerr)
        if(fms_error_handler('time_interp_list',&
           'time '//trim(terr)//' ('//date_to_string(t)//' is before range of list '//trim(tserr)//'-'//trim(teerr)//&
           '('//date_to_string(ts)//' - '//date_to_string(te)//')',&
           err_msg)) return
        deallocate(terr,tserr,teerr)
     endif
     td = te-ts
     weight = 1.0_kindl - real(((ts-t) // (period-td)), fms_ti_kind_ )
     index1 = n
     index2 = 1
 
! time falls on ending list value
  else if ( t == te ) then
    if(perthlike_behavior) then
       weight = 1.0_kindl
       index1 = n-1
       index2 = n
    else
       weight = 0.0_kindl
       index1 = n
       if (mtime == none) then
         index2 = n
       else
         index2 = 1
       endif
    endif
 
! time falls after ending list value
  else if ( t > te ) then
     if (mtime == none) then
        call time_list_error(t,terr)
        call time_list_error(ts,tserr)
        call time_list_error(te,teerr)
        if(fms_error_handler('time_interp_list',&
           'time '//trim(terr)//' ('//date_to_string(t)//' is after range of list '//trim(tserr)//'-'//trim(teerr)//&
           '('//date_to_string(ts)//' - '//date_to_string(te)//')',&
           err_msg)) return
        deallocate(terr,tserr,teerr)
     endif
     td = te-ts
     weight = real( (t-te) // (period-td), fms_ti_kind_)
     index1 = n
     index2 = 1
  endif
 
end subroutine time_interp_list_
!> }