TimeArraySeries.f90

Go to the documentation of this file.

module timearrayseriesmodule
 
  use arrayreadersmodule, only: readarray
  use blockparsermodule, only: blockparsertype
  use constantsmodule, only: linelength, undefined, stepwise, linear, &
                             lentimeseriesname, lenmodelname, dzero, done
  use mathutilmodule, only: is_close
  use inputoutputmodule, only: getunit, openfile
  use kindmodule, only: dp, i4b
  use listmodule, only: listtype
  use listnodemodule, only: listnodetype
  use simvariablesmodule, only: errmsg
  use simmodule, only: count_errors, store_error, store_error_unit
  use timearraymodule, only: timearraytype, constructtimearray, &
                             addtimearraytolist, castastimearraytype, &
                             gettimearrayfromlist
  use, intrinsic :: iso_fortran_env, only: iostat_end
 
  implicit none
  private
  public :: timearrayseriestype, constructtimearrayseries, &
            castastimearrayseriestype, gettimearrayseriesfromlist
 
  type timearrayseriestype
    ! -- Public members
    character(len=LENTIMESERIESNAME), public :: name = ''
    ! -- Private members
    integer(I4B), private :: inunit = 0
    integer(I4B), private :: iout = 0
    integer(I4B), private :: imethod = undefined
    real(dp), private :: sfac = done
    character(len=LINELENGTH), private :: datafile = ''
    logical, private :: autodeallocate = .true.
    type(listtype), pointer, private :: list => null()
    character(len=LENMODELNAME) :: modelname
    type(blockparsertype), private :: parser
 
  contains
 
    ! -- Public procedures
    procedure, public :: tas_init
    procedure, public :: getaveragevalues
    procedure, public :: getinunit
    procedure, public :: da => tas_da
    ! -- Private procedures
    procedure, private :: get_integrated_values
    procedure, private :: get_latest_preceding_node
    procedure, private :: get_values_at_time
    procedure, private :: get_surrounding_records
    procedure, private :: read_next_array
    procedure, private :: deallocatebackward
  end type timearrayseriestype
 
contains
 
  ! -- Constructor for TimeArraySeriesType
 
  !> @brief Allocate a new TimeArraySeriesType object.
  !<
  subroutine constructtimearrayseries(newTas, filename)
    ! -- dummy
    type(timearrayseriestype), pointer, intent(out) :: newtas
    character(len=*), intent(in) :: filename
    ! -- local
    logical :: lex
    ! -- formats
