Method.f90

Go to the documentation of this file.

!> @brief Particle tracking strategies
module methodmodule
 
  use kindmodule, only: dp, i4b, lgp
  use errorutilmodule, only: pstop
  use subcellmodule, only: subcelltype
  use particlemodule
  use basedismodule, only: disbasetype
  use prtfmimodule, only: prtfmitype
  use cellmodule, only: celltype
  use celldefnmodule, only: celldefntype
  use trackmodule, only: trackfilecontroltype
  use timeselectmodule, only: timeselecttype
  implicit none
 
  private
  public :: methodtype, get_iatop
 
  !> @brief Base type for particle tracking methods.
  !!
  !! The PRT tracking algorithm invokes a "tracking method" for each
  !! domain. A domain can be a model, cell in a model, or subcell in
  !! a cell. Tracking proceeds recursively, delegating to a possibly
  !! arbitrary number of subdomains (currently, only the three above
  !! are recognized). A tracking method is responsible for advancing
  !! a particle through a domain, delegating to subdomains as needed
  !! depending on cell geometry (implementing the strategy pattern).
  !<
  type, abstract :: methodtype
    character(len=40), pointer, public :: type ! method name
    logical(LGP), public :: delegates ! whether the method delegates
    type(prtfmitype), pointer, public :: fmi => null() !< ptr to fmi
    class(celltype), pointer, public :: cell => null() ! ptr to a cell
    class(subcelltype), pointer, public :: subcell => null() ! ptr to a subcell
    type(trackfilecontroltype), pointer, public :: trackfilectl => null() ! ptr to track file control
    type(timeselecttype), pointer, public :: tracktimes => null() ! ptr to user-defined tracking times
    integer(I4B), dimension(:), pointer, contiguous, public :: izone => null() !< pointer to zone numbers
    real(dp), dimension(:), pointer, contiguous, public :: flowja => null() !< pointer to intercell flows
    real(dp), dimension(:), pointer, contiguous, public :: porosity => null() !< pointer to aquifer porosity
    real(dp), dimension(:), pointer, contiguous, public :: retfactor => null() !< pointer to retardation factor
  contains
    ! Implemented in all subtypes
    procedure(apply), deferred :: apply
    procedure(destroy), deferred :: destroy
    ! Overridden in subtypes that delegate
    procedure :: pass
    procedure :: load
    ! Implemented here
    procedure :: init
    procedure :: save
    procedure :: track
    procedure :: try_pass
    procedure :: update
  end type methodtype
 
  abstract interface
    subroutine apply(this, particle, tmax)
      import dp
      import methodtype
      import particletype
      class(methodtype), intent(inout) :: this
      type(particletype), pointer, intent(inout) :: particle
      real(DP), intent(in) :: tmax
    end subroutine apply
    subroutine destroy(this)
      import methodtype
      class(methodtype), intent(inout) :: this
    end subroutine destroy
  end interface
 
