tsp-mvt.f90

Go to the documentation of this file.

! -- Groundwater Transport Mover Module
! -- This module is responsible for sending mass from providers into
! -- receiver qmfrommvr arrays and writing a mover transport budget
 
module tspmvtmodule
 
  use kindmodule, only: dp, i4b, lgp
  use constantsmodule, only: linelength, maxcharlen, dzero, lenpakloc, &
                             dnodata, lenpackagename, tabcenter, &
                             lenmodelname, lenvarname
 
  use simmodule, only: store_error
  use basedismodule, only: disbasetype
  use numericalpackagemodule, only: numericalpackagetype
  use tspfmimodule, only: tspfmitype
  use budgetmodule, only: budgettype, budget_cr
  use budgetobjectmodule, only: budgetobjecttype, budgetobject_cr
  use tablemodule, only: tabletype, table_cr
 
  implicit none
 
  private
  public :: tspmvttype
  public :: mvt_cr
 
  !> @brief Type to store mover-transport flow and concentration for one budget term
  !<í
  type :: mvrtermtype
    logical(LGP) :: prov_is_m1 !< true when provider is in model1, false when in model 2
    real(dp), dimension(:), pointer, contiguous :: qty => null() !< concentration or temperature of provider
  end type mvrtermtype
 
  type, extends(numericalpackagetype) :: tspmvttype
 
    character(len=LENMODELNAME) :: gwfmodelname1 = '' !< name of model 1
    character(len=LENMODELNAME) :: gwfmodelname2 = '' !< name of model 2 (set to modelname 1 for single model MVT)
    integer(I4B), pointer :: maxpackages !< max number of packages
    integer(I4B), pointer :: ibudgetout => null() !< unit number for budget output file
    integer(I4B), pointer :: ibudcsv => null() !< unit number for csv budget output file
    real(dp), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy
    type(tspfmitype), pointer :: fmi1 => null() !< pointer to fmi object for model 1
    type(tspfmitype), pointer :: fmi2 => null() !< pointer to fmi object for model 2 (set to fmi1 for single model)
    type(budgettype), pointer :: budget => null() !< mover transport budget object (used to write balance table)
    type(budgetobjecttype), pointer :: budobj => null() !< budget container (used to write binary file)
    type(budgetobjecttype), pointer :: mvrbudobj => null() !< pointer to the water mover budget object
    character(len=LENPACKAGENAME), &
      dimension(:), pointer, contiguous :: paknames => null() !< array of package names
    character(len=LENVARNAME) :: depvartype = ''
    !
    ! -- table objects
    type(tabletype), pointer :: outputtab => null()
    !
    ! -- mover transport terms
    type(mvrtermtype), dimension(:), pointer, contiguous :: mvrterm => null() !< qty and flow per budget term
 
  contains
    procedure :: mvt_init
    procedure :: mvt_df
    procedure :: mvt_ar
    procedure :: mvt_rp
    procedure :: mvt_fc
    procedure :: mvt_fill_mvrterm
    procedure :: mvt_update_qmfrommvr
    procedure :: mvt_cc
    procedure :: mvt_bd
    procedure :: mvt_ot_saveflow
    procedure :: mvt_ot_printflow
    procedure :: mvt_ot_bdsummary
    procedure :: mvt_da
    procedure :: allocate_scalars
    procedure :: read_options
    procedure :: mvt_setup_budobj
    procedure :: mvt_fill_budobj
    procedure :: mvt_scan_mvrbudobj
    procedure :: set_pointer_mvrbudobj
    procedure :: set_fmi
    procedure, private :: mvt_setup_outputtab
    procedure, private :: mvt_print_outputtab
  end type tspmvttype
 
