prtmodule Module Reference

Data Types
type	prtmodeltype
	Particle tracking (PRT) model. More...

Functions/Subroutines
subroutine, public	prt_cr (filename, id, modelname)
	Create a new particle tracking model object. More...
subroutine	prt_df (this)
	Define packages. More...
subroutine	prt_ar (this)
	Allocate and read. More...
subroutine	prt_rp (this)
	Read and prepare (calls package read and prepare routines) More...
subroutine	prt_ad (this)
	Time step advance (calls package advance subroutines) More...
subroutine	prt_cq (this, icnvg, isuppress_output)
	Calculate intercell flow (flowja) More...
subroutine	prt_cq_budterms (this)
	Calculate particle mass budget terms. More...
subroutine	prt_bd (this, icnvg, isuppress_output)
	Calculate flows and budget. More...
subroutine	prt_ot (this)
	Print and/or save model output. More...
subroutine	prt_ot_flow (this, icbcfl, ibudfl, icbcun)
	Save flows. More...
subroutine	prt_ot_saveflow (this, nja, flowja, icbcfl, icbcun)
	Save intercell flows. More...
subroutine	prt_ot_printflow (this, ibudfl, flowja)
	Print intercell flows. More...
subroutine	prt_ot_dv (this, idvsave, idvprint, ipflag)
	Print dependent variables. More...
subroutine	prt_ot_bdsummary (this, ibudfl, ipflag)
	Print budget summary. More...
subroutine	prt_da (this)
	Deallocate. More...
subroutine	allocate_scalars (this, modelname)
	Allocate memory for scalars. More...
subroutine	allocate_arrays (this)
	Allocate arrays. More...
subroutine	package_create (this, filtyp, ipakid, ipaknum, pakname, mempath, inunit, iout)
	Create boundary condition packages for this model. More...
subroutine	ftype_check (this, indis)
	Check to make sure required input files have been specified. More...
subroutine	prt_solve (this)
	Solve the model. More...
subroutine	create_bndpkgs (this, bndpkgs, pkgtypes, pkgnames, mempaths, inunits)
	Source package info and begin to process. More...
subroutine	create_packages (this)
	Source package info and begin to process. More...
subroutine	create_exg_prp (this)
	Create an exchange PRP package for particles entering this model from other model. More...
subroutine	log_namfile_options (this, found)
	Write model namfile options to list file. More...

Variables
integer(i4b), parameter	nbditems = 2
character(len=lenbudtxt), dimension(nbditems)	budtxt
integer(i4b), parameter, public	prt_nbasepkg = 50
	PRT base package array descriptors. More...
character(len=lenpackagetype), dimension(prt_nbasepkg), public	prt_basepkg
integer(i4b), parameter, public	prt_nmultipkg = 50
	PRT multi package array descriptors. More...
character(len=lenpackagetype), dimension(prt_nmultipkg), public	prt_multipkg
integer(i4b), parameter	niunit_prt = PRT_NBASEPKG + PRT_NMULTIPKG

Function/Subroutine Documentation

allocate_arrays()

subroutine prtmodule::allocate_arrays ( class(prtmodeltype)  this )

private

Definition at line 920 of file prt.f90.

    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

allocate_scalars()

subroutine prtmodule::allocate_scalars	(	class(prtmodeltype)	this,
		character(len=*), intent(in)	modelname
	)

Definition at line 891 of file prt.f90.

    ! 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

create_bndpkgs()

subroutine prtmodule::create_bndpkgs	(	class(prtmodeltype)	this,
		integer(i4b), dimension(:), intent(inout), allocatable	bndpkgs,
		type(characterstringtype), dimension(:), intent(inout), pointer, contiguous	pkgtypes,
		type(characterstringtype), dimension(:), intent(inout), pointer, contiguous	pkgnames,
		type(characterstringtype), dimension(:), intent(inout), pointer, contiguous	mempaths,
		integer(i4b), dimension(:), intent(inout), pointer, contiguous	inunits
	)

Definition at line 1099 of file prt.f90.

    ! 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
 

create_exg_prp()

subroutine prtmodule::create_exg_prp

(

class(prtmodeltype)

this

)

