TimeSeries.f90

Go to the documentation of this file.

module timeseriesmodule
 
  use kindmodule, only: dp, i4b
  use blockparsermodule, only: blockparsertype
  use constantsmodule, only: linelength, undefined, stepwise, linear, &
                             linearend, lentimeseriesname, lenhugeline, &
                             dzero, done, dnodata
  use mathutilmodule, only: is_close
  use inputoutputmodule, only: getunit, openfile, parseline, upcase
  use listmodule, only: listtype
  use listnodemodule, only: listnodetype
  use simvariablesmodule, only: errmsg
  use simmodule, only: count_errors, store_error, &
                       store_error_unit
  use timeseriesrecordmodule, only: timeseriesrecordtype, &
                                    constructtimeseriesrecord, &
                                    castastimeseriesrecordtype, &
                                    addtimeseriesrecordtolist
 
  private
  public :: timeseriestype, timeseriesfiletype, constructtimeseriesfile, &
            timeseriescontainertype, addtimeseriesfiletolist, &
            gettimeseriesfilefromlist, castastimeseriesfileclass, &
            sametimeseries
 
  type timeseriestype
    ! -- Public members
    integer(I4B), public :: imethod = undefined
    character(len=LENTIMESERIESNAME), public :: name = ''
    ! -- Private members
    real(dp), private :: sfac = done
    logical, public :: autodeallocate = .true.
    type(listtype), pointer, private :: list => null()
    class(timeseriesfiletype), pointer, private :: tsfile => null()
 
  contains
 
    ! -- Public procedures
    procedure, public :: addtimeseriesrecord
    procedure, public :: clear
    procedure, public :: findlatesttime
    procedure, public :: get_surrounding_records
    procedure, public :: get_surrounding_nodes
    procedure, public :: getcurrenttimeseriesrecord
    procedure, public :: getnexttimeseriesrecord
    procedure, public :: getprevioustimeseriesrecord
    procedure, public :: gettimeseriesrecord
    procedure, public :: getvalue
    procedure, public :: initializetimeseries => initialize_time_series
    procedure, public :: inserttsr
    procedure, public :: reset
    ! -- Private procedures
    procedure, private :: da => ts_da
    procedure, private :: get_average_value
    procedure, private :: get_integrated_value
    procedure, private :: get_latest_preceding_node
    procedure, private :: get_value_at_time
    procedure, private :: initialize_time_series
    procedure, private :: read_next_record
  end type timeseriestype
 
  type timeseriesfiletype
    ! -- Private members
    integer(I4B), public :: inunit = 0
    integer(I4B), public :: iout = 0
    integer(I4B), public :: ntimeseries = 0
    logical, public :: finishedreading = .false.
    character(len=LINELENGTH), public :: datafile = ''
    type(timeseriestype), dimension(:), &
      pointer, contiguous, public :: timeseries => null()
    type(blockparsertype), pointer, public :: parser
 
  contains
 
    ! -- Public procedures
    procedure, public :: count
    procedure, public :: initializetsfile
    procedure, public :: gettimeseries
    procedure, public :: da => tsf_da
    ! -- Private procedures
    procedure, private :: read_tsfile_line
  end type timeseriesfiletype
 
  type timeseriescontainertype
    ! -- Public members
    type(timeseriestype), pointer, public :: timeseries => null()
  end type timeseriescontainertype
 