contains
 
  subroutine mvt_cr(mvt, name_model, inunit, iout, fmi1, eqnsclfac, &
                    depvartype, gwfmodelname1, gwfmodelname2, fmi2)
    type(tspmvttype), pointer :: mvt
    character(len=*), intent(in) :: name_model
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    type(tspfmitype), intent(in), pointer :: fmi1
    real(dp), intent(in), pointer :: eqnsclfac !< governing equation scale factor
    character(len=LENVARNAME), intent(in) :: depvartype !< dependent variable type ('concentration' or 'temperature')
    character(len=*), intent(in), optional :: gwfmodelname1
    character(len=*), intent(in), optional :: gwfmodelname2
    type(tspfmitype), intent(in), pointer, optional :: fmi2
 
    allocate (mvt)
    call mvt%mvt_init(name_model, inunit, iout, fmi1, eqnsclfac, &
                      depvartype, gwfmodelname1, gwfmodelname2, fmi2)
 
  end subroutine mvt_cr
 
  !> @brief Create a new mover transport object
  !<
  subroutine mvt_init(this, name_model, inunit, iout, fmi1, eqnsclfac, &
                      depvartype, gwfmodelname1, gwfmodelname2, fmi2)
    ! -- dummy
    class(tspmvttype) :: this
    character(len=*), intent(in) :: name_model
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    type(tspfmitype), intent(in), pointer :: fmi1
    real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor
    character(len=LENVARNAME), intent(in) :: depvartype !< dependent variable type ('concentration' or 'temperature')
    character(len=*), intent(in), optional :: gwfmodelname1
    character(len=*), intent(in), optional :: gwfmodelname2
    type(tspfmitype), intent(in), pointer, optional :: fmi2
    !
    ! -- Create name and memory path
    call this%set_names(1, name_model, 'MVT', 'MVT')
    !
    ! -- Allocate scalars
    call this%allocate_scalars()
    !
    this%inunit = inunit
    this%iout = iout
    !
    ! -- Assume that this MVT is owned by a GWT Model
    this%fmi1 => fmi1
    this%fmi2 => fmi1
    !
    ! -- Set pointers
    if (present(fmi2)) then
      this%fmi2 => fmi2
    end if
    !
    ! -- Set model names
    if (present(gwfmodelname1)) then
      this%gwfmodelname1 = gwfmodelname1
    end if
    if (present(gwfmodelname2)) then
      this%gwfmodelname2 = gwfmodelname2
    end if
    !
    ! -- Create the budget object
    call budgetobject_cr(this%budobj, 'TRANSPORT MOVER')
    !
    ! -- Store pointer to governing equation scale factor
    this%eqnsclfac => eqnsclfac
    !
    ! -- Store pointer to labels associated with the current model so that the
    !    package has access to the corresponding dependent variable type
    this%depvartype = depvartype
 
  end subroutine mvt_init
 
  !> @brief Define mover transport object
  !<
  subroutine mvt_df(this, dis)
    ! -- dummy
    class(tspmvttype) :: this
    class(disbasetype), pointer, intent(in) :: dis
    ! -- formats
    character(len=*), parameter :: fmtmvt = &
      "(1x,/1x,'MVT -- MOVER TRANSPORT PACKAGE, VERSION 1, 4/15/2020', &
      &' INPUT READ FROM UNIT ', i0, //)"
    !
    ! -- Set pointer to dis
    this%dis => dis
    !
    ! -- Print a message identifying the MVT package.
    write (this%iout, fmtmvt) this%inunit
    !
    ! -- Initialize block parser
    call this%parser%Initialize(this%inunit, this%iout)
    !
    ! -- Initialize the budget table writer
    call budget_cr(this%budget, this%memoryPath)
    !
    ! -- Read mvt options
    call this%read_options()
  end subroutine mvt_df
 
  !> @ brief Set pointer to mvrbudobj
  !!
  !! Store a pointer to mvrbudobj, which contains the simulated water
  !! mover flows from either a gwf model MVR package or from a gwf-gwf
  !! exchange MVR package.
  !<
  subroutine set_pointer_mvrbudobj(this, mvrbudobj)
    class(tspmvttype) :: this
    type(budgetobjecttype), intent(in), target :: mvrbudobj
    this%mvrbudobj => mvrbudobj
  end subroutine set_pointer_mvrbudobj
 
  !> @brief Allocate and read mover-for-transport information
  !<
  subroutine mvt_ar(this)
    ! -- dummy
    class(tspmvttype) :: this
    !
    ! -- Setup the output table
    call this%mvt_setup_outputtab()
  end subroutine mvt_ar
 
  !> @brief Read and prepare mover transport object
  !<
  subroutine mvt_rp(this)
    ! -- modules
    use tdismodule, only: kper, kstp
    ! -- dummy
    class(tspmvttype) :: this
    !
    ! -- At this point, the mvrbudobj is available to set up the mvt budobj
    if (kper * kstp == 1) then
      !
      ! -- If mvt is for a single model then point to fmi1
      !cdl todo: this needs to be called from GwtGwtExg somehow for the 2 model case
      if (associated(this%fmi1, this%fmi2)) then
        call this%set_pointer_mvrbudobj(this%fmi1%mvrbudobj)
      end if
      !
      ! -- Set up the mvt budobject
      call this%mvt_scan_mvrbudobj()
      call this%mvt_setup_budobj()
      !
      ! -- Define the budget object to be the size of maxpackages
      call this%budget%budget_df(this%maxpackages, 'TRANSPORT MOVER', bddim='M')
      call this%budget%set_ibudcsv(this%ibudcsv)
    end if
  end subroutine mvt_rp
 
  !> @brief Calculate coefficients and fill amat and rhs
  !!
  !! The mvt package adds the mass flow rate to the provider qmfrommvr array.
  !! The advanced packages know enough to subtract any mass that is leaving, so
  !! the mvt just adds mass coming in from elsewhere.  Because the movers
  !! change by stress period, their solute effects must be added to the right-
  !! hand side of the transport matrix equations.
  !<
  subroutine mvt_fc(this, cnew1, cnew2)
    class(tspmvttype) :: this
    real(DP), intent(in), dimension(:), contiguous, target :: cnew1
    real(DP), intent(in), dimension(:), contiguous, target :: cnew2
 
    call this%mvt_fill_mvrterm(cnew1, cnew2)
    call this%mvt_update_qmfrommvr()
 
  end subroutine mvt_fc
 
  !> @brief Fill mvrterm with provider concentration and flow for each budget term
  !<
  subroutine mvt_fill_mvrterm(this, cnew1, cnew2)
    class(tspmvttype) :: this
    real(DP), intent(in), dimension(:), contiguous, target :: cnew1
    real(DP), intent(in), dimension(:), contiguous, target :: cnew2
    ! local
    integer(I4B) :: i, n
    integer(I4B) :: id1, nlist
    integer(I4B) :: ipr
    integer(I4B) :: igwtnode
    integer(I4B) :: nbudterm
    real(DP) :: cp
    real(DP), dimension(:), pointer :: concpak
    real(DP), dimension(:), contiguous, pointer :: cnew_pr !, pointer to provider dep. var.
    type(tspfmitype), pointer :: fmi_pr !< pointer to provider model fmi package
 
    nbudterm = this%mvrbudobj%nbudterm
    do i = 1, nbudterm
      nlist = this%mvrbudobj%budterm(i)%nlist
      if (nlist == 0) cycle
 
      ! fetch fmi package
      call this%set_fmi(i, fmi_pr, .true.)
      if (.not. associated(fmi_pr)) cycle
 
      cnew_pr => cnew1
      if (associated(fmi_pr, this%fmi2)) then
        cnew_pr => cnew2
      end if
 
      call fmi_pr%get_package_index(this%mvrbudobj%budterm(i)%text2id1, ipr)
 
      ! use package concentration for APT
      if (fmi_pr%iatp(ipr) /= 0) then
        concpak => fmi_pr%datp(ipr)%concpack
      end if
 
      do n = 1, nlist
 
        ! lak/sfr/maw/uzf id1 (provider)
        id1 = this%mvrbudobj%budterm(i)%id1(n)
        cp = dzero
        if (fmi_pr%iatp(ipr) /= 0) then
          ! concentration from APT
          cp = concpak(id1)
        else
          igwtnode = fmi_pr%gwfpackages(ipr)%nodelist(id1)
          cp = cnew_pr(igwtnode)
        end if
 
        this%mvrterm(i)%qty(n) = cp
      end do
    end do
  end subroutine mvt_fill_mvrterm
 
  !> @brief Assign mover transport mass into each receiver's qmfrommvr
  !<
  subroutine mvt_update_qmfrommvr(this)
    class(tspmvttype) :: this
    ! local
    integer(I4B) :: i, n
    integer(I4B) :: id2, nlist
    integer(I4B) :: irc
    integer(I4B) :: nbudterm
    type(tspfmitype), pointer :: fmi_rc !< pointer to receiver model fmi package
 
    nbudterm = this%mvrbudobj%nbudterm
    do i = 1, nbudterm
      nlist = this%mvrbudobj%budterm(i)%nlist
      if (nlist == 0) cycle
 
      call this%set_fmi(i, fmi_rc, .false.)
      if (.not. associated(fmi_rc)) cycle
 
      call fmi_rc%get_package_index(this%mvrbudobj%budterm(i)%text2id2, irc)
 
      do n = 1, nlist
 
        ! lak/sfr/maw/uzf id2 (receiver)
        id2 = this%mvrbudobj%budterm(i)%id2(n)
 
        ! add the mover rate times the provider concentration into the receiver;
        ! accumulated since multiple providers can move water into the same receiver
        if (fmi_rc%iatp(irc) /= 0) then
          fmi_rc%datp(irc)%qmfrommvr(id2) = fmi_rc%datp(irc)%qmfrommvr(id2) - &
                                            this%mvrbudobj%budterm(i)%flow(n) * &
                                            this%mvrterm(i)%qty(n) * &
                                            this%eqnsclfac
        end if
      end do
    end do
  end subroutine mvt_update_qmfrommvr
 
  subroutine set_fmi(this, ibudterm, fmi, set_provider)
    class(tspmvttype) :: this
    integer(I4B), intent(in) :: ibudterm
    type(tspfmitype), pointer :: fmi
    logical(LGP) :: set_provider !< set provider fmi when true, receiver fmi otherwise
    ! local
    character(len=LENMODELNAME) :: model_name
 
    fmi => null()
    if (this%gwfmodelname1 == '' .and. this%gwfmodelname2 == '') then
      fmi => this%fmi1 ! fmi2 equals fmi1 here
      return
    end if
 
    ! model name for provider is text1id1, receiver is text1id2
    if (set_provider) then
      model_name = this%mvrbudobj%budterm(ibudterm)%text1id1
    else
      model_name = this%mvrbudobj%budterm(ibudterm)%text1id2
    end if
 
    if (model_name == this%gwfmodelname1) then
      ! model 1 is the provider/receiver
      fmi => this%fmi1
    else if (model_name == this%gwfmodelname2) then
      ! model 2 is the provider
      fmi => this%fmi2
    end if
 
  end subroutine set_fmi
 
  !> @brief Extra convergence check for mover
  !<
  subroutine mvt_cc(this, kiter, iend, icnvgmod, cpak, dpak)
    ! -- dummy
    class(tspmvttype) :: this
    integer(I4B), intent(in) :: kiter
    integer(I4B), intent(in) :: iend
    integer(I4B), intent(in) :: icnvgmod
    character(len=LENPAKLOC), intent(inout) :: cpak
    real(DP), intent(inout) :: dpak
    ! -- formats
    character(len=*), parameter :: fmtmvrcnvg = &
      "(/,1x,'MOVER PACKAGE REQUIRES AT LEAST TWO OUTER ITERATIONS. CONVERGE &
      &FLAG HAS BEEN RESET TO FALSE.')"
    !
    ! -- If there are active movers, then at least 2 outers required
    if (associated(this%mvrbudobj)) then
      if (icnvgmod == 1 .and. kiter == 1) then
        dpak = dnodata
        cpak = trim(this%packName)
        write (this%iout, fmtmvrcnvg)
      end if
    end if
  end subroutine mvt_cc
 
  !> @brief Write mover terms to listing file
  !<
  subroutine mvt_bd(this)
    ! -- dummy
    class(tspmvttype) :: this
    !
    ! -- Fill the budget object
    call this%mvt_fill_budobj()
  end subroutine mvt_bd
 
  !> @brief Write mover budget terms
  !<
  subroutine mvt_ot_saveflow(this, icbcfl, ibudfl)
    ! -- modules
    use tdismodule, only: kstp, kper, delt, pertim, totim
    ! -- dummy
    class(tspmvttype) :: this
    integer(I4B), intent(in) :: icbcfl
    integer(I4B), intent(in) :: ibudfl
    ! -- locals
    integer(I4B) :: ibinun
    !
    ! -- Save the mover flows from the budobj to a mover binary file
    ibinun = 0
    if (this%ibudgetout /= 0) then
      ibinun = this%ibudgetout
    end if
    if (icbcfl == 0) ibinun = 0
    if (ibinun > 0) then
      call this%budobj%save_flows(this%dis, ibinun, kstp, kper, delt, &
                                  pertim, totim, this%iout)
    end if
  end subroutine mvt_ot_saveflow
 
  !> @brief Print mover flow table
  !<
  subroutine mvt_ot_printflow(this, icbcfl, ibudfl)
    ! -- dummy
    class(tspmvttype) :: this
    integer(I4B), intent(in) :: icbcfl
    integer(I4B), intent(in) :: ibudfl
    !
    ! -- Print the mover flow table
    if (ibudfl /= 0 .and. this%iprflow /= 0) then
      call this%mvt_print_outputtab()
    end if
  end subroutine mvt_ot_printflow
 
  !> @brief Write mover budget to listing file
  !<
  subroutine mvt_ot_bdsummary(this, ibudfl)
    ! -- modules
    use tdismodule, only: kstp, kper, delt, totim
    ! -- dummy
    class(tspmvttype) :: this
    integer(I4B), intent(in) :: ibudfl
    ! -- locals
    integer(I4B) :: i, j, n
    real(DP), allocatable, dimension(:) :: ratin, ratout
    !
    ! -- Allocate and initialize ratin/ratout
    allocate (ratin(this%maxpackages), ratout(this%maxpackages))
    do j = 1, this%maxpackages
      ratin(j) = dzero
      ratout(j) = dzero
    end do
    !
    ! -- Accumulate the rates
    do i = 1, this%maxpackages
      do j = 1, this%budobj%nbudterm
        do n = 1, this%budobj%budterm(j)%nlist
          !
          ! -- Provider is inflow to mover
          if (this%paknames(i) == this%budobj%budterm(j)%text2id1) then
            ratin(i) = ratin(i) + this%budobj%budterm(j)%flow(n)
          end if
          !
          ! -- Receiver is outflow from mover
          if (this%paknames(i) == this%budobj%budterm(j)%text2id2) then
            ratout(i) = ratout(i) + this%budobj%budterm(j)%flow(n)
          end if
        end do
      end do
    end do
    !
    ! -- Send rates to budget object
    call this%budget%reset()
    do j = 1, this%maxpackages
      call this%budget%addentry(ratin(j), ratout(j), delt, this%paknames(j))
    end do
    !
    ! -- Write the budget
    call this%budget%finalize_step(delt)
    if (ibudfl /= 0) then
      call this%budget%budget_ot(kstp, kper, this%iout)
    end if
    !
    ! -- Write budget csv
    call this%budget%writecsv(totim)
    !
    ! -- Deallocate
    deallocate (ratin, ratout)
    !
    ! -- Output mvr budget
    !    Not using budobj write_table here because it would result
    !    in a table that has one entry.  A custom table looks
    !    better here with a row for each package.
    !call this%budobj%write_budtable(kstp, kper, this%iout)
  end subroutine mvt_ot_bdsummary
 
  !> @ brief Deallocate memory
  !!
  !!  Method to deallocate memory for the package.
  !<
  subroutine mvt_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy
    class(tspmvttype) :: this
    ! -- local
    integer(I4B) :: i
    !
    ! -- Deallocate arrays if package was active
    if (this%inunit > 0) then
      !
      ! -- mvrterm
      if (associated(this%mvrterm)) then
        do i = 1, size(this%mvrterm)
          deallocate (this%mvrterm(i)%qty)
        end do
        deallocate (this%mvrterm)
      end if
      !
      ! -- Character array
      deallocate (this%paknames)
      !
      ! -- Budget object
      call this%budget%budget_da()
      deallocate (this%budget)
      !
      ! -- Budobj
      call this%budobj%budgetobject_da()
      deallocate (this%budobj)
      nullify (this%budobj)
      !
      ! -- Output table object
      if (associated(this%outputtab)) then
        call this%outputtab%table_da()
        deallocate (this%outputtab)
        nullify (this%outputtab)
      end if
    end if
    !
    ! -- Scalars
    this%fmi1 => null()
    this%fmi1 => null()
    this%mvrbudobj => null()
    call mem_deallocate(this%maxpackages)
    call mem_deallocate(this%ibudgetout)
    call mem_deallocate(this%ibudcsv)
    !
    ! -- Deallocate scalars in NumericalPackageType
    call this%NumericalPackageType%da()
  end subroutine mvt_da
 
  !> @ brief Allocate scalar variables for package
  !!
  !!  Method to allocate scalar variables for the MVT package.
  !<
  subroutine allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate, mem_setptr
    ! -- dummy
    class(tspmvttype) :: this
    !
    ! -- Allocate scalars in NumericalPackageType
    call this%NumericalPackageType%allocate_scalars()
    !
    ! -- Allocate
    call mem_allocate(this%maxpackages, 'MAXPACKAGES', this%memoryPath)
    call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath)
    call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath)
    !
    ! -- Initialize
    this%maxpackages = 0
    this%ibudgetout = 0
    this%ibudcsv = 0
  end subroutine allocate_scalars
 
  !> @brief Read mover-for-transport options block
  !<
  subroutine read_options(this)
    ! -- modules
    use openspecmodule, only: access, form
    use inputoutputmodule, only: assign_iounit, openfile
    ! -- dummy
    class(tspmvttype) :: this
    ! -- local
    character(len=LINELENGTH) :: errmsg, keyword
    character(len=MAXCHARLEN) :: fname
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
    ! -- formats
    character(len=*), parameter :: fmtflow = &
      "(4x, a, 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, &
      &/4x, 'OPENED ON UNIT: ', I0)"
    character(len=*), parameter :: fmtflow2 = &
      &"(4x, 'FLOWS WILL BE SAVED TO BUDGET FILE')"
    !
    ! -- Get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr, blockrequired=.false., &
                              supportopenclose=.true.)
    !
    ! -- Parse options block if detected
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING MVT OPTIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('SAVE_FLOWS')
          this%ipakcb = -1
          write (this%iout, fmtflow2)
        case ('PRINT_INPUT')
          this%iprpak = 1
          write (this%iout, '(4x,a)') 'MVT INPUT WILL BE PRINTED.'
        case ('PRINT_FLOWS')
          this%iprflow = 1
          write (this%iout, '(4x,a)') &
            'MVT FLOWS WILL BE PRINTED TO LISTING FILE.'
        case ('BUDGET')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            call assign_iounit(this%ibudgetout, this%inunit, "BUDGET fileout")
            call openfile(this%ibudgetout, this%iout, fname, 'DATA(BINARY)', &
                          form, access, 'REPLACE')
            write (this%iout, fmtflow) 'MVT', 'BUDGET', trim(adjustl(fname)), &
              this%ibudgetout
          else
            call store_error('Optional BUDGET keyword must &
                             &be followed by FILEOUT')
          end if
        case ('BUDGETCSV')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            call assign_iounit(this%ibudcsv, this%inunit, "BUDGETCSV fileout")
            call openfile(this%ibudcsv, this%iout, fname, 'CSV', &
                          filstat_opt='REPLACE')
            write (this%iout, fmtflow) 'MVT', 'BUDGET CSV', &
              trim(adjustl(fname)), this%ibudcsv
          else
            call store_error('Optional BUDGETCSV keyword must be followed by &
              &FILEOUT')
          end if
        case default
          write (errmsg, '(a,a)') 'Unknown MVT option: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF MVT OPTIONS'
    end if
  end subroutine read_options
 
  !> @brief Set up the budget object that stores all the mvr flows
  !<
  subroutine mvt_setup_budobj(this)
    ! -- modules
    use constantsmodule, only: lenbudtxt
    ! -- dummy
    class(tspmvttype) :: this
    ! -- local
    integer(I4B) :: nbudterm
    integer(I4B) :: ncv
    integer(I4B) :: maxlist
    integer(I4B) :: i
    integer(I4B) :: naux
    character(len=LENMODELNAME) :: modelname1, modelname2
    character(len=LENPACKAGENAME) :: packagename1, packagename2
    character(len=LENBUDTXT) :: text
    !
    ! -- Assign terms to set up the mover budget object
    nbudterm = this%mvrbudobj%nbudterm
    ncv = 0
    naux = 0
    if (this%depvartype == 'CONCENTRATION') then
      text = '        MVT-FLOW'
    else
      text = '        MVE-FLOW'
    end if
    !
    ! -- Set up budobj
    call this%budobj%budgetobject_df(ncv, nbudterm, 0, 0, bddim_opt='M')
    !
    ! -- Allocate mvrterm
    allocate (this%mvrterm(nbudterm))
    !
    ! -- Go through the water mover budget terms and set up the transport
    !    mover budget terms
    do i = 1, nbudterm
      modelname1 = this%mvrbudobj%budterm(i)%text1id1
      packagename1 = this%mvrbudobj%budterm(i)%text2id1
      modelname2 = this%mvrbudobj%budterm(i)%text1id2
      packagename2 = this%mvrbudobj%budterm(i)%text2id2
      maxlist = this%mvrbudobj%budterm(i)%maxlist
      call this%budobj%budterm(i)%initialize(text, &
                                             modelname1, &
                                             packagename1, &
                                             modelname2, &
                                             packagename2, &
                                             maxlist, .false., .false., &
                                             naux)
      this%mvrterm(i)%prov_is_m1 = (this%gwfmodelname1 == modelname1)
      allocate (this%mvrterm(i)%qty(maxlist))
      this%mvrterm(i)%qty = dnodata
    end do
  end subroutine mvt_setup_budobj
 
  !> @brief Copy mover-for-transport flow terms into this%budobj
  !<
  subroutine mvt_fill_budobj(this)
    ! -- dummy
    class(tspmvttype) :: this
    ! -- local
    type(tspfmitype), pointer :: fmi_rc
    integer(I4B) :: nbudterm
    integer(I4B) :: nlist
    integer(I4B) :: irc
    integer(I4B) :: i
    integer(I4B) :: j
    integer(I4B) :: n1, n2
    real(DP) :: cp
    real(DP) :: q
    real(DP) :: rate
    !
    ! -- Go through the water mover budget terms and set up the transport
    !    mover budget terms
    nbudterm = this%mvrbudobj%nbudterm
    do i = 1, nbudterm
      nlist = this%mvrbudobj%budterm(i)%nlist
      call this%set_fmi(i, fmi_rc, .false.)
      if (.not. associated(fmi_rc)) cycle
 
      call fmi_rc%get_package_index(this%mvrbudobj%budterm(i)%text2id2, irc)
      call this%budobj%budterm(i)%reset(nlist)
      do j = 1, nlist
        n1 = this%mvrbudobj%budterm(i)%id1(j)
        n2 = this%mvrbudobj%budterm(i)%id2(j)
        q = this%mvrbudobj%budterm(i)%flow(j)
        cp = this%mvrterm(i)%qty(j)
 
        ! calculate solute mover rate
        rate = dzero
        if (fmi_rc%iatp(irc) /= 0) then
          rate = -q * cp * this%eqnsclfac
        end if
        !
        ! -- Add the rate to the budterm
        call this%budobj%budterm(i)%update_term(n1, n2, rate)
      end do
    end do
    !
    ! --Terms are filled, now accumulate them for this time step
    call this%budobj%accumulate_terms()
  end subroutine mvt_fill_budobj
 
  !> @brief Determine max number of packages in use
  !!
  !! Scan through the gwf water mover budget object and determine the maximum
  !! number of packages and unique package names
  !<
  subroutine mvt_scan_mvrbudobj(this)
    class(tspmvttype) :: this
    integer(I4B) :: nbudterm
    integer(I4B) :: maxpackages
    integer(I4B) :: i, j
    integer(I4B) :: ipos
    logical :: found
    !
    ! -- Calculate maxpackages, which is the the square of nbudterm
    nbudterm = this%mvrbudobj%nbudterm
    do i = 1, nbudterm
      if (i * i == nbudterm) then
        maxpackages = i
        exit
      end if
    end do
    this%maxpackages = maxpackages
    !
    ! -- Allocate paknames
    if (associated(this%paknames)) deallocate (this%paknames)
    allocate (this%paknames(this%maxpackages))
    do i = 1, this%maxpackages
      this%paknames(i) = ''
    end do
    !
    ! -- Scan through mvrbudobj and create unique paknames
    ipos = 1
    do i = 1, nbudterm
      found = .false.
      do j = 1, ipos
        if (this%mvrbudobj%budterm(i)%text2id1 == this%paknames(j)) then
          found = .true.
          exit
        end if
      end do
      if (.not. found) then
        this%paknames(ipos) = this%mvrbudobj%budterm(i)%text2id1
        ipos = ipos + 1
      end if
    end do
  end subroutine mvt_scan_mvrbudobj
 
  !> @brief Set up the mover-for-transport output table
  !<
  subroutine mvt_setup_outputtab(this)
    ! -- dummy
    class(tspmvttype), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: text
    integer(I4B) :: ntabcol
    integer(I4B) :: maxrow
    integer(I4B) :: ilen
    !
    ! -- Allocate and initialize the output table
    if (this%iprflow /= 0) then
      !
      ! -- Dimension table
      ntabcol = 7
      maxrow = 0
      !
      ! -- Initialize the output table object
      title = 'TRANSPORT MOVER PACKAGE ('//trim(this%packName)// &
              ') FLOW RATES'
      call table_cr(this%outputtab, this%packName, title)
      call this%outputtab%table_df(maxrow, ntabcol, this%iout, &
                                   transient=.true.)
      text = 'NUMBER'
      call this%outputtab%initialize_column(text, 10, alignment=tabcenter)
      text = 'PROVIDER LOCATION'
      ilen = lenmodelname + lenpackagename + 1
      call this%outputtab%initialize_column(text, ilen)
      text = 'PROVIDER ID'
      call this%outputtab%initialize_column(text, 10)
      text = 'PROVIDER FLOW RATE'
      call this%outputtab%initialize_column(text, 10)
      text = 'PROVIDER TRANSPORT RATE'
      call this%outputtab%initialize_column(text, 10)
      text = 'RECEIVER LOCATION'
      ilen = lenmodelname + lenpackagename + 1
      call this%outputtab%initialize_column(text, ilen)
      text = 'RECEIVER ID'
      call this%outputtab%initialize_column(text, 10)
      !
    end if
  end subroutine mvt_setup_outputtab
 
  !> @brief Set up mover-for-transport output table
  !<
  subroutine mvt_print_outputtab(this)
    ! -- module
    use tdismodule, only: kstp, kper
    ! -- dummy
    class(tspmvttype), intent(inout) :: this
    ! -- local
    character(len=LINELENGTH) :: title
    character(len=LENMODELNAME + LENPACKAGENAME + 1) :: cloc1, cloc2
    integer(I4B) :: i
    integer(I4B) :: n
    integer(I4B) :: inum
    integer(I4B) :: ntabrows
    integer(I4B) :: nlist
    !
    ! -- Determine number of table rows
    ntabrows = 0
    do i = 1, this%budobj%nbudterm
      nlist = this%budobj%budterm(i)%nlist
      ntabrows = ntabrows + nlist
    end do
    !
    ! -- Set table kstp and kper
    call this%outputtab%set_kstpkper(kstp, kper)
    !
    ! -- Add terms and print the table
    title = 'TRANSPORT MOVER PACKAGE ('//trim(this%packName)// &
            ') FLOW RATES'
    call this%outputtab%set_title(title)
    call this%outputtab%set_maxbound(ntabrows)
    !
    ! -- Process each table row
    inum = 1
    do i = 1, this%budobj%nbudterm
      nlist = this%budobj%budterm(i)%nlist
      do n = 1, nlist
        cloc1 = trim(adjustl(this%budobj%budterm(i)%text1id1))//' '// &
                trim(adjustl(this%budobj%budterm(i)%text2id1))
        cloc2 = trim(adjustl(this%budobj%budterm(i)%text1id2))//' '// &
                trim(adjustl(this%budobj%budterm(i)%text2id2))
        call this%outputtab%add_term(inum)
        call this%outputtab%add_term(cloc1)
        call this%outputtab%add_term(this%budobj%budterm(i)%id1(n))
        call this%outputtab%add_term(-this%mvrbudobj%budterm(i)%flow(n))
        call this%outputtab%add_term(this%budobj%budterm(i)%flow(n))
        call this%outputtab%add_term(cloc2)
        call this%outputtab%add_term(this%budobj%budterm(i)%id2(n))
        inum = inum + 1
      end do
    end do
  end subroutine mvt_print_outputtab
 
end module tspmvtmodule