20 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
23 real(dp),
dimension(:),
pointer,
contiguous :: porosity => null()
24 real(dp),
dimension(:),
pointer,
contiguous :: alh => null()
25 real(dp),
dimension(:),
pointer,
contiguous :: alv => null()
26 real(dp),
dimension(:),
pointer,
contiguous :: ath1 => null()
27 real(dp),
dimension(:),
pointer,
contiguous :: ath2 => null()
28 real(dp),
dimension(:),
pointer,
contiguous :: atv => null()
29 real(dp),
dimension(:),
pointer,
contiguous :: ktw => null()
30 real(dp),
dimension(:),
pointer,
contiguous :: kts => null()
31 integer(I4B),
pointer :: idisp => null()
32 integer(I4B),
pointer :: ialh => null()
33 integer(I4B),
pointer :: ialv => null()
34 integer(I4B),
pointer :: iath1 => null()
35 integer(I4B),
pointer :: iath2 => null()
36 integer(I4B),
pointer :: iatv => null()
37 integer(I4B),
pointer :: ixt3doff => null()
38 integer(I4B),
pointer :: ixt3drhs => null()
39 integer(I4B),
pointer :: iktw => null()
40 integer(I4B),
pointer :: ikts => null()
41 integer(I4B),
pointer :: ixt3d => null()
43 real(dp),
dimension(:),
pointer,
contiguous :: dispcoef => null()
44 integer(I4B),
pointer :: id22 => null()
45 integer(I4B),
pointer :: id33 => null()
46 real(dp),
dimension(:),
pointer,
contiguous :: d11 => null()
47 real(dp),
dimension(:),
pointer,
contiguous :: d22 => null()
48 real(dp),
dimension(:),
pointer,
contiguous :: d33 => null()
49 real(dp),
dimension(:),
pointer,
contiguous :: angle1 => null()
50 real(dp),
dimension(:),
pointer,
contiguous :: angle2 => null()
51 real(dp),
dimension(:),
pointer,
contiguous :: angle3 => null()
52 integer(I4B),
pointer :: iangle1 => null()
53 integer(I4B),
pointer :: iangle2 => null()
54 integer(I4B),
pointer :: iangle3 => null()
55 real(dp),
pointer :: eqnsclfac => null()
84 subroutine cnd_cr(cndobj, name_model, input_mempath, inunit, iout, fmi, &
91 character(len=*),
intent(in) :: name_model
92 character(len=*),
intent(in) :: input_mempath
93 integer(I4B),
intent(in) :: inunit
94 integer(I4B),
intent(in) :: iout
96 real(dp),
intent(in),
pointer :: eqnsclfac
99 character(len=*),
parameter :: fmtcnd = &
100 "(1x,/1x,'CND-- THERMAL CONDUCTION AND DISPERSION PACKAGE, VERSION 1, ', &
101 &'01/01/2024, INPUT READ FROM MEMPATH ', A, //)"
107 call cndobj%set_names(1, name_model,
'CND',
'CND', input_mempath)
110 call cndobj%allocate_scalars()
113 cndobj%inunit = inunit
116 cndobj%eqnsclfac => eqnsclfac
117 if (
present(gwecommon))
then
118 cndobj%gwecommon => gwecommon
121 if (cndobj%inunit > 0)
then
124 if (cndobj%iout > 0)
then
125 write (cndobj%iout, fmtcnd) input_mempath
151 if (
present(cndoptions))
then
152 this%ixt3d = cndoptions%ixt3d
155 call this%allocate_arrays(this%dis%nodes)
159 call this%source_options()
160 call this%allocate_arrays(this%dis%nodes)
163 call this%source_griddata()
167 if (this%ixt3d > 0)
then
168 call xt3d_cr(this%xt3d, this%name_model, this%inunit, this%iout, &
170 this%xt3d%ixt3d = this%ixt3d
171 call this%xt3d%xt3d_df(dis)
188 integer(I4B),
intent(in) :: moffset
192 if (this%ixt3d > 0)
call this%xt3d%xt3d_ac(moffset, sparse)
202 subroutine cnd_mc(this, moffset, matrix_sln)
207 integer(I4B),
intent(in) :: moffset
211 if (this%ixt3d > 0)
call this%xt3d%xt3d_mc(moffset, matrix_sln)
221 subroutine cnd_ar(this, ibound, porosity)
225 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
226 real(DP),
dimension(:),
pointer,
contiguous :: porosity
229 character(len=*),
parameter :: fmtcnd = &
230 "(1x,/1x,'CND-- THERMAL CONDUCTION AND DISPERSION PACKAGE, VERSION 1, ', &
231 &'5/01/2023, INPUT READ FROM UNIT ', i0, //)"
234 this%ibound => ibound
235 this%porosity => porosity
253 if (this%ixt3d > 0)
then
254 call this%xt3d%xt3d_ar(this%fmi%ibdgwfsat0, this%d11, this%id33, &
255 this%d33, this%fmi%gwfsat, this%id22, this%d22, &
256 this%iangle1, this%iangle2, this%iangle3, &
257 this%angle1, this%angle2, this%angle3)
262 call this%calcdispellipse()
265 if (this%fmi%iflowsupdated == 1)
then
266 if (this%ixt3d == 0)
then
267 call this%calcdispcoef()
268 else if (this%ixt3d > 0)
then
269 call this%xt3d%xt3d_fcpc(this%dis%nodes, .true.)
281 subroutine cnd_fc(this, kiter, nodes, nja, matrix_sln, idxglo, rhs, cnew)
284 integer(I4B) :: kiter
285 integer(I4B),
intent(in) :: nodes
286 integer(I4B),
intent(in) :: nja
288 integer(I4B),
intent(in),
dimension(nja) :: idxglo
289 real(DP),
intent(inout),
dimension(nodes) :: rhs
290 real(DP),
intent(inout),
dimension(nodes) :: cnew
292 integer(I4B) :: n, m, idiag, idiagm, ipos, isympos, isymcon
295 if (this%ixt3d > 0)
then
296 call this%xt3d%xt3d_fc(kiter, matrix_sln, idxglo, rhs, cnew)
299 if (this%fmi%ibdgwfsat0(n) == 0) cycle
300 idiag = this%dis%con%ia(n)
301 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
302 if (this%dis%con%mask(ipos) == 0) cycle
303 m = this%dis%con%ja(ipos)
305 if (this%fmi%ibdgwfsat0(m) == 0) cycle
306 isympos = this%dis%con%jas(ipos)
307 dnm = this%dispcoef(isympos)
310 call matrix_sln%add_value_pos(idxglo(ipos), dnm)
311 call matrix_sln%add_value_pos(idxglo(idiag), -dnm)
314 idiagm = this%dis%con%ia(m)
315 isymcon = this%dis%con%isym(ipos)
316 call matrix_sln%add_value_pos(idxglo(isymcon), dnm)
317 call matrix_sln%add_value_pos(idxglo(idiagm), -dnm)
334 real(DP),
intent(inout),
dimension(:) :: cnew
335 real(DP),
intent(inout),
dimension(:) :: flowja
337 integer(I4B) :: n, m, ipos, isympos
341 if (this%ixt3d > 0)
then
342 call this%xt3d%xt3d_flowja(cnew, flowja)
344 do n = 1, this%dis%nodes
345 if (this%fmi%ibdgwfsat0(n) == 0) cycle
346 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
347 m = this%dis%con%ja(ipos)
348 if (this%fmi%ibdgwfsat0(m) == 0) cycle
349 isympos = this%dis%con%jas(ipos)
350 dnm = this%dispcoef(isympos)
352 qnm = dnm * (cnew(m) - cnew(n))
353 flowja(ipos) = flowja(ipos) + qnm
374 call this%NumericalPackageType%allocate_scalars()
383 call mem_allocate(this%ixt3doff,
'IXT3DOFF', this%memoryPath)
384 call mem_allocate(this%ixt3drhs,
'IXT3DRHS', this%memoryPath)
388 call mem_allocate(this%iangle1,
'IANGLE1', this%memoryPath)
389 call mem_allocate(this%iangle2,
'IANGLE2', this%memoryPath)
390 call mem_allocate(this%iangle3,
'IANGLE3', this%memoryPath)
426 integer(I4B),
intent(in) :: nodes
429 call mem_allocate(this%alh, nodes,
'ALH', trim(this%memoryPath))
430 call mem_allocate(this%alv, nodes,
'ALV', trim(this%memoryPath))
431 call mem_allocate(this%ath1, nodes,
'ATH1', trim(this%memoryPath))
432 call mem_allocate(this%ath2, nodes,
'ATH2', trim(this%memoryPath))
433 call mem_allocate(this%atv, nodes,
'ATV', trim(this%memoryPath))
434 call mem_allocate(this%d11, nodes,
'D11', trim(this%memoryPath))
435 call mem_allocate(this%d22, nodes,
'D22', trim(this%memoryPath))
436 call mem_allocate(this%d33, nodes,
'D33', trim(this%memoryPath))
437 call mem_allocate(this%angle1, nodes,
'ANGLE1', trim(this%memoryPath))
438 call mem_allocate(this%angle2, nodes,
'ANGLE2', trim(this%memoryPath))
439 call mem_allocate(this%angle3, nodes,
'ANGLE3', trim(this%memoryPath))
440 call mem_allocate(this%ktw, nodes,
'KTW', trim(this%memoryPath))
441 call mem_allocate(this%kts, nodes,
'KTS', trim(this%memoryPath))
444 if (this%ixt3d == 0)
then
445 call mem_allocate(this%dispcoef, this%dis%njas,
'DISPCOEF', &
446 trim(this%memoryPath))
448 call mem_allocate(this%dispcoef, 0,
'DISPCOEF', trim(this%memoryPath))
472 if (this%inunit /= 0)
then
487 if (this%ixt3d > 0)
call this%xt3d%xt3d_da()
491 if (this%ixt3d > 0)
deallocate (this%xt3d)
492 nullify (this%gwecommon)
513 call this%NumericalPackageType%da()
526 write (this%iout,
'(1x,a)')
'Setting CND Options'
527 write (this%iout,
'(4x,a,i0)')
'XT3D formulation [0=INACTIVE, 1=ACTIVE, &
528 &3=ACTIVE RHS] set to: ', this%ixt3d
529 write (this%iout,
'(1x,a,/)')
'End Setting CND Options'
546 call mem_set_value(this%ixt3doff,
'XT3D_OFF', this%input_mempath, &
548 call mem_set_value(this%ixt3drhs,
'XT3D_RHS', this%input_mempath, &
552 if (found%xt3d_off) this%ixt3d = 0
553 if (found%xt3d_rhs) this%ixt3d = 2
556 if (this%iout > 0)
then
557 call this%log_options(found)
571 write (this%iout,
'(1x,a)')
'Setting CND Griddata'
574 write (this%iout,
'(4x,a)')
'ALH set from input file'
578 write (this%iout,
'(4x,a)')
'ALV set from input file'
582 write (this%iout,
'(4x,a)')
'ATH1 set from input file'
586 write (this%iout,
'(4x,a)')
'ATH2 set from input file'
590 write (this%iout,
'(4x,a)')
'ATV set from input file'
594 write (this%iout,
'(4x,a)')
'KTW set from input file'
598 write (this%iout,
'(4x,a)')
'KTS set from input file'
601 write (this%iout,
'(1x,a,/)')
'End Setting CND Griddata'
619 character(len=LINELENGTH) :: errmsg
621 integer(I4B),
dimension(:),
pointer,
contiguous :: map
626 if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser
629 call mem_set_value(this%alh,
'ALH', this%input_mempath, map, found%alh)
630 call mem_set_value(this%alv,
'ALV', this%input_mempath, map, found%alv)
631 call mem_set_value(this%ath1,
'ATH1', this%input_mempath, map, found%ath1)
632 call mem_set_value(this%ath2,
'ATH2', this%input_mempath, map, found%ath2)
633 call mem_set_value(this%atv,
'ATV', this%input_mempath, map, found%atv)
634 call mem_set_value(this%ktw,
'KTW', this%input_mempath, map, found%ktw)
635 call mem_set_value(this%kts,
'KTS', this%input_mempath, map, found%kts)
638 if (found%alh) this%ialh = 1
639 if (found%alv) this%ialv = 1
640 if (found%ath1) this%iath1 = 1
641 if (found%ath2) this%iath2 = 1
642 if (found%atv) this%iatv = 1
643 if (found%ktw) this%iktw = 1
644 if (found%kts) this%ikts = 1
647 if (found%alh) this%idisp = this%idisp + 1
648 if (found%alv) this%idisp = this%idisp + 1
649 if (found%ath1) this%idisp = this%idisp + 1
650 if (found%ath2) this%idisp = this%idisp + 1
653 if (this%idisp > 0)
then
654 if (.not. (found%alh .and. found%ath1))
then
655 write (errmsg,
'(1x,a)') &
656 'if dispersivities are specified then ALH and ATH1 are required.'
660 if (.not. found%alv) &
662 'ALH', trim(this%memoryPath))
664 if (.not. found%ath2) &
666 'ATH1', trim(this%memoryPath))
668 if (.not. found%atv) &
670 'ATH2', trim(this%memoryPath))
680 if (this%iout > 0)
then
681 call this%log_griddata(found)
695 integer(I4B) :: nodes, n
696 real(DP) :: q, qx, qy, qz
697 real(DP) :: alh, alv, ath1, ath2, atv, a
698 real(DP) :: al, at1, at2
699 real(DP) :: qzoqsquared
705 nodes =
size(this%d11)
712 this%angle1(n) =
dzero
713 this%angle2(n) =
dzero
714 this%angle3(n) =
dzero
715 if (this%fmi%ibdgwfsat0(n) == 0) cycle
722 qx = this%fmi%gwfspdis(1, n)
723 qy = this%fmi%gwfspdis(2, n)
724 qz = this%fmi%gwfspdis(3, n)
725 q = qx**2 + qy**2 + qz**2
726 if (q >
dzero) q = sqrt(q)
734 if (this%idisp > 0)
then
744 if (this%iktw > 0) ktbulk = ktbulk + this%porosity(n) * this%ktw(n) * &
746 if (this%ikts > 0) ktbulk = ktbulk + (
done - this%porosity(n)) * this%kts(n)
754 dstar = ktbulk / (this%gwecommon%gwecpw * this%gwecommon%gwerhow)
761 qzoqsquared = (qz / q)**2
762 al = alh * (
done - qzoqsquared) + alv * qzoqsquared
763 at1 = ath1 * (
done - qzoqsquared) + atv * qzoqsquared
764 at2 = ath2 * (
done - qzoqsquared) + atv * qzoqsquared
768 qsw = q * this%fmi%gwfsat(n) * this%eqnsclfac
769 this%d11(n) = al * qsw + ktbulk
770 this%d22(n) = at1 * qsw + ktbulk
771 this%d33(n) = at2 * qsw + ktbulk
774 if (this%idisp > 0)
then
777 this%angle3(n) =
dzero
781 if (q >
dzero) a = qz / q
782 this%angle2(n) = asin(a)
785 a = q * cos(this%angle2(n))
795 elseif (a >=
done)
then
796 this%angle1(n) =
dzero
798 this%angle1(n) = acos(a)
817 integer(I4B) :: nodes, n, m, idiag, ipos
818 real(DP) :: clnm, clmn, dn, dm
819 real(DP) :: vg1, vg2, vg3
820 integer(I4B) :: ihc, isympos
821 integer(I4B) :: iavgmeth
822 real(DP) :: satn, satm, topn, topm, botn, botm
823 real(DP) :: hwva, cond, cn, cm, denom
824 real(DP) :: anm, amn, thksatn, thksatm
830 nodes =
size(this%d11)
832 if (this%fmi%ibdgwfsat0(n) == 0) cycle
833 idiag = this%dis%con%ia(n)
834 do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
837 m = this%dis%con%ja(ipos)
841 isympos = this%dis%con%jas(ipos)
842 this%dispcoef(isympos) =
dzero
843 if (this%fmi%ibdgwfsat0(m) == 0) cycle
846 hwva = this%dis%con%hwva(isympos)
847 clnm = this%dis%con%cl1(isympos)
848 clmn = this%dis%con%cl2(isympos)
849 ihc = this%dis%con%ihc(isympos)
850 topn = this%dis%top(n)
851 topm = this%dis%top(m)
852 botn = this%dis%bot(n)
853 botm = this%dis%bot(m)
856 satn = this%fmi%ibdgwfsat0(n)
857 satm = this%fmi%ibdgwfsat0(m)
862 call this%dis%connection_normal(n, m, ihc, vg1, vg2, vg3, ipos)
863 dn =
hyeff(this%d11(n), this%d22(n), this%d33(n), &
864 this%angle1(n), this%angle2(n), this%angle3(n), &
865 vg1, vg2, vg3, iavgmeth)
866 dm =
hyeff(this%d11(m), this%d22(m), this%d33(m), &
867 this%angle1(m), this%angle2(m), this%angle3(m), &
868 vg1, vg2, vg3, iavgmeth)
873 clnm = satn * (topn - botn) *
dhalf
874 clmn = satm * (topm - botm) *
dhalf
878 if (satn ==
dzero)
then
880 else if (n > m .and. satn <
done)
then
885 if (satm ==
dzero)
then
887 else if (m > n .and. satm <
done)
then
903 thksatn = (topn - botn) * satn
904 thksatm = (topm - botm) * satm
921 if (clnm >
dzero) cn = dn * anm / clnm
923 if (clmn >
dzero) cm = dm * amn / clmn
925 if (denom >
dzero)
then
926 cond = cn * cm / denom
932 this%dispcoef(isympos) = cond
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dpi
real constant
real(dp), parameter dzero
real constant zero
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
subroutine cnd_df(this, dis, cndOptions)
Define CND object.
subroutine cnd_ad(this)
Advance method for the package.
subroutine source_griddata(this)
Update CND simulation data from input mempath.
subroutine cnd_fc(this, kiter, nodes, nja, matrix_sln, idxglo, rhs, cnew)
Fill coefficient method for package.
subroutine, public cnd_cr(cndobj, name_model, input_mempath, inunit, iout, fmi, eqnsclfac, gwecommon)
Create a new CND object.
subroutine cnd_ar(this, ibound, porosity)
Allocate and read method for package.
subroutine calcdispellipse(this)
Calculate dispersion coefficients.
subroutine calcdispcoef(this)
Calculate dispersion coefficients.
subroutine cnd_ac(this, moffset, sparse)
Add connections to CND.
subroutine cnd_da(this)
@ brief Deallocate memory
subroutine allocate_arrays(this, nodes)
@ brief Allocate arrays for package
subroutine allocate_scalars(this)
@ brief Allocate scalar variables for package
subroutine source_options(this)
Update simulation mempath options.
subroutine cnd_mc(this, moffset, matrix_sln)
Map CND connections.
subroutine log_griddata(this, found)
Write dimensions to list file.
subroutine cnd_cq(this, cnew, flowja)
@ brief Calculate flows for package
subroutine log_options(this, found)
Write user options to list file.
This module contains stateless conductance functions.
real(dp) function, public staggered_thkfrac(top, bot, sat, topc, botc)
Calculate the thickness fraction for staggered grids.
General-purpose hydrogeologic functions.
real(dp) function, public hyeff(k11, k22, k33, ang1, ang2, ang3, vg1, vg2, vg3, iavgmeth)
Calculate the effective horizontal hydraulic conductivity from an ellipse using a specified direction...
This module defines variable data types.
subroutine, public memorylist_remove(component, subcomponent, context)
This module contains the base numerical package type.
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
This module contains simulation variables.
character(len=linelength) idm_context
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number
subroutine, public xt3d_cr(xt3dobj, name_model, inunit, iout, ldispopt)
Create a new xt3d object.
data structure (and helpers) for passing cnd option data