contains
 
  ! -- non-type-bound procedures
 
  !> @brief Construct time series file
  !<
  subroutine constructtimeseriesfile(newTimeSeriesFile)
    ! -- dummy
    type(timeseriesfiletype), pointer, intent(inout) :: newtimeseriesfile
    !
    allocate (newtimeseriesfile)
    allocate (newtimeseriesfile%parser)
  end subroutine constructtimeseriesfile
 
  !> @brief Cast an unlimited polymorphic object as class(TimeSeriesFileType)
  !<
  function castastimeseriesfiletype(obj) result(res)
    ! -- dummy
    class(*), pointer, intent(inout) :: obj
    ! -- return
    type(timeseriesfiletype), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    type is (timeseriesfiletype)
      res => obj
    end select
  end function castastimeseriesfiletype
 
  !> @brief Cast an unlimited polymorphic object as class(TimeSeriesFileType)
  !<
  function castastimeseriesfileclass(obj) result(res)
    ! -- dummy
    class(*), pointer, intent(inout) :: obj
    ! -- return
    type(timeseriesfiletype), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    class is (timeseriesfiletype)
      res => obj
    end select
  end function castastimeseriesfileclass
 
  !> @brief Add time series file to list
  !<
  subroutine addtimeseriesfiletolist(list, tsfile)
    ! -- dummy
    type(listtype), intent(inout) :: list
    class(timeseriesfiletype), pointer, intent(inout) :: tsfile
    ! -- local
    class(*), pointer :: obj => null()
    !
    obj => tsfile
    call list%Add(obj)
  end subroutine addtimeseriesfiletolist
 
  !> @brief Get time series from list
  !<
  function gettimeseriesfilefromlist(list, idx) result(res)
    ! -- dummy
    type(listtype), intent(inout) :: list
    integer(I4B), intent(in) :: idx
    ! -- return
    type(timeseriesfiletype), pointer :: res
    ! -- local
    class(*), pointer :: obj => null()
    !
    obj => list%GetItem(idx)
    res => castastimeseriesfiletype(obj)
    !
    if (.not. associated(res)) then
      res => castastimeseriesfileclass(obj)
    end if
  end function gettimeseriesfilefromlist
 
  !> @brief Compare two time series; if they are identical, return true
  !<
  function sametimeseries(ts1, ts2) result(same)
    ! -- dummy
    type(timeseriestype), intent(in) :: ts1
    type(timeseriestype), intent(in) :: ts2
    ! -- return
    logical :: same
    ! -- local
    integer :: i, n1, n2
    type(timeseriesrecordtype), pointer :: tsr1, tsr2
    !
    same = .false.
    n1 = ts1%list%Count()
    n2 = ts2%list%Count()
    if (n1 /= n2) return
    !
    call ts1%Reset()
    call ts2%Reset()
    !
    do i = 1, n1
      tsr1 => ts1%GetNextTimeSeriesRecord()
      tsr2 => ts2%GetNextTimeSeriesRecord()
      if (tsr1%tsrTime /= tsr2%tsrTime) return
      if (tsr1%tsrValue /= tsr2%tsrValue) return
    end do
    !
    same = .true.
  end function sametimeseries
 
  ! Type-bound procedures of TimeSeriesType
 
  !> @brief Get time series value
  !!
  !! If iMethod is STEPWISE or LINEAR:
  !!     Return a time-weighted average value for a specified time span.
  !! If iMethod is LINEAREND:
  !!     Return value at time1. Time0 argument is ignored.
  !! Units: (ts-value-unit)
  !<
  function getvalue(this, time0, time1, extendToEndOfSimulation)
    ! -- return
    real(dp) :: getvalue
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(dp), intent(in) :: time0
    real(dp), intent(in) :: time1
    logical, intent(in), optional :: extendtoendofsimulation
    ! -- local
    logical :: extend
    !
    if (present(extendtoendofsimulation)) then
      extend = extendtoendofsimulation
    else
      extend = .false.
    end if
    !
    select case (this%iMethod)
    case (stepwise, linear)
      getvalue = this%get_average_value(time0, time1, extend)
    case (linearend)
      getvalue = this%get_value_at_time(time1, extend)
    end select
  end function getvalue
 
  !> @brief Initialize time series
  !!
  !! Open time-series file and read options and first time-series record.
  !<
  subroutine initialize_time_series(this, tsfile, name, autoDeallocate)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    class(timeseriesfiletype), target :: tsfile
    character(len=*), intent(in) :: name
    logical, intent(in), optional :: autoDeallocate
    ! -- local
    character(len=LENTIMESERIESNAME) :: tsNameTemp
    !
    ! -- Assign the time-series tsfile, name, and autoDeallocate
    this%tsfile => tsfile
    ! Store time-series name as all caps
    tsnametemp = name
    call upcase(tsnametemp)
    this%Name = tsnametemp
    !
    this%iMethod = undefined
    !
    if (present(autodeallocate)) this%autoDeallocate = autodeallocate
    !
    ! -- allocate the list
    allocate (this%list)
    !
    ! -- ensure that NAME has been specified
    if (this%Name == '') then
      errmsg = 'Name not specified for time series.'
      call store_error(errmsg, terminate=.true.)
    end if
  end subroutine initialize_time_series
 
  !> @brief Get surrounding records
  !<
  subroutine get_surrounding_records(this, time, tsrecEarlier, tsrecLater)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(DP), intent(in) :: time
    type(timeseriesrecordtype), pointer, intent(inout) :: tsrecEarlier
    type(timeseriesrecordtype), pointer, intent(inout) :: tsrecLater
    ! -- local
    real(DP) :: time0, time1
    type(listnodetype), pointer :: currNode => null()
    type(listnodetype), pointer :: tsNode0 => null()
    type(listnodetype), pointer :: tsNode1 => null()
    type(timeseriesrecordtype), pointer :: tsr => null(), tsrec0 => null()
    type(timeseriesrecordtype), pointer :: tsrec1 => null()
    class(*), pointer :: obj => null()
    !
    tsrecearlier => null()
    tsreclater => null()
    !
    if (associated(this%list%firstNode)) then
      currnode => this%list%firstNode
    end if
    !
    ! -- If the next node is earlier than time of interest, advance along
    !    linked list until the next node is later than time of interest.
    do
      if (associated(currnode)) then
        if (associated(currnode%nextNode)) then
          obj => currnode%nextNode%GetItem()
          tsr => castastimeseriesrecordtype(obj)
          if (tsr%tsrTime < time .and. .not. is_close(tsr%tsrTime, time)) then
            currnode => currnode%nextNode
          else
            exit
          end if
        else
          ! -- read another record
          if (.not. this%read_next_record()) exit
        end if
      else
        exit
      end if
    end do
    !
    if (associated(currnode)) then
      !
      ! -- find earlier record
      tsnode0 => currnode
      obj => tsnode0%GetItem()
      tsrec0 => castastimeseriesrecordtype(obj)
      time0 = tsrec0%tsrTime
      do while (time0 > time)
        if (associated(tsnode0%prevNode)) then
          tsnode0 => tsnode0%prevNode
          obj => tsnode0%GetItem()
          tsrec0 => castastimeseriesrecordtype(obj)
          time0 = tsrec0%tsrTime
        else
          exit
        end if
      end do
      !
      ! -- find later record
      tsnode1 => currnode
      obj => tsnode1%GetItem()
      tsrec1 => castastimeseriesrecordtype(obj)
      time1 = tsrec1%tsrTime
      do while (time1 < time .and. .not. is_close(time1, time))
        if (associated(tsnode1%nextNode)) then
          tsnode1 => tsnode1%nextNode
          obj => tsnode1%GetItem()
          tsrec1 => castastimeseriesrecordtype(obj)
          time1 = tsrec1%tsrTime
        else
          ! -- get next record
          if (.not. this%read_next_record()) then
            ! -- end of file reached, so exit loop
            exit
          end if
        end if
      end do
      !
    end if
    !
    if (time0 < time .or. is_close(time0, time)) tsrecearlier => tsrec0
    if (time1 > time .or. is_close(time1, time)) tsreclater => tsrec1
  end subroutine get_surrounding_records
 
  !> @brief Get surrounding nodes
  !!
  !! This subroutine is for working with time series already entirely stored
  !! in memory -- it does not read data from a file.
  !<
  subroutine get_surrounding_nodes(this, time, nodeEarlier, nodeLater)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(DP), intent(in) :: time
    type(listnodetype), pointer, intent(inout) :: nodeEarlier
    type(listnodetype), pointer, intent(inout) :: nodeLater
    ! -- local
    real(DP) :: time0, time1
    type(listnodetype), pointer :: currNode => null()
    type(listnodetype), pointer :: tsNode0 => null()
    type(listnodetype), pointer :: tsNode1 => null()
    type(timeseriesrecordtype), pointer :: tsr => null(), tsrec0 => null()
    type(timeseriesrecordtype), pointer :: tsrec1 => null()
    type(timeseriesrecordtype), pointer :: tsrecEarlier
    type(timeseriesrecordtype), pointer :: tsrecLater
    class(*), pointer :: obj => null()
    !
    tsrecearlier => null()
    tsreclater => null()
    nodeearlier => null()
    nodelater => null()
    !
    if (associated(this%list%firstNode)) then
      currnode => this%list%firstNode
    end if
    !
    ! -- If the next node is earlier than time of interest, advance along
    !    linked list until the next node is later than time of interest.
    do
      if (associated(currnode)) then
        if (associated(currnode%nextNode)) then
          obj => currnode%nextNode%GetItem()
          tsr => castastimeseriesrecordtype(obj)
          if (tsr%tsrTime < time .and. .not. is_close(tsr%tsrTime, time)) then
            currnode => currnode%nextNode
          else
            exit
          end if
        else
          exit
        end if
      else
        exit
      end if
    end do
    !
    if (associated(currnode)) then
      !
      ! -- find earlier record
      tsnode0 => currnode
      obj => tsnode0%GetItem()
      tsrec0 => castastimeseriesrecordtype(obj)
      time0 = tsrec0%tsrTime
      do while (time0 > time)
        if (associated(tsnode0%prevNode)) then
          tsnode0 => tsnode0%prevNode
          obj => tsnode0%GetItem()
          tsrec0 => castastimeseriesrecordtype(obj)
          time0 = tsrec0%tsrTime
        else
          exit
        end if
      end do
      !
      ! -- find later record
      tsnode1 => currnode
      obj => tsnode1%GetItem()
      tsrec1 => castastimeseriesrecordtype(obj)
      time1 = tsrec1%tsrTime
      do while (time1 < time .and. .not. is_close(time1, time))
        if (associated(tsnode1%nextNode)) then
          tsnode1 => tsnode1%nextNode
          obj => tsnode1%GetItem()
          tsrec1 => castastimeseriesrecordtype(obj)
          time1 = tsrec1%tsrTime
        else
          exit
        end if
      end do
      !
    end if
    !
    if (time0 < time .or. is_close(time0, time)) then
      tsrecearlier => tsrec0
      nodeearlier => tsnode0
    end if
    if (time1 > time .or. is_close(time1, time)) then
      tsreclater => tsrec1
      nodelater => tsnode1
    end if
  end subroutine get_surrounding_nodes
 
  !> @brief Read next record
  !!
  !! Read next time-series record from input file
  !<
  logical function read_next_record(this)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    !
    ! -- If we have already encountered the end of the TIMESERIES block, do not try to read any further
    if (this%tsfile%finishedReading) then
      read_next_record = .false.
      return
    end if
    !
    read_next_record = this%tsfile%read_tsfile_line()
    if (.not. read_next_record) then
      this%tsfile%finishedReading = .true.
    end if
  end function read_next_record
 
  !> @brief Get value for a time
  !!
  !! Return a value for a specified time, same units as time-series values
  !<
  function get_value_at_time(this, time, extendToEndOfSimulation)
    ! -- return
    real(dp) :: get_value_at_time
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(dp), intent(in) :: time ! time of interest
    logical, intent(in) :: extendtoendofsimulation
    ! -- local
    integer(I4B) :: ierr
    real(dp) :: ratio, time0, time1, timediff, timediffi, val0, val1, &
                valdiff
    type(timeseriesrecordtype), pointer :: tsrearlier => null()
    type(timeseriesrecordtype), pointer :: tsrlater => null()
    ! -- formats
