prt.f90

Go to the documentation of this file.

module prtmodule
  use kindmodule, only: dp, i4b, lgp
  use errorutilmodule, only: pstop
  use inputoutputmodule, only: parseline, upcase, lowcase
  use constantsmodule, only: lenftype, lenmempath, dzero, done, &
                             lenpakloc, lenpackagetype, lenbudtxt, mnormal, &
                             linelength, lenauxname
  use versionmodule, only: write_listfile_header
  use explicitmodelmodule, only: explicitmodeltype
  use basemodelmodule, only: basemodeltype
  use bndmodule, only: bndtype, addbndtolist, getbndfromlist
  use dismodule, only: distype, dis_cr
  use disvmodule, only: disvtype, disv_cr
  use disumodule, only: disutype, disu_cr
  use prtprpmodule, only: prtprptype
  use prtfmimodule, only: prtfmitype
  use prtmipmodule, only: prtmiptype
  use prtocmodule, only: prtoctype
  use budgetmodule, only: budgettype
  use listmodule, only: listtype
  use particlemodule, only: particletype, create_particle, active, term_unreleased
  use particleeventsmodule, only: particleeventdispatchertype, handle_event
  use particletracksmodule, only: particletrackstype, &
                                  particletrackfiletype, &
                                  write_particle_event
  use simmodule, only: count_errors, store_error, store_error_filename
  use memorymanagermodule, only: mem_allocate
  use methodmodule, only: methodtype, level_feature
  use hashtablemodule, only: hashtabletype, hash_table_cr, hash_table_da
  use arrayhandlersmodule, only: expandarray
 
  implicit none
 
  private
  public :: prt_cr
  public :: prtmodeltype
  public :: prt_nbasepkg, prt_nmultipkg
  public :: prt_basepkg, prt_multipkg
 
  integer(I4B), parameter :: nbditems = 2
  character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt
  data budtxt/'         STORAGE', '     TERMINATION'/
 
  !> @brief Particle tracking (PRT) model
  type, extends(explicitmodeltype) :: prtmodeltype
    type(prtfmitype), pointer :: fmi => null() ! flow model interface
    type(prtmiptype), pointer :: mip => null() ! model input package
    type(prtoctype), pointer :: oc => null() ! output control package
    type(budgettype), pointer :: budget => null() ! budget object
    class(methodtype), pointer :: method => null() ! tracking method
    type(particleeventdispatchertype), pointer :: events => null() ! event dispatcher
    class(particletrackstype), pointer :: tracks ! track output manager
    integer(I4B), pointer :: infmi => null() ! unit number FMI
    integer(I4B), pointer :: inmip => null() ! unit number MIP
    integer(I4B), pointer :: inmvt => null() ! unit number MVT
    integer(I4B), pointer :: inmst => null() ! unit number MST
    integer(I4B), pointer :: inadv => null() ! unit number ADV
    integer(I4B), pointer :: indsp => null() ! unit number DSP
    integer(I4B), pointer :: inssm => null() ! unit number SSM
    integer(I4B), pointer :: inoc => null() ! unit number OC
    integer(I4B), pointer :: nprp => null() ! number of PRP packages in the model
    real(dp), dimension(:), pointer, contiguous :: masssto => null() !< particle mass storage in cells, new value
    real(dp), dimension(:), pointer, contiguous :: massstoold => null() !< particle mass storage in cells, old value
    real(dp), dimension(:), pointer, contiguous :: ratesto => null() !< particle mass storage rate in cells
    real(dp), dimension(:), pointer, contiguous :: masstrm => null() !< particle mass terminating in cells, new value
    real(dp), dimension(:), pointer, contiguous :: ratetrm => null() !< particle mass termination rate in cells
    type(hashtabletype), pointer :: trm_ids => null() !< terminated particle ids
  contains
    ! Override BaseModelType procs
    procedure :: model_df => prt_df
    procedure :: model_ar => prt_ar
    procedure :: model_rp => prt_rp
    procedure :: model_ad => prt_ad
    procedure :: model_cq => prt_cq
    procedure :: model_bd => prt_bd
    procedure :: model_ot => prt_ot
    procedure :: model_da => prt_da
    procedure :: model_solve => prt_solve
 
    ! Private utilities
    procedure :: allocate_scalars
    procedure :: allocate_arrays
    procedure, private :: package_create
    procedure, private :: ftype_check
    procedure, private :: prt_ot_flow
    procedure, private :: prt_ot_saveflow
    procedure, private :: prt_ot_printflow
    procedure, private :: prt_ot_dv
    procedure, private :: prt_ot_bdsummary
    procedure, private :: prt_cq_budterms
    procedure, private :: create_packages
    procedure, private :: create_bndpkgs
    procedure, private :: log_namfile_options
 
  end type prtmodeltype
 
  !> @brief PRT base package array descriptors
  !!
  !! PRT6 model base package types.  Only listed packages are candidates
  !! for input and these will be loaded in the order specified.
  !<
  integer(I4B), parameter :: prt_nbasepkg = 50
  character(len=LENPACKAGETYPE), dimension(PRT_NBASEPKG) :: prt_basepkg
  data prt_basepkg/'DIS6 ', 'DISV6', 'DISU6', 'IC6  ', 'MST6 ', & !  5
                  &'ADV6 ', 'DSP6 ', 'SSM6 ', 'MIP6 ', 'CNC6 ', & ! 10
                  &'OC6  ', '     ', 'FMI6 ', '     ', 'IST6 ', & ! 15
                  &'LKT6 ', 'SFT6 ', 'MWT6 ', 'UZT6 ', 'MVT6 ', & ! 20
                  &'API6 ', '     ', '     ', '     ', '     ', & ! 25
                  25*'     '/ ! 50
 
  !> @brief PRT multi package array descriptors
  !!
  !! PRT6 model multi-instance package types.  Only listed packages are
  !! candidates for input and these will be loaded in the order specified.
  !<
  integer(I4B), parameter :: prt_nmultipkg = 50
  character(len=LENPACKAGETYPE), dimension(PRT_NMULTIPKG) :: prt_multipkg
  data prt_multipkg/'PRP6 ', '     ', '     ', '     ', '     ', & !  5
                   &45*'     '/ ! 50
 
  ! size of supported model package arrays
  integer(I4B), parameter :: niunit_prt = prt_nbasepkg + prt_nmultipkg
 
