55 character(len=LENBUDTXT),
dimension(4) ::
budtxt = & !< text labels for budget terms
60 character(len=LENBUDTXT),
dimension(6) ::
comptxt = & !< text labels for compaction terms
75 character(len=LENLISTLABEL),
pointer :: listlabel => null()
76 character(len=LENMEMPATH),
pointer :: stomempath => null()
78 character(len=LENBOUNDNAME),
dimension(:), &
79 pointer,
contiguous :: boundname => null()
80 character(len=LENAUXNAME),
dimension(:), &
81 pointer,
contiguous :: auxname => null()
83 logical,
pointer :: lhead_based => null()
85 integer(I4B),
pointer :: istounit => null()
86 integer(I4B),
pointer :: istrainib => null()
87 integer(I4B),
pointer :: istrainsk => null()
88 integer(I4B),
pointer :: ioutcomp => null()
89 integer(I4B),
pointer :: ioutcompi => null()
90 integer(I4B),
pointer :: ioutcompe => null()
91 integer(I4B),
pointer :: ioutcompib => null()
92 integer(I4B),
pointer :: ioutcomps => null()
93 integer(I4B),
pointer :: ioutzdisp => null()
94 integer(I4B),
pointer :: ipakcsv => null()
95 integer(I4B),
pointer :: iupdatematprop => null()
96 integer(I4B),
pointer :: istoragec => null()
97 integer(I4B),
pointer :: icellf => null()
98 integer(I4B),
pointer :: ispecified_pcs => null()
99 integer(I4B),
pointer :: ispecified_dbh => null()
100 integer(I4B),
pointer :: inamedbound => null()
101 integer(I4B),
pointer :: iconvchk => null()
102 integer(I4B),
pointer :: naux => null()
103 integer(I4B),
pointer :: ninterbeds => null()
104 integer(I4B),
pointer :: maxsig0 => null()
105 integer(I4B),
pointer :: nbound => null()
106 integer(I4B),
pointer :: iscloc => null()
107 integer(I4B),
pointer :: iauxmultcol => null()
108 integer(I4B),
pointer :: ndelaycells => null()
109 integer(I4B),
pointer :: ndelaybeds => null()
110 integer(I4B),
pointer :: initialized => null()
111 integer(I4B),
pointer :: ieslag => null()
112 integer(I4B),
pointer :: ipch => null()
113 integer(I4B),
pointer :: iupdatestress => null()
115 real(dp),
pointer :: epsilon => null()
116 real(dp),
pointer :: cc_crit => null()
117 real(dp),
pointer :: gammaw => null()
118 real(dp),
pointer :: beta => null()
119 real(dp),
pointer :: brg => null()
120 real(dp),
pointer :: satomega => null()
122 integer(I4B),
pointer :: gwfiss => null()
123 integer(I4B),
pointer :: gwfiss0 => null()
125 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound => null()
126 integer(I4B),
dimension(:),
pointer,
contiguous :: stoiconv => null()
128 real(dp),
dimension(:),
pointer,
contiguous :: stoss => null()
129 real(dp),
dimension(:),
pointer,
contiguous :: buff => null()
130 real(dp),
dimension(:),
pointer,
contiguous :: buffusr => null()
131 integer,
dimension(:),
pointer,
contiguous :: nodelist => null()
132 integer,
dimension(:),
pointer,
contiguous :: unodelist => null()
135 real(dp),
dimension(:),
pointer,
contiguous :: sgm => null()
136 real(dp),
dimension(:),
pointer,
contiguous :: sgs => null()
137 real(dp),
dimension(:),
pointer,
contiguous :: cg_ske_cr => null()
138 real(dp),
dimension(:),
pointer,
contiguous :: cg_gs => null()
139 real(dp),
dimension(:),
pointer,
contiguous :: cg_es => null()
140 real(dp),
dimension(:),
pointer,
contiguous :: cg_es0 => null()
141 real(dp),
dimension(:),
pointer,
contiguous :: cg_pcs => null()
142 real(dp),
dimension(:),
pointer,
contiguous :: cg_comp => null()
143 real(dp),
dimension(:),
pointer,
contiguous :: cg_tcomp => null()
144 real(dp),
dimension(:),
pointer,
contiguous :: cg_stor => null()
145 real(dp),
dimension(:),
pointer,
contiguous :: cg_ske => null()
146 real(dp),
dimension(:),
pointer,
contiguous :: cg_sk => null()
147 real(dp),
dimension(:),
pointer,
contiguous :: cg_thickini => null()
148 real(dp),
dimension(:),
pointer,
contiguous :: cg_thetaini => null()
149 real(dp),
dimension(:),
pointer,
contiguous :: cg_thick => null()
150 real(dp),
dimension(:),
pointer,
contiguous :: cg_thick0 => null()
151 real(dp),
dimension(:),
pointer,
contiguous :: cg_theta => null()
152 real(dp),
dimension(:),
pointer,
contiguous :: cg_theta0 => null()
155 real(dp),
dimension(:),
pointer,
contiguous :: cell_wcstor => null()
156 real(dp),
dimension(:),
pointer,
contiguous :: cell_thick => null()
159 integer(I4B),
dimension(:),
pointer,
contiguous :: idelay => null()
160 integer(I4B),
dimension(:),
pointer,
contiguous :: ielastic => null()
161 integer(I4B),
dimension(:),
pointer,
contiguous :: iconvert => null()
162 real(dp),
dimension(:),
pointer,
contiguous :: ci => null()
163 real(dp),
dimension(:),
pointer,
contiguous :: rci => null()
164 real(dp),
dimension(:),
pointer,
contiguous :: pcs => null()
165 real(dp),
dimension(:),
pointer,
contiguous :: rnb => null()
166 real(dp),
dimension(:),
pointer,
contiguous :: kv => null()
167 real(dp),
dimension(:),
pointer,
contiguous :: h0 => null()
168 real(dp),
dimension(:),
pointer,
contiguous :: comp => null()
169 real(dp),
dimension(:),
pointer,
contiguous :: tcomp => null()
170 real(dp),
dimension(:),
pointer,
contiguous :: tcompi => null()
171 real(dp),
dimension(:),
pointer,
contiguous :: tcompe => null()
172 real(dp),
dimension(:),
pointer,
contiguous :: storagee => null()
173 real(dp),
dimension(:),
pointer,
contiguous :: storagei => null()
174 real(dp),
dimension(:),
pointer,
contiguous :: ske => null()
175 real(dp),
dimension(:),
pointer,
contiguous :: sk => null()
176 real(dp),
dimension(:),
pointer,
contiguous :: thickini => null()
177 real(dp),
dimension(:),
pointer,
contiguous :: thetaini => null()
178 real(dp),
dimension(:),
pointer,
contiguous :: thick => null()
179 real(dp),
dimension(:),
pointer,
contiguous :: thick0 => null()
180 real(dp),
dimension(:),
pointer,
contiguous :: theta => null()
181 real(dp),
dimension(:),
pointer,
contiguous :: theta0 => null()
182 real(dp),
dimension(:, :),
pointer,
contiguous :: auxvar => null()
185 integer(I4B),
dimension(:),
pointer,
contiguous :: idb_nconv_count => null()
186 integer(I4B),
dimension(:, :),
pointer,
contiguous :: idbconvert => null()
187 real(dp),
dimension(:),
pointer,
contiguous :: dbdhmax => null()
188 real(dp),
dimension(:, :),
pointer,
contiguous :: dbz => null()
189 real(dp),
dimension(:, :),
pointer,
contiguous :: dbrelz => null()
190 real(dp),
dimension(:, :),
pointer,
contiguous :: dbh => null()
191 real(dp),
dimension(:, :),
pointer,
contiguous :: dbh0 => null()
192 real(dp),
dimension(:, :),
pointer,
contiguous :: dbgeo => null()
193 real(dp),
dimension(:, :),
pointer,
contiguous :: dbes => null()
194 real(dp),
dimension(:, :),
pointer,
contiguous :: dbes0 => null()
195 real(dp),
dimension(:, :),
pointer,
contiguous :: dbpcs => null()
196 real(dp),
dimension(:),
pointer,
contiguous :: dbflowtop => null()
197 real(dp),
dimension(:),
pointer,
contiguous :: dbflowbot => null()
198 real(dp),
dimension(:, :),
pointer,
contiguous :: dbdzini => null()
199 real(dp),
dimension(:, :),
pointer,
contiguous :: dbthetaini => null()
200 real(dp),
dimension(:, :),
pointer,
contiguous :: dbdz => null()
201 real(dp),
dimension(:, :),
pointer,
contiguous :: dbdz0 => null()
202 real(dp),
dimension(:, :),
pointer,
contiguous :: dbtheta => null()
203 real(dp),
dimension(:, :),
pointer,
contiguous :: dbtheta0 => null()
204 real(dp),
dimension(:, :),
pointer,
contiguous :: dbcomp => null()
205 real(dp),
dimension(:, :),
pointer,
contiguous :: dbtcomp => null()
208 real(dp),
dimension(:),
pointer,
contiguous :: dbal => null()
209 real(dp),
dimension(:),
pointer,
contiguous :: dbad => null()
210 real(dp),
dimension(:),
pointer,
contiguous :: dbau => null()
211 real(dp),
dimension(:),
pointer,
contiguous :: dbrhs => null()
212 real(dp),
dimension(:),
pointer,
contiguous :: dbdh => null()
213 real(dp),
dimension(:),
pointer,
contiguous :: dbaw => null()
216 integer(I4B),
dimension(:),
pointer,
contiguous :: nodelistsig0 => null()
217 real(dp),
dimension(:),
pointer,
contiguous :: sig0 => null()
223 integer(I4B),
pointer :: inobspkg => null()
321 subroutine csub_cr(csubobj, name_model, istounit, stoPckName, inunit, iout)
324 character(len=*),
intent(in) :: name_model
325 integer(I4B),
intent(in) :: inunit
326 integer(I4B),
intent(in) :: istounit
327 character(len=*),
intent(in) :: stopckname
328 integer(I4B),
intent(in) :: iout
335 call csubobj%set_names(1, name_model,
'CSUB',
'CSUB')
338 call csubobj%csub_allocate_scalars()
344 csubobj%istounit = istounit
345 csubobj%inunit = inunit
349 call csubobj%parser%Initialize(csubobj%inunit, csubobj%iout)
368 integer(I4B),
dimension(:),
pointer,
contiguous :: ibound
370 logical :: isfound, endOfBlock
371 character(len=:),
allocatable :: line
372 character(len=LINELENGTH) :: keyword
373 character(len=20) :: cellid
375 integer(I4B) :: istheta
378 integer(I4B) :: idelay
381 integer(I4B) :: istart
382 integer(I4B) :: istop
385 integer(I4B) :: istoerr
389 real(DP) :: cg_ske_cr
393 character(len=*),
parameter :: fmtcsub = &
394 "(1x,/1x,'CSUB -- COMPACTION PACKAGE, VERSION 1, 12/15/2019', &
395 &' INPUT READ FROM UNIT ', i0, //)"
398 write (this%iout, fmtcsub) this%inunit
402 this%ibound => ibound
408 call obs_cr(this%obs, this%inobspkg)
411 call this%read_options()
415 call this%tsmanager%tsmanager_df()
418 call this%read_dimensions()
421 call this%obs%obs_ar()
425 call this%parser%StoreErrorUnit()
429 call this%csub_allocate_arrays()
438 call this%parser%GetBlock(
'GRIDDATA', isfound, ierr)
441 call this%parser%GetNextLine(endofblock)
443 call this%parser%GetStringCaps(keyword)
444 call this%parser%GetRemainingLine(line)
446 select case (keyword)
448 call this%dis%read_grid_array(line, lloc, istart, istop, &
449 this%iout, this%parser%iuactive, &
450 this%cg_ske_cr,
'CG_SKE_CR')
453 call this%dis%read_grid_array(line, lloc, istart, istop, &
454 this%iout, this%parser%iuactive, &
455 this%cg_thetaini,
'CG_THETA')
458 call this%dis%read_grid_array(line, lloc, istart, istop, &
459 this%iout, this%parser%iuactive, &
463 call this%dis%read_grid_array(line, lloc, istart, istop, &
464 this%iout, this%parser%iuactive, &
468 write (
errmsg,
'(a,1x,a,a)') &
469 "Unknown GRIDDATA tag '", trim(keyword),
"'."
474 call store_error(
'Required GRIDDATA block not found.')
479 write (
errmsg,
'(a)')
'CG_SKE GRIDDATA must be specified.'
482 if (istheta == 0)
then
483 write (
errmsg,
'(a)')
'CG_THETA GRIDDATA must be specified.'
489 do node = 1, this%dis%nodes
490 this%sgm(node) = 1.7d0
494 do node = 1, this%dis%nodes
495 this%sgs(node) = 2.0d0
503 do node = 1, this%dis%nodes
504 call this%dis%noder_to_string(node, cellid)
505 cg_ske_cr = this%cg_ske_cr(node)
506 theta = this%cg_thetaini(node)
509 if (cg_ske_cr < dzero)
then
510 write (
errmsg,
'(a,g0,a,1x,a,1x,a,a)') &
511 'Coarse-grained material CG_SKE_CR (', cg_ske_cr,
') is less', &
512 'than zero in cell', trim(adjustl(cellid)),
'.'
516 if (this%stoss(node) /= dzero)
then
521 if (theta > done .or. theta < dzero)
then
522 write (
errmsg,
'(a,g0,a,1x,a,1x,a,a)') &
523 'Coarse-grained material THETA (', theta,
') is less', &
524 'than zero or greater than 1 in cell', trim(adjustl(cellid)),
'.'
530 if (istoerr /= 0)
then
531 write (
errmsg,
'(a,3(1x,a))') &
532 'Specific storage values in the storage (STO) package must', &
533 'be zero in all active cells when using the', &
534 trim(adjustl(this%packName)), &
540 if (this%ninterbeds > 0)
then
541 call this%csub_read_packagedata()
545 do node = 1, this%dis%nodes
546 top = this%dis%top(node)
547 bot = this%dis%bot(node)
548 this%cg_thickini(node) = top - bot
549 this%cell_thick(node) = top - bot
553 do ib = 1, this%ninterbeds
554 node = this%nodelist(ib)
555 idelay = this%idelay(ib)
556 if (idelay == 0)
then
557 v = this%thickini(ib)
559 v = this%rnb(ib) * this%thickini(ib)
561 this%cg_thickini(node) = this%cg_thickini(node) - v
565 do node = 1, this%dis%nodes
566 thick = this%cg_thickini(node)
567 if (thick < dzero)
then
568 call this%dis%noder_to_string(node, cellid)
569 write (
errmsg,
'(a,g0,a,1x,a,a)') &
570 'Aquifer thickness is less than zero (', &
571 thick,
') in cell', trim(adjustl(cellid)),
'.'
578 call this%parser%StoreErrorUnit()
584 if (this%iupdatematprop /= 0)
then
585 do node = 1, this%dis%nodes
586 this%cg_thick(node) = this%cg_thickini(node)
587 this%cg_theta(node) = this%cg_thetaini(node)
609 character(len=LINELENGTH) :: keyword
610 character(len=:),
allocatable :: line
611 character(len=MAXCHARLEN) :: fname
612 character(len=LENAUXNAME),
dimension(:),
allocatable :: caux
614 logical :: endOfBlock
617 integer(I4B) :: istart
618 integer(I4B) :: istop
620 integer(I4B) :: inobs
622 integer(I4B) :: ieslag
623 integer(I4B) :: isetgamma
625 character(len=*),
parameter :: fmtts = &
626 &
"(4x,'TIME-SERIES DATA WILL BE READ FROM FILE: ',a)"
627 character(len=*),
parameter :: fmtflow = &
628 &
"(4x,'FLOWS WILL BE SAVED TO FILE: ',a,/4x,'OPENED ON UNIT: ',I7)"
629 character(len=*),
parameter :: fmtflow2 = &
630 &
"(4x,'FLOWS WILL BE SAVED TO BUDGET FILE SPECIFIED IN OUTPUT CONTROL')"
631 character(len=*),
parameter :: fmtssessv = &
632 &
"(4x,'USING SSE AND SSV INSTEAD OF CR AND CC.')"
633 character(len=*),
parameter :: fmtoffset = &
634 &
"(4x,'INITIAL_STRESS TREATED AS AN OFFSET.')"
635 character(len=*),
parameter :: fmtopt = &
637 character(len=*),
parameter :: fmtopti = &
639 character(len=*),
parameter :: fmtoptr = &
641 character(len=*),
parameter :: fmtfileout = &
642 "(4x,'CSUB ',1x,a,1x,' WILL BE SAVED TO FILE: ',a,/4x,&
643 &'OPENED ON UNIT: ',I7)"
651 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, blockrequired=.false., &
652 supportopenclose=.true.)
656 write (this%iout,
'(1x,a)')
'PROCESSING CSUB OPTIONS'
658 call this%parser%GetNextLine(endofblock)
662 call this%parser%GetStringCaps(keyword)
663 select case (keyword)
664 case (
'AUX',
'AUXILIARY')
665 call this%parser%GetRemainingLine(line)
667 call urdaux(this%naux, this%parser%iuactive, this%iout, lloc, &
668 istart, istop, caux, line, this%packName)
670 'AUXNAME', this%memoryPath)
672 this%auxname(n) = caux(n)
677 write (this%iout, fmtflow2)
680 write (this%iout,
'(4x,a)') &
681 'LISTS OF '//trim(adjustl(this%packName))//
' CELLS WILL BE PRINTED.'
684 write (this%iout,
'(4x,a)') trim(adjustl(this%packName))// &
685 ' FLOWS WILL BE PRINTED TO LISTING FILE.'
688 write (this%iout,
'(4x,a)') trim(adjustl(this%packName))// &
689 ' BOUNDARIES HAVE NAMES IN LAST COLUMN.'
691 call this%parser%GetStringCaps(keyword)
692 if (trim(adjustl(keyword)) /=
'FILEIN')
then
693 errmsg =
'TS6 keyword must be followed by "FILEIN" '// &
696 call this%parser%StoreErrorUnit()
698 call this%parser%GetString(fname)
699 write (this%iout, fmtts) trim(fname)
700 call this%TsManager%add_tsfile(fname, this%inunit)
702 call this%parser%GetStringCaps(keyword)
703 if (trim(adjustl(keyword)) /=
'FILEIN')
then
704 errmsg =
'OBS6 keyword must be followed by "FILEIN" '// &
708 if (this%obs%active)
then
709 errmsg =
'Multiple OBS6 keywords detected in OPTIONS block. '// &
710 'Only one OBS6 entry allowed for a package.'
713 this%obs%active = .true.
714 call this%parser%GetString(this%obs%inputFilename)
716 call openfile(inobs, this%iout, this%obs%inputFilename,
'OBS')
717 this%obs%inUnitObs = inobs
718 this%inobspkg = inobs
720 call this%obs%obs_df(this%iout, this%packName, this%filtyp, this%dis)
721 call this%csub_df_obs()
725 this%gammaw = this%parser%GetDouble()
728 this%beta = this%parser%GetDouble()
732 this%lhead_based = .true.
733 case (
'INITIAL_PRECONSOLIDATION_HEAD')
736 this%ndelaycells = this%parser%GetInteger()
740 case (
'COMPRESSION_INDICES')
744 case (
'UPDATE_MATERIAL_PROPERTIES')
745 this%iupdatematprop = 1
748 case (
'CELL_FRACTION')
752 case (
'SPECIFIED_INITIAL_INTERBED_STATE')
753 this%ispecified_pcs = 1
754 this%ispecified_dbh = 1
757 case (
'SPECIFIED_INITIAL_PRECONSOLIDATION_STRESS')
758 this%ispecified_pcs = 1
761 case (
'SPECIFIED_INITIAL_DELAY_HEAD')
762 this%ispecified_dbh = 1
765 case (
'EFFECTIVE_STRESS_LAG')
769 case (
'STRAIN_CSV_INTERBED')
770 call this%parser%GetStringCaps(keyword)
771 if (keyword ==
'FILEOUT')
then
772 call this%parser%GetString(fname)
774 call openfile(this%istrainib, this%iout, fname,
'CSV_OUTPUT', &
775 filstat_opt=
'REPLACE', mode_opt=
mnormal)
776 write (this%iout, fmtfileout) &
777 'INTERBED STRAIN CSV', fname, this%istrainib
779 errmsg =
'Optional STRAIN_CSV_INTERBED keyword must be '// &
780 'followed by FILEOUT.'
783 case (
'STRAIN_CSV_COARSE')
784 call this%parser%GetStringCaps(keyword)
785 if (keyword ==
'FILEOUT')
then
786 call this%parser%GetString(fname)
788 call openfile(this%istrainsk, this%iout, fname,
'CSV_OUTPUT', &
789 filstat_opt=
'REPLACE', mode_opt=
mnormal)
790 write (this%iout, fmtfileout) &
791 'COARSE STRAIN CSV', fname, this%istrainsk
793 errmsg =
'Optional STRAIN_CSV_COARSE keyword must be '// &
794 'followed by fileout.'
800 call this%parser%GetStringCaps(keyword)
801 if (keyword ==
'FILEOUT')
then
802 call this%parser%GetString(fname)
804 call openfile(this%ioutcomp, this%iout, fname,
'DATA(BINARY)', &
806 write (this%iout, fmtfileout) &
807 'COMPACTION', fname, this%ioutcomp
809 errmsg =
'Optional COMPACTION keyword must be '// &
810 'followed by FILEOUT.'
813 case (
'COMPACTION_INELASTIC')
814 call this%parser%GetStringCaps(keyword)
815 if (keyword ==
'FILEOUT')
then
816 call this%parser%GetString(fname)
818 call openfile(this%ioutcompi, this%iout, fname, &
821 write (this%iout, fmtfileout) &
822 'COMPACTION_INELASTIC', fname, this%ioutcompi
824 errmsg =
'Optional COMPACTION_INELASTIC keyword must be '// &
825 'followed by fileout.'
828 case (
'COMPACTION_ELASTIC')
829 call this%parser%GetStringCaps(keyword)
830 if (keyword ==
'FILEOUT')
then
831 call this%parser%GetString(fname)
833 call openfile(this%ioutcompe, this%iout, fname, &
836 write (this%iout, fmtfileout) &
837 'COMPACTION_ELASTIC', fname, this%ioutcompe
839 errmsg =
'Optional COMPACTION_ELASTIC keyword must be '// &
840 'followed by FILEOUT.'
843 case (
'COMPACTION_INTERBED')
844 call this%parser%GetStringCaps(keyword)
845 if (keyword ==
'FILEOUT')
then
846 call this%parser%GetString(fname)
848 call openfile(this%ioutcompib, this%iout, fname, &
851 write (this%iout, fmtfileout) &
852 'COMPACTION_INTERBED', fname, this%ioutcompib
854 errmsg =
'Optional COMPACTION_INTERBED keyword must be '// &
855 'followed by FILEOUT.'
858 case (
'COMPACTION_COARSE')
859 call this%parser%GetStringCaps(keyword)
860 if (keyword ==
'FILEOUT')
then
861 call this%parser%GetString(fname)
863 call openfile(this%ioutcomps, this%iout, fname, &
866 write (this%iout, fmtfileout) &
867 'COMPACTION_COARSE', fname, this%ioutcomps
869 errmsg =
'Optional COMPACTION_COARSE keyword must be '// &
870 'followed by FILEOUT.'
875 case (
'ZDISPLACEMENT')
876 call this%parser%GetStringCaps(keyword)
877 if (keyword ==
'FILEOUT')
then
878 call this%parser%GetString(fname)
880 call openfile(this%ioutzdisp, this%iout, fname, &
883 write (this%iout, fmtfileout) &
884 'ZDISPLACEMENT', fname, this%ioutzdisp
886 errmsg =
'Optional ZDISPLACEMENT keyword must be '// &
887 'followed by FILEOUT.'
891 case (
'PACKAGE_CONVERGENCE')
892 call this%parser%GetStringCaps(keyword)
893 if (keyword ==
'FILEOUT')
then
894 call this%parser%GetString(fname)
896 call openfile(this%ipakcsv, this%iout, fname,
'CSV', &
897 filstat_opt=
'REPLACE', mode_opt=
mnormal)
898 write (this%iout, fmtfileout) &
899 'PACKAGE_CONVERGENCE', fname, this%ipakcsv
901 call store_error(
'Optional PACKAGE_CONVERGENCE keyword must be '// &
902 'followed by FILEOUT.')
909 case (
'DEV_NO_FINAL_CHECK')
910 call this%parser%DevOpt()
912 write (this%iout,
'(4x,a)') &
913 'A FINAL CONVERGENCE CHECK OF THE CHANGE IN DELAY INTERBED '// &
914 'HEADS AND FLOWS WILL NOT BE MADE'
919 write (
errmsg,
'(a,3(1x,a),a)') &
920 'Unknown', trim(adjustl(this%packName)),
"option '", &
925 write (this%iout,
'(1x,a)') &
926 'END OF '//trim(adjustl(this%packName))//
' OPTIONS'
930 write (this%iout,
'(//2(1X,A))') trim(adjustl(this%packName)), &
932 write (this%iout, fmtopti)
'NUMBER OF DELAY CELLS =', &
934 if (this%lhead_based .EQV. .true.)
then
935 write (this%iout,
'(4x,a)') &
936 'HEAD-BASED FORMULATION'
938 write (this%iout,
'(4x,a)') &
939 'EFFECTIVE-STRESS FORMULATION'
941 if (this%istoragec == 0)
then
942 write (this%iout,
'(4x,a,1(/,6x,a))') &
943 'COMPRESSION INDICES WILL BE SPECIFIED INSTEAD OF ELASTIC AND', &
944 'INELASTIC SPECIFIC STORAGE COEFFICIENTS'
946 write (this%iout,
'(4x,a,1(/,6x,a))') &
947 'ELASTIC AND INELASTIC SPECIFIC STORAGE COEFFICIENTS WILL BE ', &
950 if (this%iupdatematprop /= 1)
then
951 write (this%iout,
'(4x,a,1(/,6x,a))') &
952 'THICKNESS AND VOID RATIO WILL NOT BE ADJUSTED DURING THE', &
955 write (this%iout,
'(4x,a)') &
956 'THICKNESS AND VOID RATIO WILL BE ADJUSTED DURING THE SIMULATION'
958 if (this%icellf /= 1)
then
959 write (this%iout,
'(4x,a)') &
960 'INTERBED THICKNESS WILL BE SPECIFIED AS A THICKNESS'
962 write (this%iout,
'(4x,a,1(/,6x,a))') &
963 'INTERBED THICKNESS WILL BE SPECIFIED AS A AS A CELL FRACTION'
965 if (this%ispecified_pcs /= 1)
then
966 if (this%ipch /= 0)
then
967 write (this%iout,
'(4x,a,1(/,6x,a))') &
968 'PRECONSOLIDATION HEAD WILL BE SPECIFIED RELATIVE TO INITIAL', &
971 write (this%iout,
'(4x,a,1(/,6x,a))') &
972 'PRECONSOLIDATION STRESS WILL BE SPECIFIED RELATIVE TO INITIAL', &
976 if (this%ipch /= 0)
then
977 write (this%iout,
'(4x,a,1(/,6x,a))') &
978 'PRECONSOLIDATION HEAD WILL BE SPECIFIED AS ABSOLUTE VALUES', &
979 'INSTEAD OF RELATIVE TO INITIAL HEAD CONDITIONS'
981 write (this%iout,
'(4x,a,1(/,6x,a))') &
982 'PRECONSOLIDATION STRESS WILL BE SPECIFIED AS ABSOLUTE VALUES', &
983 'INSTEAD OF RELATIVE TO INITIAL STRESS CONDITIONS'
986 if (this%ispecified_dbh /= 1)
then
987 write (this%iout,
'(4x,a,1(/,6x,a))') &
988 'DELAY INTERBED HEADS WILL BE SPECIFIED RELATIVE TO INITIAL ', &
991 write (this%iout,
'(4x,a,1(/,6x,a))') &
992 'DELAY INTERBED HEADS WILL BE SPECIFIED AS ABSOLUTE VALUES INSTEAD', &
993 'OF RELATIVE TO INITIAL GWF HEADS'
997 if (this%lhead_based .EQV. .false.)
then
998 if (ieslag /= 0)
then
999 write (this%iout,
'(4x,a,1(/,6x,a))') &
1000 'SPECIFIC STORAGE VALUES WILL BE CALCULATED USING THE EFFECTIVE', &
1001 'STRESS FROM THE PREVIOUS TIME STEP'
1003 write (this%iout,
'(4x,a,1(/,6x,a))') &
1004 'SPECIFIC STORAGE VALUES WILL BE CALCULATED USING THE CURRENT', &
1008 if (ieslag /= 0)
then
1010 write (this%iout,
'(4x,a,2(/,6x,a))') &
1011 'EFFECTIVE_STRESS_LAG HAS BEEN SPECIFIED BUT HAS NO EFFECT WHEN', &
1012 'USING THE HEAD-BASED FORMULATION (HEAD_BASED HAS BEEN SPECIFIED', &
1013 'IN THE OPTIONS BLOCK)'
1016 this%ieslag = ieslag
1021 this%brg = this%gammaw * this%beta
1023 write (this%iout, fmtoptr)
'GAMMAW =', this%gammaw
1024 write (this%iout, fmtoptr)
'BETA =', this%beta
1025 write (this%iout, fmtoptr)
'GAMMAW * BETA =', this%brg
1029 call this%parser%StoreErrorUnit()
1049 character(len=LENBOUNDNAME) :: keyword
1050 integer(I4B) :: ierr
1051 logical :: isfound, endOfBlock
1055 this%ninterbeds = -1
1058 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
1059 supportopenclose=.true.)
1063 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%packName))// &
1066 call this%parser%GetNextLine(endofblock)
1067 if (endofblock)
exit
1068 call this%parser%GetStringCaps(keyword)
1069 select case (keyword)
1071 this%ninterbeds = this%parser%GetInteger()
1072 write (this%iout,
'(4x,a,i0)')
'NINTERBEDS = ', this%ninterbeds
1074 this%maxsig0 = this%parser%GetInteger()
1075 write (this%iout,
'(4x,a,i0)')
'MAXSIG0 = ', this%maxsig0
1077 write (
errmsg,
'(a,3(1x,a),a)') &
1078 'Unknown', trim(this%packName),
"dimension '", trim(keyword),
"'."
1082 write (this%iout,
'(1x,a)') &
1083 'END OF '//trim(adjustl(this%packName))//
' DIMENSIONS'
1085 call store_error(
'Required dimensions block not found.')
1089 if (this%ninterbeds < 0)
then
1091 'NINTERBEDS was not specified or was specified incorrectly.'
1097 call this%parser%StoreErrorUnit()
1102 call this%define_listlabel()
1121 call this%NumericalPackageType%allocate_scalars()
1128 call mem_allocate(this%istounit,
'ISTOUNIT', this%memoryPath)
1129 call mem_allocate(this%inobspkg,
'INOBSPKG', this%memoryPath)
1130 call mem_allocate(this%ninterbeds,
'NINTERBEDS', this%memoryPath)
1131 call mem_allocate(this%maxsig0,
'MAXSIG0', this%memoryPath)
1132 call mem_allocate(this%nbound,
'NBOUND', this%memoryPath)
1133 call mem_allocate(this%iscloc,
'ISCLOC', this%memoryPath)
1134 call mem_allocate(this%iauxmultcol,
'IAUXMULTCOL', this%memoryPath)
1135 call mem_allocate(this%ndelaycells,
'NDELAYCELLS', this%memoryPath)
1136 call mem_allocate(this%ndelaybeds,
'NDELAYBEDS', this%memoryPath)
1137 call mem_allocate(this%initialized,
'INITIALIZED', this%memoryPath)
1138 call mem_allocate(this%ieslag,
'IESLAG', this%memoryPath)
1140 call mem_allocate(this%lhead_based,
'LHEAD_BASED', this%memoryPath)
1141 call mem_allocate(this%iupdatestress,
'IUPDATESTRESS', this%memoryPath)
1142 call mem_allocate(this%ispecified_pcs,
'ISPECIFIED_PCS', this%memoryPath)
1143 call mem_allocate(this%ispecified_dbh,
'ISPECIFIED_DBH', this%memoryPath)
1144 call mem_allocate(this%inamedbound,
'INAMEDBOUND', this%memoryPath)
1145 call mem_allocate(this%iconvchk,
'ICONVCHK', this%memoryPath)
1147 call mem_allocate(this%istoragec,
'ISTORAGEC', this%memoryPath)
1148 call mem_allocate(this%istrainib,
'ISTRAINIB', this%memoryPath)
1149 call mem_allocate(this%istrainsk,
'ISTRAINSK', this%memoryPath)
1150 call mem_allocate(this%ioutcomp,
'IOUTCOMP', this%memoryPath)
1151 call mem_allocate(this%ioutcompi,
'IOUTCOMPI', this%memoryPath)
1152 call mem_allocate(this%ioutcompe,
'IOUTCOMPE', this%memoryPath)
1153 call mem_allocate(this%ioutcompib,
'IOUTCOMPIB', this%memoryPath)
1154 call mem_allocate(this%ioutcomps,
'IOUTCOMPS', this%memoryPath)
1155 call mem_allocate(this%ioutzdisp,
'IOUTZDISP', this%memoryPath)
1156 call mem_allocate(this%ipakcsv,
'IPAKCSV', this%memoryPath)
1157 call mem_allocate(this%iupdatematprop,
'IUPDATEMATPROP', this%memoryPath)
1158 call mem_allocate(this%epsilon,
'EPSILON', this%memoryPath)
1159 call mem_allocate(this%cc_crit,
'CC_CRIT', this%memoryPath)
1160 call mem_allocate(this%gammaw,
'GAMMAW', this%memoryPath)
1163 call mem_allocate(this%satomega,
'SATOMEGA', this%memoryPath)
1164 call mem_allocate(this%icellf,
'ICELLF', this%memoryPath)
1165 call mem_allocate(this%gwfiss0,
'GWFISS0', this%memoryPath)
1168 allocate (this%TsManager)
1180 this%iauxmultcol = 0
1181 this%ndelaycells = 19
1183 this%initialized = 0
1186 this%lhead_based = .false.
1187 this%iupdatestress = 1
1188 this%ispecified_pcs = 0
1189 this%ispecified_dbh = 0
1190 this%inamedbound = 0
1203 this%iupdatematprop = 0
1204 this%epsilon =
dzero
1207 this%beta = 4.6512e-10_dp
1208 this%brg = this%gammaw * this%beta
1211 if (this%inewton /= 0)
then
1212 this%satomega =
dem6
1215 this%satomega =
dzero
1238 integer(I4B) :: iblen
1239 integer(I4B) :: ilen
1240 integer(I4B) :: naux
1243 if (this%ioutcomp == 0 .and. this%ioutcompi == 0 .and. &
1244 this%ioutcompe == 0 .and. this%ioutcompib == 0 .and. &
1245 this%ioutcomps == 0 .and. this%ioutzdisp == 0)
then
1246 call mem_allocate(this%buff, 1,
'BUFF', trim(this%memoryPath))
1248 call mem_allocate(this%buff, this%dis%nodes,
'BUFF', trim(this%memoryPath))
1250 if (this%ioutcomp == 0 .and. this%ioutzdisp == 0)
then
1251 call mem_allocate(this%buffusr, 1,
'BUFFUSR', trim(this%memoryPath))
1253 call mem_allocate(this%buffusr, this%dis%nodesuser,
'BUFFUSR', &
1254 trim(this%memoryPath))
1256 call mem_allocate(this%sgm, this%dis%nodes,
'SGM', trim(this%memoryPath))
1257 call mem_allocate(this%sgs, this%dis%nodes,
'SGS', trim(this%memoryPath))
1258 call mem_allocate(this%cg_ske_cr, this%dis%nodes,
'CG_SKE_CR', &
1259 trim(this%memoryPath))
1260 call mem_allocate(this%cg_es, this%dis%nodes,
'CG_ES', &
1261 trim(this%memoryPath))
1262 call mem_allocate(this%cg_es0, this%dis%nodes,
'CG_ES0', &
1263 trim(this%memoryPath))
1264 call mem_allocate(this%cg_pcs, this%dis%nodes,
'CG_PCS', &
1265 trim(this%memoryPath))
1266 call mem_allocate(this%cg_comp, this%dis%nodes,
'CG_COMP', &
1267 trim(this%memoryPath))
1268 call mem_allocate(this%cg_tcomp, this%dis%nodes,
'CG_TCOMP', &
1269 trim(this%memoryPath))
1270 call mem_allocate(this%cg_stor, this%dis%nodes,
'CG_STOR', &
1271 trim(this%memoryPath))
1272 call mem_allocate(this%cg_ske, this%dis%nodes,
'CG_SKE', &
1273 trim(this%memoryPath))
1274 call mem_allocate(this%cg_sk, this%dis%nodes,
'CG_SK', &
1275 trim(this%memoryPath))
1276 call mem_allocate(this%cg_thickini, this%dis%nodes,
'CG_THICKINI', &
1277 trim(this%memoryPath))
1278 call mem_allocate(this%cg_thetaini, this%dis%nodes,
'CG_THETAINI', &
1279 trim(this%memoryPath))
1280 if (this%iupdatematprop == 0)
then
1281 call mem_setptr(this%cg_thick,
'CG_THICKINI', trim(this%memoryPath))
1282 call mem_setptr(this%cg_thick0,
'CG_THICKINI', trim(this%memoryPath))
1283 call mem_setptr(this%cg_theta,
'CG_THETAINI', trim(this%memoryPath))
1284 call mem_setptr(this%cg_theta0,
'CG_THETAINI', trim(this%memoryPath))
1286 call mem_allocate(this%cg_thick, this%dis%nodes,
'CG_THICK', &
1287 trim(this%memoryPath))
1288 call mem_allocate(this%cg_thick0, this%dis%nodes,
'CG_THICK0', &
1289 trim(this%memoryPath))
1290 call mem_allocate(this%cg_theta, this%dis%nodes,
'CG_THETA', &
1291 trim(this%memoryPath))
1292 call mem_allocate(this%cg_theta0, this%dis%nodes,
'CG_THETA0', &
1293 trim(this%memoryPath))
1297 call mem_allocate(this%cell_wcstor, this%dis%nodes,
'CELL_WCSTOR', &
1298 trim(this%memoryPath))
1299 call mem_allocate(this%cell_thick, this%dis%nodes,
'CELL_THICK', &
1300 trim(this%memoryPath))
1304 if (this%ninterbeds > 0)
then
1305 iblen = this%ninterbeds
1308 if (this%naux > 0)
then
1311 call mem_allocate(this%auxvar, naux, iblen,
'AUXVAR', this%memoryPath)
1314 this%auxvar(j, n) =
dzero
1317 call mem_allocate(this%unodelist, iblen,
'UNODELIST', trim(this%memoryPath))
1318 call mem_allocate(this%nodelist, iblen,
'NODELIST', trim(this%memoryPath))
1319 call mem_allocate(this%cg_gs, this%dis%nodes,
'CG_GS', trim(this%memoryPath))
1320 call mem_allocate(this%pcs, iblen,
'PCS', trim(this%memoryPath))
1321 call mem_allocate(this%rnb, iblen,
'RNB', trim(this%memoryPath))
1322 call mem_allocate(this%kv, iblen,
'KV', trim(this%memoryPath))
1323 call mem_allocate(this%h0, iblen,
'H0', trim(this%memoryPath))
1324 call mem_allocate(this%ci, iblen,
'CI', trim(this%memoryPath))
1325 call mem_allocate(this%rci, iblen,
'RCI', trim(this%memoryPath))
1326 call mem_allocate(this%idelay, iblen,
'IDELAY', trim(this%memoryPath))
1327 call mem_allocate(this%ielastic, iblen,
'IELASTIC', trim(this%memoryPath))
1328 call mem_allocate(this%iconvert, iblen,
'ICONVERT', trim(this%memoryPath))
1329 call mem_allocate(this%comp, iblen,
'COMP', trim(this%memoryPath))
1330 call mem_allocate(this%tcomp, iblen,
'TCOMP', trim(this%memoryPath))
1331 call mem_allocate(this%tcompi, iblen,
'TCOMPI', trim(this%memoryPath))
1332 call mem_allocate(this%tcompe, iblen,
'TCOMPE', trim(this%memoryPath))
1333 call mem_allocate(this%storagee, iblen,
'STORAGEE', trim(this%memoryPath))
1334 call mem_allocate(this%storagei, iblen,
'STORAGEI', trim(this%memoryPath))
1335 call mem_allocate(this%ske, iblen,
'SKE', trim(this%memoryPath))
1336 call mem_allocate(this%sk, iblen,
'SK', trim(this%memoryPath))
1337 call mem_allocate(this%thickini, iblen,
'THICKINI', trim(this%memoryPath))
1338 call mem_allocate(this%thetaini, iblen,
'THETAINI', trim(this%memoryPath))
1339 if (this%iupdatematprop == 0)
then
1340 call mem_setptr(this%thick,
'THICKINI', trim(this%memoryPath))
1341 call mem_setptr(this%thick0,
'THICKINI', trim(this%memoryPath))
1342 call mem_setptr(this%theta,
'THETAINI', trim(this%memoryPath))
1343 call mem_setptr(this%theta0,
'THETAINI', trim(this%memoryPath))
1345 call mem_allocate(this%thick, iblen,
'THICK', trim(this%memoryPath))
1346 call mem_allocate(this%thick0, iblen,
'THICK0', trim(this%memoryPath))
1347 call mem_allocate(this%theta, iblen,
'THETA', trim(this%memoryPath))
1348 call mem_allocate(this%theta0, iblen,
'THETA0', trim(this%memoryPath))
1355 if (this%inamedbound /= 0)
then
1357 'BOUNDNAME', trim(this%memoryPath))
1360 'BOUNDNAME', trim(this%memoryPath))
1365 if (this%maxsig0 > 0)
then
1370 call mem_allocate(this%nodelistsig0, ilen,
'NODELISTSIG0', this%memoryPath)
1371 call mem_allocate(this%sig0, ilen,
'SIG0', this%memoryPath)
1374 call mem_setptr(this%gwfiss,
'ISS', trim(this%name_model))
1377 call mem_setptr(this%stoiconv,
'ICONVERT', this%stoMemPath)
1378 call mem_setptr(this%stoss,
'SS', this%stoMemPath)
1381 do n = 1, this%dis%nodes
1382 this%cg_gs(n) =
dzero
1383 this%cg_es(n) =
dzero
1384 this%cg_comp(n) =
dzero
1385 this%cg_tcomp(n) =
dzero
1386 this%cell_wcstor(n) =
dzero
1388 do n = 1, this%ninterbeds
1389 this%theta(n) =
dzero
1390 this%tcomp(n) =
dzero
1391 this%tcompi(n) =
dzero
1392 this%tcompe(n) =
dzero
1394 do n = 1, max(1, this%maxsig0)
1395 this%nodelistsig0(n) = 0
1396 this%sig0(n) =
dzero
1417 character(len=LINELENGTH) :: cellid
1418 character(len=LINELENGTH) :: title
1419 character(len=LINELENGTH) :: tag
1420 character(len=20) :: scellid
1421 character(len=10) :: text
1422 character(len=LENBOUNDNAME) :: bndName
1423 character(len=7) :: cdelay
1425 logical :: endOfBlock
1426 integer(I4B) :: ival
1430 integer(I4B) :: itmp
1431 integer(I4B) :: ierr
1432 integer(I4B) :: ndelaybeds
1433 integer(I4B) :: idelay
1434 integer(I4B) :: ntabrows
1435 integer(I4B) :: ntabcols
1441 integer,
allocatable,
dimension(:) :: nboundchk
1447 allocate (nboundchk(this%ninterbeds))
1448 do n = 1, this%ninterbeds
1452 call this%parser%GetBlock(
'PACKAGEDATA', isfound, ierr, &
1453 supportopenclose=.true.)
1457 write (this%iout,
'(/1x,a)')
'PROCESSING '//trim(adjustl(this%packName))// &
1460 call this%parser%GetNextLine(endofblock)
1461 if (endofblock)
then
1466 itmp = this%parser%GetInteger()
1469 if (itmp < 1 .or. itmp > this%ninterbeds)
then
1470 write (
errmsg,
'(a,1x,i0,2(1x,a),1x,i0,a)') &
1471 'Interbed number (', itmp,
') must be greater than 0 and ', &
1472 'less than or equal to', this%ninterbeds,
'.'
1478 nboundchk(itmp) = nboundchk(itmp) + 1
1481 call this%parser%GetCellid(this%dis%ndim, cellid)
1482 nn = this%dis%noder_from_cellid(cellid, &
1483 this%parser%iuactive, this%iout)
1484 n = this%dis%nodeu_from_cellid(cellid, &
1485 this%parser%iuactive, this%iout)
1486 top = this%dis%top(nn)
1487 bot = this%dis%bot(nn)
1492 write (
errmsg,
'(a,1x,i0,a)') &
1493 'Invalid cellid for packagedata entry', itmp,
'.'
1498 this%nodelist(itmp) = nn
1499 this%unodelist(itmp) = n
1502 call this%parser%GetStringCaps(cdelay)
1503 select case (cdelay)
1507 ndelaybeds = ndelaybeds + 1
1510 write (
errmsg,
'(a,1x,a,1x,i0,1x,a)') &
1511 'Invalid CDELAY ', trim(adjustl(cdelay)), &
1512 'for packagedata entry', itmp,
'.'
1517 this%idelay(itmp) = ival
1520 this%pcs(itmp) = this%parser%GetDouble()
1523 rval = this%parser%GetDouble()
1524 if (this%icellf == 0)
then
1525 if (rval < dzero .or. rval > baq)
then
1526 write (
errmsg,
'(a,g0,2(a,1x),g0,1x,a,1x,i0,a)') &
1527 'THICK (', rval,
') MUST BE greater than or equal to 0 ', &
1528 'and less than or equal to than', baq, &
1529 'for packagedata entry', itmp,
'.'
1533 if (rval < dzero .or. rval > done)
then
1534 write (
errmsg,
'(a,1x,a,1x,i0,a)') &
1535 'FRAC MUST BE greater than 0 and less than or equal to 1', &
1536 'for packagedata entry', itmp,
'.'
1541 this%thickini(itmp) = rval
1542 if (this%iupdatematprop /= 0)
then
1543 this%thick(itmp) = rval
1547 rval = this%parser%GetDouble()
1548 if (idelay > 0)
then
1549 if (rval < done)
then
1550 write (
errmsg,
'(a,g0,a,1x,a,1x,i0,a)') &
1551 'RNB (', rval,
') must be greater than or equal to 1', &
1552 'for packagedata entry', itmp,
'.'
1558 this%rnb(itmp) = rval
1561 rval = this%parser%GetDouble()
1562 if (rval < dzero)
then
1563 write (
errmsg,
'(2(a,1x),i0,a)') &
1564 '(SKV,CI) must be greater than or equal to 0', &
1565 'for packagedata entry', itmp,
'.'
1568 this%ci(itmp) = rval
1571 rval = this%parser%GetDouble()
1572 if (rval < dzero)
then
1573 write (
errmsg,
'(2(a,1x),i0,a)') &
1574 '(SKE,RCI) must be greater than or equal to 0', &
1575 'for packagedata entry', itmp,
'.'
1578 this%rci(itmp) = rval
1581 if (this%ci(itmp) == this%rci(itmp))
then
1582 this%ielastic(itmp) = 1
1584 this%ielastic(itmp) = 0
1588 rval = this%parser%GetDouble()
1589 this%thetaini(itmp) = rval
1590 if (this%iupdatematprop /= 0)
then
1591 this%theta(itmp) = rval
1593 if (rval <= dzero .or. rval > done)
then
1594 write (
errmsg,
'(a,1x,a,1x,i0,a)') &
1595 'THETA must be greater than 0 and less than or equal to 1', &
1596 'for packagedata entry', itmp,
'.'
1601 rval = this%parser%GetDouble()
1602 if (idelay > 0)
then
1603 if (rval <= 0.0)
then
1604 write (
errmsg,
'(a,1x,i0,a)') &
1605 'KV must be greater than 0 for packagedata entry', itmp,
'.'
1609 this%kv(itmp) = rval
1612 rval = this%parser%GetDouble()
1613 this%h0(itmp) = rval
1616 if (this%inamedbound /= 0)
then
1617 call this%parser%GetStringCaps(bndname)
1618 if (len_trim(bndname) < 1)
then
1619 write (
errmsg,
'(a,1x,i0,a)') &
1620 'BOUNDNAME must be specified for packagedata entry', itmp,
'.'
1623 this%boundname(itmp) = bndname
1628 write (this%iout,
'(1x,a)') &
1629 'END OF '//trim(adjustl(this%packName))//
' PACKAGEDATA'
1633 if (this%iprpak == 1)
then
1635 title = trim(adjustl(this%packName))//
' PACKAGE INTERBED DATA'
1638 ntabrows = this%ninterbeds
1640 if (this%inamedbound /= 0)
then
1641 ntabcols = ntabcols + 1
1645 call table_cr(this%inputtab, this%packName, title)
1646 call this%inputtab%table_df(ntabrows, ntabcols, this%iout)
1650 call this%inputtab%initialize_column(tag, 10, alignment=tableft)
1652 call this%inputtab%initialize_column(tag, 20, alignment=tabcenter)
1654 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1656 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1658 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1660 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1662 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1664 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1666 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1668 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1670 call this%inputtab%initialize_column(tag, 10, alignment=tabcenter)
1671 if (this%inamedbound /= 0)
then
1673 call this%inputtab%initialize_column(tag, lenboundname, &
1678 do ib = 1, this%ninterbeds
1679 call this%dis%noder_to_string(this%nodelist(ib), scellid)
1680 if (this%idelay(ib) == 0)
then
1685 call this%inputtab%add_term(ib)
1686 call this%inputtab%add_term(scellid)
1687 call this%inputtab%add_term(text)
1688 call this%inputtab%add_term(this%pcs(ib))
1689 call this%inputtab%add_term(this%thickini(ib))
1690 call this%inputtab%add_term(this%rnb(ib))
1691 call this%inputtab%add_term(this%ci(ib))
1692 call this%inputtab%add_term(this%rci(ib))
1693 call this%inputtab%add_term(this%thetaini(ib))
1694 if (this%idelay(ib) == 0)
then
1695 call this%inputtab%add_term(
'-')
1696 call this%inputtab%add_term(
'-')
1698 call this%inputtab%add_term(this%kv(ib))
1699 call this%inputtab%add_term(this%h0(ib))
1701 if (this%inamedbound /= 0)
then
1702 call this%inputtab%add_term(this%boundname(ib))
1709 do ib = 1, this%ninterbeds
1710 if (nboundchk(ib) == 0)
then
1711 write (
errmsg,
'(a,1x,i0,a)') &
1712 'Information for interbed', ib,
'not specified in packagedata block.'
1714 else if (nboundchk(ib) > 1)
then
1715 write (
errmsg,
'(2(a,1x,i0),a)') &
1716 'Information specified', nboundchk(ib),
'times for interbed', ib,
'.'
1720 deallocate (nboundchk)
1723 this%ndelaybeds = ndelaybeds
1726 if (ndelaybeds > 0)
then
1731 'IDB_NCONV_COUNT', trim(this%memoryPath))
1732 call mem_allocate(this%idbconvert, this%ndelaycells, ndelaybeds, &
1733 'IDBCONVERT', trim(this%memoryPath))
1735 'DBDHMAX', trim(this%memoryPath))
1736 call mem_allocate(this%dbz, this%ndelaycells, ndelaybeds, &
1737 'DBZ', trim(this%memoryPath))
1738 call mem_allocate(this%dbrelz, this%ndelaycells, ndelaybeds, &
1739 'DBRELZ', trim(this%memoryPath))
1740 call mem_allocate(this%dbh, this%ndelaycells, ndelaybeds, &
1741 'DBH', trim(this%memoryPath))
1742 call mem_allocate(this%dbh0, this%ndelaycells, ndelaybeds, &
1743 'DBH0', trim(this%memoryPath))
1744 call mem_allocate(this%dbgeo, this%ndelaycells, ndelaybeds, &
1745 'DBGEO', trim(this%memoryPath))
1746 call mem_allocate(this%dbes, this%ndelaycells, ndelaybeds, &
1747 'DBES', trim(this%memoryPath))
1748 call mem_allocate(this%dbes0, this%ndelaycells, ndelaybeds, &
1749 'DBES0', trim(this%memoryPath))
1750 call mem_allocate(this%dbpcs, this%ndelaycells, ndelaybeds, &
1751 'DBPCS', trim(this%memoryPath))
1753 'DBFLOWTOP', trim(this%memoryPath))
1755 'DBFLOWBOT', trim(this%memoryPath))
1756 call mem_allocate(this%dbdzini, this%ndelaycells, ndelaybeds, &
1757 'DBDZINI', trim(this%memoryPath))
1758 call mem_allocate(this%dbthetaini, this%ndelaycells, ndelaybeds, &
1759 'DBTHETAINI', trim(this%memoryPath))
1760 call mem_allocate(this%dbcomp, this%ndelaycells, ndelaybeds, &
1761 'DBCOMP', trim(this%memoryPath))
1762 call mem_allocate(this%dbtcomp, this%ndelaycells, ndelaybeds, &
1763 'DBTCOMP', trim(this%memoryPath))
1766 if (this%iupdatematprop == 0)
then
1767 call mem_setptr(this%dbdz,
'DBDZINI', trim(this%memoryPath))
1768 call mem_setptr(this%dbdz0,
'DBDZINI', trim(this%memoryPath))
1769 call mem_setptr(this%dbtheta,
'DBTHETAINI', trim(this%memoryPath))
1770 call mem_setptr(this%dbtheta0,
'DBTHETAINI', trim(this%memoryPath))
1772 call mem_allocate(this%dbdz, this%ndelaycells, ndelaybeds, &
1773 'DBDZ', trim(this%memoryPath))
1774 call mem_allocate(this%dbdz0, this%ndelaycells, ndelaybeds, &
1775 'DBDZ0', trim(this%memoryPath))
1776 call mem_allocate(this%dbtheta, this%ndelaycells, ndelaybeds, &
1777 'DBTHETA', trim(this%memoryPath))
1778 call mem_allocate(this%dbtheta0, this%ndelaycells, ndelaybeds, &
1779 'DBTHETA0', trim(this%memoryPath))
1784 'DBAL', trim(this%memoryPath))
1786 'DBAD', trim(this%memoryPath))
1788 'DBAU', trim(this%memoryPath))
1790 'DBRHS', trim(this%memoryPath))
1792 'DBDH', trim(this%memoryPath))
1794 'DBAW', trim(this%memoryPath))
1798 this%idb_nconv_count(n) = 0
1802 do ib = 1, this%ninterbeds
1803 idelay = this%idelay(ib)
1804 if (idelay == 0)
then
1809 do n = 1, this%ndelaycells
1810 rval = this%thickini(ib) / real(this%ndelaycells, dp)
1811 this%dbdzini(n, idelay) = rval
1812 this%dbh(n, idelay) = this%h0(ib)
1813 this%dbh0(n, idelay) = this%h0(ib)
1814 this%dbthetaini(n, idelay) = this%thetaini(ib)
1815 this%dbgeo(n, idelay) = dzero
1816 this%dbes(n, idelay) = dzero
1817 this%dbes0(n, idelay) = dzero
1818 this%dbpcs(n, idelay) = this%pcs(ib)
1819 this%dbcomp(n, idelay) = dzero
1820 this%dbtcomp(n, idelay) = dzero
1821 if (this%iupdatematprop /= 0)
then
1822 this%dbdz(n, idelay) = this%dbdzini(n, idelay)
1823 this%dbdz0(n, idelay) = this%dbdzini(n, idelay)
1824 this%dbtheta(n, idelay) = this%theta(ib)
1825 this%dbtheta0(n, idelay) = this%theta(ib)
1830 call this%csub_delay_init_zcell(ib)
1834 do n = 1, this%ndelaycells
1835 this%dbal(n) = dzero
1836 this%dbad(n) = dzero
1837 this%dbau(n) = dzero
1838 this%dbrhs(n) = dzero
1839 this%dbdh(n) = dzero
1840 this%dbaw(n) = dzero
1847 if (ndelaybeds > 0)
then
1848 q = mod(real(this%ndelaycells, dp), dtwo)
1849 if (q == dzero)
then
1850 write (
errmsg,
'(a,i0,a,1x,a)') &
1851 'NDELAYCELLS (', this%ndelaycells,
') must be an', &
1852 'odd number when using the effective stress formulation.'
1872 character(len=LINELENGTH) :: title
1873 character(len=LINELENGTH) :: tag
1874 character(len=LINELENGTH) :: msg
1875 character(len=10) :: ctype
1876 character(len=20) :: cellid
1877 character(len=10) :: cflag
1882 integer(I4B) :: node
1884 integer(I4B) :: idelay
1885 integer(I4B) :: iexceed
1886 integer(I4B),
parameter :: ncells = 20
1887 integer(I4B) :: nlen
1888 integer(I4B) :: ntabrows
1889 integer(I4B) :: ntabcols
1890 integer(I4B) :: ipos
1895 integer(I4B),
dimension(:),
allocatable :: imap_sel
1896 integer(I4B),
dimension(:),
allocatable :: locs
1897 real(DP),
dimension(:),
allocatable :: pctcomp_arr
1900 allocate (locs(this%dis%ndim))
1903 if (this%ninterbeds > 0)
then
1904 nlen = min(ncells, this%ninterbeds)
1905 allocate (imap_sel(nlen))
1906 allocate (pctcomp_arr(this%ninterbeds))
1908 do ib = 1, this%ninterbeds
1909 idelay = this%idelay(ib)
1910 b0 = this%thickini(ib)
1911 strain = this%tcomp(ib) / b0
1913 pctcomp_arr(ib) = pctcomp
1914 if (pctcomp >=
done)
then
1915 iexceed = iexceed + 1
1918 call selectn(imap_sel, pctcomp_arr, reverse=.true.)
1921 i0 = max(1, this%ninterbeds - ncells + 1)
1922 i1 = this%ninterbeds
1924 if (iexceed /= 0)
then
1925 write (msg,
'(1x,a,1x,i0,1x,a,1x,i0,1x,a)') &
1926 'LARGEST', (i1 - i0 + 1),
'OF', this%ninterbeds, &
1927 'INTERBED STRAIN VALUES SHOWN'
1932 title = trim(adjustl(this%packName))//
' PACKAGE INTERBED STRAIN SUMMARY'
1939 call table_cr(this%outputtab, this%packName, title)
1940 call this%outputtab%table_df(ntabrows, ntabcols, this%iout)
1943 tag =
'INTERBED NUMBER'
1944 call this%outputtab%initialize_column(tag, 10, alignment=
tabcenter)
1945 tag =
'INTERBED TYPE'
1946 call this%outputtab%initialize_column(tag, 10, alignment=
tabcenter)
1948 call this%outputtab%initialize_column(tag, 20, alignment=
tableft)
1949 tag =
'INITIAL THICKNESS'
1950 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
1951 tag =
'FINAL THICKNESS'
1952 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
1953 tag =
'TOTAL COMPACTION'
1954 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
1955 tag =
'FINAL STRAIN'
1956 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
1957 tag =
'PERCENT COMPACTION'
1958 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
1960 call this%outputtab%initialize_column(tag, 10, alignment=
tabcenter)
1965 idelay = this%idelay(ib)
1966 b0 = this%thickini(ib)
1967 b1 = this%csub_calc_interbed_thickness(ib)
1968 if (idelay == 0)
then
1972 b0 = b0 * this%rnb(ib)
1974 strain = this%tcomp(ib) / b0
1976 if (pctcomp >= 5.0_dp)
then
1978 else if (pctcomp >=
done)
then
1983 node = this%nodelist(ib)
1984 call this%dis%noder_to_string(node, cellid)
1987 call this%outputtab%add_term(ib)
1988 call this%outputtab%add_term(ctype)
1989 call this%outputtab%add_term(cellid)
1990 call this%outputtab%add_term(b0)
1991 call this%outputtab%add_term(b1)
1992 call this%outputtab%add_term(this%tcomp(ib))
1993 call this%outputtab%add_term(strain)
1994 call this%outputtab%add_term(pctcomp)
1995 call this%outputtab%add_term(cflag)
1997 write (this%iout,
'(/1X,A,1X,I0,1X,A,1X,I0,1X,A,/1X,A,/1X,A)') &
1998 'PERCENT COMPACTION IS GREATER THAN OR EQUAL TO 1 PERCENT IN', &
1999 iexceed,
'OF', this%ninterbeds,
'INTERBED(S).', &
2000 'USE THE STRAIN_CSV_INTERBED OPTION TO OUTPUT A CSV '// &
2001 'FILE WITH PERCENT COMPACTION ',
'VALUES FOR ALL INTERBEDS.'
2003 msg =
'PERCENT COMPACTION WAS LESS THAN 1 PERCENT IN ALL INTERBEDS'
2004 write (this%iout,
'(/1X,A)') trim(adjustl(msg))
2008 if (this%istrainib /= 0)
then
2011 ntabrows = this%ninterbeds
2013 if (this%dis%ndim > 1)
then
2014 ntabcols = ntabcols + 1
2016 ntabcols = ntabcols + this%dis%ndim
2019 call table_cr(this%outputtab, this%packName,
'')
2020 call this%outputtab%table_df(ntabrows, ntabcols, this%istrainib, &
2021 lineseparator=.false., separator=
',')
2024 tag =
'INTERBED_NUMBER'
2025 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2026 tag =
'INTERBED_TYPE'
2027 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2029 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2030 if (this%dis%ndim == 2)
then
2032 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2034 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2037 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2039 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2041 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2043 tag =
'INITIAL_THICKNESS'
2044 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2045 tag =
'FINAL_THICKNESS'
2046 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2047 tag =
'TOTAL_COMPACTION'
2048 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2049 tag =
'TOTAL_STRAIN'
2050 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2051 tag =
'PERCENT_COMPACTION'
2052 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2055 do ib = 1, this%ninterbeds
2056 idelay = this%idelay(ib)
2057 b0 = this%thickini(ib)
2058 b1 = this%csub_calc_interbed_thickness(ib)
2059 if (idelay == 0)
then
2063 b0 = b0 * this%rnb(ib)
2065 strain = this%tcomp(ib) / b0
2067 node = this%nodelist(ib)
2068 call this%dis%noder_to_array(node, locs)
2071 call this%outputtab%add_term(ib)
2072 call this%outputtab%add_term(ctype)
2073 if (this%dis%ndim > 1)
then
2074 call this%outputtab%add_term(this%dis%get_nodeuser(node))
2076 do ipos = 1, this%dis%ndim
2077 call this%outputtab%add_term(locs(ipos))
2079 call this%outputtab%add_term(b0)
2080 call this%outputtab%add_term(b1)
2081 call this%outputtab%add_term(this%tcomp(ib))
2082 call this%outputtab%add_term(strain)
2083 call this%outputtab%add_term(pctcomp)
2088 deallocate (imap_sel)
2089 deallocate (pctcomp_arr)
2093 nlen = min(ncells, this%dis%nodes)
2094 allocate (imap_sel(nlen))
2095 allocate (pctcomp_arr(this%dis%nodes))
2097 do node = 1, this%dis%nodes
2099 if (this%cg_thickini(node) >
dzero)
then
2100 strain = this%cg_tcomp(node) / this%cg_thickini(node)
2103 pctcomp_arr(node) = pctcomp
2104 if (pctcomp >=
done)
then
2105 iexceed = iexceed + 1
2108 call selectn(imap_sel, pctcomp_arr, reverse=.true.)
2111 i0 = max(1, this%dis%nodes - ncells + 1)
2114 if (iexceed /= 0)
then
2115 write (msg,
'(a,1x,i0,1x,a,1x,i0,1x,a)') &
2116 'LARGEST ', (i1 - i0 + 1),
'OF', this%dis%nodes, &
2117 'CELL COARSE-GRAINED VALUES SHOWN'
2121 title = trim(adjustl(this%packName))// &
2122 ' PACKAGE COARSE-GRAINED STRAIN SUMMARY'
2129 call table_cr(this%outputtab, this%packName, title)
2130 call this%outputtab%table_df(ntabrows, ntabcols, this%iout)
2134 call this%outputtab%initialize_column(tag, 20, alignment=
tableft)
2135 tag =
'INITIAL THICKNESS'
2136 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
2137 tag =
'FINAL THICKNESS'
2138 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
2139 tag =
'TOTAL COMPACTION'
2140 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
2141 tag =
'FINAL STRAIN'
2142 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
2143 tag =
'PERCENT COMPACTION'
2144 call this%outputtab%initialize_column(tag, 12, alignment=
tabcenter)
2146 call this%outputtab%initialize_column(tag, 10, alignment=
tabcenter)
2150 if (this%cg_thickini(node) >
dzero)
then
2151 strain = this%cg_tcomp(node) / this%cg_thickini(node)
2156 if (pctcomp >= 5.0_dp)
then
2158 else if (pctcomp >=
done)
then
2163 call this%dis%noder_to_string(node, cellid)
2166 call this%outputtab%add_term(cellid)
2167 call this%outputtab%add_term(this%cg_thickini(node))
2168 call this%outputtab%add_term(this%cg_thick(node))
2169 call this%outputtab%add_term(this%cg_tcomp(node))
2170 call this%outputtab%add_term(strain)
2171 call this%outputtab%add_term(pctcomp)
2172 call this%outputtab%add_term(cflag)
2174 write (this%iout,
'(/1X,A,1X,I0,1X,A,1X,I0,1X,A,/1X,A,/1X,A)') &
2175 'COARSE-GRAINED STORAGE PERCENT COMPACTION IS GREATER THAN OR '// &
2176 'EQUAL TO 1 PERCENT IN', iexceed,
'OF', this%dis%nodes,
'CELL(S).', &
2177 'USE THE STRAIN_CSV_COARSE OPTION TO OUTPUT A CSV '// &
2178 'FILE WITH PERCENT COMPACTION ',
'VALUES FOR ALL CELLS.'
2180 msg =
'COARSE-GRAINED STORAGE PERCENT COMPACTION WAS LESS THAN '// &
2181 '1 PERCENT IN ALL CELLS '
2182 write (this%iout,
'(/1X,A)') trim(adjustl(msg))
2186 if (this%istrainsk /= 0)
then
2189 ntabrows = this%dis%nodes
2191 if (this%dis%ndim > 1)
then
2192 ntabcols = ntabcols + 1
2194 ntabcols = ntabcols + this%dis%ndim
2197 call table_cr(this%outputtab, this%packName,
'')
2198 call this%outputtab%table_df(ntabrows, ntabcols, this%istrainsk, &
2199 lineseparator=.false., separator=
',')
2203 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2204 if (this%dis%ndim == 2)
then
2206 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2208 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2211 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2213 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2215 call this%outputtab%initialize_column(tag, 10, alignment=
tabright)
2217 tag =
'INITIAL_THICKNESS'
2218 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2219 tag =
'FINAL_THICKNESS'
2220 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2221 tag =
'TOTAL_COMPACTION'
2222 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2223 tag =
'TOTAL_STRAIN'
2224 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2225 tag =
'PERCENT_COMPACTION'
2226 call this%outputtab%initialize_column(tag, 20, alignment=
tabright)
2229 do node = 1, this%dis%nodes
2230 if (this%cg_thickini(node) >
dzero)
then
2231 strain = this%cg_tcomp(node) / this%cg_thickini(node)
2236 call this%dis%noder_to_array(node, locs)
2239 if (this%dis%ndim > 1)
then
2240 call this%outputtab%add_term(this%dis%get_nodeuser(node))
2242 do ipos = 1, this%dis%ndim
2243 call this%outputtab%add_term(locs(ipos))
2245 call this%outputtab%add_term(this%cg_thickini(node))
2246 call this%outputtab%add_term(this%cg_thick(node))
2247 call this%outputtab%add_term(this%cg_tcomp(node))
2248 call this%outputtab%add_term(strain)
2249 call this%outputtab%add_term(pctcomp)
2255 if (this%ndelaybeds > 0)
then
2256 if (this%idb_nconv_count(2) > 0)
then
2257 write (
warnmsg,
'(a,1x,a,1x,i0,1x,a,1x,a)') &
2258 'Delay interbed cell heads were less than the top of the interbed', &
2259 'cell in', this%idb_nconv_count(2),
'interbed cells in ', &
2260 'non-convertible GWF cells for at least one time step during '// &
2267 deallocate (imap_sel)
2269 deallocate (pctcomp_arr)
2287 if (this%inunit > 0)
then
2309 if (this%iupdatematprop == 0)
then
2310 nullify (this%cg_thick)
2311 nullify (this%cg_thick0)
2312 nullify (this%cg_theta)
2313 nullify (this%cg_theta0)
2328 call mem_deallocate(this%boundname,
'BOUNDNAME', this%memoryPath)
2345 if (this%iupdatematprop == 0)
then
2346 nullify (this%thick)
2347 nullify (this%thick0)
2348 nullify (this%theta)
2349 nullify (this%theta0)
2360 if (this%ndelaybeds > 0)
then
2361 if (this%iupdatematprop == 0)
then
2363 nullify (this%dbdz0)
2364 nullify (this%dbtheta)
2365 nullify (this%dbtheta0)
2404 nullify (this%gwfiss)
2407 nullify (this%stoiconv)
2408 nullify (this%stoss)
2411 if (this%iprpak > 0)
then
2412 call this%inputtab%table_da()
2413 deallocate (this%inputtab)
2414 nullify (this%inputtab)
2418 if (this%istrainib > 0 .or. this%istrainsk > 0)
then
2419 call this%outputtab%table_da()
2420 deallocate (this%outputtab)
2421 nullify (this%outputtab)
2426 if (this%ipakcsv > 0)
then
2427 call this%pakcsvtab%table_da()
2428 deallocate (this%pakcsvtab)
2429 nullify (this%pakcsvtab)
2433 call mem_deallocate(this%listlabel,
'LISTLABEL', this%memoryPath)
2477 if (this%inunit > 0)
then
2478 call this%obs%obs_da()
2479 call this%TsManager%da()
2482 deallocate (this%obs)
2487 deallocate (this%TsManager)
2488 nullify (this%TsManager)
2492 call this%NumericalPackageType%da()
2513 character(len=LINELENGTH) :: line
2514 character(len=LINELENGTH) :: title
2515 character(len=LINELENGTH) :: text
2516 character(len=20) :: cellid
2518 logical :: endOfBlock
2520 integer(I4B) :: ierr
2521 integer(I4B) :: node
2522 integer(I4B) :: nlist
2523 real(DP),
pointer :: bndElem => null()
2525 character(len=*),
parameter :: fmtblkerr = &
2526 &
"('Looking for BEGIN PERIOD iper. Found ',a,' instead.')"
2527 character(len=*),
parameter :: fmtlsp = &
2528 &
"(1X,/1X,'REUSING ',a,'S FROM LAST STRESS PERIOD')"
2531 if (this%inunit == 0)
return
2534 if (this%ionper <
kper)
then
2537 call this%parser%GetBlock(
'PERIOD', isfound, ierr, &
2538 supportopenclose=.true., &
2539 blockrequired=.false.)
2543 call this%read_check_ionper()
2549 this%ionper =
nper + 1
2552 call this%parser%GetCurrentLine(line)
2553 write (
errmsg, fmtblkerr) adjustl(trim(line))
2560 if (this%ionper ==
kper)
then
2563 if (this%iprpak /= 0)
then
2566 title =
'CSUB'//
' PACKAGE ('// &
2567 trim(adjustl(this%packName))//
') DATA FOR PERIOD'
2568 write (title,
'(a,1x,i6)') trim(adjustl(title)),
kper
2569 call table_cr(this%inputtab, this%packName, title)
2570 call this%inputtab%table_df(1, 2, this%iout, finalize=.false.)
2572 call this%inputtab%initialize_column(text, 20)
2574 call this%inputtab%initialize_column(text, 15, alignment=tableft)
2581 call this%TsManager%Reset(this%packName)
2585 call this%parser%GetNextLine(endofblock)
2588 if (endofblock)
then
2596 if (nlist > this%maxsig0)
then
2597 write (
errmsg,
'(a,i0,a,i0,a)') &
2598 'The number of stress period entries (', nlist, &
2599 ') exceeds the maximum number of stress period entries (', &
2606 call this%parser%GetCellid(this%dis%ndim, cellid)
2607 node = this%dis%noder_from_cellid(cellid, &
2608 this%parser%iuactive, this%iout)
2612 write (
errmsg,
'(a,2(1x,a))') &
2613 'CELLID', cellid,
'is not in the active model domain.'
2617 this%nodelistsig0(nlist) = node
2620 call this%parser%GetString(text)
2622 bndelem => this%sig0(nlist)
2624 this%packName,
'BND', &
2625 this%tsManager, this%iprpak, &
2629 if (this%iprpak /= 0)
then
2630 call this%dis%noder_to_string(node, cellid)
2631 call this%inputtab%add_term(cellid)
2632 call this%inputtab%add_term(bndelem)
2640 if (this%iprpak /= 0)
then
2641 call this%inputtab%finalize_table()
2646 write (this%iout, fmtlsp) trim(this%filtyp)
2651 call this%parser%StoreErrorUnit()
2655 call this%csub_rp_obs()
2674 integer(I4B),
intent(in) :: nodes
2675 real(DP),
dimension(nodes),
intent(in) :: hnew
2679 integer(I4B) :: idelay
2680 integer(I4B) :: node
2687 if (this%ninterbeds > 0)
then
2689 if (this%gwfiss /= 0)
then
2690 write (
errmsg,
'(a,i0,a,1x,a,1x,a,1x,i0,1x,a)') &
2691 'Only the first and last (',
nper,
')', &
2692 'stress period can be steady if interbeds are simulated.', &
2693 'Stress period',
kper,
'has been defined to be steady state.'
2700 if (this%initialized == 0)
then
2701 if (this%gwfiss == 0)
then
2702 call this%csub_set_initial_state(nodes, hnew)
2710 this%cg_comp(node) =
dzero
2711 this%cg_es0(node) = this%cg_es(node)
2712 if (this%iupdatematprop /= 0)
then
2713 this%cg_thick0(node) = this%cg_thick(node)
2714 this%cg_theta0(node) = this%cg_theta(node)
2719 do ib = 1, this%ninterbeds
2720 idelay = this%idelay(ib)
2723 this%comp(ib) =
dzero
2724 node = this%nodelist(ib)
2725 if (this%initialized /= 0)
then
2726 es = this%cg_es(node)
2732 if (this%iupdatematprop /= 0)
then
2733 this%thick0(ib) = this%thick(ib)
2734 this%theta0(ib) = this%theta(ib)
2738 if (idelay /= 0)
then
2742 if (this%gwfiss0 /= 0)
then
2743 node = this%nodelist(ib)
2745 do n = 1, this%ndelaycells
2746 this%dbh(n, idelay) = h
2752 do n = 1, this%ndelaycells
2754 if (this%initialized /= 0)
then
2755 if (this%dbes(n, idelay) > this%dbpcs(n, idelay))
then
2756 this%dbpcs(n, idelay) = this%dbes(n, idelay)
2759 this%dbh0(n, idelay) = this%dbh(n, idelay)
2760 this%dbes0(n, idelay) = this%dbes(n, idelay)
2761 if (this%iupdatematprop /= 0)
then
2762 this%dbdz0(n, idelay) = this%dbdz(n, idelay)
2763 this%dbtheta0(n, idelay) = this%dbtheta(n, idelay)
2770 this%gwfiss0 = this%gwfiss
2773 call this%TsManager%ad()
2778 call this%obs%obs_ad()
2789 subroutine csub_fc(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
2794 integer(I4B),
intent(in) :: kiter
2795 real(DP),
intent(in),
dimension(:) :: hold
2796 real(DP),
intent(in),
dimension(:) :: hnew
2798 integer(I4B),
intent(in),
dimension(:) :: idxglo
2799 real(DP),
intent(inout),
dimension(:) :: rhs
2802 integer(I4B) :: node
2803 integer(I4B) :: idiag
2804 integer(I4B) :: idelay
2812 call this%csub_cg_calc_stress(this%dis%nodes, hnew)
2815 if (this%gwfiss == 0)
then
2821 do node = 1, this%dis%nodes
2822 idiag = this%dis%con%ia(node)
2823 area = this%dis%get_area(node)
2826 if (this%ibound(node) < 1) cycle
2829 if (this%iupdatematprop /= 0)
then
2830 if (this%ieslag == 0)
then
2833 call this%csub_cg_calc_comp(node, hnew(node), hold(node), comp)
2834 this%cg_comp(node) = comp
2837 call this%csub_cg_update(node)
2842 call this%csub_cg_fc(node, tled, area, hnew(node), hold(node), &
2846 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2847 rhs(node) = rhs(node) + rhsterm
2851 if (this%brg /=
dzero)
then
2852 call this%csub_cg_wcomp_fc(node, tled, area, hnew(node), hold(node), &
2857 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2858 rhs(node) = rhs(node) + rhsterm
2863 if (this%ninterbeds /= 0)
then
2867 do ib = 1, this%ninterbeds
2868 node = this%nodelist(ib)
2869 idelay = this%idelay(ib)
2870 idiag = this%dis%con%ia(node)
2871 area = this%dis%get_area(node)
2872 call this%csub_interbed_fc(ib, node, area, hnew(node), hold(node), &
2874 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2875 rhs(node) = rhs(node) + rhsterm
2879 call this%csub_nodelay_wcomp_fc(ib, node, tled, area, &
2880 hnew(node), hold(node), &
2884 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2885 rhs(node) = rhs(node) + rhsterm
2893 call this%parser%StoreErrorUnit()
2909 subroutine csub_fn(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
2914 integer(I4B),
intent(in) :: kiter
2915 real(DP),
intent(in),
dimension(:) :: hold
2916 real(DP),
intent(in),
dimension(:) :: hnew
2918 integer(I4B),
intent(in),
dimension(:) :: idxglo
2919 real(DP),
intent(inout),
dimension(:) :: rhs
2921 integer(I4B) :: idelay
2922 integer(I4B) :: node
2923 integer(I4B) :: idiag
2931 if (this%gwfiss == 0)
then
2935 do node = 1, this%dis%nodes
2936 idiag = this%dis%con%ia(node)
2937 area = this%dis%get_area(node)
2940 if (this%ibound(node) < 1) cycle
2943 call this%csub_cg_fn(node, tled, area, &
2944 hnew(node), hcof, rhsterm)
2948 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2949 rhs(node) = rhs(node) + rhsterm
2953 if (this%brg /=
dzero)
then
2954 call this%csub_cg_wcomp_fn(node, tled, area, hnew(node), hold(node), &
2959 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2960 rhs(node) = rhs(node) + rhsterm
2965 if (this%ninterbeds /= 0)
then
2969 do ib = 1, this%ninterbeds
2970 idelay = this%idelay(ib)
2971 node = this%nodelist(ib)
2974 if (this%ibound(node) < 1) cycle
2977 idiag = this%dis%con%ia(node)
2978 area = this%dis%get_area(node)
2979 call this%csub_interbed_fn(ib, node, hnew(node), hold(node), &
2983 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2984 rhs(node) = rhs(node) + rhsterm
2987 if (this%brg /=
dzero .and. idelay == 0)
then
2988 call this%csub_nodelay_wcomp_fn(ib, node, tled, area, &
2989 hnew(node), hold(node), &
2993 call matrix_sln%add_value_pos(idxglo(idiag), hcof)
2994 rhs(node) = rhs(node) + rhsterm
3017 subroutine csub_cc(this, innertot, kiter, iend, icnvgmod, nodes, &
3018 hnew, hold, cpak, ipak, dpak)
3023 integer(I4B),
intent(in) :: innertot
3024 integer(I4B),
intent(in) :: kiter
3025 integer(I4B),
intent(in) :: iend
3026 integer(I4B),
intent(in) :: icnvgmod
3027 integer(I4B),
intent(in) :: nodes
3028 real(DP),
dimension(nodes),
intent(in) :: hnew
3029 real(DP),
dimension(nodes),
intent(in) :: hold
3030 character(len=LENPAKLOC),
intent(inout) :: cpak
3031 integer(I4B),
intent(inout) :: ipak
3032 real(DP),
intent(inout) :: dpak
3034 character(len=LINELENGTH) :: tag
3035 character(len=LENPAKLOC) :: cloc
3036 integer(I4B) :: icheck
3037 integer(I4B) :: ipakfail
3038 integer(I4B) :: ntabrows
3039 integer(I4B) :: ntabcols
3041 integer(I4B) :: node
3042 integer(I4B) :: idelay
3043 integer(I4B) :: locdhmax
3044 integer(I4B) :: locrmax
3045 integer(I4B) :: ifirst
3051 real(DP) :: hcellold
3065 icheck = this%iconvchk
3075 if (this%gwfiss /= 0)
then
3081 if (this%ipakcsv == 0)
then
3082 if (icnvgmod == 0)
then
3088 if (.not.
associated(this%pakcsvtab))
then
3095 call table_cr(this%pakcsvtab, this%packName,
'')
3096 call this%pakcsvtab%table_df(ntabrows, ntabcols, this%ipakcsv, &
3097 lineseparator=.false., separator=
',', &
3101 tag =
'total_inner_iterations'
3102 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
3104 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
3106 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
3108 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
3110 call this%pakcsvtab%initialize_column(tag, 10, alignment=
tableft)
3112 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
3114 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
3116 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
3117 tag =
'dstoragemax_loc'
3118 call this%pakcsvtab%initialize_column(tag, 15, alignment=
tableft)
3124 if (icheck /= 0)
then
3130 final_check:
do ib = 1, this%ninterbeds
3131 idelay = this%idelay(ib)
3132 node = this%nodelist(ib)
3135 if (idelay == 0) cycle
3138 if (this%ibound(node) < 1) cycle
3141 dh = this%dbdhmax(idelay)
3146 area = this%dis%get_area(node)
3148 hcellold = hold(node)
3151 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
3154 call this%csub_delay_calc_dstor(ib, hcell, stoe, stoi)
3155 v1 = (stoe + stoi) * area * this%rnb(ib) * tled
3158 call this%csub_delay_calc_wcomp(ib, dwc)
3159 v1 = v1 + dwc * area * this%rnb(ib)
3162 call this%csub_delay_fc(ib, hcof, rhs)
3163 v2 = (-hcof * hcell - rhs) * area * this%rnb(ib)
3170 df = df *
delt / area
3173 if (ifirst == 1)
then
3180 if (abs(dh) > abs(dhmax))
then
3184 if (abs(df) > abs(rmax))
then
3193 if (abs(dhmax) > abs(dpak))
then
3196 write (cloc,
"(a,'-',a)") trim(this%packName),
'head'
3201 if (abs(rmax) > abs(dpak))
then
3204 write (cloc,
"(a,'-',a)") trim(this%packName),
'storage'
3209 if (this%ipakcsv /= 0)
then
3212 call this%pakcsvtab%add_term(innertot)
3213 call this%pakcsvtab%add_term(
totim)
3214 call this%pakcsvtab%add_term(
kper)
3215 call this%pakcsvtab%add_term(
kstp)
3216 call this%pakcsvtab%add_term(kiter)
3217 call this%pakcsvtab%add_term(dhmax)
3218 call this%pakcsvtab%add_term(locdhmax)
3219 call this%pakcsvtab%add_term(rmax)
3220 call this%pakcsvtab%add_term(locrmax)
3224 call this%pakcsvtab%finalize_table()
3242 subroutine csub_cq(this, nodes, hnew, hold, isuppress_output, flowja)
3248 integer(I4B),
intent(in) :: nodes
3249 real(DP),
intent(in),
dimension(nodes) :: hnew
3250 real(DP),
intent(in),
dimension(nodes) :: hold
3251 integer(I4B),
intent(in) :: isuppress_output
3252 real(DP),
dimension(:),
contiguous,
intent(inout) :: flowja
3255 integer(I4B) :: idelay
3256 integer(I4B) :: ielastic
3257 integer(I4B) :: iconvert
3258 integer(I4B) :: node
3261 integer(I4B) :: idiag
3288 integer(I4B) :: iprobslocal
3298 do node = 1, this%dis%nodes
3299 idiag = this%dis%con%ia(node)
3300 area = this%dis%get_area(node)
3304 if (this%gwfiss == 0)
then
3310 if (this%ibound(node) > 0 .and. this%cg_thickini(node) >
dzero)
then
3313 call this%csub_cg_fc(node, tled, area, hnew(node), hold(node), &
3315 rrate = hcof * hnew(node) - rhs
3318 call this%csub_cg_calc_comp(node, hnew(node), hold(node), comp)
3321 call this%csub_cg_wcomp_fc(node, tled, area, hnew(node), hold(node), &
3323 rratewc = hcof * hnew(node) - rhs
3329 this%cg_stor(node) = rrate
3330 this%cell_wcstor(node) = rratewc
3331 this%cell_thick(node) = this%cg_thick(node)
3334 this%cg_comp(node) = comp
3338 if (isuppress_output == 0)
then
3342 if (this%iupdatematprop /= 0)
then
3343 call this%csub_cg_update(node)
3347 this%cg_tcomp(node) = this%cg_tcomp(node) + comp
3351 flowja(idiag) = flowja(idiag) + rrate
3352 flowja(idiag) = flowja(idiag) + rratewc
3358 if (this%ndelaybeds > 0)
then
3359 this%idb_nconv_count(1) = 0
3366 do ib = 1, this%ninterbeds
3368 idelay = this%idelay(ib)
3369 ielastic = this%ielastic(ib)
3373 if (idelay == 0)
then
3377 b = this%thick(ib) * this%rnb(ib)
3381 node = this%nodelist(ib)
3382 idiag = this%dis%con%ia(node)
3383 area = this%dis%get_area(node)
3386 this%cell_thick(node) = this%cell_thick(node) + b
3389 if (this%gwfiss == 0)
then
3397 if (this%ibound(node) < 1) cycle
3400 if (idelay == 0)
then
3401 iconvert = this%iconvert(ib)
3405 call this%csub_nodelay_calc_comp(ib, hnew(node), hold(node), comp, &
3409 es = this%cg_es(node)
3411 es0 = this%cg_es0(node)
3414 if (ielastic > 0 .or. iconvert == 0)
then
3417 stoi = -pcs * rho2 + (rho2 * es)
3418 stoe = pcs * rho1 - (rho1 * es0)
3424 this%storagee(ib) = stoe * tledm
3425 this%storagei(ib) = stoi * tledm
3428 this%comp(ib) = comp
3431 if (isuppress_output == 0)
then
3434 if (this%iupdatematprop /= 0)
then
3435 call this%csub_nodelay_update(ib)
3439 this%tcomp(ib) = this%tcomp(ib) + comp
3440 this%tcompe(ib) = this%tcompe(ib) + compe
3441 this%tcompi(ib) = this%tcompi(ib) + compi
3450 call this%csub_calc_sat(node, h, h0, snnew, snold)
3453 call this%csub_delay_calc_dstor(ib, h, stoe, stoi)
3454 this%storagee(ib) = stoe * area * this%rnb(ib) * tledm
3455 this%storagei(ib) = stoi * area * this%rnb(ib) * tledm
3458 q = this%csub_calc_delay_flow(ib, 1, h) * area * this%rnb(ib)
3459 this%dbflowtop(idelay) = q
3460 nn = this%ndelaycells
3461 q = this%csub_calc_delay_flow(ib, nn, h) * area * this%rnb(ib)
3462 this%dbflowbot(idelay) = q
3465 if (isuppress_output == 0)
then
3468 call this%csub_delay_calc_comp(ib, h, h0, comp, compi, compe)
3472 if (this%iupdatematprop /= 0)
then
3473 call this%csub_delay_update(ib)
3477 this%tcomp(ib) = this%tcomp(ib) + comp
3478 this%tcompi(ib) = this%tcompi(ib) + compi
3479 this%tcompe(ib) = this%tcompe(ib) + compe
3482 do n = 1, this%ndelaycells
3483 this%dbtcomp(n, idelay) = this%dbtcomp(n, idelay) + &
3484 this%dbcomp(n, idelay)
3489 call this%csub_delay_head_check(ib)
3496 if (idelay == 0)
then
3497 call this%csub_nodelay_wcomp_fc(ib, node, tledm, area, &
3498 hnew(node), hold(node), hcof, rhs)
3499 rratewc = hcof * hnew(node) - rhs
3503 call this%csub_delay_calc_wcomp(ib, q)
3504 rratewc = q * area * this%rnb(ib)
3506 this%cell_wcstor(node) = this%cell_wcstor(node) + rratewc
3509 flowja(idiag) = flowja(idiag) + rratewc
3511 this%storagee(ib) =
dzero
3512 this%storagei(ib) =
dzero
3513 if (idelay /= 0)
then
3514 this%dbflowtop(idelay) =
dzero
3515 this%dbflowbot(idelay) =
dzero
3520 flowja(idiag) = flowja(idiag) + this%storagee(ib)
3521 flowja(idiag) = flowja(idiag) + this%storagei(ib)
3525 if (this%iupdatematprop /= 0)
then
3527 call this%parser%StoreErrorUnit()
3545 subroutine csub_bd(this, isuppress_output, model_budget)
3552 integer(I4B),
intent(in) :: isuppress_output
3553 type(
budgettype),
intent(inout) :: model_budget
3560 call model_budget%addentry(rin, rout,
delt,
budtxt(1), &
3561 isuppress_output,
' CSUB')
3562 if (this%ninterbeds > 0)
then
3566 call model_budget%addentry(rin, rout,
delt,
budtxt(2), &
3567 isuppress_output,
' CSUB')
3571 call model_budget%addentry(rin, rout,
delt,
budtxt(3), &
3572 isuppress_output,
' CSUB')
3575 call model_budget%addentry(rin, rout,
delt,
budtxt(4), &
3576 isuppress_output,
' CSUB')
3588 integer(I4B),
intent(in) :: icbcfl
3589 integer(I4B),
intent(in) :: icbcun
3591 character(len=1) :: cdatafmp =
' '
3592 character(len=1) :: editdesc =
' '
3593 integer(I4B) :: ibinun
3594 integer(I4B) :: iprint
3595 integer(I4B) :: nvaluesp
3596 integer(I4B) :: nwidthp
3598 integer(I4B) :: node
3599 integer(I4B) :: naux
3605 if (this%ipakcb < 0)
then
3607 elseif (this%ipakcb == 0)
then
3610 ibinun = this%ipakcb
3612 if (icbcfl == 0) ibinun = 0
3615 if (ibinun /= 0)
then
3620 call this%dis%record_array(this%cg_stor, this%iout, iprint, -ibinun, &
3621 budtxt(1), cdatafmp, nvaluesp, &
3622 nwidthp, editdesc, dinact)
3623 if (this%ninterbeds > 0)
then
3627 call this%dis%record_srcdst_list_header(
budtxt(2), &
3637 do ib = 1, this%ninterbeds
3638 q = this%storagee(ib)
3639 node = this%nodelist(ib)
3640 call this%dis%record_mf6_list_entry(ibinun, node, node, q, naux, &
3645 call this%dis%record_srcdst_list_header(
budtxt(3), &
3655 do ib = 1, this%ninterbeds
3656 q = this%storagei(ib)
3657 node = this%nodelist(ib)
3658 call this%dis%record_mf6_list_entry(ibinun, node, node, q, naux, &
3664 call this%dis%record_array(this%cell_wcstor, this%iout, iprint, -ibinun, &
3665 budtxt(4), cdatafmp, nvaluesp, &
3666 nwidthp, editdesc, dinact)
3682 integer(I4B),
intent(in) :: idvfl
3683 integer(I4B),
intent(in) :: idvprint
3685 character(len=1) :: cdatafmp =
' '
3686 character(len=1) :: editdesc =
' '
3687 integer(I4B) :: ibinun
3688 integer(I4B) :: iprint
3689 integer(I4B) :: nvaluesp
3690 integer(I4B) :: nwidthp
3692 integer(I4B) :: node
3693 integer(I4B) :: nodem
3694 integer(I4B) :: nodeu
3697 integer(I4B) :: idx_conn
3699 integer(I4B) :: ncpl
3700 integer(I4B) :: nlay
3703 real(DP) :: va_scale
3705 character(len=*),
parameter :: fmtnconv = &
3706 "(/4x, 'DELAY INTERBED CELL HEADS IN ', i0, ' INTERBEDS IN', &
3707 &' NON-CONVERTIBLE GWF CELLS WERE LESS THAN THE TOP OF THE INTERBED CELL')"
3712 if (this%ioutcomp /= 0 .or. this%ioutzdisp /= 0)
then
3717 if (idvfl == 0) ibinun = 0
3720 if (ibinun /= 0)
then
3725 do node = 1, this%dis%nodes
3726 this%buff(node) = this%cg_tcomp(node)
3728 do ib = 1, this%ninterbeds
3729 node = this%nodelist(ib)
3730 this%buff(node) = this%buff(node) + this%tcomp(ib)
3734 if (this%ioutcomp /= 0)
then
3735 ibinun = this%ioutcomp
3736 call this%dis%record_array(this%buff, this%iout, iprint, ibinun, &
3737 comptxt(1), cdatafmp, nvaluesp, &
3738 nwidthp, editdesc, dinact)
3742 if (this%ioutzdisp /= 0)
then
3743 ibinun = this%ioutzdisp
3746 do nodeu = 1, this%dis%nodesuser
3747 this%buffusr(nodeu) =
dzero
3751 do node = 1, this%dis%nodes
3752 nodeu = this%dis%get_nodeuser(node)
3753 this%buffusr(nodeu) = this%buff(node)
3757 ncpl = this%dis%get_ncpl()
3760 if (this%dis%ndim == 1)
then
3761 do node = this%dis%nodes, 1, -1
3762 do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
3765 nodem = this%dis%con%ja(ii)
3766 idx_conn = this%dis%con%jas(ii)
3769 ihc = this%dis%con%ihc(idx_conn)
3773 if (node < nodem)
then
3774 va_scale = this%dis%get_area_factor(node, idx_conn)
3775 this%buffusr(node) = this%buffusr(node) + &
3776 va_scale * this%buffusr(nodem)
3783 nlay = this%dis%nodesuser / ncpl
3784 do k = nlay - 1, 1, -1
3786 node = (k - 1) * ncpl + i
3787 nodem = k * ncpl + i
3788 this%buffusr(node) = this%buffusr(node) + this%buffusr(nodem)
3794 do nodeu = 1, this%dis%nodesuser
3795 node = this%dis%get_nodenumber_idx1(nodeu, 1)
3797 this%buff(node) = this%buffusr(nodeu)
3802 call this%dis%record_array(this%buff, this%iout, iprint, ibinun, &
3803 comptxt(6), cdatafmp, nvaluesp, &
3804 nwidthp, editdesc, dinact)
3810 if (this%ioutcompi /= 0)
then
3811 ibinun = this%ioutcompi
3815 if (idvfl == 0) ibinun = 0
3818 if (ibinun /= 0)
then
3823 do node = 1, this%dis%nodes
3824 this%buff(node) =
dzero
3826 do ib = 1, this%ninterbeds
3827 node = this%nodelist(ib)
3828 this%buff(node) = this%buff(node) + this%tcompi(ib)
3832 call this%dis%record_array(this%buff, this%iout, iprint, ibinun, &
3833 comptxt(2), cdatafmp, nvaluesp, &
3834 nwidthp, editdesc, dinact)
3838 if (this%ioutcompe /= 0)
then
3839 ibinun = this%ioutcompe
3843 if (idvfl == 0) ibinun = 0
3846 if (ibinun /= 0)
then
3851 do node = 1, this%dis%nodes
3852 this%buff(node) =
dzero
3854 do ib = 1, this%ninterbeds
3855 node = this%nodelist(ib)
3856 this%buff(node) = this%buff(node) + this%tcompe(ib)
3860 call this%dis%record_array(this%buff, this%iout, iprint, ibinun, &
3861 comptxt(3), cdatafmp, nvaluesp, &
3862 nwidthp, editdesc, dinact)
3866 if (this%ioutcompib /= 0)
then
3867 ibinun = this%ioutcompib
3871 if (idvfl == 0) ibinun = 0
3874 if (ibinun /= 0)
then
3879 do node = 1, this%dis%nodes
3880 this%buff(node) =
dzero
3882 do ib = 1, this%ninterbeds
3883 node = this%nodelist(ib)
3884 this%buff(node) = this%buff(node) + this%tcompe(ib) + this%tcompi(ib)
3888 call this%dis%record_array(this%buff, this%iout, iprint, ibinun, &
3889 comptxt(4), cdatafmp, nvaluesp, &
3890 nwidthp, editdesc, dinact)
3894 if (this%ioutcomps /= 0)
then
3895 ibinun = this%ioutcomps
3899 if (idvfl == 0) ibinun = 0
3902 if (ibinun /= 0)
then
3907 do node = 1, this%dis%nodes
3908 this%buff(node) = this%cg_tcomp(node)
3912 call this%dis%record_array(this%buff, this%iout, iprint, ibinun, &
3913 comptxt(5), cdatafmp, nvaluesp, &
3914 nwidthp, editdesc, dinact)
3919 if (this%gwfiss == 0)
then
3920 call this%csub_cg_chk_stress()
3927 if (this%ndelaybeds > 0)
then
3928 if (this%idb_nconv_count(1) > this%idb_nconv_count(2))
then
3929 this%idb_nconv_count(2) = this%idb_nconv_count(1)
3931 if (this%idb_nconv_count(1) > 0)
then
3932 write (this%iout, fmtnconv) this%idb_nconv_count(1)
3950 integer(I4B),
intent(in) :: nodes
3951 real(DP),
dimension(nodes),
intent(in) :: hnew
3953 integer(I4B) :: node
3957 integer(I4B) :: idx_conn
3962 real(DP) :: va_scale
3971 if (this%iupdatestress /= 0)
then
3972 do node = 1, this%dis%nodes
3977 top = this%dis%top(node)
3978 bot = this%dis%bot(node)
3982 if (this%ibound(node) /= 0)
then
3992 if (hcell < top)
then
3993 gs = (top - hbar) * this%sgm(node) + (hbar - bot) * this%sgs(node)
3995 gs = thick * this%sgs(node)
3999 this%cg_gs(node) = gs
4003 do nn = 1, this%nbound
4004 node = this%nodelistsig0(nn)
4005 sadd = this%sig0(nn)
4006 this%cg_gs(node) = this%cg_gs(node) + sadd
4010 do node = 1, this%dis%nodes
4013 gs = this%cg_gs(node)
4017 do ii = this%dis%con%ia(node) + 1, this%dis%con%ia(node + 1) - 1
4020 m = this%dis%con%ja(ii)
4021 idx_conn = this%dis%con%jas(ii)
4024 if (this%dis%con%ihc(idx_conn) == 0)
then
4030 if (this%dis%ndim /= 1)
then
4031 gs = gs + this%cg_gs(m)
4035 va_scale = this%dis%get_area_factor(node, idx_conn)
4036 gs_conn = this%cg_gs(m)
4037 gs = gs + (gs_conn * va_scale)
4045 this%cg_gs(node) = gs
4051 do node = 1, this%dis%nodes
4052 top = this%dis%top(node)
4053 bot = this%dis%bot(node)
4054 if (this%ibound(node) /= 0)
then
4067 es = this%cg_gs(node) - phead
4068 this%cg_es(node) = es
4088 character(len=20) :: cellid
4089 integer(I4B) :: ierr
4090 integer(I4B) :: node
4104 do node = 1, this%dis%nodes
4105 if (this%ibound(node) < 1) cycle
4106 bot = this%dis%bot(node)
4107 gs = this%cg_gs(node)
4108 es = this%cg_es(node)
4110 if (this%ibound(node) /= 0)
then
4114 if (this%lhead_based .EQV. .false.)
then
4117 call this%dis%noder_to_string(node, cellid)
4118 write (
errmsg,
'(a,g0,a,1x,a,1x,a,4(g0,a))') &
4119 'Small to negative effective stress (', es,
') in cell', &
4120 trim(adjustl(cellid)),
'. (', es,
' = ', this%cg_gs(node), &
4121 ' - (', hcell,
' - ', bot,
').'
4129 write (
errmsg,
'(a,1x,i0,3(1x,a))') &
4130 'Solution: small to negative effective stress values in', ierr, &
4131 'cells can be eliminated by increasing storage values and/or ', &
4132 'adding/modifying stress boundaries to prevent water-levels from', &
4133 'exceeding the top of the model.'
4135 call this%parser%StoreErrorUnit()
4152 integer(I4B),
intent(in) :: i
4159 comp = this%tcomp(i) + this%comp(i)
4160 if (abs(comp) >
dzero)
then
4161 thick = this%thickini(i)
4162 theta = this%thetaini(i)
4163 call this%csub_adj_matprop(comp, thick, theta)
4164 if (thick <=
dzero)
then
4165 write (
errmsg,
'(a,1x,i0,1x,a,g0,a)') &
4166 'Adjusted thickness for no-delay interbed', i, &
4167 'is less than or equal to 0 (', thick,
').'
4170 if (theta <=
dzero)
then
4171 write (
errmsg,
'(a,1x,i0,1x,a,g0,a)') &
4172 'Adjusted theta for no-delay interbed', i, &
4173 'is less than or equal to 0 (', theta,
').'
4176 this%thick(i) = thick
4177 this%theta(i) = theta
4202 integer(I4B),
intent(in) :: ib
4203 real(DP),
intent(in) :: hcell
4204 real(DP),
intent(in) :: hcellold
4205 real(DP),
intent(inout) :: rho1
4206 real(DP),
intent(inout) :: rho2
4207 real(DP),
intent(inout) :: rhs
4208 real(DP),
intent(in),
optional :: argtled
4210 integer(I4B) :: node
4220 real(DP) :: sto_fac0
4230 if (
present(argtled))
then
4235 node = this%nodelist(ib)
4236 area = this%dis%get_area(node)
4237 bot = this%dis%bot(node)
4238 top = this%dis%top(node)
4239 thick = this%thickini(ib)
4245 this%iconvert(ib) = 0
4248 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
4249 if (this%lhead_based .EQV. .true.)
then
4253 znode = this%csub_calc_znode(top, bot, hbar)
4254 es = this%cg_es(node)
4255 es0 = this%cg_es0(node)
4256 theta = this%thetaini(ib)
4261 call this%csub_calc_sfacts(node, bot, znode, theta, es, es0, f)
4263 sto_fac = tled * snnew * thick * f
4264 sto_fac0 = tled * snold * thick * f
4267 rho1 = this%rci(ib) * sto_fac0
4268 rho2 = this%rci(ib) * sto_fac
4269 if (this%cg_es(node) > this%pcs(ib))
then
4270 this%iconvert(ib) = 1
4271 rho2 = this%ci(ib) * sto_fac
4275 rcorr = rho2 * (hcell - hbar)
4278 if (this%ielastic(ib) /= 0)
then
4279 rhs = rho1 * this%cg_es0(node) - &
4280 rho2 * (this%cg_gs(node) + bot) - &
4283 rhs = -rho2 * (this%cg_gs(node) + bot) + &
4284 (this%pcs(ib) * (rho2 - rho1)) + &
4285 (rho1 * this%cg_es0(node)) - &
4311 integer(I4B),
intent(in) :: ib
4312 real(DP),
intent(in) :: hcell
4313 real(DP),
intent(in) :: hcellold
4314 real(DP),
intent(inout) :: comp
4315 real(DP),
intent(inout) :: rho1
4316 real(DP),
intent(inout) :: rho2
4318 integer(I4B) :: node
4326 node = this%nodelist(ib)
4328 es = this%cg_es(node)
4329 es0 = this%cg_es0(node)
4333 call this%csub_nodelay_fc(ib, hcell, hcellold, rho1, rho2, rhs, argtled=tled)
4336 if (this%ielastic(ib) /= 0)
then
4337 comp = rho2 * es - rho1 * es0
4339 comp = -pcs * (rho2 - rho1) - (rho1 * es0) + (rho2 * es)
4357 integer(I4B),
intent(in) :: nodes
4358 real(DP),
dimension(nodes),
intent(in) :: hnew
4360 character(len=LINELENGTH) :: title
4361 character(len=LINELENGTH) :: tag
4362 character(len=20) :: cellid
4364 integer(I4B) :: node
4366 integer(I4B) :: idelay
4367 integer(I4B) :: ntabrows
4368 integer(I4B) :: ntabcols
4374 real(DP) :: void_ratio
4384 call this%csub_cg_calc_stress(nodes, hnew)
4389 this%cg_es0(node) = this%cg_es(node)
4393 do ib = 1, this%ninterbeds
4394 idelay = this%idelay(ib)
4395 node = this%nodelist(ib)
4396 top = this%dis%top(node)
4397 bot = this%dis%bot(node)
4401 if (this%ispecified_pcs == 0)
then
4403 if (this%ipch /= 0)
then
4404 pcs = this%cg_es(node) - pcs0
4406 pcs = this%cg_es(node) + pcs0
4410 if (this%ipch /= 0)
then
4411 pcs = this%cg_gs(node) - (pcs0 - bot)
4413 if (pcs < this%cg_es(node))
then
4414 pcs = this%cg_es(node)
4420 if (idelay /= 0)
then
4421 dzhalf =
dhalf * this%dbdzini(1, idelay)
4426 do n = 1, this%ndelaycells
4427 if (this%ispecified_dbh == 0)
then
4428 this%dbh(n, idelay) = hcell + this%dbh(n, idelay)
4430 this%dbh(n, idelay) = hcell
4432 this%dbh0(n, idelay) = this%dbh(n, idelay)
4436 call this%csub_delay_calc_stress(ib, hcell)
4440 do n = 1, this%ndelaycells
4441 zbot = this%dbz(n, idelay) - dzhalf
4444 dbpcs = pcs - (zbot - bot) * (this%sgs(node) -
done)
4445 this%dbpcs(n, idelay) = dbpcs
4448 this%dbes0(n, idelay) = this%dbes(n, idelay)
4455 top = this%dis%top(node)
4456 bot = this%dis%bot(node)
4459 if (this%istoragec == 1)
then
4462 if (this%lhead_based .EQV. .true.)
then
4468 void_ratio = this%csub_calc_void_ratio(this%cg_theta(node))
4469 es = this%cg_es(node)
4476 znode = this%csub_calc_znode(top, bot, hbar)
4477 fact = this%csub_calc_adjes(node, es, bot, znode)
4478 fact = fact * (
done + void_ratio)
4485 this%cg_ske_cr(node) = this%cg_ske_cr(node) * fact
4488 if (fact <=
dzero)
then
4489 call this%dis%noder_to_string(node, cellid)
4490 write (
errmsg,
'(a,1x,a,a)') &
4491 'Negative recompression index calculated for cell', &
4492 trim(adjustl(cellid)),
'.'
4498 do ib = 1, this%ninterbeds
4499 idelay = this%idelay(ib)
4500 node = this%nodelist(ib)
4501 top = this%dis%top(node)
4502 bot = this%dis%bot(node)
4505 if (this%istoragec == 1)
then
4508 if (this%lhead_based .EQV. .true.)
then
4514 void_ratio = this%csub_calc_void_ratio(this%theta(ib))
4515 es = this%cg_es(node)
4522 znode = this%csub_calc_znode(top, bot, hbar)
4523 fact = this%csub_calc_adjes(node, es, bot, znode)
4524 fact = fact * (
done + void_ratio)
4531 this%ci(ib) = this%ci(ib) * fact
4532 this%rci(ib) = this%rci(ib) * fact
4535 if (fact <=
dzero)
then
4536 call this%dis%noder_to_string(node, cellid)
4537 write (
errmsg,
'(a,1x,i0,2(1x,a),a)') &
4538 'Negative compression indices calculated for interbed', ib, &
4539 'in cell', trim(adjustl(cellid)),
'.'
4545 if (this%iprpak == 1)
then
4547 title = trim(adjustl(this%packName))// &
4548 ' PACKAGE CALCULATED INITIAL INTERBED STRESSES AT THE CELL BOTTOM'
4551 ntabrows = this%ninterbeds
4553 if (this%inamedbound /= 0)
then
4554 ntabcols = ntabcols + 1
4558 call table_cr(this%inputtab, this%packName, title)
4559 call this%inputtab%table_df(ntabrows, ntabcols, this%iout)
4562 tag =
'INTERBED NUMBER'
4563 call this%inputtab%initialize_column(tag, 10, alignment=
tableft)
4565 call this%inputtab%initialize_column(tag, 20)
4566 tag =
'GEOSTATIC STRESS'
4567 call this%inputtab%initialize_column(tag, 16)
4568 tag =
'EFFECTIVE STRESS'
4569 call this%inputtab%initialize_column(tag, 16)
4570 tag =
'PRECONSOLIDATION STRESS'
4571 call this%inputtab%initialize_column(tag, 16)
4572 if (this%inamedbound /= 0)
then
4574 call this%inputtab%initialize_column(tag,
lenboundname, &
4579 do ib = 1, this%ninterbeds
4580 node = this%nodelist(ib)
4581 call this%dis%noder_to_string(node, cellid)
4584 call this%inputtab%add_term(ib)
4585 call this%inputtab%add_term(cellid)
4586 call this%inputtab%add_term(this%cg_gs(node))
4587 call this%inputtab%add_term(this%cg_es(node))
4588 call this%inputtab%add_term(this%pcs(ib))
4589 if (this%inamedbound /= 0)
then
4590 call this%inputtab%add_term(this%boundname(ib))
4597 title = trim(adjustl(this%packName))// &
4598 ' PACKAGE CALCULATED INITIAL DELAY INTERBED STRESSES'
4602 do ib = 1, this%ninterbeds
4603 idelay = this%idelay(ib)
4604 if (idelay /= 0)
then
4605 ntabrows = ntabrows + this%ndelaycells
4609 if (this%inamedbound /= 0)
then
4610 ntabcols = ntabcols + 1
4614 call table_cr(this%inputtab, this%packName, title)
4615 call this%inputtab%table_df(ntabrows, ntabcols, this%iout)
4618 tag =
'INTERBED NUMBER'
4619 call this%inputtab%initialize_column(tag, 10, alignment=
tableft)
4621 call this%inputtab%initialize_column(tag, 20)
4623 call this%inputtab%initialize_column(tag, 10, alignment=
tableft)
4624 tag =
'GEOSTATIC STRESS'
4625 call this%inputtab%initialize_column(tag, 16)
4626 tag =
'EFFECTIVE STRESS'
4627 call this%inputtab%initialize_column(tag, 16)
4628 tag =
'PRECONSOLIDATION STRESS'
4629 call this%inputtab%initialize_column(tag, 16)
4630 if (this%inamedbound /= 0)
then
4632 call this%inputtab%initialize_column(tag,
lenboundname, &
4637 do ib = 1, this%ninterbeds
4638 idelay = this%idelay(ib)
4639 if (idelay /= 0)
then
4640 node = this%nodelist(ib)
4641 call this%dis%noder_to_string(node, cellid)
4644 do n = 1, this%ndelaycells
4646 call this%inputtab%add_term(ib)
4647 call this%inputtab%add_term(cellid)
4649 call this%inputtab%add_term(
' ')
4650 call this%inputtab%add_term(
' ')
4652 call this%inputtab%add_term(n)
4653 call this%inputtab%add_term(this%dbgeo(n, idelay))
4654 call this%inputtab%add_term(this%dbes(n, idelay))
4655 call this%inputtab%add_term(this%dbpcs(n, idelay))
4656 if (this%inamedbound /= 0)
then
4658 call this%inputtab%add_term(this%boundname(ib))
4660 call this%inputtab%add_term(
' ')
4668 if (this%istoragec == 1)
then
4669 if (this%lhead_based .EQV. .false.)
then
4671 title = trim(adjustl(this%packName))// &
4672 ' PACKAGE COMPRESSION INDICES'
4675 ntabrows = this%ninterbeds
4677 if (this%inamedbound /= 0)
then
4678 ntabcols = ntabcols + 1
4682 call table_cr(this%inputtab, this%packName, title)
4683 call this%inputtab%table_df(ntabrows, ntabcols, this%iout)
4686 tag =
'INTERBED NUMBER'
4687 call this%inputtab%initialize_column(tag, 10, alignment=
tableft)
4689 call this%inputtab%initialize_column(tag, 20)
4691 call this%inputtab%initialize_column(tag, 16)
4693 call this%inputtab%initialize_column(tag, 16)
4694 if (this%inamedbound /= 0)
then
4696 call this%inputtab%initialize_column(tag,
lenboundname, &
4701 do ib = 1, this%ninterbeds
4703 node = this%nodelist(ib)
4704 call this%dis%noder_to_string(node, cellid)
4707 call this%inputtab%add_term(ib)
4708 call this%inputtab%add_term(cellid)
4709 call this%inputtab%add_term(this%ci(ib) * fact)
4710 call this%inputtab%add_term(this%rci(ib) * fact)
4711 if (this%inamedbound /= 0)
then
4712 call this%inputtab%add_term(this%boundname(ib))
4721 call this%parser%StoreErrorUnit()
4725 this%initialized = 1
4728 if (this%lhead_based .EQV. .true.)
then
4729 this%iupdatestress = 0
4745 subroutine csub_cg_fc(this, node, tled, area, hcell, hcellold, hcof, rhs)
4748 integer(I4B),
intent(in) :: node
4749 real(DP),
intent(in) :: tled
4750 real(DP),
intent(in) :: area
4751 real(DP),
intent(in) :: hcell
4752 real(DP),
intent(in) :: hcellold
4753 real(DP),
intent(inout) :: hcof
4754 real(DP),
intent(inout) :: rhs
4770 top = this%dis%top(node)
4771 bot = this%dis%bot(node)
4772 tthk = this%cg_thickini(node)
4775 if (tthk >
dzero)
then
4778 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
4784 call this%csub_cg_calc_sske(node, sske, hcell)
4785 rho1 = sske * area * tthk * tled
4788 this%cg_ske(node) = sske * tthk * snold
4789 this%cg_sk(node) = sske * tthk * snnew
4792 hcof = -rho1 * snnew
4793 rhs = rho1 * snold * this%cg_es0(node) - &
4794 rho1 * snnew * (this%cg_gs(node) + bot)
4797 rhs = rhs - rho1 * snnew * (hcell - hbar)
4817 integer(I4B),
intent(in) :: node
4818 real(DP),
intent(in) :: tled
4819 real(DP),
intent(in) :: area
4820 real(DP),
intent(in) :: hcell
4821 real(DP),
intent(inout) :: hcof
4822 real(DP),
intent(inout) :: rhs
4831 real(DP) :: hbarderv
4840 top = this%dis%top(node)
4841 bot = this%dis%bot(node)
4842 tthk = this%cg_thickini(node)
4845 if (tthk >
dzero)
then
4848 call this%csub_calc_sat(node, hcell, top, snnew, snold)
4851 satderv = this%csub_calc_sat_derivative(node, hcell)
4860 call this%csub_cg_calc_sske(node, sske, hcell)
4861 rho1 = sske * area * tthk * tled
4864 hcof = rho1 * snnew * (
done - hbarderv) + &
4865 rho1 * (this%cg_gs(node) - hbar + bot) * satderv
4868 if (this%ieslag /= 0)
then
4869 hcof = hcof - rho1 * this%cg_es0(node) * satderv
4892 integer(I4B),
intent(in) :: ib
4893 integer(I4B),
intent(in) :: node
4894 real(DP),
intent(in) :: area
4895 real(DP),
intent(in) :: hcell
4896 real(DP),
intent(in) :: hcellold
4897 real(DP),
intent(inout) :: hcof
4898 real(DP),
intent(inout) :: rhs
4917 if (this%ibound(node) > 0)
then
4918 if (this%idelay(ib) == 0)
then
4921 if (this%iupdatematprop /= 0)
then
4922 if (this%ieslag == 0)
then
4925 call this%csub_nodelay_calc_comp(ib, hcell, hcellold, comp, &
4927 this%comp(ib) = comp
4930 call this%csub_nodelay_update(ib)
4935 call this%csub_nodelay_fc(ib, hcell, hcellold, rho1, hcof, rhs)
4940 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
4943 if (this%iupdatematprop /= 0)
then
4944 if (this%ieslag == 0)
then
4947 call this%csub_delay_calc_comp(ib, hcell, hcellold, &
4949 this%comp(ib) = comp
4952 call this%csub_delay_update(ib)
4957 call this%csub_delay_sln(ib, hcell)
4958 call this%csub_delay_fc(ib, hcof, rhs)
4959 f = area * this%rnb(ib)
4983 integer(I4B),
intent(in) :: ib
4984 integer(I4B),
intent(in) :: node
4985 real(DP),
intent(in) :: hcell
4986 real(DP),
intent(in) :: hcellold
4987 real(DP),
intent(inout) :: hcof
4988 real(DP),
intent(inout) :: rhs
4990 integer(I4B) :: idelay
5002 real(DP) :: hbarderv
5012 idelay = this%idelay(ib)
5013 top = this%dis%top(node)
5014 bot = this%dis%bot(node)
5017 if (this%ibound(node) > 0)
then
5019 tthk = this%thickini(ib)
5022 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
5025 if (idelay == 0)
then
5031 satderv = this%csub_calc_sat_derivative(node, hcell)
5040 call this%csub_nodelay_fc(ib, hcell, hcellold, rho1, rho2, rhsn)
5043 hcofn = rho2 * (
done - hbarderv) * snnew + &
5044 rho2 * (this%cg_gs(node) - hbar + bot) * satderv
5045 if (this%ielastic(ib) == 0)
then
5046 hcofn = hcofn - rho2 * this%pcs(ib) * satderv
5050 if (this%ieslag /= 0)
then
5051 if (this%ielastic(ib) /= 0)
then
5052 hcofn = hcofn - rho1 * this%cg_es0(node) * satderv
5054 hcofn = hcofn - rho1 * (this%pcs(ib) - this%cg_es0(node)) * satderv
5074 integer(I4B),
intent(in) :: n
5075 real(DP),
intent(inout) :: sske
5076 real(DP),
intent(in) :: hcell
5092 if (this%lhead_based .EQV. .true.)
then
5098 top = this%dis%top(n)
5099 bot = this%dis%bot(n)
5105 znode = this%csub_calc_znode(top, bot, hbar)
5109 es0 = this%cg_es0(n)
5110 theta = this%cg_thetaini(n)
5115 call this%csub_calc_sfacts(n, bot, znode, theta, es, es0, f)
5117 sske = f * this%cg_ske_cr(n)
5133 integer(I4B),
intent(in) :: node
5134 real(DP),
intent(in) :: hcell
5135 real(DP),
intent(in) :: hcellold
5136 real(DP),
intent(inout) :: comp
5148 call this%csub_cg_fc(node, tled, area, hcell, hcellold, hcof, rhs)
5151 comp = hcof * hcell - rhs
5165 integer(I4B),
intent(in) :: node
5167 character(len=20) :: cellid
5173 comp = this%cg_tcomp(node) + this%cg_comp(node)
5174 call this%dis%noder_to_string(node, cellid)
5175 if (abs(comp) >
dzero)
then
5176 thick = this%cg_thickini(node)
5177 theta = this%cg_thetaini(node)
5178 call this%csub_adj_matprop(comp, thick, theta)
5179 if (thick <=
dzero)
then
5180 write (
errmsg,
'(a,1x,a,1x,a,g0,a)') &
5181 'Adjusted thickness for cell', trim(adjustl(cellid)), &
5182 'is less than or equal to 0 (', thick,
').'
5185 if (theta <=
dzero)
then
5186 write (
errmsg,
'(a,1x,a,1x,a,g0,a)') &
5187 'Adjusted theta for cell', trim(adjustl(cellid)), &
5188 'is less than or equal to 0 (', theta,
').'
5191 this%cg_thick(node) = thick
5192 this%cg_theta(node) = theta
5213 integer(I4B),
intent(in) :: node
5214 real(DP),
intent(in) :: tled
5215 real(DP),
intent(in) :: area
5216 real(DP),
intent(in) :: hcell
5217 real(DP),
intent(in) :: hcellold
5218 real(DP),
intent(inout) :: hcof
5219 real(DP),
intent(inout) :: rhs
5235 top = this%dis%top(node)
5236 bot = this%dis%bot(node)
5237 tthk = this%cg_thick(node)
5238 tthk0 = this%cg_thick0(node)
5241 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
5244 wc0 = this%brg * area * tthk0 * this%cg_theta0(node) * tled
5245 wc = this%brg * area * tthk * this%cg_theta(node) * tled
5251 rhs = -wc0 * snold * hcellold
5270 integer(I4B),
intent(in) :: node
5271 real(DP),
intent(in) :: tled
5272 real(DP),
intent(in) :: area
5273 real(DP),
intent(in) :: hcell
5274 real(DP),
intent(in) :: hcellold
5275 real(DP),
intent(inout) :: hcof
5276 real(DP),
intent(inout) :: rhs
5292 top = this%dis%top(node)
5293 bot = this%dis%bot(node)
5294 tthk = this%cg_thick(node)
5297 satderv = this%csub_calc_sat_derivative(node, hcell)
5300 f = this%brg * area * tled
5303 wc = f * tthk * this%cg_theta(node)
5306 hcof = -wc * hcell * satderv
5309 if (this%ieslag /= 0)
then
5310 tthk0 = this%cg_thick0(node)
5311 wc0 = f * tthk0 * this%cg_theta0(node)
5312 hcof = hcof + wc * hcellold * satderv
5333 hcell, hcellold, hcof, rhs)
5336 integer(I4B),
intent(in) :: ib
5337 integer(I4B),
intent(in) :: node
5338 real(DP),
intent(in) :: tled
5339 real(DP),
intent(in) :: area
5340 real(DP),
intent(in) :: hcell
5341 real(DP),
intent(in) :: hcellold
5342 real(DP),
intent(inout) :: hcof
5343 real(DP),
intent(inout) :: rhs
5358 top = this%dis%top(node)
5359 bot = this%dis%bot(node)
5362 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
5365 f = this%brg * area * tled
5366 wc0 = f * this%theta0(ib) * this%thick0(ib)
5367 wc = f * this%theta(ib) * this%thick(ib)
5369 rhs = -wc0 * snold * hcellold
5386 hcell, hcellold, hcof, rhs)
5389 integer(I4B),
intent(in) :: ib
5390 integer(I4B),
intent(in) :: node
5391 real(DP),
intent(in) :: tled
5392 real(DP),
intent(in) :: area
5393 real(DP),
intent(in) :: hcell
5394 real(DP),
intent(in) :: hcellold
5395 real(DP),
intent(inout) :: hcof
5396 real(DP),
intent(inout) :: rhs
5410 top = this%dis%top(node)
5411 bot = this%dis%bot(node)
5414 f = this%brg * area * tled
5417 satderv = this%csub_calc_sat_derivative(node, hcell)
5420 wc = f * this%theta(ib) * this%thick(ib)
5423 hcof = -wc * hcell * satderv
5426 if (this%ieslag /= 0)
then
5427 wc0 = f * this%theta0(ib) * this%thick0(ib)
5428 hcof = hcof + wc0 * hcellold * satderv
5447 real(dp),
intent(in) :: theta
5449 real(dp) :: void_ratio
5451 void_ratio = theta / (
done - theta)
5466 real(dp),
intent(in) :: void_ratio
5471 theta = void_ratio / (
done + void_ratio)
5486 integer(I4B),
intent(in) :: ib
5488 integer(I4B) :: idelay
5492 idelay = this%idelay(ib)
5493 thick = this%thick(ib)
5494 if (idelay /= 0)
then
5495 thick = thick * this%rnb(ib)
5515 real(dp),
intent(in) :: top
5516 real(dp),
intent(in) :: bottom
5517 real(dp),
intent(in) :: zbar
5523 if (zbar > top)
then
5528 znode =
dhalf * (v + bottom)
5545 integer(I4B),
intent(in) :: node
5546 real(dp),
intent(in) :: es0
5547 real(dp),
intent(in) :: z0
5548 real(dp),
intent(in) :: z
5553 es = es0 - (z - z0) * (this%sgs(node) -
done)
5569 integer(I4B),
intent(in) :: ib
5571 integer(I4B) :: iviolate
5572 integer(I4B) :: idelay
5573 integer(I4B) :: node
5582 idelay = this%idelay(ib)
5583 node = this%nodelist(ib)
5586 idelaycells:
do n = 1, this%ndelaycells
5587 z = this%dbz(n, idelay)
5588 h = this%dbh(n, idelay)
5589 dzhalf =
dhalf * this%dbdzini(1, idelay)
5592 if (this%stoiconv(node) == 0)
then
5595 this%idb_nconv_count(1) = this%idb_nconv_count(1) + 1
5601 if (iviolate > 0)
then
5622 integer(I4B),
intent(in) :: node
5623 real(DP),
intent(in) :: hcell
5624 real(DP),
intent(in) :: hcellold
5625 real(DP),
intent(inout) :: snnew
5626 real(DP),
intent(inout) :: snold
5632 if (this%stoiconv(node) /= 0)
then
5633 top = this%dis%top(node)
5634 bot = this%dis%bot(node)
5641 if (this%ieslag /= 0)
then
5659 integer(I4B),
intent(in) :: node
5660 real(dp),
intent(in) :: hcell
5666 if (this%stoiconv(node) /= 0)
then
5667 top = this%dis%top(node)
5668 bot = this%dis%bot(node)
5690 integer(I4B),
intent(in) :: node
5691 real(DP),
intent(in) :: bot
5692 real(DP),
intent(in) :: znode
5693 real(DP),
intent(in) :: theta
5694 real(DP),
intent(in) :: es
5695 real(DP),
intent(in) :: es0
5696 real(DP),
intent(inout) :: fact
5699 real(DP) :: void_ratio
5704 if (this%ieslag /= 0)
then
5711 void_ratio = this%csub_calc_void_ratio(theta)
5712 denom = this%csub_calc_adjes(node, esv, bot, znode)
5713 denom = denom * (
done + void_ratio)
5714 if (denom /=
dzero)
then
5733 real(DP),
intent(in) :: comp
5734 real(DP),
intent(inout) :: thick
5735 real(DP),
intent(inout) :: theta
5738 real(DP) :: void_ratio
5742 void_ratio = this%csub_calc_void_ratio(theta)
5745 if (thick >
dzero) strain = -comp / thick
5748 void_ratio = void_ratio + strain * (
done + void_ratio)
5749 theta = this%csub_calc_theta(void_ratio)
5750 thick = thick - comp
5766 integer(I4B),
intent(in) :: ib
5767 real(DP),
intent(in) :: hcell
5768 logical,
intent(in),
optional :: update
5774 integer(I4B) :: icnvg
5775 integer(I4B) :: iter
5776 integer(I4B) :: idelay
5783 if (
present(update))
then
5790 call this%csub_delay_calc_stress(ib, hcell)
5794 call this%parser%StoreErrorUnit()
5798 if (this%thickini(ib) >
dzero)
then
5801 idelay = this%idelay(ib)
5806 call this%csub_delay_assemble(ib, hcell)
5810 this%dbal, this%dbad, this%dbau, &
5811 this%dbrhs, this%dbdh, this%dbaw)
5815 do n = 1, this%ndelaycells
5816 dh = this%dbdh(n) - this%dbh(n, idelay)
5817 if (abs(dh) > abs(dhmax))
then
5820 this%dbdhmax(idelay) = dhmax
5824 this%dbh(n, idelay) = this%dbdh(n)
5828 call this%csub_delay_calc_stress(ib, hcell)
5831 if (abs(dhmax) < dclose)
then
5833 else if (iter /= 1)
then
5834 if (abs(dhmax) - abs(dhmax0) <
dprec)
then
5838 if (icnvg == 1)
then
5862 integer(I4B),
intent(in) :: ib
5865 integer(I4B) :: node
5866 integer(I4B) :: idelay
5878 idelay = this%idelay(ib)
5879 node = this%nodelist(ib)
5880 b = this%thickini(ib)
5881 bot = this%dis%bot(node)
5887 znode = this%csub_calc_znode(top, bot, hbar)
5888 dz =
dhalf * this%dbdzini(1, idelay)
5895 do n = 1, this%ndelaycells
5898 this%dbz(n, idelay) = z
5902 if (abs(zr) < dz)
then
5905 this%dbrelz(n, idelay) = zr
5923 integer(I4B),
intent(in) :: ib
5924 real(DP),
intent(in) :: hcell
5927 integer(I4B) :: idelay
5928 integer(I4B) :: node
5944 idelay = this%idelay(ib)
5945 node = this%nodelist(ib)
5946 sigma = this%cg_gs(node)
5947 topaq = this%dis%top(node)
5948 botaq = this%dis%bot(node)
5949 dzhalf =
dhalf * this%dbdzini(1, idelay)
5950 top = this%dbz(1, idelay) + dzhalf
5956 sgm = this%sgm(node)
5957 sgs = this%sgs(node)
5958 if (hcell < top)
then
5959 sadd = ((top - hbar) * sgm) + ((hbar - botaq) * sgs)
5961 sadd = (top - botaq) * sgs
5963 sigma = sigma - sadd
5966 do n = 1, this%ndelaycells
5967 h = this%dbh(n, idelay)
5970 z = this%dbz(n, idelay)
5979 sadd = ((top - hbar) * sgm) + ((hbar - bot) * sgs)
5981 sadd = (top - bot) * sgs
5983 sigma = sigma + sadd
5985 this%dbgeo(n, idelay) = sigma
5986 this%dbes(n, idelay) = sigma - phead
6006 integer(I4B),
intent(in) :: ib
6007 integer(I4B),
intent(in) :: n
6008 real(DP),
intent(in) :: hcell
6009 real(DP),
intent(inout) :: ssk
6010 real(DP),
intent(inout) :: sske
6012 integer(I4B) :: idelay
6013 integer(I4B) :: ielastic
6014 integer(I4B) :: node
6017 real(DP) :: hbarcell
6036 idelay = this%idelay(ib)
6037 ielastic = this%ielastic(ib)
6040 if (this%lhead_based .EQV. .true.)
then
6046 node = this%nodelist(ib)
6047 theta = this%dbthetaini(n, idelay)
6050 topcell = this%dis%top(node)
6051 botcell = this%dis%bot(node)
6058 zcell = this%csub_calc_znode(topcell, botcell, hbarcell)
6061 zcenter = zcell + this%dbrelz(n, idelay)
6062 dzhalf =
dhalf * this%dbdzini(1, idelay)
6063 top = zcenter + dzhalf
6064 bot = zcenter - dzhalf
6065 h = this%dbh(n, idelay)
6072 znode = this%csub_calc_znode(top, bot, hbar)
6076 zbot = this%dbz(n, idelay) - dzhalf
6079 es = this%dbes(n, idelay)
6080 es0 = this%dbes0(n, idelay)
6085 call this%csub_calc_sfacts(node, zbot, znode, theta, es, es0, f)
6087 this%idbconvert(n, idelay) = 0
6088 sske = f * this%rci(ib)
6089 ssk = f * this%rci(ib)
6090 if (ielastic == 0)
then
6091 if (this%dbes(n, idelay) > this%dbpcs(n, idelay))
then
6092 this%idbconvert(n, idelay) = 1
6093 ssk = f * this%ci(ib)
6111 integer(I4B),
intent(in) :: ib
6112 real(DP),
intent(in) :: hcell
6121 do n = 1, this%ndelaycells
6124 if (this%inewton == 0)
then
6125 call this%csub_delay_assemble_fc(ib, n, hcell, aii, au, al, r)
6127 call this%csub_delay_assemble_fn(ib, n, hcell, aii, au, al, r)
6153 integer(I4B),
intent(in) :: ib
6154 integer(I4B),
intent(in) :: n
6155 real(DP),
intent(in) :: hcell
6156 real(DP),
intent(inout) :: aii
6157 real(DP),
intent(inout) :: au
6158 real(DP),
intent(inout) :: al
6159 real(DP),
intent(inout) :: r
6161 integer(I4B) :: node
6162 integer(I4B) :: idelay
6163 integer(I4B) :: ielastic
6199 idelay = this%idelay(ib)
6200 ielastic = this%ielastic(ib)
6201 node = this%nodelist(ib)
6202 dzini = this%dbdzini(1, idelay)
6203 dzhalf =
dhalf * dzini
6205 c = this%kv(ib) / dzini
6213 if (n == 1 .or. n == this%ndelaycells)
then
6224 if (n < this%ndelaycells)
then
6229 z = this%dbz(n, idelay)
6232 h = this%dbh(n, idelay)
6233 h0 = this%dbh0(n, idelay)
6234 dz = this%dbdz(n, idelay)
6235 dz0 = this%dbdz0(n, idelay)
6236 theta = this%dbtheta(n, idelay)
6237 theta0 = this%dbtheta0(n, idelay)
6243 call this%csub_delay_calc_sat(node, idelay, n, h, h0, dsn, dsn0)
6246 call this%csub_delay_calc_ssksske(ib, n, hcell, ssk, sske)
6249 smult = dzini * tled
6250 gs = this%dbgeo(n, idelay)
6251 es0 = this%dbes0(n, idelay)
6252 pcs = this%dbpcs(n, idelay)
6253 aii = aii - smult * dsn * ssk
6254 if (ielastic /= 0)
then
6256 (dsn * ssk * (gs + zbot) - dsn0 * sske * es0)
6259 (dsn * ssk * (gs + zbot - pcs) + dsn0 * sske * (pcs - es0))
6263 r = r + smult * dsn * ssk * (h - hbar)
6266 wcf = this%brg * tled
6267 wc = dz * wcf * theta
6268 wc0 = dz0 * wcf * theta0
6269 aii = aii - dsn * wc
6270 r = r - dsn0 * wc0 * h0
6288 integer(I4B),
intent(in) :: ib
6289 integer(I4B),
intent(in) :: n
6290 real(DP),
intent(in) :: hcell
6291 real(DP),
intent(inout) :: aii
6292 real(DP),
intent(inout) :: au
6293 real(DP),
intent(inout) :: al
6294 real(DP),
intent(inout) :: r
6296 integer(I4B) :: node
6297 integer(I4B) :: idelay
6298 integer(I4B) :: ielastic
6324 real(DP) :: hbarderv
6340 idelay = this%idelay(ib)
6341 ielastic = this%ielastic(ib)
6342 node = this%nodelist(ib)
6343 dzini = this%dbdzini(1, idelay)
6344 dzhalf =
dhalf * dzini
6346 c = this%kv(ib) / dzini
6354 if (n == 1 .or. n == this%ndelaycells)
then
6365 if (n < this%ndelaycells)
then
6370 z = this%dbz(n, idelay)
6373 h = this%dbh(n, idelay)
6374 h0 = this%dbh0(n, idelay)
6375 dz = this%dbdz(n, idelay)
6376 dz0 = this%dbdz0(n, idelay)
6377 theta = this%dbtheta(n, idelay)
6378 theta0 = this%dbtheta0(n, idelay)
6387 call this%csub_delay_calc_sat(node, idelay, n, h, h0, dsn, dsn0)
6390 dsnderv = this%csub_delay_calc_sat_derivative(node, idelay, n, hcell)
6393 call this%csub_delay_calc_ssksske(ib, n, hcell, ssk, sske)
6396 smult = dzini * tled
6397 gs = this%dbgeo(n, idelay)
6398 es0 = this%dbes0(n, idelay)
6399 pcs = this%dbpcs(n, idelay)
6400 if (ielastic /= 0)
then
6401 qsto = smult * (dsn * ssk * (gs - hbar + zbot) - dsn0 * sske * es0)
6402 stoderv = -smult * dsn * ssk * hbarderv + &
6403 smult * ssk * (gs - hbar + zbot) * dsnderv
6405 qsto = smult * (dsn * ssk * (gs - hbar + zbot - pcs) + &
6406 dsn0 * sske * (pcs - es0))
6407 stoderv = -smult * dsn * ssk * hbarderv + &
6408 smult * ssk * (gs - hbar + zbot - pcs) * dsnderv
6412 if (this%ieslag /= 0)
then
6413 if (ielastic /= 0)
then
6414 stoderv = stoderv - smult * sske * es0 * dsnderv
6416 stoderv = stoderv + smult * sske * (pcs - es0) * dsnderv
6422 r = r - qsto + stoderv * h
6425 wcf = this%brg * tled
6426 wc = dz * wcf * theta
6427 wc0 = dz0 * wcf * theta0
6428 qwc = dsn0 * wc0 * h0 - dsn * wc * h
6429 wcderv = -dsn * wc - wc * h * dsnderv
6432 if (this%ieslag /= 0)
then
6433 wcderv = wcderv + wc0 * h0 * dsnderv
6438 r = r - qwc + wcderv * h
6457 integer(I4B),
intent(in) :: node
6458 integer(I4B),
intent(in) :: idelay
6459 integer(I4B),
intent(in) :: n
6460 real(DP),
intent(in) :: hcell
6461 real(DP),
intent(in) :: hcellold
6462 real(DP),
intent(inout) :: snnew
6463 real(DP),
intent(inout) :: snold
6470 if (this%stoiconv(node) /= 0)
then
6471 dzhalf =
dhalf * this%dbdzini(n, idelay)
6472 top = this%dbz(n, idelay) + dzhalf
6473 bot = this%dbz(n, idelay) - dzhalf
6480 if (this%ieslag /= 0)
then
6499 integer(I4B),
intent(in) :: node
6500 integer(I4B),
intent(in) :: idelay
6501 integer(I4B),
intent(in) :: n
6502 real(dp),
intent(in) :: hcell
6509 if (this%stoiconv(node) /= 0)
then
6510 dzhalf =
dhalf * this%dbdzini(n, idelay)
6511 top = this%dbz(n, idelay) + dzhalf
6512 bot = this%dbz(n, idelay) - dzhalf
6533 integer(I4B),
intent(in) :: ib
6534 real(DP),
intent(in) :: hcell
6535 real(DP),
intent(inout) :: stoe
6536 real(DP),
intent(inout) :: stoi
6538 integer(I4B) :: idelay
6539 integer(I4B) :: ielastic
6540 integer(I4B) :: node
6559 idelay = this%idelay(ib)
6560 ielastic = this%ielastic(ib)
6561 node = this%nodelist(ib)
6568 if (this%thickini(ib) >
dzero)
then
6569 fmult = this%dbdzini(1, idelay)
6570 dzhalf =
dhalf * this%dbdzini(1, idelay)
6571 do n = 1, this%ndelaycells
6572 call this%csub_delay_calc_ssksske(ib, n, hcell, ssk, sske)
6573 z = this%dbz(n, idelay)
6575 h = this%dbh(n, idelay)
6576 h0 = this%dbh0(n, idelay)
6577 call this%csub_delay_calc_sat(node, idelay, n, h, h0, dsn, dsn0)
6579 if (ielastic /= 0)
then
6580 v1 = dsn * ssk * (this%dbgeo(n, idelay) - hbar + zbot) - &
6581 dsn0 * sske * this%dbes0(n, idelay)
6584 v1 = dsn * ssk * (this%dbgeo(n, idelay) - hbar + zbot - &
6585 this%dbpcs(n, idelay))
6586 v2 = dsn0 * sske * (this%dbpcs(n, idelay) - this%dbes0(n, idelay))
6590 if (this%idbconvert(n, idelay) /= 0)
then
6591 stoi = stoi + v1 * fmult
6592 stoe = stoe + v2 * fmult
6594 stoe = stoe + (v1 + v2) * fmult
6598 ske = ske + sske * fmult
6599 sk = sk + ssk * fmult
6623 integer(I4B),
intent(in) :: ib
6624 real(DP),
intent(inout) :: dwc
6626 integer(I4B) :: idelay
6627 integer(I4B) :: node
6644 if (this%thickini(ib) >
dzero)
then
6645 idelay = this%idelay(ib)
6646 node = this%nodelist(ib)
6648 do n = 1, this%ndelaycells
6649 h = this%dbh(n, idelay)
6650 h0 = this%dbh0(n, idelay)
6651 dz = this%dbdz(n, idelay)
6652 dz0 = this%dbdz0(n, idelay)
6653 call this%csub_delay_calc_sat(node, idelay, n, h, h0, dsn, dsn0)
6654 wc = dz * this%brg * this%dbtheta(n, idelay)
6655 wc0 = dz0 * this%brg * this%dbtheta0(n, idelay)
6656 v = dsn0 * wc0 * h0 - dsn * wc * h
6657 dwc = dwc + v * tled
6677 integer(I4B),
intent(in) :: ib
6678 real(DP),
intent(in) :: hcell
6679 real(DP),
intent(in) :: hcellold
6680 real(DP),
intent(inout) :: comp
6681 real(DP),
intent(inout) :: compi
6682 real(DP),
intent(inout) :: compe
6684 integer(I4B) :: idelay
6685 integer(I4B) :: ielastic
6686 integer(I4B) :: node
6702 idelay = this%idelay(ib)
6703 ielastic = this%ielastic(ib)
6704 node = this%nodelist(ib)
6710 call this%csub_calc_sat(node, hcell, hcellold, snnew, snold)
6713 if (this%thickini(ib) >
dzero)
then
6714 fmult = this%dbdzini(1, idelay)
6715 do n = 1, this%ndelaycells
6716 h = this%dbh(n, idelay)
6717 h0 = this%dbh0(n, idelay)
6718 call this%csub_delay_calc_sat(node, idelay, n, h, h0, dsn, dsn0)
6719 call this%csub_delay_calc_ssksske(ib, n, hcell, ssk, sske)
6720 if (ielastic /= 0)
then
6721 v1 = dsn * ssk * this%dbes(n, idelay) - sske * this%dbes0(n, idelay)
6724 v1 = dsn * ssk * (this%dbes(n, idelay) - this%dbpcs(n, idelay))
6725 v2 = dsn0 * sske * (this%dbpcs(n, idelay) - this%dbes0(n, idelay))
6727 v = (v1 + v2) * fmult
6731 this%dbcomp(n, idelay) = v * snnew
6734 if (this%idbconvert(n, idelay) /= 0)
then
6735 compi = compi + v1 * fmult
6736 compe = compe + v2 * fmult
6738 compe = compe + (v1 + v2) * fmult
6744 comp = comp * this%rnb(ib)
6745 compi = compi * this%rnb(ib)
6746 compe = compe * this%rnb(ib)
6760 integer(I4B),
intent(in) :: ib
6762 integer(I4B) :: idelay
6771 idelay = this%idelay(ib)
6777 do n = 1, this%ndelaycells
6780 comp = this%dbtcomp(n, idelay) + this%dbcomp(n, idelay)
6784 comp = comp / this%rnb(ib)
6787 if (abs(comp) >
dzero)
then
6788 thick = this%dbdzini(n, idelay)
6789 theta = this%dbthetaini(n, idelay)
6790 call this%csub_adj_matprop(comp, thick, theta)
6791 if (thick <=
dzero)
then
6792 write (
errmsg,
'(2(a,i0),a,g0,a)') &
6793 'Adjusted thickness for delay interbed (', ib, &
6794 ') cell (', n,
') is less than or equal to 0 (', thick,
').'
6797 if (theta <=
dzero)
then
6798 write (
errmsg,
'(2(a,i0),a,g0,a)') &
6799 'Adjusted theta for delay interbed (', ib, &
6800 ') cell (', n,
'is less than or equal to 0 (', theta,
').'
6803 this%dbdz(n, idelay) = thick
6804 this%dbtheta(n, idelay) = theta
6805 tthick = tthick + thick
6806 wtheta = wtheta + thick * theta
6808 thick = this%dbdz(n, idelay)
6809 theta = this%dbtheta(n, idelay)
6810 tthick = tthick + thick
6811 wtheta = wtheta + thick * theta
6817 if (tthick >
dzero)
then
6818 wtheta = wtheta / tthick
6823 this%thick(ib) = tthick
6824 this%theta(ib) = wtheta
6843 integer(I4B),
intent(in) :: ib
6844 real(DP),
intent(inout) :: hcof
6845 real(DP),
intent(inout) :: rhs
6847 integer(I4B) :: idelay
6852 idelay = this%idelay(ib)
6855 if (this%thickini(ib) >
dzero)
then
6857 c1 =
dtwo * this%kv(ib) / this%dbdzini(1, idelay)
6858 rhs = -c1 * this%dbh(1, idelay)
6860 this%kv(ib) / this%dbdzini(this%ndelaycells, idelay)
6861 rhs = rhs - c2 * this%dbh(this%ndelaycells, idelay)
6879 integer(I4B),
intent(in) :: ib
6880 integer(I4B),
intent(in) :: n
6881 real(dp),
intent(in) :: hcell
6883 integer(I4B) :: idelay
6888 idelay = this%idelay(ib)
6889 c =
dtwo * this%kv(ib) / this%dbdzini(n, idelay)
6890 q = c * (hcell - this%dbh(n, idelay))
6925 integer(I4B) :: indx
6929 call this%obs%StoreObsType(
'csub', .true., indx)
6934 call this%obs%StoreObsType(
'inelastic-csub', .true., indx)
6939 call this%obs%StoreObsType(
'elastic-csub', .true., indx)
6944 call this%obs%StoreObsType(
'coarse-csub', .false., indx)
6949 call this%obs%StoreObsType(
'csub-cell', .true., indx)
6954 call this%obs%StoreObsType(
'wcomp-csub-cell', .false., indx)
6959 call this%obs%StoreObsType(
'ske', .true., indx)
6964 call this%obs%StoreObsType(
'sk', .true., indx)
6969 call this%obs%StoreObsType(
'ske-cell', .true., indx)
6974 call this%obs%StoreObsType(
'sk-cell', .true., indx)
6979 call this%obs%StoreObsType(
'gstress-cell', .false., indx)
6984 call this%obs%StoreObsType(
'estress-cell', .false., indx)
6989 call this%obs%StoreObsType(
'interbed-compaction', .true., indx)
6994 call this%obs%StoreObsType(
'inelastic-compaction', .true., indx)
6999 call this%obs%StoreObsType(
'elastic-compaction', .true., indx)
7004 call this%obs%StoreObsType(
'coarse-compaction', .false., indx)
7009 call this%obs%StoreObsType(
'inelastic-compaction-cell', .true., indx)
7014 call this%obs%StoreObsType(
'elastic-compaction-cell', .true., indx)
7019 call this%obs%StoreObsType(
'compaction-cell', .true., indx)
7024 call this%obs%StoreObsType(
'thickness', .true., indx)
7029 call this%obs%StoreObsType(
'coarse-thickness', .false., indx)
7034 call this%obs%StoreObsType(
'thickness-cell', .false., indx)
7039 call this%obs%StoreObsType(
'theta', .true., indx)
7044 call this%obs%StoreObsType(
'coarse-theta', .false., indx)
7049 call this%obs%StoreObsType(
'theta-cell', .true., indx)
7054 call this%obs%StoreObsType(
'preconstress-cell', .false., indx)
7059 call this%obs%StoreObsType(
'delay-preconstress', .false., indx)
7064 call this%obs%StoreObsType(
'delay-head', .false., indx)
7069 call this%obs%StoreObsType(
'delay-gstress', .false., indx)
7074 call this%obs%StoreObsType(
'delay-estress', .false., indx)
7079 call this%obs%StoreObsType(
'delay-compaction', .false., indx)
7084 call this%obs%StoreObsType(
'delay-thickness', .false., indx)
7089 call this%obs%StoreObsType(
'delay-theta', .false., indx)
7094 call this%obs%StoreObsType(
'delay-flowtop', .true., indx)
7099 call this%obs%StoreObsType(
'delay-flowbot', .true., indx)
7118 integer(I4B) :: idelay
7119 integer(I4B) :: ncol
7120 integer(I4B) :: node
7126 if (this%obs%npakobs > 0)
then
7127 call this%obs%obs_bd_clear()
7128 do i = 1, this%obs%npakobs
7129 obsrv => this%obs%pakobs(i)%obsrv
7130 if (obsrv%BndFound)
then
7131 if (obsrv%ObsTypeId ==
'SKE' .or. &
7132 obsrv%ObsTypeId ==
'SK' .or. &
7133 obsrv%ObsTypeId ==
'SKE-CELL' .or. &
7134 obsrv%ObsTypeId ==
'SK-CELL' .or. &
7135 obsrv%ObsTypeId ==
'DELAY-HEAD' .or. &
7136 obsrv%ObsTypeId ==
'DELAY-PRECONSTRESS' .or. &
7137 obsrv%ObsTypeId ==
'DELAY-GSTRESS' .or. &
7138 obsrv%ObsTypeId ==
'DELAY-ESTRESS' .or. &
7139 obsrv%ObsTypeId ==
'PRECONSTRESS-CELL')
then
7140 if (this%gwfiss /= 0)
then
7141 call this%obs%SaveOneSimval(obsrv,
dnodata)
7144 do j = 1, obsrv%indxbnds_count
7145 n = obsrv%indxbnds(j)
7146 select case (obsrv%ObsTypeId)
7167 case (
'DELAY-HEAD',
'DELAY-PRECONSTRESS', &
7168 'DELAY-GSTRESS',
'DELAY-ESTRESS')
7169 if (n > this%ndelaycells)
then
7170 r = real(n - 1, dp) / real(this%ndelaycells, dp)
7171 idelay = int(floor(r)) + 1
7172 ncol = n - int(floor(r)) * this%ndelaycells
7177 select case (obsrv%ObsTypeId)
7179 v = this%dbh(ncol, idelay)
7180 case (
'DELAY-PRECONSTRESS')
7181 v = this%dbpcs(ncol, idelay)
7182 case (
'DELAY-GSTRESS')
7183 v = this%dbgeo(ncol, idelay)
7184 case (
'DELAY-ESTRESS')
7185 v = this%dbes(ncol, idelay)
7187 case (
'PRECONSTRESS-CELL')
7190 errmsg =
"Unrecognized observation type '"// &
7191 trim(obsrv%ObsTypeId)//
"'."
7194 call this%obs%SaveOneSimval(obsrv, v)
7199 do j = 1, obsrv%indxbnds_count
7200 n = obsrv%indxbnds(j)
7201 select case (obsrv%ObsTypeId)
7203 v = this%storagee(n) + this%storagei(n)
7204 case (
'INELASTIC-CSUB')
7205 v = this%storagei(n)
7206 case (
'ELASTIC-CSUB')
7207 v = this%storagee(n)
7208 case (
'COARSE-CSUB')
7210 case (
'WCOMP-CSUB-CELL')
7211 v = this%cell_wcstor(n)
7218 v = this%storagee(n) + this%storagei(n)
7222 case (
'COARSE-THETA')
7223 v = this%cg_theta(n)
7228 f = this%cg_thick(n) / this%cell_thick(n)
7229 v = f * this%cg_theta(n)
7231 node = this%nodelist(n)
7232 f = this%csub_calc_interbed_thickness(n) / this%cell_thick(node)
7233 v = f * this%theta(n)
7235 case (
'GSTRESS-CELL')
7237 case (
'ESTRESS-CELL')
7239 case (
'INTERBED-COMPACTION')
7241 case (
'INELASTIC-COMPACTION')
7243 case (
'ELASTIC-COMPACTION')
7245 case (
'COARSE-COMPACTION')
7246 v = this%cg_tcomp(n)
7247 case (
'INELASTIC-COMPACTION-CELL')
7253 case (
'ELASTIC-COMPACTION-CELL')
7257 v = this%cg_tcomp(n)
7261 case (
'COMPACTION-CELL')
7265 v = this%cg_tcomp(n)
7270 idelay = this%idelay(n)
7272 if (idelay /= 0)
then
7275 case (
'COARSE-THICKNESS')
7276 v = this%cg_thick(n)
7277 case (
'THICKNESS-CELL')
7278 v = this%cell_thick(n)
7279 case (
'DELAY-COMPACTION',
'DELAY-THICKNESS', &
7281 if (n > this%ndelaycells)
then
7282 r = real(n, dp) / real(this%ndelaycells, dp)
7283 idelay = int(floor(r)) + 1
7284 ncol = mod(n, this%ndelaycells)
7289 select case (obsrv%ObsTypeId)
7290 case (
'DELAY-COMPACTION')
7291 v = this%dbtcomp(ncol, idelay)
7292 case (
'DELAY-THICKNESS')
7293 v = this%dbdz(ncol, idelay)
7294 case (
'DELAY-THETA')
7295 v = this%dbtheta(ncol, idelay)
7297 case (
'DELAY-FLOWTOP')
7298 idelay = this%idelay(n)
7299 v = this%dbflowtop(idelay)
7300 case (
'DELAY-FLOWBOT')
7301 idelay = this%idelay(n)
7302 v = this%dbflowbot(idelay)
7304 errmsg =
"Unrecognized observation type: '"// &
7305 trim(obsrv%ObsTypeId)//
"'."
7308 call this%obs%SaveOneSimval(obsrv, v)
7312 call this%obs%SaveOneSimval(obsrv,
dnodata)
7318 call this%parser%StoreErrorUnit()
7337 character(len=LENBOUNDNAME) :: bname
7342 integer(I4B) :: idelay
7345 if (.not. this%csub_obs_supported())
then
7353 do i = 1, this%obs%npakobs
7354 obsrv => this%obs%pakobs(i)%obsrv
7357 obsrv%BndFound = .false.
7359 bname = obsrv%FeatureName
7360 if (bname /=
'')
then
7365 do j = 1, this%ninterbeds
7366 if (this%boundname(j) == bname)
then
7367 obsrv%BndFound = .true.
7368 obsrv%CurrentTimeStepEndValue =
dzero
7369 call obsrv%AddObsIndex(j)
7374 else if (obsrv%ObsTypeId ==
'GSTRESS-CELL' .or. &
7375 obsrv%ObsTypeId ==
'ESTRESS-CELL' .or. &
7376 obsrv%ObsTypeId ==
'THICKNESS-CELL' .or. &
7377 obsrv%ObsTypeId ==
'COARSE-CSUB' .or. &
7378 obsrv%ObsTypeId ==
'WCOMP-CSUB-CELL' .or. &
7379 obsrv%ObsTypeId ==
'COARSE-COMPACTION' .or. &
7380 obsrv%ObsTypeId ==
'COARSE-THETA' .or. &
7381 obsrv%ObsTypeId ==
'COARSE-THICKNESS')
then
7382 obsrv%BndFound = .true.
7383 obsrv%CurrentTimeStepEndValue =
dzero
7384 call obsrv%AddObsIndex(obsrv%NodeNumber)
7385 else if (obsrv%ObsTypeId ==
'DELAY-PRECONSTRESS' .or. &
7386 obsrv%ObsTypeId ==
'DELAY-HEAD' .or. &
7387 obsrv%ObsTypeId ==
'DELAY-GSTRESS' .or. &
7388 obsrv%ObsTypeId ==
'DELAY-ESTRESS' .or. &
7389 obsrv%ObsTypeId ==
'DELAY-COMPACTION' .or. &
7390 obsrv%ObsTypeId ==
'DELAY-THICKNESS' .or. &
7391 obsrv%ObsTypeId ==
'DELAY-THETA')
then
7392 if (this%ninterbeds > 0)
then
7393 n = obsrv%NodeNumber
7394 idelay = this%idelay(n)
7395 if (idelay /= 0)
then
7396 j = (idelay - 1) * this%ndelaycells + 1
7397 n2 = obsrv%NodeNumber2
7398 if (n2 < 1 .or. n2 > this%ndelaycells)
then
7399 write (
errmsg,
'(a,2(1x,a),1x,i0,1x,a,i0,a)') &
7400 trim(adjustl(obsrv%ObsTypeId)),
'interbed cell must be ', &
7401 'greater than 0 and less than or equal to', this%ndelaycells, &
7402 '(specified value is ', n2,
').'
7405 j = (idelay - 1) * this%ndelaycells + n2
7407 obsrv%BndFound = .true.
7408 call obsrv%AddObsIndex(j)
7413 else if (obsrv%ObsTypeId ==
'CSUB' .or. &
7414 obsrv%ObsTypeId ==
'INELASTIC-CSUB' .or. &
7415 obsrv%ObsTypeId ==
'ELASTIC-CSUB' .or. &
7416 obsrv%ObsTypeId ==
'SK' .or. &
7417 obsrv%ObsTypeId ==
'SKE' .or. &
7418 obsrv%ObsTypeId ==
'THICKNESS' .or. &
7419 obsrv%ObsTypeId ==
'THETA' .or. &
7420 obsrv%ObsTypeId ==
'INTERBED-COMPACTION' .or. &
7421 obsrv%ObsTypeId ==
'INELASTIC-COMPACTION' .or. &
7422 obsrv%ObsTypeId ==
'ELASTIC-COMPACTION')
then
7423 if (this%ninterbeds > 0)
then
7424 j = obsrv%NodeNumber
7425 if (j < 1 .or. j > this%ninterbeds)
then
7426 write (
errmsg,
'(a,2(1x,a),1x,i0,1x,a,i0,a)') &
7427 trim(adjustl(obsrv%ObsTypeId)),
'interbed cell must be greater', &
7428 'than 0 and less than or equal to', this%ninterbeds, &
7429 '(specified value is ', j,
').'
7432 obsrv%BndFound = .true.
7433 obsrv%CurrentTimeStepEndValue =
dzero
7434 call obsrv%AddObsIndex(j)
7437 else if (obsrv%ObsTypeId ==
'DELAY-FLOWTOP' .or. &
7438 obsrv%ObsTypeId ==
'DELAY-FLOWBOT')
then
7439 if (this%ninterbeds > 0)
then
7440 j = obsrv%NodeNumber
7441 if (j < 1 .or. j > this%ninterbeds)
then
7442 write (
errmsg,
'(a,2(1x,a),1x,i0,1x,a,i0,a)') &
7443 trim(adjustl(obsrv%ObsTypeId)), &
7444 'interbed cell must be greater ', &
7445 'than 0 and less than or equal to', this%ninterbeds, &
7446 '(specified value is ', j,
').'
7449 idelay = this%idelay(j)
7450 if (idelay /= 0)
then
7451 obsrv%BndFound = .true.
7452 obsrv%CurrentTimeStepEndValue =
dzero
7453 call obsrv%AddObsIndex(j)
7461 if (obsrv%ObsTypeId ==
'CSUB-CELL' .or. &
7462 obsrv%ObsTypeId ==
'SKE-CELL' .or. &
7463 obsrv%ObsTypeId ==
'SK-CELL' .or. &
7464 obsrv%ObsTypeId ==
'THETA-CELL' .or. &
7465 obsrv%ObsTypeId ==
'INELASTIC-COMPACTION-CELL' .or. &
7466 obsrv%ObsTypeId ==
'ELASTIC-COMPACTION-CELL' .or. &
7467 obsrv%ObsTypeId ==
'COMPACTION-CELL')
then
7468 if (.NOT. obsrv%BndFound)
then
7469 obsrv%BndFound = .true.
7470 obsrv%CurrentTimeStepEndValue =
dzero
7471 call obsrv%AddObsIndex(obsrv%NodeNumber)
7474 jloop:
do j = 1, this%ninterbeds
7475 if (this%nodelist(j) == obsrv%NodeNumber)
then
7476 obsrv%BndFound = .true.
7477 obsrv%CurrentTimeStepEndValue =
dzero
7478 call obsrv%AddObsIndex(j)
7507 integer(I4B),
intent(in) :: inunitobs
7508 integer(I4B),
intent(in) :: iout
7512 integer(I4B) :: icol, istart, istop
7513 character(len=LINELENGTH) :: string
7514 character(len=LENBOUNDNAME) :: bndname
7515 logical :: flag_string
7518 string = obsrv%IDstring
7526 if (obsrv%ObsTypeId ==
'CSUB' .or. &
7527 obsrv%ObsTypeId ==
'INELASTIC-CSUB' .or. &
7528 obsrv%ObsTypeId ==
'ELASTIC-CSUB' .or. &
7529 obsrv%ObsTypeId ==
'SK' .or. &
7530 obsrv%ObsTypeId ==
'SKE' .or. &
7531 obsrv%ObsTypeId ==
'THETA' .or. &
7532 obsrv%ObsTypeId ==
'THICKNESS' .or. &
7533 obsrv%ObsTypeId ==
'INTERBED-COMPACTION' .or. &
7534 obsrv%ObsTypeId ==
'INELASTIC-COMPACTION' .or. &
7535 obsrv%ObsTypeId ==
'ELASTIC-COMPACTION' .or. &
7536 obsrv%ObsTypeId ==
'DELAY-HEAD' .or. &
7537 obsrv%ObsTypeId ==
'DELAY-GSTRESS' .or. &
7538 obsrv%ObsTypeId ==
'DELAY-ESTRESS' .or. &
7539 obsrv%ObsTypeId ==
'DELAY-PRECONSTRESS' .or. &
7540 obsrv%ObsTypeId ==
'DELAY-COMPACTION' .or. &
7541 obsrv%ObsTypeId ==
'DELAY-THICKNESS' .or. &
7542 obsrv%ObsTypeId ==
'DELAY-THETA' .or. &
7543 obsrv%ObsTypeId ==
'DELAY-FLOWTOP' .or. &
7544 obsrv%ObsTypeId ==
'DELAY-FLOWBOT')
then
7547 nn1 = dis%noder_from_string(icol, istart, istop, inunitobs, &
7548 iout, string, flag_string)
7551 obsrv%FeatureName = bndname
7553 if (obsrv%ObsTypeId ==
'DELAY-HEAD' .or. &
7554 obsrv%ObsTypeId ==
'DELAY-GSTRESS' .or. &
7555 obsrv%ObsTypeId ==
'DELAY-ESTRESS' .or. &
7556 obsrv%ObsTypeId ==
'DELAY-PRECONSTRESS' .or. &
7557 obsrv%ObsTypeId ==
'DELAY-COMPACTION' .or. &
7558 obsrv%ObsTypeId ==
'DELAY-THICKNESS' .or. &
7559 obsrv%ObsTypeId ==
'DELAY-THETA')
then
7562 obsrv%FeatureName = bndname
7566 obsrv%NodeNumber2 = nn2
7572 obsrv%NodeNumber = nn1
7589 this%listlabel = trim(this%filtyp)//
' NO.'
7590 if (this%dis%ndim == 3)
then
7591 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
7592 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'ROW'
7593 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'COL'
7594 elseif (this%dis%ndim == 2)
then
7595 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'LAYER'
7596 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'CELL2D'
7598 write (this%listlabel,
'(a, a7)') trim(this%listlabel),
'NODE'
7600 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'SIG0'
7601 if (this%inamedbound == 1)
then
7602 write (this%listlabel,
'(a, a16)') trim(this%listlabel),
'BOUNDARY NAME'
This module contains block parser methods.
This module contains the BudgetModule.
subroutine, public rate_accumulator(flow, rin, rout)
@ brief Rate accumulator subroutine
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
@ tabcenter
centered table column
@ tabright
right justified table column
@ tableft
left justified table column
@ mnormal
normal output mode
real(dp), parameter dem20
real constant 1e-20
@ tabucstring
upper case string table data
@ tabstring
string table data
@ tabinteger
integer table data
integer(i4b), parameter lenpackagename
maximum length of the package name
real(dp), parameter dp9
real constant 9/10
real(dp), parameter dem10
real constant 1e-10
real(dp), parameter dem7
real constant 1e-7
real(dp), parameter dem8
real constant 1e-8
integer(i4b), parameter namedboundflag
named bound flag
real(dp), parameter dnodata
real no data constant
real(dp), parameter dhnoflo
real no flow constant
integer(i4b), parameter lenlistlabel
maximum length of a llist label
real(dp), parameter dhundred
real constant 100
integer(i4b), parameter lenpakloc
maximum length of a package location
real(dp), parameter dem1
real constant 1e-1
real(dp), parameter dhalf
real constant 1/2
integer(i4b), parameter lenftype
maximum length of a package type (DIS, WEL, OC, etc.)
real(dp), parameter dgravity
real constant gravitational acceleration (m/(s s))
integer(i4b), parameter lenauxname
maximum length of a aux variable
integer(i4b), parameter lenboundname
maximum length of a bound name
real(dp), parameter dem4
real constant 1e-4
real(dp), parameter dem6
real constant 1e-6
real(dp), parameter dzero
real constant zero
real(dp), parameter dten
real constant 10
real(dp), parameter dprec
real constant machine precision
integer(i4b), parameter maxcharlen
maximum length of char string
real(dp), parameter dem15
real constant 1e-15
real(dp), parameter dtwo
real constant 2
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter dthree
real constant 3
real(dp), parameter done
real constant 1
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
This module contains the CSUB package methods.
subroutine csub_nodelay_wcomp_fn(this, ib, node, tled, area, hcell, hcellold, hcof, rhs)
@ brief Formulate no-delay interbed water compressibility coefficients
subroutine csub_read_packagedata(this)
@ brief Read packagedata for package
real(dp) function csub_calc_delay_flow(this, ib, n, hcell)
Calculate the flow from delay interbed top or bottom.
subroutine csub_cg_wcomp_fc(this, node, tled, area, hcell, hcellold, hcof, rhs)
@ brief Formulate coarse-grained water compressibility coefficients
subroutine csub_calc_sfacts(this, node, bot, znode, theta, es, es0, fact)
Calculate specific storage coefficient factor.
subroutine csub_read_dimensions(this)
@ brief Read dimensions for package
subroutine csub_delay_assemble_fn(this, ib, n, hcell, aii, au, al, r)
Assemble delay interbed Newton-Raphson formulation coefficients.
subroutine csub_ar(this, dis, ibound)
@ brief Allocate and read method for package
subroutine csub_nodelay_wcomp_fc(this, ib, node, tled, area, hcell, hcellold, hcof, rhs)
@ brief Formulate no-delay interbed water compressibility coefficients
real(dp) function csub_calc_sat_derivative(this, node, hcell)
Calculate the saturation derivative.
character(len=lenbudtxt), dimension(4) budtxt
subroutine read_options(this)
@ brief Read options for package
subroutine csub_cg_calc_comp(this, node, hcell, hcellold, comp)
@ brief Calculate coarse-grained compaction in a cell
real(dp) function csub_calc_adjes(this, node, es0, z0, z)
Calculate the effective stress at elevation z.
subroutine csub_cg_wcomp_fn(this, node, tled, area, hcell, hcellold, hcof, rhs)
@ brief Formulate coarse-grained water compressibility coefficients
subroutine csub_interbed_fc(this, ib, node, area, hcell, hcellold, hcof, rhs)
@ brief Formulate the coefficients for a interbed
subroutine csub_delay_fc(this, ib, hcof, rhs)
Calculate delay interbed contribution to the cell.
subroutine csub_delay_update(this, ib)
Update delay interbed material properties.
subroutine csub_delay_init_zcell(this, ib)
Calculate delay interbed znode and z relative to interbed center.
subroutine csub_nodelay_update(this, i)
@ brief Update no-delay material properties
subroutine csub_allocate_arrays(this)
@ brief Allocate package arrays
subroutine csub_delay_calc_ssksske(this, ib, n, hcell, ssk, sske)
Calculate delay interbed cell storage coefficients.
subroutine csub_adj_matprop(this, comp, thick, theta)
Calculate new material properties.
subroutine csub_cg_calc_sske(this, n, sske, hcell)
@ brief Calculate Sske for a cell
real(dp) function csub_calc_void_ratio(this, theta)
Calculate the void ratio.
subroutine csub_fc(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
@ brief Fill A and r for the package
subroutine csub_calc_sat(this, node, hcell, hcellold, snnew, snold)
Calculate cell saturation.
real(dp) function csub_calc_theta(this, void_ratio)
Calculate the porosity.
subroutine csub_cc(this, innertot, kiter, iend, icnvgmod, nodes, hnew, hold, cpak, ipak, dpak)
@ brief Final convergence check
subroutine csub_delay_calc_wcomp(this, ib, dwc)
Calculate delay interbed water compressibility.
subroutine csub_delay_calc_sat(this, node, idelay, n, hcell, hcellold, snnew, snold)
Calculate delay interbed saturation.
real(dp) function csub_calc_znode(this, top, bottom, zbar)
Calculate the cell node.
subroutine csub_delay_calc_comp(this, ib, hcell, hcellold, comp, compi, compe)
Calculate delay interbed compaction.
subroutine csub_delay_calc_stress(this, ib, hcell)
Calculate delay interbed stress values.
subroutine csub_nodelay_calc_comp(this, ib, hcell, hcellold, comp, rho1, rho2)
@ brief Calculate no-delay interbed compaction
subroutine csub_set_initial_state(this, nodes, hnew)
@ brief Set initial states for the package
subroutine csub_cg_calc_stress(this, nodes, hnew)
@ brief Calculate the stress for model cells
real(dp) function csub_calc_interbed_thickness(this, ib)
Calculate the interbed thickness.
real(dp), parameter dlog10es
derivative of the log of effective stress
subroutine csub_delay_assemble_fc(this, ib, n, hcell, aii, au, al, r)
Assemble delay interbed standard formulation coefficients.
subroutine csub_interbed_fn(this, ib, node, hcell, hcellold, hcof, rhs)
@ brief Formulate the coefficients for a interbed
subroutine csub_rp_obs(this)
Read and prepare the observations.
subroutine csub_rp(this)
@ brief Read and prepare stress period data for package
subroutine csub_nodelay_fc(this, ib, hcell, hcellold, rho1, rho2, rhs, argtled)
@ brief Calculate no-delay interbed storage coefficients
subroutine csub_ad(this, nodes, hnew)
@ brief Advance the package
subroutine csub_bd_obs(this)
Set the observations for this time step.
subroutine csub_cg_update(this, node)
@ brief Update coarse-grained material properties
subroutine csub_delay_assemble(this, ib, hcell)
Assemble delay interbed coefficients.
subroutine csub_bd(this, isuppress_output, model_budget)
@ brief Model budget calculation for package
subroutine define_listlabel(this)
@ brief Define the list label for the package
subroutine, public csub_cr(csubobj, name_model, istounit, stoPckName, inunit, iout)
@ brief Create a new package object
subroutine csub_ot_dv(this, idvfl, idvprint)
@ brief Save and print dependent values for package
real(dp) function csub_delay_calc_sat_derivative(this, node, idelay, n, hcell)
Calculate the delay interbed cell saturation derivative.
subroutine csub_da(this)
@ brief Deallocate package memory
subroutine csub_save_model_flows(this, icbcfl, icbcun)
@ brief Save model flows for package
subroutine csub_cg_fn(this, node, tled, area, hcell, hcof, rhs)
@ brief Formulate coarse-grained Newton-Raphson terms
subroutine csub_delay_head_check(this, ib)
Check delay interbed head.
subroutine csub_delay_sln(this, ib, hcell, update)
Solve delay interbed continuity equation.
subroutine csub_fp(this)
@ brief Final processing for package
subroutine csub_process_obsid(obsrv, dis, inunitobs, iout)
Process the observation IDs for the package.
subroutine csub_fn(this, kiter, hold, hnew, matrix_sln, idxglo, rhs)
@ brief Fill Newton-Raphson terms in A and r for the package
logical function csub_obs_supported(this)
Determine if observations are supported.
character(len=lenbudtxt), dimension(6) comptxt
subroutine csub_delay_calc_dstor(this, ib, hcell, stoe, stoi)
Calculate delay interbed storage change.
subroutine csub_cg_chk_stress(this)
@ brief Check effective stress values
subroutine csub_cg_fc(this, node, tled, area, hcell, hcellold, hcof, rhs)
@ brief Formulate the coefficients for coarse-grained materials
subroutine csub_cq(this, nodes, hnew, hold, isuppress_output, flowja)
@ brief Calculate flows for package
subroutine csub_allocate_scalars(this)
@ brief Allocate scalars
subroutine csub_df_obs(this)
Define the observation types available in the package.
subroutine, public ims_misc_thomas(n, tl, td, tu, b, x, w)
Tridiagonal solve using the Thomas algorithm.
This module defines variable data types.
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
Store and issue logging messages to output units.
subroutine, public write_message(text, iunit, fmt, skipbefore, skipafter, advance)
Write a message to an output unit.
This module contains the base numerical package type.
This module contains the derived types ObserveType and ObsDataType.
This module contains the derived type ObsType.
subroutine, public obs_cr(obs, inobs)
@ brief Create a new ObsType object
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.
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
This module contains simulation variables.
character(len=maxcharlen) errmsg
error message string
character(len=maxcharlen) warnmsg
warning message string
real(dp) function squadraticsaturation(top, bot, x, eps)
@ brief sQuadraticSaturation
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function
real(dp) function squadratic0spderivative(x, xi, tomega)
@ brief sQuadratic0spDerivative
real(dp) function squadratic0sp(x, xi, tomega)
@ brief sQuadratic0sp
subroutine, public selectn(indx, v, reverse)
subroutine, public table_cr(this, name, title)
real(dp), pointer, public totim
time relative to start of simulation
integer(i4b), pointer, public kstp
current time step number
integer(i4b), pointer, public kper
current stress period number
real(dp), pointer, public delt
length of the current time step
integer(i4b), pointer, public nper
number of stress period
type(timeserieslinktype) function, pointer, public gettimeserieslinkfromlist(list, indx)
Get time series link from a list.
subroutine, public read_value_or_time_series_adv(textInput, ii, jj, bndElem, pkgName, auxOrBnd, tsManager, iprpak, varName)
Call this subroutine from advanced packages to define timeseries link for a variable (varName).
subroutine, public tsmanager_cr(this, iout, removeTsLinksOnCompletion, extendTsToEndOfSimulation)
Create the tsmanager.
Derived type for the Budget object.
A generic heterogeneous doubly-linked list.