10  format('Error: Time-array-series file "', a, '" does not exist.')
    !
    ! -- Allocate a new object of type TimeArraySeriesType
    allocate (newtas)
    allocate (newtas%list)
    !
    ! -- Ensure that input file exists
    inquire (file=filename, exist=lex)
    if (.not. lex) then
      write (errmsg, 10) trim(filename)
      call store_error(errmsg, terminate=.true.)
    end if
    newtas%datafile = filename
  end subroutine constructtimearrayseries
 
  ! -- Public procedures
 
  !> @brief Initialize the time array series
  !<
  subroutine tas_init(this, fname, modelname, iout, tasname, autoDeallocate)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    character(len=*), intent(in) :: fname
    character(len=*), intent(in) :: modelname
    integer(I4B), intent(in) :: iout
    character(len=*), intent(inout) :: tasname
    logical, optional, intent(in) :: autoDeallocate
    ! -- local
    integer(I4B) :: istatus
    integer(I4B) :: ierr
    integer(I4B) :: inunit
    character(len=40) :: keyword, keyvalue
    logical :: found, continueread, endOfBlock
    !
    ! -- initialize some variables
    if (present(autodeallocate)) this%autoDeallocate = autodeallocate
    this%dataFile = fname
    allocate (this%list)
    !
    ! -- assign members
    this%modelname = modelname
    this%iout = iout
    !
    ! -- open time-array series input file
    inunit = getunit()
    this%inunit = inunit
    call openfile(inunit, 0, fname, 'TAS6')
    !
    ! -- initialize block parser
    call this%parser%Initialize(this%inunit, this%iout)
    !
    ! -- read ATTRIBUTES block
    continueread = .false.
    ierr = 0
    !
    ! -- get BEGIN line of ATTRIBUTES block
    call this%parser%GetBlock('ATTRIBUTES', found, ierr, &
                              supportopenclose=.true.)
    if (.not. found) then
      errmsg = 'Error: Attributes block not found in file: '// &
               trim(fname)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- parse ATTRIBUTES entries
    do
      ! -- read line from input
      call this%parser%GetNextLine(endofblock)
      if (endofblock) exit
      !
      ! -- get the keyword
      call this%parser%GetStringCaps(keyword)
      !
      ! -- get the word following the keyword (the key value)
      call this%parser%GetStringCaps(keyvalue)
      select case (keyword)
      case ('NAME')
        this%Name = keyvalue
        tasname = keyvalue
      case ('METHOD')
        select case (keyvalue)
        case ('STEPWISE')
          this%iMethod = stepwise
        case ('LINEAR')
          this%iMethod = linear
        case default
          errmsg = 'Unknown interpolation method: "'//trim(keyvalue)//'"'
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      case ('AUTODEALLOCATE')
        this%autoDeallocate = (keyvalue == 'TRUE')
      case ('SFAC')
        read (keyvalue, *, iostat=istatus) this%sfac
        if (istatus /= 0) then
          errmsg = 'Error reading numeric SFAC value from "'//trim(keyvalue) &
                   //'"'
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end if
      case default
        errmsg = 'Unknown option found in ATTRIBUTES block: "'// &
                 trim(keyword)//'"'
        call store_error(errmsg)
        call this%parser%StoreErrorUnit()
      end select
    end do
    !
    ! -- ensure that NAME and METHOD have been specified
    if (this%Name == '') then
      errmsg = 'Name not specified for time array series in file: '// &
               trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    if (this%iMethod == undefined) then
      errmsg = 'Interpolation method not specified for time'// &
               ' array series in file: '//trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- handle any errors encountered so far
    if (count_errors() > 0) then
      errmsg = 'Error(s) encountered initializing time array series from file: ' &
               //trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- try to read first time array into linked list
    if (.not. this%read_next_array()) then
      errmsg = 'Error encountered reading time-array data from file: '// &
               trim(this%dataFile)
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
  end subroutine tas_init
 
  !> @brief Populate an array time-weighted average value for a specified time
  !! span
  !<
  subroutine getaveragevalues(this, nvals, values, time0, time1)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    integer(I4B), intent(in) :: nvals
    real(DP), dimension(nvals), intent(inout) :: values
    real(DP), intent(in) :: time0
    real(DP), intent(in) :: time1
    ! -- local
    integer(I4B) :: i
    real(DP) :: timediff
    !
    timediff = time1 - time0
    if (timediff > 0) then
      call this%get_integrated_values(nvals, values, time0, time1)
      do i = 1, nvals
        values(i) = values(i) / timediff
      end do
    else
      ! -- time0 and time1 are the same, so skip the integration step.
      call this%get_values_at_time(nvals, values, time0)
    end if
  end subroutine getaveragevalues
 
  !> @brief Return unit number
  !<
  function getinunit(this)
    ! -- return
    integer(I4B) :: getinunit
    ! -- dummy
    class(timearrayseriestype) :: this
    !
    getinunit = this%inunit
  end function getinunit
 
  ! -- Private procedures
 
  !> @brief Get surrounding records
  !<
  subroutine get_surrounding_records(this, time, taEarlier, taLater)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    real(DP), intent(in) :: time
    type(timearraytype), pointer, intent(inout) :: taEarlier
    type(timearraytype), pointer, intent(inout) :: taLater
    ! -- local
    real(DP) :: time0, time1
    type(listnodetype), pointer :: currNode => null()
    type(listnodetype), pointer :: node0 => null()
    type(listnodetype), pointer :: node1 => null()
    type(timearraytype), pointer :: ta => null(), ta0 => null(), ta1 => null()
    class(*), pointer :: obj
    !
    taearlier => null()
    talater => 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()
          ta => castastimearraytype(obj)
          if (ta%taTime <= time) then
            currnode => currnode%nextNode
          else
            exit
          end if
        else
          ! -- read another array
          if (.not. this%read_next_array()) exit
        end if
      else
        exit
      end if
    end do
    !
    if (associated(currnode)) then
      !
      ! -- find earlier record
      node0 => currnode
      obj => node0%GetItem()
      ta0 => castastimearraytype(obj)
      time0 = ta0%taTime
      do while (time0 > time)
        if (associated(node0%prevNode)) then
          node0 => node0%prevNode
          obj => node0%GetItem()
          ta0 => castastimearraytype(obj)
          time0 = ta0%taTime
        else
          exit
        end if
      end do
      !
      ! -- find later record
      node1 => currnode
      obj => node1%GetItem()
      ta1 => castastimearraytype(obj)
      time1 = ta1%taTime
      do while (time1 < time)
        if (associated(node1%nextNode)) then
          node1 => node1%nextNode
          obj => node1%GetItem()
          ta1 => castastimearraytype(obj)
          time1 = ta1%taTime
        else
          ! -- get next array
          if (.not. this%read_next_array()) then
            ! -- end of file reached, so exit loop
            exit
          end if
        end if
      end do
      !
    end if
    !
    if (time0 <= time) taearlier => ta0
    if (time1 >= time) talater => ta1
  end subroutine get_surrounding_records
 
  !> @brief Read next time array from input file and append to list
  !<
  logical function read_next_array(this)
    ! -- modules
    use constantsmodule, only: lenmempath
    use memorymanagermodule, only: mem_setptr
    use memoryhelpermodule, only: create_mem_path
    use simvariablesmodule, only: idm_context
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    ! -- local
    integer(I4B) :: i, ierr, istart, istat, istop, lloc, nrow, ncol, &
                    nodesperlayer
    logical :: lopen, isfound
    type(timearraytype), pointer :: ta => null()
    character(len=LENMEMPATH) :: mempath
    integer(I4B), dimension(:), contiguous, pointer :: mshape
    !
    ! -- initialize
    istart = 1
    istat = 0
    istop = 1
    lloc = 1
    nullify (mshape)
    !
    ! -- create mempath to model input context where MODEL_SHAPE is
    ! -- stored during static DIS loading
    mempath = create_mem_path(component=this%modelname, context=idm_context)
    !
    ! -- set mshape pointer
    call mem_setptr(mshape, 'MODEL_SHAPE', mempath)
    !
    ! Get dimensions for supported discretization type
    if (size(mshape) == 2) then
      nodesperlayer = mshape(2)
      nrow = 1
      ncol = mshape(2)
    else if (size(mshape) == 3) then
      nodesperlayer = mshape(2) * mshape(3)
      nrow = mshape(2)
      ncol = mshape(3)
    else
      errmsg = 'Time array series is not supported for selected &
               &discretization type.'
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
    !
    read_next_array = .false.
    inquire (unit=this%inunit, opened=lopen)
    if (lopen) then
      call constructtimearray(ta, this%modelname)
      ! -- read a time and an array from the input file
      ! -- Get a TIME block and read the time
      call this%parser%GetBlock('TIME', isfound, ierr, &
                                supportopenclose=.false.)
      if (isfound) then
        ta%taTime = this%parser%GetDouble()
        ! -- Read the array
        call readarray(this%parser%iuactive, ta%taArray, this%Name, &
                       size(mshape), ncol, nrow, 1, nodesperlayer, &
                       this%iout, 0, 0)
        !
        ! -- multiply values by sfac
        do i = 1, nodesperlayer
          ta%taArray(i) = ta%taArray(i) * this%sfac
        end do
        !
        ! -- append the new time array to the list
        call addtimearraytolist(this%list, ta)
        read_next_array = .true.
        !
        ! -- make sure block is closed
        call this%parser%terminateblock()
      end if
    end if
  end function read_next_array
 
  !> @brief Return an array of values for a specified time, same units as
  !! time-series values
  !<
  subroutine get_values_at_time(this, nvals, values, time)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    integer(I4B), intent(in) :: nvals
    real(DP), dimension(nvals), intent(inout) :: values
    real(DP), intent(in) :: time ! time of interest
    ! -- local
    integer(I4B) :: i, ierr
    real(DP) :: ratio, time0, time1, timediff, timediffi, val0, val1, &
                valdiff
    type(timearraytype), pointer :: taEarlier => null()
    type(timearraytype), pointer :: taLater => null()
    ! -- formats