10  format('Error getting value at time ', g10.3, ' for time series "', a, '"')
    !
    ierr = 0
    call this%get_surrounding_records(time, tsrearlier, tsrlater)
    if (associated(tsrearlier)) then
      if (associated(tsrlater)) then
        ! -- values are available for both earlier and later times
        if (this%iMethod == stepwise) then
          get_value_at_time = tsrearlier%tsrValue
        elseif (this%iMethod == linear .or. this%iMethod == linearend) then
          ! -- For get_value_at_time, result is the same for either
          !    linear method.
          ! -- Perform linear interpolation.
          time0 = tsrearlier%tsrTime
          time1 = tsrlater%tsrtime
          timediff = time1 - time0
          timediffi = time - time0
          if (timediff > 0) then
            ratio = timediffi / timediff
          else
            ! -- should not happen if TS does not contain duplicate times
            ratio = 0.5d0
          end if
          val0 = tsrearlier%tsrValue
          val1 = tsrlater%tsrValue
          valdiff = val1 - val0
          get_value_at_time = val0 + (ratio * valdiff)
        else
          ierr = 1
        end if
      else
        if (extendtoendofsimulation .or. is_close(tsrearlier%tsrTime, time)) then
          get_value_at_time = tsrearlier%tsrValue
        else
          ! -- Only earlier time is available, and it is not time of interest;
          !    however, if method is STEPWISE, use value for earlier time.
          if (this%iMethod == stepwise) then
            get_value_at_time = tsrearlier%tsrValue
          else
            ierr = 1
          end if
        end if
      end if
    else
      if (associated(tsrlater)) then
        if (is_close(tsrlater%tsrTime, time)) then
          get_value_at_time = tsrlater%tsrValue
        else
          ! -- only later time is available, and it is not time of interest
          ierr = 1
        end if
      else
        ! -- Neither earlier nor later time is available.
        !    This should never happen!
        ierr = 1
      end if
    end if
    !
    if (ierr > 0) then
      write (errmsg, 10) time, trim(this%Name)
      call store_error(errmsg, terminate=.true.)
    end if
  end function get_value_at_time
 
  !> @brief Get integrated value
  !!
  !! Return an integrated value for a specified time span.
  !! Units: (ts-value-unit)*time
  !<
  function get_integrated_value(this, time0, time1, extendToEndOfSimulation)
    ! -- return
    real(dp) :: get_integrated_value
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(dp), intent(in) :: time0
    real(dp), intent(in) :: time1
    logical, intent(in) :: extendtoendofsimulation
    ! -- local
    real(dp) :: area, currtime, nexttime, ratio0, ratio1, t0, t01, t1, &
                timediff, value, value0, value1, valuediff, currval, nextval
    logical :: ldone, lprocess
    type(listnodetype), pointer :: tslnodepreceding => null()
    type(listnodetype), pointer :: currnode => null(), nextnode => null()
    type(timeseriesrecordtype), pointer :: currrecord => null()
    type(timeseriesrecordtype), pointer :: nextrecord => null()
    class(*), pointer :: currobj => null(), nextobj => null()
    ! -- formats
