This subroutine consists partly of code written by David W. Pollock of the USGS for MODPATH 7. PRT's authors take responsibility for its application in this context and for any modifications or errors.
90 class(MethodSubcellPollockType),
intent(inout) :: this
91 class(SubcellRectType),
intent(in) :: subcell
92 type(ParticleType),
pointer,
intent(inout) :: particle
93 real(DP),
intent(in) :: tmax
112 integer(I4B) :: statusVX
113 integer(I4B) :: statusVY
114 integer(I4B) :: statusVZ
119 integer(I4B) :: exitFace
120 integer(I4B) :: reason
121 integer(I4B) :: tslice(2)
126 initialx = particle%x / subcell%dx
127 initialy = particle%y / subcell%dy
128 initialz = particle%z / subcell%dz
131 statusvx = calculate_dt(subcell%vx1, subcell%vx2, subcell%dx, &
132 initialx, vx, dvxdx, dtexitx)
133 statusvy = calculate_dt(subcell%vy1, subcell%vy2, subcell%dy, &
134 initialy, vy, dvydy, dtexity)
135 statusvz = calculate_dt(subcell%vz1, subcell%vz2, subcell%dz, &
136 initialz, vz, dvzdz, dtexitz)
139 if ((statusvx .eq. 3) .and. (statusvy .eq. 3) .and. (statusvz .eq. 3))
then
140 print *,
"Subcell with no exit face:", &
142 "cell", particle%idomain(2)
149 if ((statusvx .lt. 2) .or. (statusvy .lt. 2) .or. (statusvz .lt. 2))
then
152 if (vx .lt. 0d0)
then
154 else if (vx .gt. 0)
then
158 if (dtexity .lt. dtexit)
then
160 if (vy .lt. 0d0)
then
162 else if (vy .gt. 0d0)
then
167 if (dtexitz .lt. dtexit)
then
169 if (vz .lt. 0d0)
then
171 else if (vz .gt. 0d0)
then
179 texit = particle%ttrack + dtexit
186 call this%tracktimes%try_advance()
187 tslice = this%tracktimes%selection
189 if (all(tslice > 0))
then
190 do i = tslice(1), tslice(2)
191 t = this%tracktimes%times(i)
192 if (t < particle%ttrack .or. t >= texit .or. t >= tmax) cycle
194 x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
195 dt, initialx, subcell%dx, statusvx == 1)
196 y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
197 dt, initialy, subcell%dy, statusvy == 1)
198 z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
199 dt, initialz, subcell%dz, statusvz == 1)
200 particle%x = x * subcell%dx
201 particle%y = y * subcell%dy
202 particle%z = z * subcell%dz
205 call this%save(particle, reason=5)
209 if (texit .gt. tmax)
then
215 x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
216 dt, initialx, subcell%dx, statusvx == 1)
217 y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
218 dt, initialy, subcell%dy, statusvy == 1)
219 z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
220 dt, initialz, subcell%dz, statusvz == 1)
223 particle%advancing = .false.
231 if ((exitface .eq. 1) .or. (exitface .eq. 2))
then
233 y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
234 dt, initialy, subcell%dy, statusvy == 1)
235 z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, &
236 dt, initialz, subcell%dz, statusvz == 1)
237 if (exitface .eq. 2) x = 1.0d0
238 else if ((exitface .eq. 3) .or. (exitface .eq. 4))
then
239 x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, dt, &
240 initialx, subcell%dx, statusvx == 1)
242 z = new_x(vz, dvzdz, subcell%vz1, subcell%vz2, dt, &
243 initialz, subcell%dz, statusvz == 1)
244 if (exitface .eq. 4) y = 1.0d0
245 else if ((exitface .eq. 5) .or. (exitface .eq. 6))
then
246 x = new_x(vx, dvxdx, subcell%vx1, subcell%vx2, &
247 dt, initialx, subcell%dx, statusvx == 1)
248 y = new_x(vy, dvydy, subcell%vy1, subcell%vy2, &
249 dt, initialy, subcell%dy, statusvy == 1)
251 if (exitface .eq. 6) z = 1.0d0
253 print *,
"programmer error, invalid exit face", exitface
261 particle%x = x * subcell%dx
262 particle%y = y * subcell%dy
263 particle%z = z * subcell%dz
265 particle%iboundary(3) = exitface
268 if (reason >= 0)
call this%save(particle, reason=reason)
pure character(len=lenmempath) function, public get_particle_id(particle)
Return the particle's composite ID.