Definition at line 1252 of file prt.f90.

    class(PrtModelType) :: this
    ! local
    class(BndType), pointer :: packobj
    character(len=LENPACKAGENAME) :: exgprp_name
 
    exgprp_name = 'EXGPRP'
 
    call prp_create(packobj, &
                    id=0, &
                    ibcnum=0, &
                    inunit=-1, &
                    iout=this%iout, &
                    namemodel=this%name, &
                    pakname=exgprp_name, &
                    fmi=this%fmi)
    call addbndtolist(this%bndlist, packobj)

Here is the call graph for this function:

create_packages()

subroutine prtmodule::create_packages

(

class(prtmodeltype)

this

)

Definition at line 1151 of file prt.f90.

    ! 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 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)
 
    ! determine which packages we have. create
    ! dis up front as the others depend on it.
    do n = 1, size(pkgtypes)
      pkgtype = pkgtypes(n)
      pkgname = pkgnames(n)
      mempath = mempaths(n)
      inunit => inunits(n)
 
      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 methods
    call create_method_dis(this%method_dis)
    call create_method_disv(this%method_disv)
 
    ! Create non-boundary packages
    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 required input files
    call this%ftype_check(indis)
 
    ! Create boundary packages
    call this%create_bndpkgs(bndpkgs, pkgtypes, pkgnames, mempaths, inunits)
    call this%create_exg_prp()

Here is the call graph for this function:

ftype_check()

subroutine prtmodule::ftype_check	(	class(prtmodeltype)	this,
		integer(i4b), intent(in)	indis
	)

Definition at line 997 of file prt.f90.

    ! 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

Here is the call graph for this function:

log_namfile_options()

subroutine prtmodule::log_namfile_options ( class(prtmodeltype)  this, type(prtnamparamfoundtype), intent(in)  found  )

private

Definition at line 1272 of file prt.f90.

    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:'

package_create()

subroutine prtmodule::package_create	(	class(prtmodeltype)	this,
		character(len=*), intent(in)	filtyp,
		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
	)

Definition at line 949 of file prt.f90.

    ! modules
    use constantsmodule, only: linelength
    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, this%fmi, mempath)
    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)

Here is the call graph for this function:

prt_ad()

subroutine prtmodule::prt_ad

(

class(prtmodeltype)

this

)

Definition at line 361 of file prt.f90.

    ! 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

Here is the call graph for this function:

prt_ar()

subroutine prtmodule::prt_ar

(

class(prtmodeltype)

this

)

(1) allocates and reads packages part of this model, (2) allocates memory for arrays part of this model object

Definition at line 245 of file prt.f90.

    ! modules
    use constantsmodule, only: dhnoflo
    use prtprpmodule, only: prtprptype
    use prtmipmodule, only: prtmiptype
    ! 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.)
 
    ! Initialize and select the tracking method based on discretization
    select type (dis => this%dis)
    type is (distype)
      call this%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 => this%method_dis
    type is (disvtype)
      call this%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 => this%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

Here is the call graph for this function:

prt_bd()

subroutine prtmodule::prt_bd	(	class(prtmodeltype)	this,
		integer(i4b), intent(in)	icnvg,
		integer(i4b), intent(in)	isuppress_output
	)

(1) Calculate intercell flows (flowja) (2) Calculate package contributions to model budget

Definition at line 529 of file prt.f90.

    ! 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

Here is the call graph for this function:

prt_cq()

subroutine prtmodule::prt_cq	(	class(prtmodeltype)	this,
		integer(i4b), intent(in)	icnvg,
		integer(i4b), intent(in)	isuppress_output
	)

Definition at line 404 of file prt.f90.

    ! 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)

Here is the call graph for this function:

prt_cq_budterms()

subroutine prtmodule::prt_cq_budterms

(

class(prtmodeltype)

this

)

Definition at line 450 of file prt.f90.

    ! 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)
 

Here is the call graph for this function:

prt_cr()

subroutine, public prtmodule::prt_cr	(	character(len=*), intent(in)	filename,
		integer(i4b), intent(in)	id,
		character(len=*), intent(in)	modelname
	)

Definition at line 132 of file prt.f90.

    ! 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
 

Here is the call graph for this function:

Here is the caller graph for this function:

prt_da()

subroutine prtmodule::prt_da

(

class(prtmodeltype)

this

)

Definition at line 826 of file prt.f90.

    ! modules
    use memorymanagermodule, only: mem_deallocate
    use memorymanagerextmodule, only: memorystore_remove
    use simvariablesmodule, only: idm_context
    ! 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 this%method_dis%deallocate()
    deallocate (this%method_dis)
    call this%method_disv%deallocate()
    deallocate (this%method_disv)
 
    ! 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()

