refactor(prt): use logical arrays (#2476)

Use logical arrays in PRT as enabled by #2473. This involved a few DFN tweaks to work with IDM. Also some cleanup.

Author: wpbonelli

doc/mf6io/mf6ivar/dfn/prt-prp.dfn

@@ -45,6 +45,7 @@ reader urword
 optional true
 longname whether to use local z coordinates
 description indicates that ``zrpt'' defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell.  If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user.
+mf6internal localz
 
 block options
 name extend_tracking
@@ -53,6 +54,7 @@ reader urword
 optional true
 longname whether to extend tracking beyond the end of the simulation
 description indicates that particles should be tracked beyond the end of the simulation's final time step (using that time step's flows) until particles terminate or reach a specified stop time.  By default, particles are terminated at the end of the simulation's final time step.
+mf6internal extend
 
 block options
 name track_filerecord
@@ -262,7 +264,7 @@ reader urword
 optional false
 longname force ternary tracking method
 description force use of the ternary tracking method regardless of cell type in DISV grids.
-mf6internal ifrctrn
+mf6internal frctrn
 
 block options
 name release_time_tolerance

src/Idm/prt-prpidm.f90

@@ -16,8 +16,8 @@ module PrtPrpInputModule
     logical :: iprpak = .false.
     logical :: iexmeth = .false.
     logical :: extol = .false.
-    logical :: local_z = .false.
-    logical :: extend_tracking = .false.
+    logical :: localz = .false.
+    logical :: extend = .false.
     logical :: track_filerecord = .false.
     logical :: track = .false.
     logical :: fileout = .false.
@@ -37,7 +37,7 @@ module PrtPrpInputModule
     logical :: release_timesfn = .false.
     logical :: timesfile = .false.
     logical :: idrymeth = .false.
-    logical :: ifrctrn = .false.
+    logical :: frctrn = .false.
     logical :: rttol = .false.
     logical :: rtfreq = .false.
     logical :: ichkmeth = .false.
@@ -141,13 +141,13 @@ module PrtPrpInputModule
     )
 
   type(InputParamDefinitionType), parameter :: &
-    prtprp_local_z = InputParamDefinitionType &
+    prtprp_localz = InputParamDefinitionType &
     ( &
     'PRT', & ! component
     'PRP', & ! subcomponent
     'OPTIONS', & ! block
     'LOCAL_Z', & ! tag name
-    'LOCAL_Z', & ! fortran variable
+    'LOCALZ', & ! fortran variable
     'KEYWORD', & ! type
     '', & ! shape
     'whether to use local z coordinates', & ! longname
@@ -159,13 +159,13 @@ module PrtPrpInputModule
     )
 
   type(InputParamDefinitionType), parameter :: &
-    prtprp_extend_tracking = InputParamDefinitionType &
+    prtprp_extend = InputParamDefinitionType &
     ( &
     'PRT', & ! component
     'PRP', & ! subcomponent
     'OPTIONS', & ! block
     'EXTEND_TRACKING', & ! tag name
-    'EXTEND_TRACKING', & ! fortran variable
+    'EXTEND', & ! fortran variable
     'KEYWORD', & ! type
     '', & ! shape
     'whether to extend tracking beyond the end of the simulation', & ! longname
@@ -519,13 +519,13 @@ module PrtPrpInputModule
     )
 
   type(InputParamDefinitionType), parameter :: &
