25 real(dp),
dimension(:),
pointer :: conc => null()
26 integer(I4B),
dimension(:),
pointer :: icbund => null()
30 integer(I4B),
pointer :: thermivisc => null()
31 integer(I4B),
pointer :: idxtmpr => null()
32 integer(I4B),
pointer :: ioutvisc => null()
33 integer(I4B),
pointer :: iconcset => null()
34 integer(I4B),
pointer :: ireadelev => null()
35 integer(I4B),
dimension(:),
pointer,
contiguous :: ivisc => null()
36 real(dp),
pointer :: viscref => null()
37 real(dp),
dimension(:),
pointer,
contiguous :: visc => null()
38 real(dp),
dimension(:),
pointer,
contiguous :: elev => null()
39 integer(I4B),
dimension(:),
pointer :: ibound => null()
41 integer(I4B),
pointer :: nviscspecies => null()
42 real(dp),
dimension(:),
pointer,
contiguous :: dviscdc => null()
43 real(dp),
dimension(:),
pointer,
contiguous :: cviscref => null()
44 real(dp),
dimension(:),
pointer,
contiguous :: ctemp => null()
45 character(len=LENMODELNAME),
dimension(:),
allocatable :: cmodelname
46 character(len=LENAUXNAME),
dimension(:),
allocatable :: cauxspeciesname
47 character(len=LENAUXNAME) :: name_temp_spec =
'TEMPERATURE'
50 real(dp),
pointer :: a2 => null()
51 real(dp),
pointer :: a3 => null()
52 real(dp),
pointer :: a4 => null()
56 real(dp),
dimension(:),
pointer,
contiguous :: k11 => null()
57 real(dp),
dimension(:),
pointer,
contiguous :: k22 => null()
58 real(dp),
dimension(:),
pointer,
contiguous :: k33 => null()
59 real(dp),
dimension(:),
pointer,
contiguous :: k11input => null()
60 real(dp),
dimension(:),
pointer,
contiguous :: k22input => null()
61 real(dp),
dimension(:),
pointer,
contiguous :: k33input => null()
62 integer(I4B),
pointer :: kchangeper => null()
63 integer(I4B),
pointer :: kchangestp => null()
64 integer(I4B),
dimension(:),
pointer,
contiguous :: nodekchange => null()
97 function calc_visc(ivisc, viscref, dviscdc, cviscref, conc, &
98 a2, a3, a4)
result(visc)
100 integer(I4B),
dimension(:),
intent(in) :: ivisc
101 real(dp),
intent(in) :: viscref
102 real(dp),
dimension(:),
intent(in) :: dviscdc
103 real(dp),
dimension(:),
intent(in) :: cviscref
104 real(dp),
dimension(:),
intent(in) :: conc
105 real(dp),
intent(in) :: a2, a3, a4
109 integer(I4B) :: nviscspec
114 nviscspec =
size(dviscdc)
118 if (ivisc(i) == 1)
then
119 visc = visc + dviscdc(i) * (conc(i) - cviscref(i))
121 expon = -1 * a3 * ((conc(i) - cviscref(i)) / &
122 ((conc(i) + a4) * (cviscref(i) + a4)))
123 mu_t = viscref * a2**expon
129 visc = (visc - viscref) + mu_t
142 subroutine vsc_cr(vscobj, name_model, inunit, iout)
145 character(len=*),
intent(in) :: name_model
146 integer(I4B),
intent(in) :: inunit
147 integer(I4B),
intent(in) :: iout
153 call vscobj%set_names(1, name_model,
'VSC',
'VSC')
156 call vscobj%allocate_scalars()
159 vscobj%inunit = inunit
163 call vscobj%parser%Initialize(vscobj%inunit, vscobj%iout)
179 character(len=*),
parameter :: fmtvsc = &
180 "(1x,/1x,'VSC -- Viscosity Package, version 1, 11/15/2022', &
181 &' input read from unit ', i0, //)"
184 write (this%iout, fmtvsc) this%inunit
189 if (.not.
present(vsc_input))
then
192 call this%read_options()
195 call this%read_dimensions()
198 call this%set_options(vsc_input)
199 this%nviscspecies = vsc_input%nviscspecies
203 call this%allocate_arrays(dis%nodes)
205 if (.not.
present(vsc_input))
then
208 call this%read_packagedata()
211 call this%set_packagedata(vsc_input)
226 integer(I4B),
dimension(:),
pointer :: ibound
229 this%ibound => ibound
232 call this%set_npf_pointers()
255 class(
bndtype),
pointer :: packobj
258 select case (packobj%filtyp)
262 select type (packobj)
264 call packobj%bnd_activate_viscosity()
269 select type (packobj)
271 call packobj%bnd_activate_viscosity()
276 select type (packobj)
278 call packobj%bnd_activate_viscosity()
283 select type (packobj)
285 call packobj%lak_activate_viscosity()
290 select type (packobj)
292 call packobj%sfr_activate_viscosity()
297 select type (packobj)
299 call packobj%maw_activate_viscosity()
319 character(len=LENMEMPATH) :: npfMemoryPath
324 call mem_setptr(this%k11,
'K11', npfmemorypath)
325 call mem_setptr(this%k22,
'K22', npfmemorypath)
326 call mem_setptr(this%k33,
'K33', npfmemorypath)
327 call mem_setptr(this%k11input,
'K11INPUT', npfmemorypath)
328 call mem_setptr(this%k22input,
'K22INPUT', npfmemorypath)
329 call mem_setptr(this%k33input,
'K33INPUT', npfmemorypath)
330 call mem_setptr(this%kchangeper,
'KCHANGEPER', npfmemorypath)
331 call mem_setptr(this%kchangestp,
'KCHANGESTP', npfmemorypath)
332 call mem_setptr(this%nodekchange,
'NODEKCHANGE', npfmemorypath)
348 character(len=LINELENGTH) :: errmsg
351 character(len=*),
parameter :: fmtc = &
352 "('Viscosity Package does not have a concentration set &
353 &for species ',i0,'. One or more model names may be specified &
354 &incorrectly in the PACKAGEDATA block or a GWF-GWT exchange may need &
359 do i = 1, this%nviscspecies
360 if (.not.
associated(this%modelconc(i)%conc))
then
361 write (errmsg, fmtc) i
366 call this%parser%StoreErrorUnit()
384 call this%vsc_calcvisc()
400 class(
bndtype),
pointer :: packobj
401 real(DP),
intent(in),
dimension(:) :: hnew
404 integer(I4B) :: n, locvisc, locelev
405 integer(I4B),
dimension(:),
allocatable :: locconc
410 allocate (locconc(this%nviscspecies))
414 do n = 1, packobj%naux
415 if (packobj%auxname(n) ==
'VISCOSITY')
then
417 else if (packobj%auxname(n) ==
'ELEVATION')
then
423 do i = 1, this%nviscspecies
425 do j = 1, packobj%naux
426 if (this%cauxspeciesname(i) == packobj%auxname(j))
then
431 if (locconc(i) == 0)
then
439 select case (packobj%filtyp)
440 case (
'GHB',
'DRN',
'RIV')
444 locelev, locvisc, locconc, this%dviscdc, &
445 this%cviscref, this%ivisc, this%a2, this%a3, &
452 call vsc_ad_lak(packobj, this%visc, this%viscref, this%elev, locvisc, &
453 locconc, this%dviscdc, this%cviscref, this%ivisc, &
454 this%a2, this%a3, this%a4, this%ctemp)
460 call vsc_ad_sfr(packobj, this%visc, this%viscref, this%elev, locvisc, &
461 locconc, this%dviscdc, this%cviscref, this%ivisc, &
462 this%a2, this%a3, this%a4, this%ctemp)
466 call vsc_ad_maw(packobj, this%visc, this%viscref, this%elev, locvisc, &
467 locconc, this%dviscdc, this%cviscref, this%ivisc, &
468 this%a2, this%a3, this%a4, this%ctemp)
490 locvisc, locconc, dviscdc, cviscref, &
491 ivisc, a2, a3, a4, ctemp)
497 class(
bndtype),
pointer :: packobj
499 real(DP),
intent(in),
dimension(:) :: hnew
500 real(DP),
intent(in),
dimension(:) :: visc
501 real(DP),
intent(in) :: a2, a3, a4
502 real(DP),
intent(in) :: viscref
503 integer(I4B),
intent(in) :: locelev
504 integer(I4B),
intent(in) :: locvisc
505 integer(I4B),
dimension(:),
intent(in) :: locconc
506 integer(I4B),
dimension(:),
intent(in) :: ivisc
507 real(DP),
dimension(:),
intent(in) :: dviscdc
508 real(DP),
dimension(:),
intent(in) :: cviscref
509 real(DP),
dimension(:),
intent(inout) :: ctemp
516 do n = 1, packobj%nbound
517 node = packobj%nodelist(n)
520 if (packobj%ibound(node) <= 0) cycle
524 cviscref, ctemp, ivisc, a2, a3, a4, &
528 select case (packobj%filtyp)
530 select type (packobj)
533 packobj%condinput(n))
536 select type (packobj)
539 packobj%condinput(n))
542 select type (packobj)
545 packobj%condinput(n))
549 packobj%condinput(n))
563 subroutine vsc_ad_sfr(packobj, visc, viscref, elev, locvisc, locconc, &
564 dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
568 class(
bndtype),
pointer :: packobj
570 real(DP),
intent(in) :: viscref
571 real(DP),
intent(in) :: a2, a3, a4
572 integer(I4B),
intent(in) :: locvisc
573 integer(I4B),
dimension(:),
intent(in) :: locconc
574 integer(I4B),
dimension(:),
intent(in) :: ivisc
575 real(DP),
dimension(:),
intent(in) :: visc
576 real(DP),
dimension(:),
intent(in) :: elev
577 real(DP),
dimension(:),
intent(in) :: dviscdc
578 real(DP),
dimension(:),
intent(in) :: cviscref
579 real(DP),
dimension(:),
intent(inout) :: ctemp
586 select type (packobj)
588 do n = 1, packobj%nbound
591 node = packobj%nodelist(n)
594 if (packobj%ibound(node) <= 0) cycle
600 cviscref, ctemp, ivisc, a2, a3, a4, &
620 subroutine vsc_ad_lak(packobj, visc, viscref, elev, locvisc, locconc, &
621 dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
625 class(
bndtype),
pointer :: packobj
627 real(DP),
intent(in) :: viscref
628 real(DP),
intent(in) :: a2, a3, a4
629 integer(I4B),
intent(in) :: locvisc
630 integer(I4B),
dimension(:),
intent(in) :: locconc
631 integer(I4B),
dimension(:),
intent(in) :: ivisc
632 real(DP),
dimension(:),
intent(in) :: visc
633 real(DP),
dimension(:),
intent(in) :: elev
634 real(DP),
dimension(:),
intent(in) :: dviscdc
635 real(DP),
dimension(:),
intent(in) :: cviscref
636 real(DP),
dimension(:),
intent(inout) :: ctemp
643 select type (packobj)
645 do n = 1, packobj%nbound
648 node = packobj%nodelist(n)
651 if (packobj%ibound(node) <= 0) cycle
657 cviscref, ctemp, ivisc, a2, a3, a4, &
677 subroutine vsc_ad_maw(packobj, visc, viscref, elev, locvisc, locconc, &
678 dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
682 class(
bndtype),
pointer :: packobj
684 real(DP),
intent(in) :: viscref
685 real(DP),
intent(in) :: a2, a3, a4
686 integer(I4B),
intent(in) :: locvisc
687 integer(I4B),
dimension(:),
intent(in) :: locconc
688 integer(I4B),
dimension(:),
intent(in) :: ivisc
689 real(DP),
dimension(:),
intent(in) :: visc
690 real(DP),
dimension(:),
intent(in) :: elev
691 real(DP),
dimension(:),
intent(in) :: dviscdc
692 real(DP),
dimension(:),
intent(in) :: cviscref
693 real(DP),
dimension(:),
intent(inout) :: ctemp
700 select type (packobj)
702 do n = 1, packobj%nbound
705 node = packobj%nodelist(n)
708 if (packobj%ibound(node) <= 0) cycle
714 cviscref, ctemp, ivisc, a2, a3, a4, &
736 real(dp),
intent(in) :: viscref
737 real(dp),
intent(in) :: bndvisc
738 real(dp),
intent(in) :: spcfdcond
741 real(dp) :: updatedcond
746 updatedcond = vscratio * spcfdcond
756 real(dp),
intent(in) :: viscref
757 real(dp),
intent(in) :: bndvisc
759 real(dp) :: viscratio
761 viscratio = viscref / bndvisc
776 ctemp, ivisc, a2, a3, a4, auxvar)
result(viscbnd)
779 integer(I4B),
intent(in) :: n
780 integer(I4B),
intent(in) :: locvisc
781 real(dp),
intent(in) :: a2, a3, a4
782 integer(I4B),
dimension(:),
intent(in) :: ivisc
783 integer(I4B),
dimension(:),
intent(in) :: locconc
784 real(dp),
intent(in) :: viscref
785 real(dp),
dimension(:),
intent(in) :: dviscdc
786 real(dp),
dimension(:),
intent(in) :: cviscref
787 real(dp),
dimension(:),
intent(inout) :: ctemp
788 real(dp),
dimension(:, :),
intent(in) :: auxvar
795 if (locvisc > 0)
then
797 viscbnd = auxvar(locvisc, n)
798 else if (locconc(1) > 0)
then
800 do i = 1,
size(locconc)
802 if (locconc(i) > 0)
then
803 ctemp(i) = auxvar(locconc(i), n)
806 viscbnd =
calc_visc(ivisc, viscref, dviscdc, cviscref, ctemp, a2, a3, a4)
828 integer(I4B),
intent(in) :: n, m
829 real(DP),
intent(in) :: gwhdn, gwhdm
830 real(DP),
intent(inout) :: viscratio
832 integer(I4B) :: cellid
835 if (gwhdm > gwhdn)
then
837 else if (gwhdn >= gwhdm)
then
840 call this%calc_q_visc(cellid, viscratio)
855 integer(I4B),
intent(in) :: cellid
857 real(DP),
intent(inout) :: viscratio
862 visc = this%visc(cellid)
883 real(DP) :: viscratio
888 do n = 1, this%dis%nodes
889 call this%calc_q_visc(n, viscratio)
890 this%k11(n) = this%k11input(n) * viscratio
891 this%k22(n) = this%k22input(n) * viscratio
892 this%k33(n) = this%k33input(n) * viscratio
893 this%nodekchange(n) = 1
897 call this%vsc_set_changed_at(
kper,
kstp)
911 integer(I4B),
intent(in) :: kper
912 integer(I4B),
intent(in) :: kstp
914 this%kchangeper = kper
915 this%kchangestp = kstp
928 integer(I4B),
intent(in) :: idvfl
930 character(len=1) :: cdatafmp =
' ', editdesc =
' '
931 integer(I4B) :: ibinun
932 integer(I4B) :: iprint
933 integer(I4B) :: nvaluesp
934 integer(I4B) :: nwidthp
938 if (this%ioutvisc /= 0)
then
943 if (idvfl == 0) ibinun = 0
946 if (ibinun /= 0)
then
951 if (this%ioutvisc /= 0)
then
952 ibinun = this%ioutvisc
953 call this%dis%record_array(this%visc, this%iout, iprint, ibinun, &
954 ' VISCOSITY', cdatafmp, nvaluesp, &
955 nwidthp, editdesc, dinact)
972 if (this%inunit > 0)
then
978 deallocate (this%cmodelname)
979 deallocate (this%cauxspeciesname)
980 deallocate (this%modelconc)
999 nullify (this%k11input)
1000 nullify (this%k22input)
1001 nullify (this%k33input)
1002 nullify (this%kchangeper)
1003 nullify (this%kchangestp)
1004 nullify (this%nodekchange)
1007 call this%NumericalPackageType%da()
1022 character(len=LINELENGTH) :: errmsg, keyword
1023 integer(I4B) :: ierr
1024 logical :: isfound, endOfBlock
1027 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
1028 supportopenclose=.true.)
1032 write (this%iout,
'(/1x,a)')
'Processing VSC DIMENSIONS block'
1034 call this%parser%GetNextLine(endofblock)
1035 if (endofblock)
exit
1036 call this%parser%GetStringCaps(keyword)
1037 select case (keyword)
1038 case (
'NVISCSPECIES')
1039 this%nviscspecies = this%parser%GetInteger()
1040 write (this%iout,
'(4x,a,i0)')
'NVISCSPECIES = ', this%nviscspecies
1042 write (errmsg,
'(a,a)') &
1043 'Unknown VSC dimension: ', trim(keyword)
1045 call this%parser%StoreErrorUnit()
1048 write (this%iout,
'(1x,a)')
'End of VSC DIMENSIONS block'
1050 call store_error(
'Required VSC DIMENSIONS block not found.')
1051 call this%parser%StoreErrorUnit()
1055 if (this%nviscspecies < 1)
then
1056 call store_error(
'NVISCSPECIES must be greater than zero.')
1057 call this%parser%StoreErrorUnit()
1072 character(len=LINELENGTH) :: errmsg
1073 character(len=LINELENGTH) :: line
1074 integer(I4B) :: ierr
1075 integer(I4B) :: iviscspec
1076 logical :: isfound, endOfBlock
1077 logical :: blockrequired
1078 integer(I4B),
dimension(:),
allocatable :: itemp
1079 character(len=10) :: c10
1080 character(len=16) :: c16
1082 character(len=*),
parameter :: fmterr = &
1083 "('Invalid value for IRHOSPEC (',i0,') detected in VSC Package. &
1084 &IRHOSPEC must be > 0 and <= NVISCSPECIES, and duplicate values &
1085 &are not allowed.')"
1088 allocate (itemp(this%nviscspecies))
1092 blockrequired = .true.
1093 call this%parser%GetBlock(
'PACKAGEDATA', isfound, ierr, &
1094 blockrequired=blockrequired, &
1095 supportopenclose=.true.)
1099 write (this%iout,
'(1x,a)')
'Procesing VSC PACKAGEDATA block'
1101 call this%parser%GetNextLine(endofblock)
1102 if (endofblock)
exit
1103 iviscspec = this%parser%GetInteger()
1104 if (iviscspec < 1 .or. iviscspec > this%nviscspecies)
then
1105 write (errmsg, fmterr) iviscspec
1108 if (itemp(iviscspec) /= 0)
then
1109 write (errmsg, fmterr) iviscspec
1112 itemp(iviscspec) = 1
1114 this%dviscdc(iviscspec) = this%parser%GetDouble()
1115 this%cviscref(iviscspec) = this%parser%GetDouble()
1116 call this%parser%GetStringCaps(this%cmodelname(iviscspec))
1117 call this%parser%GetStringCaps(this%cauxspeciesname(iviscspec))
1119 if (this%cauxspeciesname(iviscspec) == this%name_temp_spec)
then
1120 if (this%idxtmpr > 0)
then
1121 write (errmsg,
'(a)')
'More than one species in VSC input identified &
1122 &as '//trim(this%name_temp_spec)//
'. Only one species may be &
1123 &designated to represent temperature.'
1126 this%idxtmpr = iviscspec
1127 if (this%thermivisc == 2)
then
1128 this%ivisc(iviscspec) = 2
1134 call store_error(
'Required VSC PACKAGEDATA block not found.')
1135 call this%parser%StoreErrorUnit()
1140 call this%parser%StoreErrorUnit()
1144 write (this%iout,
'(/,1x,a)')
'Summary of species information in VSC Package'
1145 write (this%iout,
'(1a11,5a17)') &
1146 'Species',
'DVISCDC',
'CVISCREF',
'Model',
'AUXSPECIESNAME'
1147 do iviscspec = 1, this%nviscspecies
1148 write (c10,
'(i0)') iviscspec
1149 line =
' '//adjustr(c10)
1151 write (c16,
'(g15.6)') this%dviscdc(iviscspec)
1152 line = trim(line)//
' '//adjustr(c16)
1153 write (c16,
'(g15.6)') this%cviscref(iviscspec)
1154 line = trim(line)//
' '//adjustr(c16)
1155 write (c16,
'(a)') this%cmodelname(iviscspec)
1156 line = trim(line)//
' '//adjustr(c16)
1157 write (c16,
'(a)') this%cauxspeciesname(iviscspec)
1158 line = trim(line)//
' '//adjustr(c16)
1159 write (this%iout,
'(a)') trim(line)
1165 write (this%iout,
'(/,1x,a)')
'End of VSC PACKAGEDATA block'
1178 integer(I4B) :: ispec
1180 do ispec = 1, this%nviscspecies
1181 this%dviscdc(ispec) = input_data%dviscdc(ispec)
1182 this%cviscref(ispec) = input_data%cviscref(ispec)
1183 this%cmodelname(ispec) = input_data%cmodelname(ispec)
1184 this%cauxspeciesname(ispec) = input_data%cauxspeciesname(ispec)
1205 do n = 1, this%dis%nodes
1206 do i = 1, this%nviscspecies
1207 if (this%modelconc(i)%icbund(n) == 0)
then
1210 this%ctemp(i) = this%modelconc(i)%conc(n)
1214 this%visc(n) =
calc_visc(this%ivisc, this%viscref, this%dviscdc, &
1215 this%cviscref, this%ctemp, this%a2, &
1235 call this%NumericalPackageType%allocate_scalars()
1238 call mem_allocate(this%thermivisc,
'THERMIVISC', this%memoryPath)
1239 call mem_allocate(this%idxtmpr,
'IDXTMPR', this%memoryPath)
1240 call mem_allocate(this%ioutvisc,
'IOUTVISC', this%memoryPath)
1241 call mem_allocate(this%ireadelev,
'IREADELEV', this%memoryPath)
1242 call mem_allocate(this%iconcset,
'ICONCSET', this%memoryPath)
1243 call mem_allocate(this%viscref,
'VISCREF', this%memoryPath)
1248 call mem_allocate(this%nviscspecies,
'NVISCSPECIES', this%memoryPath)
1261 this%nviscspecies = 0
1274 integer(I4B),
intent(in) :: nodes
1279 call mem_allocate(this%visc, nodes,
'VISC', this%memoryPath)
1280 call mem_allocate(this%ivisc, this%nviscspecies,
'IVISC', this%memoryPath)
1281 call mem_allocate(this%dviscdc, this%nviscspecies,
'DRHODC', &
1283 call mem_allocate(this%cviscref, this%nviscspecies,
'CRHOREF', &
1285 call mem_allocate(this%ctemp, this%nviscspecies,
'CTEMP', this%memoryPath)
1286 allocate (this%cmodelname(this%nviscspecies))
1287 allocate (this%cauxspeciesname(this%nviscspecies))
1288 allocate (this%modelconc(this%nviscspecies))
1292 this%visc(i) = this%viscref
1296 do i = 1, this%nviscspecies
1298 this%dviscdc(i) =
dzero
1299 this%cviscref(i) =
dzero
1300 this%ctemp(i) =
dzero
1301 this%cmodelname(i) =
''
1302 this%cauxspeciesname(i) =
''
1320 character(len=LINELENGTH) :: warnmsg, errmsg, keyword, keyword2
1321 character(len=MAXCHARLEN) :: fname
1322 integer(I4B) :: ierr
1323 logical :: isfound, endOfBlock
1325 character(len=*),
parameter :: fmtfileout = &
1326 "(1x, 'VSC', 1x, a, 1x, 'Will be saved to file: ', &
1327 &a, /4x, 'opened on unit: ', I7)"
1328 character(len=*),
parameter :: fmtlinear = &
1329 "(/,1x,'Viscosity will vary linearly with temperature &
1331 character(len=*),
parameter :: fmtnonlinear = &
1332 "(/,1x,'Viscosity will vary non-linearly with temperature &
1336 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
1337 supportopenclose=.true., blockrequired=.false.)
1341 write (this%iout,
'(1x,a)')
'Processing VSC OPTIONS block'
1343 call this%parser%GetNextLine(endofblock)
1344 if (endofblock)
exit
1345 call this%parser%GetStringCaps(keyword)
1346 select case (keyword)
1348 this%viscref = this%parser%GetDouble()
1349 write (this%iout,
'(4x,a,1pg15.6)') &
1350 'Reference viscosity has been set to: ', &
1353 call this%parser%GetStringCaps(keyword)
1354 if (keyword ==
'FILEOUT')
then
1355 call this%parser%GetString(fname)
1357 call openfile(this%ioutvisc, this%iout, fname,
'DATA(BINARY)', &
1359 write (this%iout, fmtfileout) &
1360 'VISCOSITY', fname, this%ioutvisc
1362 errmsg =
'Optional VISCOSITY keyword must be '// &
1363 'followed by FILEOUT'
1366 case (
'TEMPERATURE_SPECIES_NAME')
1367 call this%parser%GetStringCaps(this%name_temp_spec)
1368 write (this%iout,
'(4x, a)')
'Temperature species name set to: '// &
1369 trim(this%name_temp_spec)
1370 case (
'THERMAL_FORMULATION')
1371 call this%parser%GetStringCaps(keyword2)
1372 if (trim(adjustl(keyword2)) ==
'LINEAR') this%thermivisc = 1
1373 if (trim(adjustl(keyword2)) ==
'NONLINEAR') this%thermivisc = 2
1374 select case (this%thermivisc)
1376 write (this%iout, fmtlinear)
1378 write (this%iout, fmtnonlinear)
1381 this%a2 = this%parser%GetDouble()
1382 if (this%thermivisc == 2)
then
1383 write (this%iout,
'(4x,a,1pg15.6)') &
1384 'A2 in nonlinear viscosity formulation has been set to: ', &
1387 write (this%iout,
'(4x,a,/,4x,a,/,4x,a)')
'THERMAL_A2 specified by user &
1388 &in VSC Package input file. LINEAR viscosity ',
'formulation also &
1389 &specified. THERMAL_A2 will not affect ',
'viscosity calculations.'
1392 this%a3 = this%parser%GetDouble()
1393 if (this%thermivisc == 2)
then
1394 write (this%iout,
'(4x,a,1pg15.6)') &
1395 'A3 in nonlinear viscosity formulation has been set to: ', &
1398 write (this%iout,
'(4x,a,/,4x,a,/,4x,a)')
'THERMAL_A3 specified by user &
1399 &in VSC Package input file. LINEAR viscosity ',
'formulation also &
1400 &specified. THERMAL_A3 will not affect ',
'viscosity calculations.'
1403 this%a4 = this%parser%GetDouble()
1404 if (this%thermivisc == 2)
then
1405 write (this%iout,
'(4x,a,1pg15.6)') &
1406 'A4 in nonlinear viscosity formulation has been set to: ', &
1409 write (this%iout,
'(4x,a,/,4x,a,/,4x,a)')
'THERMAL_A4 specified by user &
1410 &in VSC Package input file. LINEAR viscosity ',
'formulation also &
1411 &specified. THERMAL_A4 will not affect ',
'viscosity calculations.'
1414 write (errmsg,
'(a,a)')
'Unknown VSC option: ', &
1417 call this%parser%StoreErrorUnit()
1421 if (this%thermivisc == 1)
then
1422 if (this%a2 == 0.0)
then
1423 write (errmsg,
'(a)')
'LINEAR option selected for varying &
1424 &viscosity with temperature, but A1, a surrogate for &
1425 &dVISC/dT, set equal to 0.0'
1429 if (this%thermivisc > 1)
then
1430 if (this%a2 == 0)
then
1431 write (warnmsg,
'(a)')
'NONLINEAR option selected for &
1432 &varying viscosity with temperature, but A2 set equal to &
1433 &zero which may lead to unintended values for viscosity'
1436 if (this%a3 == 0)
then
1437 write (warnmsg,
'(a)')
'NONLINEAR option selected for &
1438 &varying viscosity with temperature,, but A3 set equal to &
1439 &zero which may lead to unintended values for viscosity'
1442 if (this%a4 == 0)
then
1443 write (warnmsg,
'(a)')
'NONLINEAR option selected for &
1444 &varying viscosity with temperature, BUT A4 SET EQUAL TO &
1445 &zero which may lead to unintended values for viscosity'
1451 write (this%iout,
'(/,1x,a)')
'end of VSC options block'
1464 this%viscref = input_data%viscref
1480 character(len=LENMODELNAME),
intent(in) :: modelname
1481 real(DP),
dimension(:),
pointer :: conc
1482 integer(I4B),
dimension(:),
pointer :: icbund
1483 integer(I4B),
optional,
intent(in) :: istmpr
1490 do i = 1, this%nviscspecies
1491 if (this%cmodelname(i) == modelname)
then
1492 this%modelconc(i)%conc => conc
1493 this%modelconc(i)%icbund => icbund
This module contains the base boundary package.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dep3
real constant 1000
integer(i4b), parameter lenmodelname
maximum length of the model name
real(dp), parameter dhnoflo
real no flow constant
real(dp), parameter dhalf
real constant 1/2
integer(i4b), parameter lenauxname
maximum length of a aux variable
real(dp), parameter dzero
real constant zero
real(dp), parameter dten
real constant 10
integer(i4b), parameter maxcharlen
maximum length of char string
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
subroutine vsc_da(this)
@ brief Deallocate viscosity package memory
subroutine read_dimensions(this)
@ brief Read dimensions
subroutine allocate_arrays(this, nodes)
@ brief Allocate arrays
subroutine vsc_set_changed_at(this, kper, kstp)
Mark K changes as having occurred at (kper, kstp)
subroutine update_k_with_vsc(this)
Appled the viscosity ratio (mu_o/mu) to the hydraulic conductivity.
subroutine vsc_ad_standard_bnd(packobj, hnew, visc, viscref, locelev, locvisc, locconc, dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
advance ghb while accounting for viscosity
subroutine allocate_scalars(this)
@ brief Allocate scalars
subroutine vsc_ad_lak(packobj, visc, viscref, elev, locvisc, locconc, dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
Update lak-related viscosity ratios.
subroutine calc_q_visc(this, cellid, viscratio)
Account for viscosity in the aquiferhorizontal flow barriers.
subroutine, public vsc_cr(vscobj, name_model, inunit, iout)
@ brief Create a new package object
real(dp) function calc_vsc_ratio(viscref, bndvisc)
calculate and return the viscosity ratio
subroutine vsc_calcvisc(this)
Calculate fluid viscosity.
subroutine vsc_ad(this)
@ brief Advance the viscosity package
subroutine vsc_ot_dv(this, idvfl)
@ brief Output viscosity package dependent-variable terms.
real(dp) function calc_bnd_viscosity(n, locvisc, locconc, viscref, dviscdc, cviscref, ctemp, ivisc, a2, a3, a4, auxvar)
@ brief Calculate the boundary viscosity
real(dp) function update_bnd_cond(bndvisc, viscref, spcfdcond)
Apply viscosity to the conductance term.
subroutine vsc_ar_bnd(this, packobj)
Activate viscosity in advanced packages.
subroutine read_options(this)
@ brief Read Options block
subroutine set_packagedata(this, input_data)
Sets package data instead of reading from file.
subroutine vsc_df(this, dis, vsc_input)
@ brief Define viscosity package options and dimensions
subroutine vsc_ad_sfr(packobj, visc, viscref, elev, locvisc, locconc, dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
Update sfr-related viscosity ratios.
subroutine vsc_ar(this, ibound)
@ brief Allocate and read method for viscosity package
subroutine set_options(this, input_data)
Sets options as opposed to reading them from a file.
subroutine vsc_rp(this)
@ brief Read new period data in viscosity package
subroutine vsc_ad_bnd(this, packobj, hnew)
Advance the boundary packages when viscosity is active.
subroutine read_packagedata(this)
@ brief Read data for package
subroutine get_visc_ratio(this, n, m, gwhdn, gwhdm, viscratio)
Calculate the viscosity ratio.
subroutine set_npf_pointers(this)
Set pointers to NPF variables.
subroutine vsc_ad_maw(packobj, visc, viscref, elev, locvisc, locconc, dviscdc, cviscref, ivisc, a2, a3, a4, ctemp)
Update maw-related viscosity ratios.
subroutine set_concentration_pointer(this, modelname, conc, icbund, istmpr)
@ brief Set pointers to concentration(s)
real(dp) function calc_visc(ivisc, viscref, dviscdc, cviscref, conc, a2, a3, a4)
Generic function to calculate changes in fluid viscosity using a linear formulation.
This module defines variable data types.
type(listtype), public basemodellist
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
This module contains the base numerical package type.
This module contains the SFR package methods.
This module contains simulation methods.
subroutine, public store_warning(msg, substring)
Store warning message.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number