10  format('Error encountered while performing integration', &
           ' for time series "', a, '" for time interval: ', g12.5, ' to ', g12.5)
    !
    value = dzero
    ldone = .false.
    t1 = -done
    call this%get_latest_preceding_node(time0, tslnodepreceding)
    if (associated(tslnodepreceding)) then
      currnode => tslnodepreceding
      do while (.not. ldone)
        currobj => currnode%GetItem()
        currrecord => castastimeseriesrecordtype(currobj)
        currtime = currrecord%tsrTime
        if (is_close(currtime, time1)) then
          ! Current node time = time1 so should be ldone
          ldone = .true.
        elseif (currtime < time1) then
          if (.not. associated(currnode%nextNode)) then
            ! -- try to read the next record
            if (.not. this%read_next_record()) then
              if (.not. extendtoendofsimulation) then
                write (errmsg, 10) trim(this%Name), time0, time1
                call store_error(errmsg, terminate=.true.)
              end if
            end if
          end if
          !
          currval = currrecord%tsrValue
          lprocess = .false.
          if (associated(currnode%nextNode)) then
            nextnode => currnode%nextNode
            nextobj => nextnode%GetItem()
            nextrecord => castastimeseriesrecordtype(nextobj)
            nexttime = nextrecord%tsrTime
            nextval = nextrecord%tsrValue
            lprocess = .true.
          elseif (extendtoendofsimulation) then
            ! -- Last time series value extends forever, so integrate the final value over all simulation time after the end of the series
            nexttime = time1
            nextval = currval
            lprocess = .true.
          end if
          !
          if (lprocess) then
            ! -- determine lower and upper limits of time span of interest
            !    within current interval
            if (currtime > time0 .or. is_close(currtime, time0)) then
              t0 = currtime
            else
              t0 = time0
            end if
            if (nexttime < time1 .or. is_close(nexttime, time1)) then
              t1 = nexttime
            else
              t1 = time1
            end if
            ! -- find area of rectangle or trapezoid delimited by t0 and t1
            t01 = t1 - t0
            select case (this%iMethod)
            case (stepwise)
              ! -- compute area of a rectangle
              value0 = currval
              area = value0 * t01
            case (linear, linearend)
              ! -- compute area of a trapezoid
              timediff = nexttime - currtime
              ratio0 = (t0 - currtime) / timediff
              ratio1 = (t1 - currtime) / timediff
              valuediff = nextval - currval
              value0 = currval + ratio0 * valuediff
              value1 = currval + ratio1 * valuediff
              if (this%iMethod == linear) then
                area = 0.5d0 * t01 * (value0 + value1)
              elseif (this%iMethod == linearend) then
                area = dzero
                value = value1
              end if
            end select
            ! -- add area to integrated value
            value = value + area
          end if
        end if
        !
        ! -- Are we done yet?
        if (t1 > time1) then
          ldone = .true.
        elseif (is_close(t1, time1)) then
          ldone = .true.
        else
          ! -- We are not done yet
          if (.not. associated(currnode%nextNode)) then
            ! -- Not done and no more data, so try to read the next record
            if (.not. this%read_next_record()) then
              write (errmsg, 10) trim(this%Name), time0, time1
              call store_error(errmsg, terminate=.true.)
            end if
          elseif (associated(currnode%nextNode)) then
            currnode => currnode%nextNode
          end if
        end if
      end do
    end if
    !
    get_integrated_value = value
    if (this%autoDeallocate) then
      if (associated(tslnodepreceding)) then
        if (associated(tslnodepreceding%prevNode)) then
          call this%list%DeallocateBackward(tslnodepreceding%prevNode)
        end if
      end if
    end if
  end function get_integrated_value
 
  !> @brief Get average value
  !!
  !! Return a time-weighted average value for a specified time span.
  !! Units: (ts-value-unit)
  !<
  function get_average_value(this, time0, time1, extendToEndOfSimulation)
    ! -- return
    real(dp) :: get_average_value
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(dp), intent(in) :: time0
    real(dp), intent(in) :: time1
    logical, intent(in) :: extendtoendofsimulation
    ! -- local
    real(dp) :: timediff, value, valueintegrated
    !
    timediff = time1 - time0
    if (timediff > 0) then
      valueintegrated = this%get_integrated_value(time0, time1, &
                                                  extendtoendofsimulation)
      if (this%iMethod == linearend) then
        value = valueintegrated
      else
        value = valueintegrated / timediff
      end if
    else
      ! -- time0 and time1 are the same
      value = this%get_value_at_time(time0, extendtoendofsimulation)
    end if
    get_average_value = value
  end function get_average_value
 
  !> @brief Get latest preceding node
  !!
  !! Return pointer to ListNodeType object for the node representing the
  !! latest preceding time in the time series
  !<
  subroutine get_latest_preceding_node(this, time, tslNode)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    real(DP), intent(in) :: time
    type(listnodetype), pointer, intent(inout) :: tslNode
    ! -- local
    real(DP) :: time0
    type(listnodetype), pointer :: currNode => null()
    type(listnodetype), pointer :: tsNode0 => null()
    type(timeseriesrecordtype), pointer :: tsr => null()
    type(timeseriesrecordtype), pointer :: tsrec0 => null()
    class(*), pointer :: obj => null()
    !
    tslnode => null()
    if (associated(this%list%firstNode)) then
      currnode => this%list%firstNode
    else
      call store_error('probable programming error in &
                       &get_latest_preceding_node', &
                       terminate=.true.)
    end if
    !
    ! -- If the next node is earlier than time of interest, advance along
    !    linked list until the next node is later than time of interest.
    do
      if (associated(currnode)) then
        if (associated(currnode%nextNode)) then
          obj => currnode%nextNode%GetItem()
          tsr => castastimeseriesrecordtype(obj)
          if (tsr%tsrTime < time .or. is_close(tsr%tsrTime, time)) then
            currnode => currnode%nextNode
          else
            exit
          end if
        else
          ! -- read another record
          if (.not. this%read_next_record()) exit
        end if
      else
        exit
      end if
    end do
    !
    if (associated(currnode)) then
      !
      ! -- find earlier record
      tsnode0 => currnode
      obj => tsnode0%GetItem()
      tsrec0 => castastimeseriesrecordtype(obj)
      time0 = tsrec0%tsrTime
      do while (time0 > time)
        if (associated(tsnode0%prevNode)) then
          tsnode0 => tsnode0%prevNode
          obj => tsnode0%GetItem()
          tsrec0 => castastimeseriesrecordtype(obj)
          time0 = tsrec0%tsrTime
        else
          exit
        end if
      end do
    end if
    !
    if (time0 < time .or. is_close(time0, time)) tslnode => tsnode0
  end subroutine get_latest_preceding_node
 
  !> @brief Deallocate
  !<
  subroutine ts_da(this)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    !
    if (associated(this%list)) then
      call this%list%Clear(.true.)
      deallocate (this%list)
    end if
  end subroutine ts_da
 
  !> @brief Add ts record
  !<
  subroutine addtimeseriesrecord(this, tsr)
    ! -- dummy
    class(timeseriestype) :: this
    type(timeseriesrecordtype), pointer, intent(inout) :: tsr
    ! -- local
    class(*), pointer :: obj => null()
    !
    obj => tsr
    call this%list%Add(obj)
  end subroutine addtimeseriesrecord
 
  !> @brief Get current ts record
  !<
  function getcurrenttimeseriesrecord(this) result(res)
    ! -- dummy
    class(timeseriestype) :: this
    ! -- result
    type(timeseriesrecordtype), pointer :: res
    ! -- local
    class(*), pointer :: obj => null()
    !
    obj => null()
    res => null()
    obj => this%list%GetItem()
    if (associated(obj)) then
      res => castastimeseriesrecordtype(obj)
    end if
  end function getcurrenttimeseriesrecord
 
  !> @brief Get previous ts record
  !<
  function getprevioustimeseriesrecord(this) result(res)
    ! -- dummy
    class(timeseriestype) :: this
    ! -- result
    type(timeseriesrecordtype), pointer :: res
    ! -- local
    class(*), pointer :: obj => null()
    !
    obj => null()
    res => null()
    obj => this%list%GetPreviousItem()
    if (associated(obj)) then
      res => castastimeseriesrecordtype(obj)
    end if
  end function getprevioustimeseriesrecord
 
  !> @brief Get next ts record
  !<
  function getnexttimeseriesrecord(this) result(res)
    ! -- dummy
    class(timeseriestype) :: this
    ! -- result
    type(timeseriesrecordtype), pointer :: res
    ! -- local
    class(*), pointer :: obj => null()
    !
    obj => null()
    res => null()
    obj => this%list%GetNextItem()
    if (associated(obj)) then
      res => castastimeseriesrecordtype(obj)
    end if
  end function getnexttimeseriesrecord
 
  !> @brief Get ts record
  !<
  function gettimeseriesrecord(this, time, epsi) result(res)
    ! -- dummy
    class(timeseriestype) :: this
    double precision, intent(in) :: time
    double precision, intent(in) :: epsi
    ! -- result
    type(timeseriesrecordtype), pointer :: res
    ! -- local
    type(timeseriesrecordtype), pointer :: tsr
    !
    call this%list%Reset()
    res => null()
    do
      tsr => this%GetNextTimeSeriesRecord()
      if (associated(tsr)) then
        if (is_close(tsr%tsrTime, time)) then
          res => tsr
          exit
        end if
        if (tsr%tsrTime > time) exit
      else
        exit
      end if
    end do
  end function gettimeseriesrecord
 
  !> @brief Reset
  !<
  subroutine reset(this)
    ! -- dummy
    class(timeseriestype) :: this
    !
    call this%list%Reset()
  end subroutine reset
 
  !> @brief Insert a time series record
  !<
  subroutine inserttsr(this, tsr)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    type(timeseriesrecordtype), pointer, intent(inout) :: tsr
    ! -- local
    double precision :: badtime, time, time0, time1
    type(timeseriesrecordtype), pointer :: tsrEarlier, tsrLater
    type(listnodetype), pointer :: nodeEarlier, nodeLater
    class(*), pointer :: obj => null()
    !
    badtime = -9.0d30
    time0 = badtime
    time1 = badtime
    time = tsr%tsrTime
    call this%get_surrounding_nodes(time, nodeearlier, nodelater)
    !
    if (associated(nodeearlier)) then
      obj => nodeearlier%GetItem()
      tsrearlier => castastimeseriesrecordtype(obj)
      if (associated(tsrearlier)) then
        time0 = tsrearlier%tsrTime
      end if
    end if
    !
    if (associated(nodelater)) then
      obj => nodelater%GetItem()
      tsrlater => castastimeseriesrecordtype(obj)
      if (associated(tsrlater)) then
        time1 = tsrlater%tsrTime
      end if
    end if
    !
    if (time0 > badtime) then
      ! Time0 is valid
      if (time1 > badtime) then
        ! Both time0 and time1 are valid
        if (time > time0 .and. time < time1) then
          ! Insert record between two list nodes
          obj => tsr
          call this%list%InsertBefore(obj, nodelater)
        else
          ! No need to insert a time series record, but if existing record
          ! for time of interest has NODATA as tsrValue, replace tsrValue
          if (time == time0 .and. tsrearlier%tsrValue == dnodata .and. &
              tsr%tsrValue /= dnodata) then
            tsrearlier%tsrValue = tsr%tsrValue
          elseif (time == time1 .and. tsrlater%tsrValue == dnodata .and. &
                  tsr%tsrValue /= dnodata) then
            tsrlater%tsrValue = tsr%tsrValue
          end if
        end if
      else
        ! Time0 is valid and time1 is invalid. Just add tsr to the list.
        call this%AddTimeSeriesRecord(tsr)
      end if
    else
      ! Time0 is invalid, so time1 must be for first node in list
      if (time1 > badtime) then
        ! Time 1 is valid
        if (time < time1) then
          ! Insert tsr at beginning of list
          obj => tsr
          call this%list%InsertBefore(obj, nodelater)
        elseif (time == time1) then
          ! No need to insert a time series record, but if existing record
          ! for time of interest has NODATA as tsrValue, replace tsrValue
          if (tsrlater%tsrValue == dnodata .and. tsr%tsrValue /= dnodata) then
            tsrlater%tsrValue = tsr%tsrValue
          end if
        end if
      else
        ! Both time0 and time1 are invalid. Just add tsr to the list.
        call this%AddTimeSeriesRecord(tsr)
      end if
    end if
  end subroutine inserttsr
 
  !> @brief Find latest time
  !<
  function findlatesttime(this, readToEnd) result(endtime)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    logical, intent(in), optional :: readtoend
    ! -- local
    integer :: nrecords
    type(timeseriesrecordtype), pointer :: tsr
    class(*), pointer :: obj => null()
    ! -- return
    double precision :: endtime
    !
    ! -- If the caller requested the very last time in the series (readToEnd is true), check that we have first read all records
    if (present(readtoend)) then
      if (readtoend) then
        do while (this%read_next_record())
        end do
      end if
    end if
    !
    nrecords = this%list%Count()
    obj => this%list%GetItem(nrecords)
    tsr => castastimeseriesrecordtype(obj)
    endtime = tsr%tsrTime
  end function findlatesttime
 
  !> @brief Clear the list of time series records
  !<
  subroutine clear(this, destroy)
    ! -- dummy
    class(timeseriestype), intent(inout) :: this
    logical, optional, intent(in) :: destroy
    !
    call this%list%Clear(destroy)
  end subroutine clear
 
