Releases may be scheduled via a global RELEASETIME, or within a stress period via ALL, FIRST, FREQUENCY or STEPS (with optional FRACTION). If no release option is specified, a single release is conducted at the first moment of the first time step of the first stress period.
308 class(PrtPrpType) :: this
310 character(len=LINELENGTH) :: errmsg
311 integer(I4B) :: ic, icu, nps, nts, nrel, &
312 nreleasets, np, irow, icol, ilay, icpl
313 real(DP) :: x, y, z, trelease, tend, top, bot, hds
314 real(DP),
allocatable :: polyverts(:, :)
315 type(ParticleType),
pointer :: particle
322 (this%rlsfirst .and.
kstp == 1) .or. &
324 any(this%rlskstp ==
kstp) .or. &
326 this%rlstimelist))
return
328 if (this%rlstimelist)
then
329 nreleasets =
size(this%releasetimes%times)
333 nrel = this%nreleasepts * nreleasets
336 do nps = 1, this%nreleasepts
337 this%rptmass(nps) = dzero
342 if ((this%nparticles + nrel) >
size(this%particles%irpt)) &
343 call this%particles%resize( &
344 size(this%particles%irpt) + nrel, &
348 do nps = 1, this%nreleasepts
349 ic = this%rptnode(nps)
350 icu = this%dis%get_nodeuser(ic)
352 tsloop:
do nts = 1, nreleasets
353 if (this%rlstimelist)
then
354 trelease = this%releasetimes%times(nts)
356 if (trelease <
totimc .or. trelease >= tend) cycle tsloop
361 np = this%nparticles + 1
369 call this%dis%get_polyverts(ic, polyverts)
370 if (.not. point_in_polygon(x, y, polyverts))
then
371 write (errmsg,
'(a,g0,a,g0,a,i0)') &
372 'Error: release point (x=', x,
', y=', y,
') is not in cell ', icu
373 call store_error(errmsg, terminate=.false.)
374 call store_error_unit(this%inunit, terminate=.true.)
376 if (z > maxval(this%dis%top))
then
377 write (errmsg,
'(a,g0,a,g0,a,i0)') &
378 'Error: release point (z=', z,
') is above grid top ', &
380 call store_error(errmsg, terminate=.false.)
381 call store_error_unit(this%inunit, terminate=.true.)
382 else if (z < minval(this%dis%bot))
then
383 write (errmsg,
'(a,g0,a,g0,a,i0)') &
384 'Error: release point (z=', z,
') is below grid bottom ', &
386 call store_error(errmsg, terminate=.false.)
387 call store_error_unit(this%inunit, terminate=.true.)
397 call create_particle(particle)
398 if (
size(this%boundname) /= 0)
then
399 particle%name = this%boundname(nps)
404 particle%istopweaksink = this%istopweaksink
405 particle%istopzone = this%istopzone
407 select type (dis => this%dis)
409 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
411 call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
414 particle%izone = this%rptzone(ic)
417 if (this%ibound(ic) == 0)
then
420 if (this%idrape /= 0) &
421 call this%dis%highest_active(ic, this%ibound)
423 if (this%ibound(ic) == 0) &
428 if (this%localz)
then
429 top = this%fmi%dis%top(ic)
430 bot = this%fmi%dis%bot(ic)
431 hds = this%fmi%gwfhead(ic)
432 particle%z = bot + this%rptz(nps) * (hds - bot)
434 particle%z = this%rptz(nps)
436 particle%trelease = trelease
438 if (this%stoptraveltime == huge(1d0))
then
439 particle%tstop = this%stoptime
441 particle%tstop = particle%trelease + this%stoptraveltime
442 if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
444 particle%ttrack = particle%trelease
445 particle%idomain(0) = 0
446 particle%iboundary(0) = 0
447 particle%idomain(1) = 0
448 particle%iboundary(1) = 0
449 particle%idomain(2) = ic
450 particle%iboundary(2) = 0
451 particle%idomain(3) = 0
452 particle%iboundary(3) = 0
453 call this%particles%load_from_particle(particle, np)
456 this%rptmass(nps) = this%rptmass(nps) + done
real(dp), pointer, public totimc
simulation time at start of time step
integer(i4b), pointer, public kstp
current time step number
real(dp), pointer, public delt
length of the current time step
Structured grid discretization.
Vertex grid discretization.