-    prtprp_ifrctrn = InputParamDefinitionType &
+    prtprp_frctrn = InputParamDefinitionType &
     ( &
     'PRT', & ! component
     'PRP', & ! subcomponent
     'OPTIONS', & ! block
     'DEV_FORCETERNARY', & ! tag name
-    'IFRCTRN', & ! fortran variable
+    'FRCTRN', & ! fortran variable
     'KEYWORD', & ! type
     '', & ! shape
     'force ternary tracking method', & ! longname
@@ -903,8 +903,8 @@ module PrtPrpInputModule
     prtprp_iprpak, &
     prtprp_iexmeth, &
     prtprp_extol, &
-    prtprp_local_z, &
-    prtprp_extend_tracking, &
+    prtprp_localz, &
+    prtprp_extend, &
     prtprp_track_filerecord, &
     prtprp_track, &
     prtprp_fileout, &
@@ -924,7 +924,7 @@ module PrtPrpInputModule
     prtprp_release_timesfn, &
     prtprp_timesfile, &
     prtprp_idrymeth, &
-    prtprp_ifrctrn, &
+    prtprp_frctrn, &
     prtprp_rttol, &
     prtprp_rtfreq, &
     prtprp_ichkmeth, &

src/Model/ParticleTracking/prt-prp.f90

@@ -36,23 +36,18 @@ module PrtPrpModule
 
   !> @brief Particle release point (PRP) package
   type, extends(BndExtType) :: PrtPrpType
-    type(PrtFmiType), pointer :: fmi => null() !< flow model interface
-    type(ParticleStoreType), pointer :: particles => null() !< particle store
-    type(ParticleReleaseScheduleType), pointer :: schedule !< particle release schedule
-    integer(I4B), pointer :: nreleasepoints => null() !< number of release points
-    integer(I4B), pointer :: nreleasetimes => null() !< number of user-specified particle release times
-    integer(I4B), pointer :: nparticles => null() !< number of particles released
+    ! options
+    logical(LGP), pointer :: extend => null() !< extend tracking beyond simulation's end
+    logical(LGP), pointer :: frctrn => null() !< force ternary solution for quad grids
+    logical(LGP), pointer :: drape => null() !< whether to drape particle to topmost active cell
+    logical(LGP), pointer :: localz => null() !< compute z coordinates local to the release cell
     integer(I4B), pointer :: istopweaksink => null() !< weak sink option: 0 = no stop, 1 = stop
     integer(I4B), pointer :: istopzone => null() !< optional stop zone number: 0 = no stop zone
-    integer(I4B), pointer :: idrape => null() !< drape option: 0 = do not drape, 1 = drape to topmost active cell
     integer(I4B), pointer :: idrymeth => null() !< dry tracking method: 0 = drop, 1 = stop, 2 = stay
     integer(I4B), pointer :: itrkout => null() !< binary track file
     integer(I4B), pointer :: itrkhdr => null() !< track header file
     integer(I4B), pointer :: itrkcsv => null() !< CSV track file
     integer(I4B), pointer :: irlstls => null() !< release time file
-    integer(I4B), pointer :: ilocalz => null() !< compute z coordinates local to the cell
-    integer(I4B), pointer :: iextend => null() !< extend tracking beyond simulation's end
-    integer(I4B), pointer :: ifrctrn => null() !< force ternary solution for quad grids
     integer(I4B), pointer :: iexmeth => null() !< method for iterative solution of particle exit location and time in generalized Pollock's method
     integer(I4B), pointer :: ichkmeth => null() !< method for checking particle release coordinates are in the specified cells, 0 = none, 1 = eager
     integer(I4B), pointer :: icycwin => null() !< cycle detection window size
@@ -62,6 +57,13 @@ module PrtPrpModule
     real(DP), pointer :: offset => null() !< release time offset
     real(DP), pointer :: stoptime => null() !< stop time for all release points
     real(DP), pointer :: stoptraveltime => null() !< stop travel time for all points
+    !
+    type(PrtFmiType), pointer :: fmi => null() !< flow model interface
+    type(ParticleStoreType), pointer :: particles => null() !< particle store
+    type(ParticleReleaseScheduleType), pointer :: schedule => null() !< particle release schedule
+    integer(I4B), pointer :: nreleasepoints => null() !< number of release points
+    integer(I4B), pointer :: nreleasetimes => null() !< number of user-specified particle release times
+    integer(I4B), pointer :: nparticles => null() !< number of particles released
     integer(I4B), pointer, contiguous :: rptnode(:) => null() !< release point reduced nns
     integer(I4B), pointer, contiguous :: rptzone(:) => null() !< release point zone numbers
     real(DP), pointer, contiguous :: rptx(:) => null() !< release point x coordinates
@@ -151,14 +153,14 @@ subroutine prp_da(this)
     call this%BndExtType%bnd_da()
 
     ! Deallocate scalars
-    call mem_deallocate(this%ilocalz)
-    call mem_deallocate(this%iextend)
+    call mem_deallocate(this%localz)
+    call mem_deallocate(this%extend)
     call mem_deallocate(this%offset)
     call mem_deallocate(this%stoptime)
     call mem_deallocate(this%stoptraveltime)
     call mem_deallocate(this%istopweaksink)
     call mem_deallocate(this%istopzone)
-    call mem_deallocate(this%idrape)
+    call mem_deallocate(this%drape)
     call mem_deallocate(this%idrymeth)
     call mem_deallocate(this%nreleasepoints)
     call mem_deallocate(this%nreleasetimes)
@@ -167,7 +169,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%itrkhdr)
     call mem_deallocate(this%itrkcsv)
     call mem_deallocate(this%irlstls)