! Type-bound procedures of TimeSeriesFileType
 
  !> @brief Count number of time series
  !<
  function count(this)
    ! -- return
    integer(I4B) :: count
    ! -- dummy
    class(timeseriesfiletype) :: this
    !
    if (associated(this%timeSeries)) then
      count = size(this%timeSeries)
    else
      count = 0
    end if
  end function count
 
  !> @brief Get time series
  !<
  function gettimeseries(this, indx) result(res)
    ! -- dummy
    class(timeseriesfiletype) :: this
    integer(I4B), intent(in) :: indx
    ! -- return
    type(timeseriestype), pointer :: res
    !
    res => null()
    if (indx > 0 .and. indx <= this%nTimeSeries) then
      res => this%timeSeries(indx)
    end if
  end function gettimeseries
 
  !> @brief Open time-series tsfile file and read options and first record,
  !! which may contain data to define multiple time series.
  !<
  subroutine initializetsfile(this, filename, iout, autoDeallocate)
    ! -- dummy
    class(timeseriesfiletype), target, intent(inout) :: this
    character(len=*), intent(in) :: filename
    integer(I4B), intent(in) :: iout
    logical, optional, intent(in) :: autoDeallocate
    ! -- local
    integer(I4B) :: iMethod, istatus, j, nwords
    integer(I4B) :: ierr, inunit
    logical :: autoDeallocateLocal = .true.
    logical :: continueread, found, endOfBlock
    logical :: methodWasSet
    real(DP) :: sfaclocal
    character(len=40) :: keyword, keyvalue
    character(len=:), allocatable :: line
    character(len=LENTIMESERIESNAME), allocatable, dimension(:) :: words
    !
    ! -- Initialize some variables
    if (present(autodeallocate)) autodeallocatelocal = autodeallocate
    imethod = undefined
    methodwasset = .false.
    !
    ! -- Assign members
    this%iout = iout
    this%datafile = filename
    !
    ! -- Open the time-series tsfile input file
    this%inunit = getunit()
    inunit = this%inunit
    call openfile(inunit, 0, filename, 'TS6')
    !
    ! -- Initialize block parser
    call this%parser%Initialize(this%inunit, this%iout)
    !
    ! -- Read the ATTRIBUTES block and count time series
    continueread = .false.
    ierr = 0
    !
    ! -- get BEGIN line of ATTRIBUTES block
    call this%parser%GetBlock('ATTRIBUTES', found, ierr, &
                              supportopenclose=.true.)
    if (ierr /= 0) then
      ! end of file
      errmsg = 'End-of-file encountered while searching for'// &
               ' ATTRIBUTES in time-series '// &
               'input file "'//trim(this%datafile)//'"'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    elseif (.not. found) then
      errmsg = 'ATTRIBUTES block not found in time-series '// &
               'tsfile input file "'//trim(this%datafile)//'"'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- parse ATTRIBUTES entries
    do
      ! -- read a line from input
      call this%parser%GetNextLine(endofblock)
      if (endofblock) exit
      !
      ! -- get the keyword
      call this%parser%GetStringCaps(keyword)
      !
      ! support either NAME or NAMES as equivalent keywords
      if (keyword == 'NAMES') keyword = 'NAME'
      !
      if (keyword /= 'NAME' .and. keyword /= 'METHODS' .and. &
          keyword /= 'SFACS') then
        ! -- get the word following the keyword (the key value)
        call this%parser%GetStringCaps(keyvalue)
      end if
      !
      select case (keyword)
      case ('NAME')