10  format('Error getting array at time ', g10.3, &
           ' for time-array series "', a, '"')
    !
    ierr = 0
    call this%get_surrounding_records(time, taearlier, talater)
    if (associated(taearlier)) then
      if (associated(talater)) then
        ! -- values are available for both earlier and later times
        if (this%iMethod == stepwise) then
          ! -- Just populate values from elements of earlier time array
          values = taearlier%taArray
        elseif (this%iMethod == linear) then
          ! -- perform linear interpolation
          time0 = taearlier%taTime
          time1 = talater%tatime
          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
          ! -- Iterate through all elements and perform interpolation.
          do i = 1, nvals
            val0 = taearlier%taArray(i)
            val1 = talater%taArray(i)
            valdiff = val1 - val0
            values(i) = val0 + (ratio * valdiff)
          end do
        else
          ierr = 1
        end if
      else
        if (is_close(taearlier%taTime, time)) then
          values = taearlier%taArray
        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
            values = taearlier%taArray
          else
            ierr = 1
          end if
        end if
      end if
    else
      if (associated(talater)) then
        if (is_close(talater%taTime, time)) then
          values = talater%taArray
        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)
      call store_error_unit(this%inunit)
    end if
  end subroutine get_values_at_time
 
  !> @brief Populates an array with integrated values for a specified time span
  !!
  !! Units: (ts-value-unit)*time
  !<
  subroutine get_integrated_values(this, nvals, values, time0, time1)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    integer(I4B), intent(in) :: nvals
    real(DP), dimension(nvals), intent(inout) :: values
    real(DP), intent(in) :: time0
    real(DP), intent(in) :: time1
    ! -- local
    integer(I4B) :: i
    real(DP) :: area, currTime, nextTime, ratio0, ratio1, t0, &
                t01, t1, timediff, value, value0, value1, valuediff
    logical :: ldone
    type(listnodetype), pointer :: precNode => null()
    type(listnodetype), pointer :: currNode => null(), nextnode => null()
    type(timearraytype), pointer :: currRecord => null(), nextrecord => null()
    class(*), pointer :: currObj => null(), nextobj => null()
    ! -- formats