Here is the call graph for this function:

prt_df()

subroutine prtmodule::prt_df

(

class(prtmodeltype)

this

)

(1) call df routines for each package (2) set variables and pointers

Definition at line 212 of file prt.f90.

    ! 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()
 

Here is the call graph for this function:

prt_ot()

subroutine prtmodule::prt_ot

(

class(prtmodeltype)

this

)

Definition at line 564 of file prt.f90.

    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)

Here is the call graph for this function:

prt_ot_bdsummary()

subroutine prtmodule::prt_ot_bdsummary ( class(prtmodeltype)  this, integer(i4b), intent(in)  ibudfl, integer(i4b), intent(inout)  ipflag  )

private

Definition at line 796 of file prt.f90.

    ! 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)

Here is the call graph for this function:

prt_ot_dv()

subroutine prtmodule::prt_ot_dv	(	class(prtmodeltype)	this,
		integer(i4b), intent(in)	idvsave,
		integer(i4b), intent(in)	idvprint,
		integer(i4b), intent(inout)	ipflag
	)

Definition at line 775 of file prt.f90.

    ! 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)

Here is the call graph for this function:

prt_ot_flow()

subroutine prtmodule::prt_ot_flow	(	class(prtmodeltype)	this,
		integer(i4b), intent(in)	icbcfl,
		integer(i4b), intent(in)	ibudfl,
		integer(i4b), intent(in)	icbcun
	)

Definition at line 609 of file prt.f90.

    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

Here is the call graph for this function:

prt_ot_printflow()

subroutine prtmodule::prt_ot_printflow ( class(prtmodeltype)  this, integer(i4b), intent(in)  ibudfl, real(dp), dimension(:), intent(inout)  flowja  )

private

Definition at line 736 of file prt.f90.

    ! 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

prt_ot_saveflow()

subroutine prtmodule::prt_ot_saveflow	(	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
	)

Definition at line 646 of file prt.f90.

    ! 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
 

prt_rp()

subroutine prtmodule::prt_rp

(

class(prtmodeltype)

this

)

Definition at line 341 of file prt.f90.

    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

Here is the call graph for this function:

prt_solve()

subroutine prtmodule::prt_solve ( class(prtmodeltype)  this )

private

Definition at line 1024 of file prt.f90.

    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)

Here is the call graph for this function:

Variable Documentation

budtxt

character(len=lenbudtxt), dimension(nbditems) prtmodule::budtxt

private

Definition at line 43 of file prt.f90.

  character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt

nbditems

integer(i4b), parameter prtmodule::nbditems = 2

private

Definition at line 42 of file prt.f90.

  integer(I4B), parameter :: NBDITEMS = 2

niunit_prt

integer(i4b), parameter prtmodule::niunit_prt = PRT_NBASEPKG + PRT_NMULTIPKG

private

Definition at line 127 of file prt.f90.

  integer(I4B), parameter :: NIUNIT_PRT = prt_nbasepkg + prt_nmultipkg

prt_basepkg

character(len=lenpackagetype), dimension(prt_nbasepkg), public prtmodule::prt_basepkg

Definition at line 108 of file prt.f90.

  character(len=LENPACKAGETYPE), dimension(PRT_NBASEPKG) :: PRT_BASEPKG

prt_multipkg

character(len=lenpackagetype), dimension(prt_nmultipkg), public prtmodule::prt_multipkg

Definition at line 122 of file prt.f90.

  character(len=LENPACKAGETYPE), dimension(PRT_NMULTIPKG) :: PRT_MULTIPKG

prt_nbasepkg

integer(i4b), parameter, public prtmodule::prt_nbasepkg = 50

PRT6 model base package types. Only listed packages are candidates for input and these will be loaded in the order specified.

Definition at line 107 of file prt.f90.

  integer(I4B), parameter :: PRT_NBASEPKG = 50

prt_nmultipkg

integer(i4b), parameter, public prtmodule::prt_nmultipkg = 50

PRT6 model multi-instance package types. Only listed packages are candidates for input and these will be loaded in the order specified.

Definition at line 121 of file prt.f90.

  integer(I4B), parameter :: PRT_NMULTIPKG = 50