-    call mem_deallocate(this%ifrctrn)
+    call mem_deallocate(this%frctrn)
     call mem_deallocate(this%iexmeth)
     call mem_deallocate(this%ichkmeth)
     call mem_deallocate(this%icycwin)
@@ -243,14 +245,14 @@ subroutine prp_allocate_scalars(this)
     call this%BndExtType%allocate_scalars()
 
     ! Allocate scalars for this type
-    call mem_allocate(this%ilocalz, 'ILOCALZ', this%memoryPath)
-    call mem_allocate(this%iextend, 'IEXTEND', this%memoryPath)
+    call mem_allocate(this%localz, 'LOCALZ', this%memoryPath)
+    call mem_allocate(this%extend, 'EXTEND', this%memoryPath)
     call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
     call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
     call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
     call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
     call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
-    call mem_allocate(this%idrape, 'IDRAPE', this%memoryPath)
+    call mem_allocate(this%drape, 'DRAPE', this%memoryPath)
     call mem_allocate(this%idrymeth, 'IDRYMETH', this%memoryPath)
     call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
     call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
@@ -259,7 +261,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%itrkhdr, 'ITRKHDR', this%memoryPath)
     call mem_allocate(this%itrkcsv, 'ITRKCSV', this%memoryPath)
     call mem_allocate(this%irlstls, 'IRLSTLS', this%memoryPath)
-    call mem_allocate(this%ifrctrn, 'IFRCTRN', this%memoryPath)
+    call mem_allocate(this%frctrn, 'FRCTRN', this%memoryPath)
     call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
     call mem_allocate(this%ichkmeth, 'ICHKMETH', this%memoryPath)
     call mem_allocate(this%icycwin, 'ICYCWIN', this%memoryPath)
@@ -268,14 +270,14 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
 
     ! Set values
-    this%ilocalz = 0
-    this%iextend = 0
+    this%localz = .false.
+    this%extend = .false.
     this%offset = DZERO
     this%stoptime = huge(1d0)
     this%stoptraveltime = huge(1d0)
     this%istopweaksink = 0
     this%istopzone = 0
-    this%idrape = 0
+    this%drape = .false.
     this%idrymeth = 0
     this%nreleasepoints = 0
     this%nreleasetimes = 0
@@ -284,7 +286,7 @@ subroutine prp_allocate_scalars(this)
     this%itrkhdr = 0
     this%itrkcsv = 0
     this%irlstls = 0
-    this%ifrctrn = 0
+    this%frctrn = .false.
     this%iexmeth = 0
     this%ichkmeth = 1
     this%icycwin = 0
@@ -361,7 +363,7 @@ subroutine prp_ad(this)
             'Skipping negative release time (t=', t, ').'
           call store_warning(warnmsg)
           cycle
-        else if (t > totalsimtime .and. this%iextend == 0) then
+        else if (t > totalsimtime .and. .not. this%extend) then
           write (warnmsg, '(a,g0,a)') &
             'Skipping release time falling after the end of the &
             &simulation (t=', t, '). Enable EXTEND_TRACKING to &
@@ -494,7 +496,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     ! enabled, or else terminate permanently unreleased.
     if (this%ibound(ic) == 0) then
       ic_old = ic
-      if (this%idrape > 0) then
+      if (this%drape) then
         call this%dis%highest_active(ic, this%ibound)
         if (ic == ic_old .or. this%ibound(ic) == 0) then
           ! negative unreleased status signals to the
@@ -510,7 +512,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     ! Load coordinates and transform if needed
     x = this%rptx(ip)
     y = this%rpty(ip)
-    if (this%ilocalz > 0) then
+    if (this%localz) then
       top = this%fmi%dis%top(ic)
       bot = this%fmi%dis%bot(ic)
       hds = this%fmi%gwfhead(ic)
@@ -542,9 +544,9 @@ subroutine initialize_particle(this, particle, ip, trelease)
     particle%iboundary(LEVEL_FEATURE) = 0
     particle%itrdomain(LEVEL_SUBFEATURE) = 0
     particle%iboundary(LEVEL_SUBFEATURE) = 0
-    particle%ifrctrn = this%ifrctrn
+    particle%frctrn = this%frctrn
     particle%iexmeth = this%iexmeth
-    particle%iextend = this%iextend
+    particle%extend = this%extend
     particle%icycwin = this%icycwin
     particle%extol = this%extol
   end subroutine initialize_particle
@@ -694,26 +696,26 @@ subroutine prp_options(this)
                        found%istopweaksink)
     call mem_set_value(this%istopzone, 'ISTOPZONE', this%input_mempath, &
                        found%istopzone)
