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_sto (this)
	Calculate particle mass storage. 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_obs (this)
	Calculate and save observations. 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 non-allocatable members. 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	log_namfile_options (this, found)
	Write model namfile options to list file. More...

Variables
integer(i4b), parameter	nbditems = 1
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 822 of file prt.f90.

    use memorymanagermodule, only: mem_allocate
    class(PrtModelType) :: this
    integer(I4B) :: n
 
    ! -- Allocate arrays in parent type
    this%nja = this%dis%nja
    call this%NumericalModelType%allocate_arrays()
 
    ! -- Allocate and initialize 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)
    ! -- explicit model, so these must be manually allocated
    call mem_allocate(this%x, this%dis%nodes, 'X', this%memoryPath)
    call mem_allocate(this%rhs, this%dis%nodes, 'RHS', this%memoryPath)
    call mem_allocate(this%ibound, this%dis%nodes, 'IBOUND', this%memoryPath)
    do n = 1, this%dis%nodes
      this%masssto(n) = dzero
      this%massstoold(n) = dzero
      this%ratesto(n) = dzero
      this%x(n) = dzero
      this%rhs(n) = dzero
      this%ibound(n) = 1
    end do

allocate_scalars()

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

Definition at line 789 of file prt.f90.

    ! -- dummy
    class(PrtModelType) :: this
    character(len=*), intent(in) :: modelname
 
    ! -- allocate members from parent class
    call this%NumericalModelType%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)
    call mem_allocate(this%inobs, 'INOBS', this%memoryPath)
    call mem_allocate(this%nprp, 'NPRP', this%memoryPath) ! kluge?
 
    this%infmi = 0
    this%inmip = 0
    this%inmvt = 0
    this%inmst = 0
    this%inadv = 0
    this%indsp = 0
    this%inssm = 0
    this%inoc = 0
    this%inobs = 0
    this%nprp = 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 1008 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_packages()

subroutine prtmodule::create_packages

(

class(prtmodeltype)

this

)

Definition at line 1062 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 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
    use prtobsmodule, only: prt_obs_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 = ''
 
    ! -- 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 = inunit
      case ('OC6')
        this%inoc = inunit
      case ('OBS6')
        this%inobs = inunit
      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, this%infmi, this%iout)
    call oc_cr(this%oc, this%name, this%inoc, this%iout)
    call prt_obs_cr(this%obs, this%inobs)
 
    ! -- 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)

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 903 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(gwfnamparamfoundtype), intent(in)	found
	)

Definition at line 1165 of file prt.f90.

    use gwfnaminputmodule, only: gwfnamparamfoundtype
    class(PrtModelType) :: this
    type(GwfNamParamFoundType), intent(in) :: found
 
    write (this%iout, '(1x,a)') 'NAMEFILE OPTIONS:'
 
    if (found%newton) then
      write (this%iout, '(4x,a)') &
        'NEWTON-RAPHSON method enabled for the model.'
      if (found%under_relaxation) then
        write (this%iout, '(4x,a,a)') &
          'NEWTON-RAPHSON UNDER-RELAXATION based on the bottom ', &
          'elevation of the model will be applied to the model.'
      end if
    end if
 
    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 853 of file prt.f90.

    ! -- 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')
      this%nprp = this%nprp + 1
      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)
    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 325 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
 
    ! -- Copy masssto into massstoold
    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
 
    ! -- Push simulated values to preceding time/subtime step
    call this%obs%obs_ad()
    !
    ! -- 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 232 of file prt.f90.

    ! -- 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
    class(BndType), pointer :: packobj
 
    ! -- Allocate and read modules attached to model
    call this%fmi%fmi_ar(this%ibound)
    if (this%inmip > 0) call this%mip%mip_ar()
 
    ! -- set up output control
    call this%oc%oc_ar(this%masssto, this%dis, dhnoflo)
    call this%budget%set_ibudcsv(this%oc%ibudcsv)
 
    ! -- Package input files now open, so allocate and read
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      select type (packobj)
      type is (prtprptype)
        call packobj%prp_set_pointers(this%ibound, this%mip%izone, &
                                      this%trackfilectl)
      end select
      ! -- Read and allocate package
      call packobj%bnd_ar()
    end do
 
    ! -- Initialize tracking method
    select type (dis => this%dis)
    type is (distype)
      call method_dis%init( &
        fmi=this%fmi, &
        trackfilectl=this%trackfilectl, &
        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, &
        trackfilectl=this%trackfilectl, &
        izone=this%mip%izone, &
        flowja=this%flowja, &
        porosity=this%mip%porosity, &
        retfactor=this%mip%retfactor, &
        tracktimes=this%oc%tracktimes)
      this%method => method_disv
      method_disv%zeromethod = this%mip%zeromethod
    end select
 
    ! -- Initialize track output files and reporting options
    if (this%oc%itrkout > 0) &
      call this%trackfilectl%init_track_file(this%oc%itrkout)
    if (this%oc%itrkcsv > 0) &
      call this%trackfilectl%init_track_file(this%oc%itrkcsv, csv=.true.)
    call this%trackfilectl%set_track_events( &
      this%oc%trackrelease, &
      this%oc%tracktransit, &
      this%oc%tracktimestep, &
      this%oc%trackterminate, &
      this%oc%trackweaksink, &
      this%oc%trackusertime)

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 471 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()
    call rate_accumulator(this%ratesto, rin, rout)
    call this%budget%addentry(rin, rout, delt, budtxt(1), &
                              isuppress_output, '             PRT')
    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 371 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 storage
    call this%prt_cq_sto()
 
    ! -- 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_sto()