contains
 
  !> @brief Create a new particle tracking model object
  subroutine prt_cr(filename, id, modelname)
    ! modules
    use listsmodule, only: basemodellist
    use basemodelmodule, only: addbasemodeltolist
    use constantsmodule, only: linelength, lenpackagename
    use compilerversion
    use memoryhelpermodule, only: create_mem_path
    use memorymanagerextmodule, only: mem_set_value
    use simvariablesmodule, only: idm_context
    use prtnaminputmodule, only: prtnamparamfoundtype
    ! dummy
    character(len=*), intent(in) :: filename
    integer(I4B), intent(in) :: id
    character(len=*), intent(in) :: modelname
    ! local
    type(prtmodeltype), pointer :: this
    class(basemodeltype), pointer :: model
    character(len=LENMEMPATH) :: input_mempath
    character(len=LINELENGTH) :: lst_fname
    type(prtnamparamfoundtype) :: found
 
    ! Allocate a new PRT Model (this)
    allocate (this)
 
    ! Set this before any allocs in the memory manager can be done
    this%memoryPath = create_mem_path(modelname)
 
    ! Allocate event system and track output manager
    allocate (this%events)
    allocate (this%tracks)
 
    ! Allocate scalars and add model to basemodellist
    call this%allocate_scalars(modelname)
    model => this
    call addbasemodeltolist(basemodellist, model)
 
    ! Assign variables
    this%filename = filename
    this%name = modelname
    this%macronym = 'PRT'
    this%id = id
 
    ! Set input model namfile memory path
    input_mempath = create_mem_path(modelname, 'NAM', idm_context)
 
    ! Copy options from input context
    call mem_set_value(this%iprpak, 'PRINT_INPUT', input_mempath, &
                       found%print_input)
    call mem_set_value(this%iprflow, 'PRINT_FLOWS', input_mempath, &
                       found%print_flows)
    call mem_set_value(this%ipakcb, 'SAVE_FLOWS', input_mempath, &
                       found%save_flows)
 
    ! Create the list file
    call this%create_lstfile(lst_fname, filename, found%list, &
                             'PARTICLE TRACKING MODEL (PRT)')
 
    ! Activate save_flows if found
    if (found%save_flows) then
      this%ipakcb = -1
    end if
 
    ! Create model packages
    call this%create_packages()
 
    ! Create hash table for terminated particle ids
    call hash_table_cr(this%trm_ids)
 
    ! Log options
    if (this%iout > 0) then
      call this%log_namfile_options(found)
    end if
 
  end subroutine prt_cr
 
  !> @brief Define packages
  !!
  !! (1) call df routines for each package
  !! (2) set variables and pointers
  !<
  subroutine prt_df(this)
    ! modules
    use prtprpmodule, only: prtprptype
    ! dummy
    class(prtmodeltype) :: this
    ! local
    integer(I4B) :: ip
    class(bndtype), pointer :: packobj
 
    ! Define packages and utility objects
    call this%dis%dis_df()
    call this%fmi%fmi_df(this%dis, 1)
    call this%oc%oc_df()
    call this%budget%budget_df(niunit_prt, 'MASS', 'M')
 
    ! Define packages and assign iout for time series managers
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_df(this%dis%nodes, this%dis)
      packobj%TsManager%iout = this%iout
      packobj%TasManager%iout = this%iout
    end do
 
    ! Allocate model arrays
    call this%allocate_arrays()
 
  end subroutine prt_df
 
  !> @brief Allocate and read
  !!
  !! (1) allocates and reads packages part of this model,
  !! (2) allocates memory for arrays part of this model object
  !<
  subroutine prt_ar(this)
    ! modules
    use constantsmodule, only: dhnoflo
    use prtprpmodule, only: prtprptype
    use prtmipmodule, only: prtmiptype
    use methodpoolmodule, only: method_dis, method_disv
    ! dummy
    class(prtmodeltype) :: this
    ! locals
    integer(I4B) :: ip, nprp
    class(bndtype), pointer :: packobj
    class(*), pointer :: p
 
    ! Set up basic packages
    call this%fmi%fmi_ar(this%ibound)
    if (this%inmip > 0) call this%mip%mip_ar()
 
    ! Set up output control and budget
    call this%oc%oc_ar(this%dis, dhnoflo)
    call this%budget%set_ibudcsv(this%oc%ibudcsv)
 
    ! Select tracking events
    call this%tracks%select_events( &
      this%oc%trackrelease, &
      this%oc%trackfeatexit, &
      this%oc%tracktimestep, &
      this%oc%trackterminate, &
      this%oc%trackweaksink, &
      this%oc%trackusertime, &
      this%oc%tracksubfexit, &
      this%oc%trackdropped)
 
    ! Set up boundary pkgs and pkg-scoped track files
    nprp = 0
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      select type (packobj)
      type is (prtprptype)
        nprp = nprp + 1
        call packobj%prp_set_pointers(this%ibound, this%mip%izone)
        call packobj%bnd_ar()
        call packobj%bnd_ar()
        if (packobj%itrkout > 0) then
          call this%tracks%init_file( &
            packobj%itrkout, &
            iprp=nprp)
        end if
        if (packobj%itrkcsv > 0) then
          call this%tracks%init_file( &
            packobj%itrkcsv, &
            csv=.true., &
            iprp=nprp)
        end if
      class default
        call packobj%bnd_ar()
      end select
    end do
 
    ! Set up model-scoped track files
    if (this%oc%itrkout > 0) &
      call this%tracks%init_file(this%oc%itrkout)
    if (this%oc%itrkcsv > 0) &
      call this%tracks%init_file(this%oc%itrkcsv, csv=.true.)
 
    ! Set up the tracking method
    select type (dis => this%dis)
    type is (distype)
      call method_dis%init( &
        fmi=this%fmi, &
        events=this%events, &
        izone=this%mip%izone, &
        flowja=this%flowja, &
        porosity=this%mip%porosity, &
        retfactor=this%mip%retfactor, &
        tracktimes=this%oc%tracktimes)
      this%method => method_dis
    type is (disvtype)
      call method_disv%init( &
        fmi=this%fmi, &
        events=this%events, &
        izone=this%mip%izone, &
        flowja=this%flowja, &
        porosity=this%mip%porosity, &
        retfactor=this%mip%retfactor, &
        tracktimes=this%oc%tracktimes)
      this%method => method_disv
    end select
 
    ! Subscribe particle track output manager to events
    p => this%tracks
    call this%events%subscribe(write_particle_event, p)
 
    ! Set verbose tracing if requested
    if (this%oc%dump_event_trace) this%tracks%iout = 0
  end subroutine prt_ar
 
  !> @brief Read and prepare (calls package read and prepare routines)
  subroutine prt_rp(this)
    use tdismodule, only: readnewdata
    ! dummy
    class(prtmodeltype) :: this
    ! local
    class(bndtype), pointer :: packobj
    integer(I4B) :: ip
 
    ! Check with TDIS on whether or not it is time to RP
    if (.not. readnewdata) return
 
    ! Read and prepare
    if (this%inoc > 0) call this%oc%oc_rp()
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_rp()
    end do
  end subroutine prt_rp
 
  !> @brief Time step advance (calls package advance subroutines)
  subroutine prt_ad(this)
    ! modules
    use simvariablesmodule, only: isimcheck, ifailedstepretry
    ! dummy
    class(prtmodeltype) :: this
    class(bndtype), pointer :: packobj
    ! local
    integer(I4B) :: irestore
    integer(I4B) :: ip, n, i
 
    ! Reset state variable
    irestore = 0
    if (ifailedstepretry > 0) irestore = 1
 
    ! Update look-behind mass
    do n = 1, this%dis%nodes
      this%massstoold(n) = this%masssto(n)
    end do
 
    ! Advance fmi
    call this%fmi%fmi_ad()
 
    ! Advance
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ad()
      if (isimcheck > 0) then
        call packobj%bnd_ck()
      end if
    end do
    !
    ! Initialize the flowja array.  Flowja is calculated each time,
    ! even if output is suppressed.  (Flowja represents flow of particle
    ! mass and is positive into a cell.  Currently, each particle is assigned
    ! unit mass.)  Flowja is updated continually as particles are tracked
    ! over the time step and at the end of the time step.  The diagonal
    ! position of the flowja array will contain the flow residual.
    do i = 1, this%dis%nja
      this%flowja(i) = dzero
    end do
  end subroutine prt_ad
 
  !> @brief Calculate intercell flow (flowja)
  subroutine prt_cq(this, icnvg, isuppress_output)
    ! modules
    use sparsemodule, only: csr_diagsum
    use tdismodule, only: delt
    use prtprpmodule, only: prtprptype
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: icnvg
    integer(I4B), intent(in) :: isuppress_output
    ! local
    integer(I4B) :: i
    integer(I4B) :: ip
    class(bndtype), pointer :: packobj
    real(DP) :: tled
 
    ! Flowja is calculated each time, even if output is suppressed.
    ! Flowja represents flow of particle mass and is positive into a cell.
    ! Currently, each particle is assigned unit mass.
    !
    ! Reciprocal of time step size.
    tled = done / delt
    !
    ! Flowja was updated continually as particles were tracked over the
    ! time step.  At this point, flowja contains the net particle mass
    ! exchanged between cells during the time step.  To convert these to
    ! flow rates (particle mass per time), divide by the time step size.
    do i = 1, this%dis%nja
      this%flowja(i) = this%flowja(i) * tled
    end do
 
    ! Particle mass budget terms
    call this%prt_cq_budterms()
 
    ! Go through packages and call cq routines. Just a formality.
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_cq(this%masssto, this%flowja)
    end do
 
    ! Finalize calculation of flowja by adding face flows to the diagonal.
    ! This results in the flow residual being stored in the diagonal
    ! position for each cell.
    call csr_diagsum(this%dis%con%ia, this%flowja)
  end subroutine prt_cq
 
  !> @brief Calculate particle mass budget terms
  subroutine prt_cq_budterms(this)
    ! modules
    use tdismodule, only: delt
    use prtprpmodule, only: prtprptype
    ! dummy
    class(prtmodeltype) :: this
    ! local
    integer(I4B) :: ip
    class(bndtype), pointer :: packobj
    integer(I4B) :: n
    integer(I4B) :: np
    integer(I4B) :: idiag
    integer(I4B) :: iprp
    integer(I4B) :: istatus
    real(DP) :: tled
    real(DP) :: ratesto, ratetrm
    character(len=:), allocatable :: particle_id
    type(particletype), pointer :: particle
 
    call create_particle(particle)
 
    ! Reciprocal of time step size.
    tled = done / delt
 
    ! Reset mass and rate arrays
    do n = 1, this%dis%nodes
      this%masssto(n) = dzero
      this%masstrm(n) = dzero
      this%ratesto(n) = dzero
      this%ratetrm(n) = dzero
    end do
 
    ! Loop over PRP packages and assign particle mass to the
    ! appropriate budget term based on the particle status.
    iprp = 0
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      select type (packobj)
      type is (prtprptype)
        do np = 1, packobj%nparticles
          call packobj%particles%get(particle, this%id, iprp, np)
          istatus = packobj%particles%istatus(np)
          particle_id = particle%get_id()
          if (istatus == active) then
            ! calculate storage mass
            n = packobj%particles%itrdomain(np, level_feature)
            this%masssto(n) = this%masssto(n) + done ! unit mass
          else if (istatus > active) then
            if (this%trm_ids%get(particle_id) /= 0) cycle
            ! calculate terminating mass
            n = packobj%particles%itrdomain(np, level_feature)
            this%masstrm(n) = this%masstrm(n) + done ! unit mass
            call this%trm_ids%add(particle_id, 1) ! mark id terminated
          end if
        end do
      end select
    end do
 
    ! Calculate rates and update flowja
    do n = 1, this%dis%nodes
      ratesto = -(this%masssto(n) - this%massstoold(n)) * tled
      ratetrm = -this%masstrm(n) * tled
      this%ratesto(n) = ratesto
      this%ratetrm(n) = ratetrm
      idiag = this%dis%con%ia(n)
      this%flowja(idiag) = this%flowja(idiag) + ratesto
    end do
 
    call particle%destroy()
    deallocate (particle)
 
  end subroutine prt_cq_budterms
 
  !> @brief Calculate flows and budget
  !!
  !! (1) Calculate intercell flows (flowja)
  !! (2) Calculate package contributions to model budget
  !!
  !<
  subroutine prt_bd(this, icnvg, isuppress_output)
    ! modules
    use tdismodule, only: delt
    use budgetmodule, only: rate_accumulator
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: icnvg
    integer(I4B), intent(in) :: isuppress_output
    ! local
    integer(I4B) :: ip
    class(bndtype), pointer :: packobj
    real(DP) :: rin
    real(DP) :: rout
 
    ! Budget routines (start by resetting). Sole purpose of this section
    ! is to add in and outs to model budget. All ins and out for a model
    ! should be added here to this%budget. In a subsequent exchange call,
    ! exchange flows might also be added.
    call this%budget%reset()
    ! storage term
    call rate_accumulator(this%ratesto, rin, rout)
    call this%budget%addentry(rin, rout, delt, budtxt(1), &
                              isuppress_output, '             PRT')
    ! termination term
    call rate_accumulator(this%ratetrm, rin, rout)
    call this%budget%addentry(rin, rout, delt, budtxt(2), &
                              isuppress_output, '             PRT')
    ! boundary packages
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_bd(this%budget)
    end do
  end subroutine prt_bd
 
  !> @brief Print and/or save model output
  subroutine prt_ot(this)
    use tdismodule, only: tdis_ot, endofperiod
    ! dummy
    class(prtmodeltype) :: this
    ! local
    integer(I4B) :: idvsave
    integer(I4B) :: idvprint
    integer(I4B) :: icbcfl
    integer(I4B) :: icbcun
    integer(I4B) :: ibudfl
    integer(I4B) :: ipflag
 
    ! Note: particle tracking output is handled elsewhere
 
    ! Set write and print flags
    idvsave = 0
    idvprint = 0
    icbcfl = 0
    ibudfl = 0
    if (this%oc%oc_save('CONCENTRATION')) idvsave = 1
    if (this%oc%oc_print('CONCENTRATION')) idvprint = 1
    if (this%oc%oc_save('BUDGET')) icbcfl = 1
    if (this%oc%oc_print('BUDGET')) ibudfl = 1
    icbcun = this%oc%oc_save_unit('BUDGET')
 
    ! Override ibudfl and idvprint flags for nonconvergence
    ! and end of period
    ibudfl = this%oc%set_print_flag('BUDGET', 1, endofperiod)
    idvprint = this%oc%set_print_flag('CONCENTRATION', 1, endofperiod)
 
    ! Save and print flows
    call this%prt_ot_flow(icbcfl, ibudfl, icbcun)
 
    ! Save and print dependent variables
    call this%prt_ot_dv(idvsave, idvprint, ipflag)
 
    ! Print budget summaries
    call this%prt_ot_bdsummary(ibudfl, ipflag)
 
    ! Timing Output; if any dependent variables or budgets
    ! are printed, then ipflag is set to 1.
    if (ipflag == 1) call tdis_ot(this%iout)
  end subroutine prt_ot
 
  !> @brief Save flows
  subroutine prt_ot_flow(this, icbcfl, ibudfl, icbcun)
    use prtprpmodule, only: prtprptype
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: icbcfl
    integer(I4B), intent(in) :: ibudfl
    integer(I4B), intent(in) :: icbcun
    class(bndtype), pointer :: packobj
    integer(I4B) :: ip
 
    ! Save PRT flows
    call this%prt_ot_saveflow(this%dis%nja, this%flowja, icbcfl, icbcun)
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun)
    end do
 
    ! Save advanced package flows
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ot_package_flows(icbcfl=icbcfl, ibudfl=0)
    end do
 
    ! Print PRT flows
    call this%prt_ot_printflow(ibudfl, this%flowja)
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0)
    end do
 
    ! Print advanced package flows
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ot_package_flows(icbcfl=0, ibudfl=ibudfl)
    end do
  end subroutine prt_ot_flow
 
  !> @brief Save intercell flows
  subroutine prt_ot_saveflow(this, nja, flowja, icbcfl, icbcun)
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: nja
    real(DP), dimension(nja), intent(in) :: flowja
    integer(I4B), intent(in) :: icbcfl
    integer(I4B), intent(in) :: icbcun
    ! local
    integer(I4B) :: ibinun
    integer(I4B) :: naux
    real(DP), dimension(0) :: auxrow
    character(len=LENAUXNAME), dimension(0) :: auxname
    logical(LGP) :: header_written
    integer(I4B) :: i, nn
    real(DP) :: m
    integer(I4B) :: nsto, ntrm
    logical(LGP), allocatable :: msto_mask(:), mtrm_mask(:)
    integer(I4B), allocatable :: msto_nns(:), mtrm_nns(:)
    real(DP), allocatable :: msto_vals(:), mtrm_vals(:)
 
    ! Set unit number for binary output
    if (this%ipakcb < 0) then
      ibinun = icbcun
    elseif (this%ipakcb == 0) then
      ibinun = 0
    else
      ibinun = this%ipakcb
    end if
    if (icbcfl == 0) ibinun = 0
 
    ! Return if nothing to do
    if (ibinun == 0) return
 
    ! Write mass face flows
    call this%dis%record_connection_array(flowja, ibinun, this%iout)
 
    ! Write mass storage term
    naux = 0
    header_written = .false.
    msto_mask = this%masssto > dzero
    msto_vals = pack(this%masssto, msto_mask)
    msto_nns = [(i, i=1, size(this%masssto))]
    msto_nns = pack(msto_nns, msto_mask)
    nsto = size(msto_nns)
    do i = 1, nsto
      nn = msto_nns(i)
      m = msto_vals(i)
      if (.not. header_written) then
        call this%dis%record_srcdst_list_header(budtxt(1), &
                                                'PRT             ', &
                                                'PRT             ', &
                                                'PRT             ', &
                                                'STORAGE         ', &
                                                naux, auxname, ibinun, &
                                                nsto, this%iout)
        header_written = .true.
      end if
      call this%dis%record_mf6_list_entry(ibinun, nn, nn, m, &
                                          0, auxrow, &
                                          olconv2=.false.)
    end do
 
    ! Write mass termination term
    header_written = .false.
    mtrm_mask = this%masstrm > dzero
    mtrm_vals = pack(this%masstrm, mtrm_mask)
    mtrm_nns = [(i, i=1, size(this%masstrm))]
    mtrm_nns = pack(mtrm_nns, mtrm_mask)
    ntrm = size(mtrm_nns)
    do i = 1, ntrm
      nn = mtrm_nns(i)
      m = mtrm_vals(i)
      if (.not. header_written) then
        call this%dis%record_srcdst_list_header(budtxt(2), &
                                                'PRT             ', &
                                                'PRT             ', &
                                                'PRT             ', &
                                                'TERMINATION     ', &
                                                naux, auxname, ibinun, &
                                                ntrm, this%iout)
        header_written = .true.
      end if
      call this%dis%record_mf6_list_entry(ibinun, nn, nn, m, &
                                          0, auxrow, &
                                          olconv2=.false.)
    end do
 
  end subroutine prt_ot_saveflow
 
  !> @brief Print intercell flows
  subroutine prt_ot_printflow(this, ibudfl, flowja)
    ! modules
    use tdismodule, only: kper, kstp
    use constantsmodule, only: lenbigline
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: ibudfl
    real(DP), intent(inout), dimension(:) :: flowja
    ! local
    character(len=LENBIGLINE) :: line
    character(len=30) :: tempstr
    integer(I4B) :: n, ipos, m
    real(DP) :: qnm
    ! formats
    character(len=*), parameter :: fmtiprflow = &
                "(/,4x,'CALCULATED INTERCELL FLOW &
                &FOR PERIOD ', i0, ' STEP ', i0)"
 
    ! Write flowja to list file if requested
    if (ibudfl /= 0 .and. this%iprflow > 0) then
      write (this%iout, fmtiprflow) kper, kstp
      do n = 1, this%dis%nodes
        line = ''
        call this%dis%noder_to_string(n, tempstr)
        line = trim(tempstr)//':'
        do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
          m = this%dis%con%ja(ipos)
          call this%dis%noder_to_string(m, tempstr)
          line = trim(line)//' '//trim(tempstr)
          qnm = flowja(ipos)
          write (tempstr, '(1pg15.6)') qnm
          line = trim(line)//' '//trim(adjustl(tempstr))
        end do
        write (this%iout, '(a)') trim(line)
      end do
    end if
  end subroutine prt_ot_printflow
 
  !> @brief Print dependent variables
  subroutine prt_ot_dv(this, idvsave, idvprint, ipflag)
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: idvsave
    integer(I4B), intent(in) :: idvprint
    integer(I4B), intent(inout) :: ipflag
    ! local
    class(bndtype), pointer :: packobj
    integer(I4B) :: ip
 
    ! Print advanced package dependent variables
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ot_dv(idvsave, idvprint)
    end do
 
    ! save head and print head
    call this%oc%oc_ot(ipflag)
  end subroutine prt_ot_dv
 
  !> @brief Print budget summary
  subroutine prt_ot_bdsummary(this, ibudfl, ipflag)
    ! modules
    use tdismodule, only: kstp, kper, totim, delt
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: ibudfl
    integer(I4B), intent(inout) :: ipflag
    ! local
    class(bndtype), pointer :: packobj
    integer(I4B) :: ip
 
    ! Package budget summary
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_ot_bdsummary(kstp, kper, this%iout, ibudfl)
    end do
 
    ! model budget summary
    call this%budget%finalize_step(delt)
    if (ibudfl /= 0) then
      ipflag = 1
      ! model budget summary
      call this%budget%budget_ot(kstp, kper, this%iout)
    end if
 
    ! Write to budget csv
    call this%budget%writecsv(totim)
  end subroutine prt_ot_bdsummary
 
  !> @brief Deallocate
  subroutine prt_da(this)
    ! modules
    use memorymanagermodule, only: mem_deallocate
    use memorymanagerextmodule, only: memorystore_remove
    use simvariablesmodule, only: idm_context
    use methodpoolmodule, only: destroy_method_pool
    use methodcellpoolmodule, only: destroy_method_cell_pool
    use methodsubcellpoolmodule, only: destroy_method_subcell_pool
    ! dummy
    class(prtmodeltype) :: this
    ! local
    integer(I4B) :: ip
    class(bndtype), pointer :: packobj
 
    ! Deallocate idm memory
    call memorystore_remove(this%name, 'NAM', idm_context)
    call memorystore_remove(component=this%name, context=idm_context)
 
    ! Internal packages
    call this%dis%dis_da()
    call this%fmi%fmi_da()
    call this%mip%mip_da()
    call this%budget%budget_da()
    call this%oc%oc_da()
    deallocate (this%dis)
    deallocate (this%fmi)
    deallocate (this%mip)
    deallocate (this%budget)
    deallocate (this%oc)
 
    ! Method objects
    call destroy_method_subcell_pool()
    call destroy_method_cell_pool()
    call destroy_method_pool()
 
    ! Boundary packages
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_da()
      deallocate (packobj)
    end do
 
    ! Scalars
    call mem_deallocate(this%infmi)
    call mem_deallocate(this%inmip)
    call mem_deallocate(this%inadv)
    call mem_deallocate(this%indsp)
    call mem_deallocate(this%inssm)
    call mem_deallocate(this%inmst)
    call mem_deallocate(this%inmvt)
    call mem_deallocate(this%inoc)
 
    ! Arrays
    call mem_deallocate(this%masssto)
    call mem_deallocate(this%massstoold)
    call mem_deallocate(this%ratesto)
    call mem_deallocate(this%masstrm)
    call mem_deallocate(this%ratetrm)
 
    call this%tracks%destroy()
    deallocate (this%events)
    deallocate (this%tracks)
 
    call this%ExplicitModelType%model_da()
  end subroutine prt_da
 
  !> @brief Allocate memory for scalars
  subroutine allocate_scalars(this, modelname)
    ! dummy
    class(prtmodeltype) :: this
    character(len=*), intent(in) :: modelname
 
    ! allocate members from parent class
    call this%ExplicitModelType%allocate_scalars(modelname)
 
    ! allocate members that are part of model class
    call mem_allocate(this%infmi, 'INFMI', this%memoryPath)
    call mem_allocate(this%inmip, 'INMIP', this%memoryPath)
    call mem_allocate(this%inmvt, 'INMVT', this%memoryPath)
    call mem_allocate(this%inmst, 'INMST', this%memoryPath)
    call mem_allocate(this%inadv, 'INADV', this%memoryPath)
    call mem_allocate(this%indsp, 'INDSP', this%memoryPath)
    call mem_allocate(this%inssm, 'INSSM', this%memoryPath)
    call mem_allocate(this%inoc, 'INOC ', this%memoryPath)
 
    this%infmi = 0
    this%inmip = 0
    this%inmvt = 0
    this%inmst = 0
    this%inadv = 0
    this%indsp = 0
    this%inssm = 0
    this%inoc = 0
  end subroutine allocate_scalars
 
  !> @brief Allocate arrays
  subroutine allocate_arrays(this)
    use memorymanagermodule, only: mem_allocate
    class(prtmodeltype) :: this
    integer(I4B) :: n
 
    ! Allocate arrays in parent type (ibound, flowja, nja)
    call this%ExplicitModelType%allocate_arrays()
 
    ! Allocate and initialize PRT-specific arrays
    call mem_allocate(this%masssto, this%dis%nodes, &
                      'MASSSTO', this%memoryPath)
    call mem_allocate(this%massstoold, this%dis%nodes, &
                      'MASSSTOOLD', this%memoryPath)
    call mem_allocate(this%ratesto, this%dis%nodes, &
                      'RATESTO', this%memoryPath)
    call mem_allocate(this%masstrm, this%dis%nodes, &
                      'MASSTRM', this%memoryPath)
    call mem_allocate(this%ratetrm, this%dis%nodes, &
                      'RATETRM', this%memoryPath)
    do n = 1, this%dis%nodes
      this%masssto(n) = dzero
      this%massstoold(n) = dzero
      this%ratesto(n) = dzero
      this%masstrm(n) = dzero
      this%ratetrm(n) = dzero
    end do
  end subroutine allocate_arrays
 
  !> @brief Create boundary condition packages for this model
  subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
                            inunit, iout)
    ! modules
    use constantsmodule, only: linelength
    use prtprpmodule, only: prp_create
    use apimodule, only: api_create
    ! dummy
    class(prtmodeltype) :: this
    character(len=*), intent(in) :: filtyp
    character(len=LINELENGTH) :: errmsg
    integer(I4B), intent(in) :: ipakid
    integer(I4B), intent(in) :: ipaknum
    character(len=*), intent(in) :: pakname
    character(len=*), intent(in) :: mempath
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    ! local
    class(bndtype), pointer :: packobj
    class(bndtype), pointer :: packobj2
    integer(I4B) :: ip
 
    ! This part creates the package object
    select case (filtyp)
    case ('PRP6')
      call prp_create(packobj, ipakid, ipaknum, inunit, iout, &
                      this%name, pakname, mempath, this%fmi)
    case ('API6')
      call api_create(packobj, ipakid, ipaknum, inunit, iout, &
                      this%name, pakname, mempath)
    case default
      write (errmsg, *) 'Invalid package type: ', filtyp
      call store_error(errmsg, terminate=.true.)
    end select
 
    ! Packages is the bndlist that is associated with the parent model
    ! The following statement puts a pointer to this package in the ipakid
    ! position of packages.
    do ip = 1, this%bndlist%Count()
      packobj2 => getbndfromlist(this%bndlist, ip)
      if (packobj2%packName == pakname) then
        write (errmsg, '(a,a)') 'Cannot create package.  Package name  '// &
          'already exists: ', trim(pakname)
        call store_error(errmsg, terminate=.true.)
      end if
    end do
    call addbndtolist(this%bndlist, packobj)
  end subroutine package_create
 
  !> @brief Check to make sure required input files have been specified
  subroutine ftype_check(this, indis)
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), intent(in) :: indis
    ! local
    character(len=LINELENGTH) :: errmsg
 
    ! Check for DIS(u) and MIP. Stop if not present.
    if (indis == 0) then
      write (errmsg, '(1x,a)') &
        'Discretization (DIS6, DISV6, or DISU6) package not specified.'
      call store_error(errmsg)
    end if
    if (this%inmip == 0) then
      write (errmsg, '(1x,a)') &
        'Model input (MIP6) package not specified.'
      call store_error(errmsg)
    end if
 
    if (count_errors() > 0) then
      write (errmsg, '(1x,a)') 'One or more required package(s) not specified.'
      call store_error(errmsg)
      call store_error_filename(this%filename)
    end if
  end subroutine ftype_check
 
  !> @brief Solve the model
  subroutine prt_solve(this)
    use tdismodule, only: totimc, delt, endofsimulation
    use prtprpmodule, only: prtprptype
    use particlemodule, only: active, term_unreleased, term_timeout
    use particleeventmodule, only: release, terminate
    ! dummy
    class(prtmodeltype) :: this
    ! local
    integer(I4B) :: np, ip
    class(bndtype), pointer :: packobj
    type(particletype), pointer :: particle
    real(DP) :: tmax
    integer(I4B) :: iprp
 
    ! A single particle is reused in the tracking loops
    ! to avoid allocating and deallocating it each time.
    ! get() and put() retrieve and store particle state.
    call create_particle(particle)
    ! Loop over PRP packages and particles within them.
    iprp = 0
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      select type (packobj)
      type is (prtprptype)
        iprp = iprp + 1
        do np = 1, packobj%nparticles
          ! Get the particle from the store
          call packobj%particles%get(particle, this%id, iprp, np)
          ! If particle is permanently unreleased, cycle.
          ! Raise a termination event if we haven't yet.
          ! TODO: when we have generic dynamic vectors,
          ! consider terminating permanently unreleased
          ! in PRP instead of here. For now, status -8
          ! indicates the permanently unreleased event
          ! is not yet recorded, status 8 it has been.
          if (particle%istatus == (-1 * term_unreleased)) then
            call this%method%terminate(particle, status=term_unreleased)
            call packobj%particles%put(particle, np)
          end if
          if (particle%istatus > active) cycle ! Skip terminated particles
          particle%istatus = active ! Set active status in case of release
          ! If the particle was released this time step, emit a release event
          if (particle%trelease >= totimc) call this%method%release(particle)
          ! Maximum time is the end of the time step or the particle
          ! stop time, whichever comes first, unless it's the final
          ! time step and the extend option is on, in which case
          ! it's just the particle stop time.
          if (endofsimulation .and. particle%extend) then
            tmax = particle%tstop
          else
            tmax = min(totimc + delt, particle%tstop)
          end if
          ! Apply the tracking method until the maximum time.
          call this%method%apply(particle, tmax)
          ! If the particle timed out, terminate it.
          ! "Timed out" means it's still active but
          !   - it reached its stop time, or
          !   - the simulation is over.
          ! We can't detect timeout within the tracking
          ! method because the method just receives the
          ! maximum time with no context on what it is.
          ! TODO maybe think about changing that?
          if (particle%istatus <= active .and. &
              (particle%ttrack == particle%tstop .or. endofsimulation)) &
            call this%method%terminate(particle, status=term_timeout)
          ! Return the particle to the store
          call packobj%particles%put(particle, np)
        end do
      end select
    end do
    call particle%destroy()
    deallocate (particle)
  end subroutine prt_solve
 
  !> @brief Source package info and begin to process
  subroutine create_bndpkgs(this, bndpkgs, pkgtypes, pkgnames, &
                            mempaths, inunits)
    ! modules
    use constantsmodule, only: linelength, lenpackagename
    use characterstringmodule, only: characterstringtype
    ! dummy
    class(prtmodeltype) :: this
    integer(I4B), dimension(:), allocatable, intent(inout) :: bndpkgs
    type(characterstringtype), dimension(:), contiguous, &
      pointer, intent(inout) :: pkgtypes
    type(characterstringtype), dimension(:), contiguous, &
      pointer, intent(inout) :: pkgnames
    type(characterstringtype), dimension(:), contiguous, &
      pointer, intent(inout) :: mempaths
    integer(I4B), dimension(:), contiguous, &
      pointer, intent(inout) :: inunits
    ! local
    integer(I4B) :: ipakid, ipaknum
    character(len=LENFTYPE) :: pkgtype, bndptype
    character(len=LENPACKAGENAME) :: pkgname
    character(len=LENMEMPATH) :: mempath
    integer(I4B), pointer :: inunit
    integer(I4B) :: n
 
    if (allocated(bndpkgs)) then
      ! create stress packages
      ipakid = 1
      bndptype = ''
      do n = 1, size(bndpkgs)
        pkgtype = pkgtypes(bndpkgs(n))
        pkgname = pkgnames(bndpkgs(n))
        mempath = mempaths(bndpkgs(n))
        inunit => inunits(bndpkgs(n))
 
        if (bndptype /= pkgtype) then
          ipaknum = 1
          bndptype = pkgtype
        end if
 
        call this%package_create(pkgtype, ipakid, ipaknum, pkgname, mempath, &
                                 inunit, this%iout)
        ipakid = ipakid + 1
        ipaknum = ipaknum + 1
      end do
 
      ! cleanup
      deallocate (bndpkgs)
    end if
 
  end subroutine create_bndpkgs
 
  !> @brief Source package info and begin to process
  subroutine create_packages(this)
    ! modules
    use constantsmodule, only: linelength, lenpackagename
    use characterstringmodule, only: characterstringtype
    use arrayhandlersmodule, only: expandarray
    use memorymanagermodule, only: mem_setptr
    use memoryhelpermodule, only: create_mem_path
    use simvariablesmodule, only: idm_context
    use budgetmodule, only: budget_cr
    use methodpoolmodule, only: create_method_pool
    use methodcellpoolmodule, only: create_method_cell_pool
    use methodsubcellpoolmodule, only: create_method_subcell_pool
    use prtmipmodule, only: mip_cr
    use prtfmimodule, only: fmi_cr
    use prtocmodule, only: oc_cr
    ! dummy
    class(prtmodeltype) :: this
    ! local
    type(characterstringtype), dimension(:), contiguous, &
      pointer :: pkgtypes => null()
    type(characterstringtype), dimension(:), contiguous, &
      pointer :: pkgnames => null()
    type(characterstringtype), dimension(:), contiguous, &
      pointer :: mempaths => null()
    integer(I4B), dimension(:), contiguous, &
      pointer :: inunits => null()
    character(len=LENMEMPATH) :: model_mempath
    character(len=LENFTYPE) :: pkgtype
    character(len=LENPACKAGENAME) :: pkgname
    character(len=LENMEMPATH) :: mempath
    integer(I4B), pointer :: inunit
    integer(I4B), dimension(:), allocatable :: bndpkgs
    integer(I4B) :: n
    integer(I4B) :: indis = 0 ! DIS enabled flag
    character(len=LENMEMPATH) :: mempathmip = ''
    character(len=LENMEMPATH) :: mempathfmi = ''
    character(len=LENMEMPATH) :: mempathoc = ''
 
    ! set input memory paths, input/model and input/model/namfile
    model_mempath = create_mem_path(component=this%name, context=idm_context)
 
    ! set pointers to model path package info
    call mem_setptr(pkgtypes, 'PKGTYPES', model_mempath)
    call mem_setptr(pkgnames, 'PKGNAMES', model_mempath)
    call mem_setptr(mempaths, 'MEMPATHS', model_mempath)
    call mem_setptr(inunits, 'INUNITS', model_mempath)
 
    do n = 1, size(pkgtypes)
      ! attributes for this input package
      pkgtype = pkgtypes(n)
      pkgname = pkgnames(n)
      mempath = mempaths(n)
      inunit => inunits(n)
 
      ! create dis package first as it is a prerequisite for other packages
      select case (pkgtype)
      case ('DIS6')
        indis = 1
        call dis_cr(this%dis, this%name, mempath, indis, this%iout)
      case ('DISV6')
        indis = 1
        call disv_cr(this%dis, this%name, mempath, indis, this%iout)
      case ('DISU6')
        indis = 1
        call disu_cr(this%dis, this%name, mempath, indis, this%iout)
      case ('MIP6')
        this%inmip = 1
        mempathmip = mempath
      case ('FMI6')
        this%infmi = 1
        mempathfmi = mempath
      case ('OC6')
        this%inoc = 1
        mempathoc = mempath
      case ('PRP6')
        call expandarray(bndpkgs)
        bndpkgs(size(bndpkgs)) = n
      case default
        call pstop(1, "Unrecognized package type: "//pkgtype)
      end select
    end do
 
    ! Create budget manager
    call budget_cr(this%budget, this%name)
 
    ! Create tracking method pools
    call create_method_pool()
    call create_method_cell_pool()
    call create_method_subcell_pool()
 
    ! Create packages that are tied directly to model
    call mip_cr(this%mip, this%name, mempathmip, this%inmip, this%iout, this%dis)
    call fmi_cr(this%fmi, this%name, mempathfmi, this%infmi, this%iout)
    call oc_cr(this%oc, this%name, mempathoc, this%inoc, this%iout)
 
    ! Check to make sure that required ftype's have been specified
    call this%ftype_check(indis)
 
    ! Create boundary packages
    call this%create_bndpkgs(bndpkgs, pkgtypes, pkgnames, mempaths, inunits)
  end subroutine create_packages
 
  !> @brief Write model namfile options to list file
  subroutine log_namfile_options(this, found)
    use prtnaminputmodule, only: prtnamparamfoundtype
    class(prtmodeltype) :: this
    type(prtnamparamfoundtype), intent(in) :: found
 
    write (this%iout, '(1x,a)') 'NAMEFILE OPTIONS:'
 
    if (found%print_input) then
      write (this%iout, '(4x,a)') 'STRESS PACKAGE INPUT WILL BE PRINTED '// &
        'FOR ALL MODEL STRESS PACKAGES'
    end if
 
    if (found%print_flows) then
      write (this%iout, '(4x,a)') 'PACKAGE FLOWS WILL BE PRINTED '// &
        'FOR ALL MODEL PACKAGES'
    end if
 
    if (found%save_flows) then
      write (this%iout, '(4x,a)') &
        'FLOWS WILL BE SAVED TO BUDGET FILE SPECIFIED IN OUTPUT CONTROL'
    end if
 
    write (this%iout, '(1x,a)') 'END NAMEFILE OPTIONS:'
  end subroutine log_namfile_options
 
end module prtmodule