10  format('Error encountered while performing integration', &
           ' for time-array series "', a, '" for time interval: ', &
           g12.5, ' to ', g12.5)
    !
    values = dzero
    value = dzero
    ldone = .false.
    t1 = -done
    call this%get_latest_preceding_node(time0, precnode)
    if (associated(precnode)) then
      currnode => precnode
      do while (.not. ldone)
        currobj => currnode%GetItem()
        currrecord => castastimearraytype(currobj)
        currtime = currrecord%taTime
        if (currtime < time1) then
          if (.not. associated(currnode%nextNode)) then
            ! -- try to read the next array
            if (.not. this%read_next_array()) then
              write (errmsg, 10) trim(this%Name), time0, time1
              call store_error(errmsg)
              call store_error_unit(this%inunit)
            end if
          end if
          if (associated(currnode%nextNode)) then
            nextnode => currnode%nextNode
            nextobj => nextnode%GetItem()
            nextrecord => castastimearraytype(nextobj)
            nexttime = nextrecord%taTime
            ! -- determine lower and upper limits of time span of interest
            !    within current interval
            if (currtime >= time0) then
              t0 = currtime
            else
              t0 = time0
            end if
            if (nexttime <= time1) then
              t1 = nexttime
            else
              t1 = time1
            end if
            ! -- For each element, find area of rectangle
            !    or trapezoid delimited by t0 and t1.
            t01 = t1 - t0
            select case (this%iMethod)
            case (stepwise)
              do i = 1, nvals
                ! -- compute area of a rectangle
                value0 = currrecord%taArray(i)
                area = value0 * t01
                ! -- add area to integrated value
                values(i) = values(i) + area
              end do
            case (linear)
              do i = 1, nvals
                ! -- compute area of a trapezoid
                timediff = nexttime - currtime
                ratio0 = (t0 - currtime) / timediff
                ratio1 = (t1 - currtime) / timediff
                valuediff = nextrecord%taArray(i) - currrecord%taArray(i)
                value0 = currrecord%taArray(i) + ratio0 * valuediff
                value1 = currrecord%taArray(i) + ratio1 * valuediff
                area = 0.5d0 * t01 * (value0 + value1)
                ! -- add area to integrated value
                values(i) = values(i) + area
              end do
            end select
          else
            write (errmsg, 10) trim(this%Name), time0, time1
            call store_error(errmsg)
            call store_error('(Probable programming error)', terminate=.true.)
          end if
        else
          ! Current node time = time1 so should be done
          ldone = .true.
        end if
        !
        ! -- Are we done yet?
        if (t1 >= time1) then
          ldone = .true.
        else
          if (.not. associated(currnode%nextNode)) then
            ! -- try to read the next array
            if (.not. this%read_next_array()) then
              write (errmsg, 10) trim(this%Name), time0, time1
              call store_error(errmsg)
              call this%parser%StoreErrorUnit()
            end if
          end if
          if (associated(currnode%nextNode)) then
            currnode => currnode%nextNode
          else
            write (errmsg, 10) trim(this%Name), time0, time1
            call store_error(errmsg)
            call store_error('(Probable programming error)', terminate=.true.)
          end if
        end if
      end do
    end if
    !
    if (this%autoDeallocate) then
      if (associated(precnode)) then
        if (associated(precnode%prevNode)) then
          call this%DeallocateBackward(precnode%prevNode)
        end if
      end if
    end if
  end subroutine get_integrated_values
 
  !> @brief Deallocate fromNode and all previous nodes in list;
  !! reassign firstNode
  !<
  subroutine deallocatebackward(this, fromNode)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    type(listnodetype), pointer, intent(inout) :: fromNode
    !
    ! -- local
    type(listnodetype), pointer :: current => null()
    type(listnodetype), pointer :: prev => null()
    type(timearraytype), pointer :: ta => null()
    class(*), pointer :: obj => null()
    !
    if (associated(fromnode)) then
      ! -- reassign firstNode
      if (associated(fromnode%nextNode)) then
        this%list%firstNode => fromnode%nextNode
      else
        this%list%firstNode => null()
      end if
      ! -- deallocate fromNode and all previous nodes
      current => fromnode
      do while (associated(current))
        prev => current%prevNode
        obj => current%GetItem()
        ta => castastimearraytype(obj)
        ! -- Deallocate the contents of this time array,
        !    then remove it from the list
        call ta%da()
        call this%list%RemoveNode(current, .true.)
        current => prev
      end do
      fromnode => null()
    end if
  end subroutine deallocatebackward
 
  !> @brief 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(timearrayseriestype), 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 :: node0 => null()
    type(timearraytype), pointer :: ta => null()
    type(timearraytype), pointer :: ta0 => 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
    !
    continue
    ! -- 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()
          ta => castastimearraytype(obj)
          if (ta%taTime < time .or. is_close(ta%taTime, time)) then
            currnode => currnode%nextNode
          else
            exit
          end if
        else
          ! -- read another record
          if (.not. this%read_next_array()) exit
        end if
      else
        exit
      end if
    end do
    !
    if (associated(currnode)) then
      !
      ! -- find earlier record
      node0 => currnode
      obj => node0%GetItem()
      ta0 => castastimearraytype(obj)
      time0 = ta0%taTime
      do while (time0 > time)
        if (associated(node0%prevNode)) then
          node0 => node0%prevNode
          obj => node0%GetItem()
          ta0 => castastimearraytype(obj)
          time0 = ta0%taTime
        else
          exit
        end if
      end do
    end if
    !
    if (time0 <= time) tslnode => node0
  end subroutine get_latest_preceding_node
 
  !> @brief Deallocate memory
  !<
  subroutine tas_da(this)
    ! -- dummy
    class(timearrayseriestype), intent(inout) :: this
    ! -- local
    integer :: i, n
    type(timearraytype), pointer :: ta => null()
    !
    ! -- Deallocate contents of each time array in list
    n = this%list%Count()
    do i = 1, n
      ta => gettimearrayfromlist(this%list, i)
      call ta%da()
    end do
    !
    ! -- Deallocate the list of time arrays
    call this%list%Clear(.true.)
    deallocate (this%list)
  end subroutine tas_da
 
  ! -- Procedures not type-bound
 
  !> @brief Cast an unlimited polymorphic object as class(TimeArraySeriesType)
  !<
  function castastimearrayseriestype(obj) result(res)
    ! -- dummy
    class(*), pointer, intent(inout) :: obj
    ! -- return
    type(timearrayseriestype), pointer :: res
    !
    res => null()
    if (.not. associated(obj)) return
    !
    select type (obj)
    type is (timearrayseriestype)
      res => obj
    end select
  end function castastimearrayseriestype
 
  !> @brief Get time array from list
  !<
  function gettimearrayseriesfromlist(list, indx) result(res)
    ! -- dummy
    type(listtype), intent(inout) :: list
    integer, intent(in) :: indx
    ! -- return
    type(timearrayseriestype), pointer :: res
    ! -- local
    class(*), pointer :: obj
    !
    obj => list%GetItem(indx)
    res => castastimearrayseriestype(obj)
  end function gettimearrayseriesfromlist
 
end module timearrayseriesmodule