subroutine prtmodule::prt_cq_sto

(

class(prtmodeltype)

this

)

Definition at line 417 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) :: istatus
    real(DP) :: tled
    real(DP) :: rate
 
    ! -- Reciprocal of time step size.
    tled = done / delt
 
    ! -- Particle mass storage rate
    do n = 1, this%dis%nodes
      this%masssto(n) = dzero
      this%ratesto(n) = dzero
    end do
    do ip = 1, this%bndlist%Count() ! kluge note: could accumulate masssto on the fly in prt_solve instead
      packobj => getbndfromlist(this%bndlist, ip)
      select type (packobj)
      type is (prtprptype)
        do np = 1, packobj%nparticles
          istatus = packobj%particles%istatus(np)
          ! refine these conditions as necessary
          ! (status 8 is permanently unreleased)
          if ((istatus > 0) .and. (istatus /= 8)) then
            n = packobj%particles%idomain(np, 2)
            ! -- Each particle currently assigned unit mass
            this%masssto(n) = this%masssto(n) + done
          end if
        end do
      end select
    end do
    do n = 1, this%dis%nodes ! kluge note: set rate to zero and skip inactive nodes?
      rate = -(this%masssto(n) - this%massstoold(n)) * tled
      this%ratesto(n) = rate
      idiag = this%dis%con%ia(n)
      this%flowja(idiag) = this%flowja(idiag) + rate
    end do

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 122 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 gwfnaminputmodule, only: gwfnamparamfoundtype
    ! -- 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(GwfNamParamFoundType) :: 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 track control object
    allocate (this%trackfilectl)
 
    ! -- 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
 
    ! -- Log options
    if (this%iout > 0) then
      call this%log_namfile_options(found)
    end if
 
    ! -- Create model packages
    call this%create_packages()

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 720 of file prt.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    use memorymanagerextmodule, only: memorylist_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 memorylist_remove(this%name, 'NAM', idm_context)
    call memorylist_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()
    call this%obs%obs_da()
    deallocate (this%dis)
    deallocate (this%fmi)
    deallocate (this%mip)
    deallocate (this%budget)
    deallocate (this%oc)
    deallocate (this%obs)
 
    ! -- 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)
    call mem_deallocate(this%inobs)
    call mem_deallocate(this%nprp)
 
    ! -- Arrays
    call mem_deallocate(this%masssto)
    call mem_deallocate(this%massstoold)
    call mem_deallocate(this%ratesto)
 
    ! -- Track file control
    deallocate (this%trackfilectl)
 
    ! -- Parent type
    call this%NumericalModelType%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 197 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()
 
    ! -- Store information needed for observations
    call this%obs%obs_df(this%iout, this%name, 'PRT', this%dis)

