time_interp/time_interp.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup time_interp_mod time_interp_mod
  20 !> @ingroup time_interp
  21 !> @brief Computes a weight and dates/indices for linearly interpolating between two dates.
  22 !> @author Bruce Wyman
  23 !!
  24 !! A time type is converted into two consecutive dates plus
  25 !! a fraction representing the distance between the dates.
  26 !! This information can be used to interpolate between the dates.
  27 !! The dates may be expressed as years, months, or days or
  28 !! as indices in an array.
  29
  30 module time_interp_mod
  31
  32 use time_manager_mod, only: time_type, get_date, set_date, set_time, &
  33                             days_in_year, days_in_month, leap_year,  &
  34                             time_type_to_real, real_to_time_type,    &
  35                             get_calendar_type, JULIAN, GREGORIAN, NO_CALENDAR, &
  36                             operator(+), operator(-), operator(>),   &
  37                             operator(<), operator( // ), operator( / ),  &
  38                             operator(>=), operator(<=), operator( * ), &
  39                             operator(==), print_date, print_time,&
  40                             time_list_error, date_to_string
  41
  42 use          fms_mod, only: write_version_number, &
  43                             error_mesg, FATAL, stdout, stdlog, &
  44                             check_nml_error, &
  45                             fms_error_handler
  46 use          mpp_mod, only: input_nml_file
  47 use     platform_mod
  48
  49 implicit none
  50 private
  51
  52 !-----------------------------------------------------------------------
  53
  54 public :: time_interp_init, time_interp, fraction_of_year
  55
  56 !> Returns a weight and dates or indices for interpolating between two dates. The
  57 !! interface fraction_of_year is provided for backward compatibility with the
  58 !! previous version.\n
  59 !!
  60 !! Returns weight by interpolating Time between Time1 and Time2.
  61 !! i.e. weight = (Time-Time1)/(Time2-Time1)
  62 !! Time1 and Time2 may be specified by any of several different ways,
  63 !! which is the reason for multiple interfaces.\n
  64 !!
  65 !! - If Time1 and Time2 are the begining and end of the year in which
  66 !!   Time falls, use first interface.\n
  67 !!
  68 !! - If Time1 and Time2 fall on year boundaries, use second interface.\n
  69 !!
  70 !! - If Time1 and Time2 fall on month boundaries, use third.\n
  71 !!
  72 !! - If Time1 and Time2 fall on day boundaries, use fourth.\n
  73 !!
  74 !! - If Time1 and Time2 are consecutive elements of an assending list, use fifth.
  75 !!   The fifth also returns the indices of Timelist between which Time falls.\n
  76 !!
  77 !! - The sixth interface is for cyclical data. Time_beg and Time_end specify the
  78 !!   begining and end of a repeating period. In this case:<br>
  79 !!              weight = (Time_adjusted - Time1) / (Time2 - Time1)
  80 !! <br>Where:
  81 !! @code{.F90}
  82 !!              Time1 = Timelist(index1)
  83 !!              Time2 = Timelist(index2)
  84 !!              Time_adjusted = Time - N*Period
  85 !!              Period = Time_end-Time_beg
  86 !! @endcode
  87 !! N is between (Time-Time_end)/Period and (Time-Time_beg)/Period
  88 !! That is, N is the integer that results in Time_adjusted that is between Time_beg and Time_end.
  89 !!
  90 !! <br>Example usages:
  91 !! @code{.F90}
  92 !!              call time_interp( Time, weight )
  93 !!              call time_interp( Time, weight, year1, year2 )
  94 !!              call time_interp( Time, weight, year1, year2, month1, month2 )
  95 !!              call time_interp( Time, weight, year1, year2, month1, month2, day1, day2 )
  96 !!              call time_interp( Time, Timelist, weight, index1, index2 [, modtime] )
  97 !!              call time_interp( Time, Time_beg, Time_end, Timelist, weight, index1, index2
  98 !!              [,correct_leap_year_inconsistency])
  99 !! @endcode
 100 !!
 101 !!   For all routines in this module the calendar type in module
 102 !!   time_manager must be set.
 103 !!
 104 !!     The following private parameters are set by this module:
 105 !!
 106 !!          seconds per minute = 60
 107 !!          minutes per hour   = 60
 108 !!          hours   per day    = 24
 109 !!          months  per year   = 12
 110 !!
 111 !! @param Time The time at which the the weight is computed.
 112 !! @param Time_beg For cyclical interpolation: Time_beg specifies the begining time of a cycle.
 113 !! @param Time_end For cyclical interpolation: Time_end specifies the ending time of a cycle.
 114 !! @param Timelist For cyclical interpolation: Timelist is an array of times between Time_beg and Time_end.
 115 !!                 Must be monotonically increasing.
 116 !! @param index1 Timelist(index1) = The largest value of Timelist which is less than mod(Time,Time_end-Time_beg)
 117 !! @param index2 Timelist(index2) = The smallest value of Timelist which is greater than mod(Time,Time_end-Time_beg)
 118 !! @param correct_leap_year_inconsistency Turns on a kluge for an inconsistency which may occur in a special case.
 119 !!       When the modulo time period (i.e. Time_end - Time_beg) is a whole number of years
 120 !!       and is not a multiple of 4, and the calendar in use has leap years, then it is
 121 !!       likely that the interpolation will involve mapping a common year onto a leap year.
 122 !!       In this case it is often desirable, but not absolutely necessary, to use data for
 123 !!       Feb 28 of the leap year when it is mapped onto a common year.
 124 !!       To turn this on, set correct_leap_year_inconsistency=.true.
 125 !! @param weight weight = (mod(Time,Time_end-Time_beg) - Timelist(index1)) / (Timelist(index2) - Timelist(index1))
 126 !> @ingroup time_interp_mod
 127 interface time_interp
 128     module procedure time_interp_frac_r8,  time_interp_year_r8, &
 129                      time_interp_month_r8, time_interp_day_r8,  &
 130                      time_interp_list_r8,  time_interp_modulo_r8
 131     module procedure time_interp_frac_r4,  time_interp_year_r4, &
 132                      time_interp_month_r4, time_interp_day_r4,  &
 133                      time_interp_list_r4,  time_interp_modulo_r4
 134 end interface
 135
 136 !> @addtogroup time_interp_mod
 137 !> @{
 138 integer, public, parameter :: NONE=0, YEAR=1, MONTH=2, DAY=3
 139
 140 !-----------------------------------------------------------------------
 141
 142    integer, parameter ::  secmin = 60, minhour = 60, hourday = 24,  &
 143                           sechour = secmin*minhour,                  &
 144                           secday = secmin*minhour*hourday
 145
 146    integer, parameter :: monyear = 12
 147    integer, parameter :: halfday = secday/2
 148
 149    integer :: yrmod, momod, dymod
 150    logical :: mod_leapyear
 151
 152 ! Include variable "version" to be written to log file.
 153 #include<file_version.h>
 154
 155    logical :: module_is_initialized=.FALSE.
 156    logical :: perthlike_behavior=.FALSE.
 157
 158    namelist / time_interp_nml / perthlike_behavior
 159
 160 contains
 161
 162
 163  subroutine time_interp_init()
 164    integer :: ierr, io, logunit
 165
 166    if ( module_is_initialized ) return
 167
 168    read (input_nml_file, time_interp_nml, iostat=io)
 169    ierr = check_nml_error (io, 'time_interp_nml')
 170
 171    call write_version_number("TIME_INTERP_MOD", version)
 172    logunit = stdlog()
 173    write(logunit,time_interp_nml)
 174
 175    module_is_initialized = .TRUE.
 176
 177  end subroutine time_interp_init
 178
 179
 180 !> @brief Wrapper function to return the fractional time into the current year
 181 !! Always returns an r8_kind, conversion to r4 will be done implicitly if needed
 182 !> @param Time time to calculate fraction with
 183 !> @return real(kind=8) fraction of time passed in current year
 184  function fraction_of_year (Time)
 185  type(time_type), intent(in)  :: Time
 186  real(r8_kind) :: fraction_of_year
 187
 188   call time_interp ( Time, fraction_of_year )
 189
 190  end function fraction_of_year
 191
 192 !#######################################################################
 193 !> Given an array of times in ascending order and a specific time returns
 194 !! values of index1 and index2 such that the Timelist(index1)<=Time and
 195 !! Time<=Timelist(index2), and index2=index1+1
 196 !! index1=0, index2=1 or index=n, index2=n+1 are returned to indicate that
 197 !! the time is out of range
 198 subroutine bisect(Timelist,Time,index1,index2)
 199   type(time_type)  , intent(in)  :: Timelist(:)
 200   type(time_type)  , intent(in)  :: Time
 201   integer, optional, intent(out) :: index1, index2
 202
 203   integer :: i,il,iu,n,i1,i2
 204
 205   n = size(Timelist(:))
 206
 207   if (Time==Timelist(1)) then
 208      i1 = 1 ; i2 = 2
 209   else if (Time==Timelist(n)) then
 210      i1 = n ; i2 = n+1
 211   else
 212      il = 0; iu=n+1
 213      do while(iu-il > 1)
 214         i = (iu+il)/2
 215         if(Timelist(i) > Time) then
 216            iu = i
 217         else
 218            il = i
 219         endif
 220      enddo
 221      i1 = il ; i2 = il+1
 222   endif
 223
 224   if(PRESENT(index1)) index1 = i1
 225   if(PRESENT(index2)) index2 = i2
 226 end subroutine bisect
 227
 228 !#######################################################################
 229 !  private routines
 230 !#######################################################################
 231
 232  function year_midpt (yr)
 233
 234    integer, intent(in) :: yr
 235    type (time_type)    :: year_midpt, year_beg, year_end
 236
 237
 238    year_beg = set_date(yr  , 1, 1)
 239    year_end = set_date(yr+1, 1, 1)
 240
 241    year_midpt = (year_beg + year_end) / 2
 242
 243  end function year_midpt
 244
 245  function month_midpt (yr, mo)
 246
 247    integer, intent(in) :: yr, mo
 248    type (time_type)    :: month_midpt, month_beg, month_end
 249
 250 !  --- beginning of this month ---
 251    month_beg = set_date(yr, mo, 1)
 252
 253 !  --- start of next month ---
 254    if (mo < 12) then
 255       month_end = set_date(yr, mo+1, 1)
 256    else
 257       month_end = set_date(yr+1, 1, 1)
 258    endif
 259
 260    month_midpt = (month_beg + month_end) / 2
 261
 262  end function month_midpt
 263
 264 function set_modtime (Tin, modtime) result (Tout)
 265 type(time_type), intent(in) :: Tin
 266 integer, intent(in), optional :: modtime
 267 type(time_type)             :: Tout
 268 integer :: yr, mo, dy, hr, mn, se, mtime
 269
 270   if(present(modtime)) then
 271     mtime = modtime
 272   else
 273     mtime = NONE
 274   endif
 275
 276   select case (mtime)
 277     case (NONE)
 278        Tout = Tin
 279     case (YEAR)
 280        call get_date (Tin, yr, mo, dy, hr, mn, se)
 281        yr = yrmod
 282         ! correct leap year dates
 283           if (.not.mod_leapyear .and. mo == 2 .and. dy > 28) then
 284              mo = 3; dy = dy-28
 285           endif
 286        Tout = set_date (yr, mo, dy, hr, mn, se)
 287     case (MONTH)
 288        call get_date (Tin, yr, mo, dy, hr, mn, se)
 289        yr = yrmod; mo = momod
 290        Tout = set_date (yr, mo, dy, hr, mn, se)
 291     case (DAY)
 292        call get_date (Tin, yr, mo, dy, hr, mn, se)
 293        yr = yrmod; mo = momod; dy = dymod
 294        Tout = set_date (yr, mo, dy, hr, mn, se)
 295   end select
 296
 297 end function set_modtime
 298
 299 subroutine error_handler(string)
 300   character(len=*), intent(in) :: string
 301
 302   call error_mesg ('time_interp_mod', trim(string), FATAL)
 303
 304 end subroutine error_handler
 305
 306
 307 #include "time_interp_r4.fh"
 308 #include "time_interp_r8.fh"
 309
 310 end module time_interp_mod
 311
 312 !> @}
 313 ! close documentation grouping