contains
 
  subroutine init(this, fmi, cell, subcell, trackfilectl, tracktimes, &
                  izone, flowja, porosity, retfactor)
    class(methodtype), intent(inout) :: this
    type(prtfmitype), intent(in), pointer, optional :: fmi
    class(celltype), intent(in), pointer, optional :: cell
    class(subcelltype), intent(in), pointer, optional :: subcell
    type(trackfilecontroltype), intent(in), pointer, optional :: trackfilectl
    type(timeselecttype), intent(in), pointer, optional :: tracktimes
    integer(I4B), intent(in), pointer, optional :: izone(:)
    real(DP), intent(in), pointer, optional :: flowja(:)
    real(DP), intent(in), pointer, optional :: porosity(:)
    real(DP), intent(in), pointer, optional :: retfactor(:)
 
    if (present(fmi)) this%fmi => fmi
    if (present(cell)) this%cell => cell
    if (present(subcell)) this%subcell => subcell
    if (present(trackfilectl)) this%trackfilectl => trackfilectl
    if (present(tracktimes)) this%tracktimes => tracktimes
    if (present(izone)) this%izone => izone
    if (present(flowja)) this%flowja => flowja
    if (present(porosity)) this%porosity => porosity
    if (present(retfactor)) this%retfactor => retfactor
  end subroutine init
 
  !> @brief Track particle through subdomains
  recursive subroutine track(this, particle, level, tmax)
    ! dummy
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer(I4B) :: level
    real(dp), intent(in) :: tmax
    ! local
    logical(LGP) :: advancing
    integer(I4B) :: nextlevel
    class(methodtype), pointer :: submethod
 
    advancing = .true.
    nextlevel = level + 1
    do while (advancing)
      call this%load(particle, nextlevel, submethod)
      call submethod%apply(particle, tmax)
      call this%try_pass(particle, nextlevel, advancing)
    end do
  end subroutine track
 
  !> @brief Try passing the particle to the next subdomain
  subroutine try_pass(this, particle, nextlevel, advancing)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer(I4B) :: nextlevel
    logical(LGP) :: advancing
 
    ! tracking submethod marked tracking complete?
    ! reset domain boundary flag and don't advance
    if (.not. particle%advancing) then
      particle%iboundary = 0
      advancing = .false.
    else
      ! otherwise pass particle to next subdomain
      ! and if it's on a boundary, stop advancing
      call this%pass(particle)
      if (particle%iboundary(nextlevel - 1) .ne. 0) &
        advancing = .false.
    end if
  end subroutine try_pass
 
  !> @brief Load subdomain tracking method (submethod)
  subroutine load(this, particle, next_level, submethod)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer, intent(in) :: next_level
    class(methodtype), pointer, intent(inout) :: submethod
    call pstop(1, "load must be overridden")
  end subroutine load
 
  !> @brief Pass a particle to the next subdomain, internal use only
  subroutine pass(this, particle)
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    call pstop(1, "pass must be overridden")
  end subroutine pass
 
  !> @brief Save a particle's current state.
  subroutine save(this, particle, reason)
    use tdismodule, only: kper, kstp, totimc
    ! dummy
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    integer(I4B), intent(in) :: reason
    ! local
    integer(I4B) :: per, stp
 
    per = kper
    stp = kstp
 
    ! If tracking time falls exactly on a boundary between time steps,
    ! report the previous time step for this datum. This is to follow
    ! MP7's behavior, and because the particle will have been tracked
    ! up to this instant under the previous time step's conditions, so
    ! the time step we're about to start shouldn't get "credit" for it.
    if (particle%ttrack == totimc .and. (per > 1 .or. stp > 1)) then
      if (stp > 1) then
        stp = stp - 1
      else if (per > 1) then
        per = per - 1
        stp = 1
      end if
    end if
 
    ! Save the particle's state to any registered tracking output files
    call this%trackfilectl%save(particle, kper=per, &
                                kstp=stp, reason=reason)
  end subroutine save
 
  !> @brief Update particle state and check termination conditions
  !!
  !! Update the particle's properties (e.g. advancing flag, zone number,
  !! status). If any termination conditions apply, the particle's status
  !! will be set to the appropriate termination value. If any reporting
  !! conditions apply, save particle state with the proper reason code.
  !<
  subroutine update(this, particle, cell_defn)
    ! -- dummy
    class(methodtype), intent(inout) :: this
    type(particletype), pointer, intent(inout) :: particle
    type(celldefntype), pointer, intent(inout) :: cell_defn
 
    particle%izone = cell_defn%izone
    if (cell_defn%izone .ne. 0) then
      if (particle%istopzone .eq. cell_defn%izone) then
        particle%advancing = .false.
        particle%istatus = 6
        call this%save(particle, reason=3) ! reason=3: termination
        return
      end if
    end if
    if (cell_defn%inoexitface .ne. 0) then
      particle%advancing = .false.
      particle%istatus = 5
      call this%save(particle, reason=3) ! reason=3: termination
      return
    end if
    if (cell_defn%iweaksink .ne. 0) then
      if (particle%istopweaksink .ne. 0) then
        particle%advancing = .false.
        particle%istatus = 3
        call this%save(particle, reason=3) ! reason=3: termination
        return
      else
        call this%save(particle, reason=4) ! reason=4: exited weak sink
        return
      end if
    end if
  end subroutine update
 
  !> @brief Get the index corresponding to top elevation of a cell in the grid.
  !! This is -1 if the cell is in the top layer and 1 otherwise.
  function get_iatop(ncpl, icu) result(iatop)
    integer(I4B), intent(in) :: ncpl, icu
    integer(I4B) :: iatop
 
    if (icu .le. ncpl) then
      iatop = -1
    else
      iatop = 1
    end if
  end function get_iatop
 
end module methodmodule