!        line = line(istart:linelen)
        call this%parser%GetRemainingLine(line)
        call parseline(line, nwords, words, this%parser%iuactive)
        this%nTimeSeries = nwords
        ! -- Allocate the timeSeries array and initialize each
        !    time series.
        allocate (this%timeSeries(this%nTimeSeries))
        do j = 1, this%nTimeSeries
          call this%timeSeries(j)%initialize_time_series(this, words(j), &
                                                         autodeallocatelocal)
        end do
      case ('METHOD')
        methodwasset = .true.
        if (this%nTimeSeries == 0) then
          errmsg = 'Error: NAME attribute not provided before METHOD in file: ' &
                   //trim(filename)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
        select case (keyvalue)
        case ('STEPWISE')
          imethod = stepwise
        case ('LINEAR')
          imethod = linear
        case ('LINEAREND')
          imethod = linearend
        case default
          errmsg = 'Unknown interpolation method: "'//trim(keyvalue)//'"'
          call store_error(errmsg)
        end select
        do j = 1, this%nTimeSeries
          this%timeSeries(j)%iMethod = imethod
        end do
      case ('METHODS')
        methodwasset = .true.
        if (this%nTimeSeries == 0) then
          errmsg = 'Error: NAME attribute not provided before METHODS in file: ' &
                   //trim(filename)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
        call this%parser%GetRemainingLine(line)
        call parseline(line, nwords, words, this%parser%iuactive)
        if (nwords < this%nTimeSeries) then
          errmsg = 'METHODS attribute does not list a method for'// &
                   ' all time series.'
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
        do j = 1, this%nTimeSeries
          call upcase(words(j))
          select case (words(j))
          case ('STEPWISE')
            imethod = stepwise
          case ('LINEAR')
            imethod = linear
          case ('LINEAREND')
            imethod = linearend
          case default
            errmsg = 'Unknown interpolation method: "'//trim(words(j))//'"'
            call store_error(errmsg)
          end select
          this%timeSeries(j)%iMethod = imethod
        end do
      case ('SFAC')
        if (this%nTimeSeries == 0) then
          errmsg = 'NAME attribute not provided before SFAC in file: ' &
                   //trim(filename)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
        read (keyvalue, *, iostat=istatus) sfaclocal
        if (istatus /= 0) then
          errmsg = 'Error reading numeric value from: "'//trim(keyvalue)//'"'
          call store_error(errmsg)
        end if
        do j = 1, this%nTimeSeries
          this%timeSeries(j)%sfac = sfaclocal
        end do
      case ('SFACS')
        if (this%nTimeSeries == 0) then
          errmsg = 'NAME attribute not provided before SFACS in file: ' &
                   //trim(filename)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
        do j = 1, this%nTimeSeries
          sfaclocal = this%parser%GetDouble()
          this%timeSeries(j)%sfac = sfaclocal
        end do
      case ('AUTODEALLOCATE')
        do j = 1, this%nTimeSeries
          this%timeSeries(j)%autoDeallocate = (keyvalue == 'TRUE')
        end do
      case default
        errmsg = 'Unknown option found in ATTRIBUTES block: "'// &
                 trim(keyword)//'"'
        call store_error(errmsg)
        call this%parser%StoreErrorUnit()
      end select
    end do
    !
    ! -- Get TIMESERIES block
    call this%parser%GetBlock('TIMESERIES', found, ierr, &
                              supportopenclose=.true.)
    !
    ! -- Read the first line of time-series data
    if (.not. this%read_tsfile_line()) then
      errmsg = 'Error: No time-series data contained in file: '// &
               trim(this%datafile)
      call store_error(errmsg)
    end if
    !
    ! -- Ensure method was set
    if (.not. methodwasset) then
      errmsg = 'Interpolation method was not set.  METHOD or METHODS &
      &must be specified in the ATTRIBUTES block for this time series file.'
      call store_error(errmsg)
    end if
    !
    ! -- Clean up and return
    if (allocated(words)) deallocate (words)
    !
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
  end subroutine initializetsfile
 
  !> @brief Read time series file line
  logical function read_tsfile_line(this)
    ! -- dummy
    class(timeseriesfiletype), intent(inout) :: this
    ! -- local
    real(dp) :: tsrtime, tsrvalue
    integer(I4B) :: i
    logical :: endofblock
    type(timeseriesrecordtype), pointer :: tsrecord => null()
    !
    read_tsfile_line = .false.
    !
    ! -- Get an arbitrary length, non-comment, non-blank line
    !    from the input file.
    call this%parser%GetNextLine(endofblock)
    !
    ! -- Check if we've reached the end of the TIMESERIES block
    if (endofblock) then
      return
    end if
    !
    ! -- Get the time
    tsrtime = this%parser%GetDouble()
    !
    ! -- Construct a new record and append a new node to each time series
    tsloop: do i = 1, this%nTimeSeries
      tsrvalue = this%parser%GetDouble()
      if (tsrvalue == dnodata) cycle tsloop
      ! -- multiply value by sfac
      tsrvalue = tsrvalue * this%timeSeries(i)%sfac
      call constructtimeseriesrecord(tsrecord, tsrtime, tsrvalue)
      call addtimeseriesrecordtolist(this%timeSeries(i)%list, tsrecord)
    end do tsloop
    read_tsfile_line = .true.
  end function read_tsfile_line
 
  !> @brief Deallocate memory
  !<
  subroutine tsf_da(this)
    ! -- dummy
    class(timeseriesfiletype), intent(inout) :: this
    ! -- local
    integer :: i, n
    type(timeseriestype), pointer :: ts => null()
    !
    n = this%Count()
    do i = 1, n
      ts => this%GetTimeSeries(i)
      if (associated(ts)) then
        call ts%da()
!        deallocate(ts)
      end if
    end do
    !
    deallocate (this%timeSeries)
    deallocate (this%parser)
  end subroutine tsf_da
 
end module timeseriesmodule