-    call mem_set_value(this%idrape, 'DRAPE', this%input_mempath, &
+    call mem_set_value(this%drape, 'DRAPE', this%input_mempath, &
                        found%drape)
     call mem_set_value(this%idrymeth, 'IDRYMETH', this%input_mempath, &
                        drytrack_method, found%idrymeth)
     call mem_set_value(trackfile, 'TRACKFILE', this%input_mempath, &
                        found%trackfile)
     call mem_set_value(trackcsvfile, 'TRACKCSVFILE', this%input_mempath, &
                        found%trackcsvfile)
-    call mem_set_value(this%ilocalz, 'LOCAL_Z', this%input_mempath, &
-                       found%local_z)
-    call mem_set_value(this%iextend, 'EXTEND_TRACKING', this%input_mempath, &
-                       found%extend_tracking)
+    call mem_set_value(this%localz, 'LOCALZ', this%input_mempath, &
+                       found%localz)
+    call mem_set_value(this%extend, 'EXTEND', this%input_mempath, &
+                       found%extend)
     call mem_set_value(this%extol, 'EXTOL', this%input_mempath, &
                        found%extol)
     call mem_set_value(this%rttol, 'RTTOL', this%input_mempath, &
                        found%rttol)
     call mem_set_value(this%rtfreq, 'RTFREQ', this%input_mempath, &
                        found%rtfreq)
-    call mem_set_value(this%ifrctrn, 'IFRCTRN', this%input_mempath, &
-                       found%ifrctrn)
+    call mem_set_value(this%frctrn, 'FRCTRN', this%input_mempath, &
+                       found%frctrn)
     call mem_set_value(this%iexmeth, 'IEXMETH', this%input_mempath, &
                        found%iexmeth)
     call mem_set_value(this%ichkmeth, 'ICHKMETH', this%input_mempath, &
@@ -823,7 +825,7 @@ subroutine prp_log_options(this, found, trackfile, trackcsvfile)
 
     write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
 
-    if (found%ifrctrn) then
+    if (found%frctrn) then
       write (this%iout, '(4x,a)') &
         'IF DISV, TRACKING WILL USE THE TERNARY METHOD REGARDLESS OF CELL TYPE'
     end if
@@ -948,7 +950,7 @@ subroutine prp_packagedata(this)
         this%rptnode(rptno) = noder
       end if
 
-      if (this%ilocalz > 0 .and. (zrpts(n) < 0 .or. zrpts(n) > 1)) then
+      if (this%localz .and. (zrpts(n) < 0 .or. zrpts(n) > 1)) then
         call store_error('Local z coordinate must fall in the interval [0, 1]')
         cycle
       end if

src/Model/ParticleTracking/prt.f90

@@ -961,7 +961,7 @@ subroutine prt_solve(this)
           ! 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%iextend > 0) then
+          if (endofsimulation .and. particle%extend) then
             tmax = particle%tstop
           else
             tmax = min(totimc + delt, particle%tstop)
@@ -982,7 +982,7 @@ subroutine prt_solve(this)
           ! TODO maybe think about changing that?
           if (particle%istatus <= ACTIVE .and. &
               (particle%ttrack == particle%tstop .or. &
-               (endofsimulation .and. particle%iextend == 0))) &
+               (endofsimulation .and. .not. particle%extend))) &
             call this%method%terminate(particle, status=TERM_TIMEOUT)
           ! Return the particle to the store
           call packobj%particles%put(particle, np)

src/Solution/ParticleTracker/Method/MethodDisv.f90

@@ -94,7 +94,7 @@ subroutine load_disv(this, particle, next_level, submethod)
           events=this%events, &
           tracktimes=this%tracktimes)
         submethod => method_cell_ptb
-      else if (particle%ifrctrn > 0) then
+      else if (particle%frctrn) then
         ! Force the ternary method
         call method_cell_tern%init( &
           fmi=this%fmi, &

src/Solution/ParticleTracker/Method/MethodSubcellPollock.f90

@@ -150,7 +150,7 @@ subroutine track_subcell(this, subcell, particle, tmax)
     ! ideally that would be handled at a higher scope but with extended tracking
     ! tmax is not the end of the simulation, it's just a wildly high upper bound.
     if ((statusVX .eq. 2) .and. (statusVY .eq. 2) .and. (statusVZ .eq. 2) .and. &
-        particle%iextend > 0 .and. endofsimulation) then
+        particle%extend .and. endofsimulation) then
       call this%terminate(particle, status=TERM_TIMEOUT)
       return
     end if