Here is the call graph for this function:

prt_ot()

subroutine prtmodule::prt_ot

(

class(prtmodeltype)

this

)

Definition at line 500 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)
 
    ! -- Calculate and save observations
    call this%prt_ot_obs()
 
    ! -- 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 690 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 669 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  )

private

Definition at line 566 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_obs()

subroutine prtmodule::prt_ot_obs

(

class(prtmodeltype)

this

)

Definition at line 548 of file prt.f90.

    class(PrtModelType) :: this
    class(BndType), pointer :: packobj
    integer(I4B) :: ip
 
    ! -- Calculate and save observations
    call this%obs%obs_bd()
    call this%obs%obs_ot()
 
    ! -- Calculate and save package obserations
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      call packobj%bnd_bd_obs()
      call packobj%bnd_ot_obs()
    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 630 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 603 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
 
    ! -- 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
 
    ! -- Write the face flows if requested
    if (ibinun /= 0) then
      call this%dis%record_connection_array(flowja, ibinun, this%iout)
    end if

prt_rp()

subroutine prtmodule::prt_rp

(

class(prtmodeltype)

this

)

Definition at line 304 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()
      call packobj%bnd_rp_obs()
    end do

Here is the call graph for this function:

prt_solve()

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

private

Definition at line 930 of file prt.f90.

    ! -- modules
    use tdismodule, only: kper, kstp, totimc, nper, nstp, delt
    use prtprpmodule, only: prtprptype
    ! -- dummy variables
    class(PrtModelType) :: this
    ! -- local variables
    integer(I4B) :: np, ip
    class(BndType), pointer :: packobj
    type(ParticleType), pointer :: particle
    real(DP) :: tmax
    integer(I4B) :: iprp
 
    ! -- Initialize particle
    call create_particle(particle)
 
    ! -- Loop over PRP packages
    iprp = 0
    do ip = 1, this%bndlist%Count()
      packobj => getbndfromlist(this%bndlist, ip)
      select type (packobj)
      type is (prtprptype)
        ! -- Update PRP index
        iprp = iprp + 1
 
        ! -- Initialize PRP-specific track files, if enabled
        if (packobj%itrkout > 0) then
          call this%trackfilectl%init_track_file( &
            packobj%itrkout, &
            iprp=iprp)
        end if
        if (packobj%itrkcsv > 0) then
          call this%trackfilectl%init_track_file( &
            packobj%itrkcsv, &
            csv=.true., &
            iprp=iprp)
        end if
 
        ! -- Loop over particles in package
        do np = 1, packobj%nparticles
          ! -- Load particle from storage
          call particle%load_from_store(packobj%particles, &
                                        this%id, iprp, np)
 
          ! -- If particle is permanently unreleased, record its initial/terminal state
          if (particle%istatus == 8) &
            call this%method%save(particle, reason=3) ! reason=3: termination
 
          ! -- If particle is inactive or not yet to be released, cycle
          if (particle%istatus > 1) cycle
 
          ! -- If particle released this time step, record its initial state
          particle%istatus = 1
          if (particle%trelease >= totimc) &
            call this%method%save(particle, reason=0) ! reason=0: release
 
          ! -- Maximum time is end of time step unless this is the last
          !    time step in the simulation, which case it's the particle
          !    stop time.
          tmax = totimc + delt
          if (nper == kper .and. nstp(kper) == kstp) &
            tmax = particle%tstop
 
          ! -- Get and apply the tracking method
          call this%method%apply(particle, tmax)
 
          ! -- Update particle storage
          call packobj%particles%load_from_particle(particle, np)
        end do
      end select
    end do
 
    ! -- Destroy particle
    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 37 of file prt.f90.

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

nbditems

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

private

Definition at line 36 of file prt.f90.

  integer(I4B), parameter :: NBDITEMS = 1

niunit_prt

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

private

Definition at line 117 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 98 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 112 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 97 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 111 of file prt.f90.

  integer(I4B), parameter :: PRT_NMULTIPKG = 50