MODFLOW 6  version 6.7.0.dev0
USGS Modular Hydrologic Model
tsp-fmi.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b
7  use simvariablesmodule, only: errmsg
9  use basedismodule, only: disbasetype
10  use listmodule, only: listtype
16 
17  implicit none
18  private
19  public :: tspfmitype
20  public :: fmi_cr
21 
22  character(len=LENPACKAGENAME) :: text = ' GWTFMI'
23 
24  integer(I4B), parameter :: nbditems = 2
25  character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt
26  data budtxt/' FLOW-ERROR', ' FLOW-CORRECTION'/
27 
29  real(dp), dimension(:), contiguous, pointer :: concpack => null()
30  real(dp), dimension(:), contiguous, pointer :: qmfrommvr => null()
31  end type
32 
34  type(budgetobjecttype), pointer :: ptr
35  end type budobjptrarray
36 
38 
39  integer(I4B), dimension(:), pointer, contiguous :: iatp => null() !< advanced transport package applied to gwfpackages
40  integer(I4B), pointer :: iflowerr => null() !< add the flow error correction
41  real(dp), dimension(:), pointer, contiguous :: flowcorrect => null() !< mass flow correction
42  real(dp), pointer :: eqnsclfac => null() !< governing equation scale factor; =1. for solute; =rhow*cpw for energy
44  dimension(:), pointer, contiguous :: datp => null()
45  type(budobjptrarray), dimension(:), allocatable :: aptbudobj !< flow budget objects for the advanced packages
46 
47  contains
48 
49  procedure :: allocate_arrays => gwtfmi_allocate_arrays
50  procedure :: allocate_gwfpackages => gwtfmi_allocate_gwfpackages
51  procedure :: allocate_scalars => gwtfmi_allocate_scalars
52  procedure :: deallocate_gwfpackages => gwtfmi_deallocate_gwfpackages
53  procedure :: fmi_rp
54  procedure :: fmi_ad
55  procedure :: fmi_fc
56  procedure :: fmi_cq
57  procedure :: fmi_bd
58  procedure :: fmi_ot_flow
59  procedure :: fmi_da => gwtfmi_da
60  procedure :: gwfsatold
63  procedure :: read_options => gwtfmi_read_options
64  procedure :: set_aptbudobj_pointer
65  procedure :: read_packagedata => gwtfmi_read_packagedata
66  procedure :: set_active_status
67 
68  end type tspfmitype
69 
70 contains
71 
72  !> @brief Create a new FMI object
73  !<
74  subroutine fmi_cr(fmiobj, name_model, inunit, iout, eqnsclfac, depvartype)
75  ! -- dummy
76  type(tspfmitype), pointer :: fmiobj
77  character(len=*), intent(in) :: name_model
78  integer(I4B), intent(in) :: inunit
79  integer(I4B), intent(in) :: iout
80  real(dp), intent(in), pointer :: eqnsclfac !< governing equation scale factor
81  character(len=LENVARNAME), intent(in) :: depvartype
82  !
83  ! -- Create the object
84  allocate (fmiobj)
85  !
86  ! -- create name and memory path
87  call fmiobj%set_names(1, name_model, 'FMI', 'FMI')
88  fmiobj%text = text
89  !
90  ! -- Allocate scalars
91  call fmiobj%allocate_scalars()
92  !
93  ! -- Set variables
94  fmiobj%inunit = inunit
95  fmiobj%iout = iout
96  !
97  ! -- Initialize block parser
98  call fmiobj%parser%Initialize(fmiobj%inunit, fmiobj%iout)
99  !
100  ! -- Assign label based on dependent variable
101  fmiobj%depvartype = depvartype
102  !
103  ! -- Store pointer to governing equation scale factor
104  fmiobj%eqnsclfac => eqnsclfac
105  end subroutine fmi_cr
106 
107  !> @brief Read and prepare
108  !<
109  subroutine fmi_rp(this, inmvr)
110  ! -- modules
111  use tdismodule, only: kper, kstp
112  ! -- dummy
113  class(tspfmitype) :: this
114  integer(I4B), intent(in) :: inmvr
115  ! -- local
116  ! -- formats
117  !
118  ! --Check to make sure MVT Package is active if mvr flows are available.
119  ! This cannot be checked until RP because exchange doesn't set a pointer
120  ! to mvrbudobj until exg_ar().
121  if (kper * kstp == 1) then
122  if (associated(this%mvrbudobj) .and. inmvr == 0) then
123  write (errmsg, '(a)') 'GWF water mover is active but the GWT MVT &
124  &package has not been specified. activate GWT MVT package.'
125  call store_error(errmsg, terminate=.true.)
126  end if
127  if (.not. associated(this%mvrbudobj) .and. inmvr > 0) then
128  write (errmsg, '(a)') 'GWF water mover terms are not available &
129  &but the GWT MVT package has been activated. Activate GWF-GWT &
130  &exchange or specify GWFMOVER in FMI PACKAGEDATA.'
131  call store_error(errmsg, terminate=.true.)
132  end if
133  end if
134  end subroutine fmi_rp
135 
136  !> @brief Advance routine for FMI object
137  !<
138  subroutine fmi_ad(this, cnew)
139  ! -- modules
140  use constantsmodule, only: dhdry
141  ! -- dummy
142  class(tspfmitype) :: this
143  real(DP), intent(inout), dimension(:) :: cnew
144  ! -- local
145  integer(I4B) :: n
146  character(len=*), parameter :: fmtdry = &
147  &"(/1X,'WARNING: DRY CELL ENCOUNTERED AT ',a,'; RESET AS INACTIVE &
148  &WITH DRY CONCENTRATION = ', G13.5)"
149  character(len=*), parameter :: fmtrewet = &
150  &"(/1X,'DRY CELL REACTIVATED AT ', a,&
151  &' WITH STARTING CONCENTRATION =',G13.5)"
152  !
153  ! -- Set flag to indicated that flows are being updated. For the case where
154  ! flows may be reused (only when flows are read from a file) then set
155  ! the flag to zero to indicate that flows were not updated
156  this%iflowsupdated = 1
157  !
158  ! -- If reading flows from a budget file, read the next set of records
159  if (this%iubud /= 0) then
160  call this%advance_bfr()
161  end if
162  !
163  ! -- If reading heads from a head file, read the next set of records
164  if (this%iuhds /= 0) then
165  call this%advance_hfr()
166  end if
167  !
168  ! -- If mover flows are being read from file, read the next set of records
169  if (this%iumvr /= 0) then
170  call this%mvrbudobj%bfr_advance(this%dis, this%iout)
171  end if
172  !
173  ! -- If advanced package flows are being read from file, read the next set of records
174  if (this%flows_from_file .and. this%inunit /= 0) then
175  do n = 1, size(this%aptbudobj)
176  call this%aptbudobj(n)%ptr%bfr_advance(this%dis, this%iout)
177  end do
178  end if
179  !
180  ! -- set inactive transport cell status
181  if (this%idryinactive /= 0) then
182  call this%set_active_status(cnew)
183  end if
184  end subroutine fmi_ad
185 
186  !> @brief Calculate coefficients and fill matrix and rhs terms associated
187  !! with FMI object
188  !<
189  subroutine fmi_fc(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
190  ! -- dummy
191  class(tspfmitype) :: this
192  integer, intent(in) :: nodes
193  real(DP), intent(in), dimension(nodes) :: cold
194  integer(I4B), intent(in) :: nja
195  class(matrixbasetype), pointer :: matrix_sln
196  integer(I4B), intent(in), dimension(nja) :: idxglo
197  real(DP), intent(inout), dimension(nodes) :: rhs
198  ! -- local
199  integer(I4B) :: n, idiag, idiag_sln
200  real(DP) :: qcorr
201  !
202  ! -- Calculate the flow imbalance error and make a correction for it
203  if (this%iflowerr /= 0) then
204  !
205  ! -- Correct the transport solution for the flow imbalance by adding
206  ! the flow residual to the diagonal
207  do n = 1, nodes
208  idiag = this%dis%con%ia(n)
209  idiag_sln = idxglo(idiag)
210  !call matrix_sln%add_value_pos(idiag_sln, -this%gwfflowja(idiag))
211  qcorr = -this%gwfflowja(idiag) * this%eqnsclfac
212  call matrix_sln%add_value_pos(idiag_sln, qcorr)
213  end do
214  end if
215  end subroutine fmi_fc
216 
217  !> @brief Calculate flow correction
218  !!
219  !! Where there is a flow imbalance for a given cell, a correction may be
220  !! applied if selected
221  !<
222  subroutine fmi_cq(this, cnew, flowja)
223  ! -- modules
224  ! -- dummy
225  class(tspfmitype) :: this
226  real(DP), intent(in), dimension(:) :: cnew
227  real(DP), dimension(:), contiguous, intent(inout) :: flowja
228  ! -- local
229  integer(I4B) :: n
230  integer(I4B) :: idiag
231  real(DP) :: rate
232  !
233  ! -- If not adding flow error correction, return
234  if (this%iflowerr /= 0) then
235  !
236  ! -- Accumulate the flow correction term
237  do n = 1, this%dis%nodes
238  rate = dzero
239  idiag = this%dis%con%ia(n)
240  if (this%ibound(n) > 0) then
241  rate = -this%gwfflowja(idiag) * cnew(n) * this%eqnsclfac
242  end if
243  this%flowcorrect(n) = rate
244  flowja(idiag) = flowja(idiag) + rate
245  end do
246  end if
247  end subroutine fmi_cq
248 
249  !> @brief Calculate budget terms associated with FMI object
250  !<
251  subroutine fmi_bd(this, isuppress_output, model_budget)
252  ! -- modules
253  use tdismodule, only: delt
255  ! -- dummy
256  class(tspfmitype) :: this
257  integer(I4B), intent(in) :: isuppress_output
258  type(budgettype), intent(inout) :: model_budget
259  ! -- local
260  real(DP) :: rin
261  real(DP) :: rout
262  !
263  ! -- flow correction
264  if (this%iflowerr /= 0) then
265  call rate_accumulator(this%flowcorrect, rin, rout)
266  call model_budget%addentry(rin, rout, delt, budtxt(2), isuppress_output)
267  end if
268  end subroutine fmi_bd
269 
270  !> @brief Save budget terms associated with FMI object
271  !<
272  subroutine fmi_ot_flow(this, icbcfl, icbcun)
273  ! -- dummy
274  class(tspfmitype) :: this
275  integer(I4B), intent(in) :: icbcfl
276  integer(I4B), intent(in) :: icbcun
277  ! -- local
278  integer(I4B) :: ibinun
279  integer(I4B) :: iprint, nvaluesp, nwidthp
280  character(len=1) :: cdatafmp = ' ', editdesc = ' '
281  real(DP) :: dinact
282  !
283  ! -- Set unit number for binary output
284  if (this%ipakcb < 0) then
285  ibinun = icbcun
286  elseif (this%ipakcb == 0) then
287  ibinun = 0
288  else
289  ibinun = this%ipakcb
290  end if
291  if (icbcfl == 0) ibinun = 0
292  !
293  ! -- Do not save flow corrections if not active
294  if (this%iflowerr == 0) ibinun = 0
295  !
296  ! -- Record the storage rates if requested
297  if (ibinun /= 0) then
298  iprint = 0
299  dinact = dzero
300  !
301  ! -- flow correction
302  call this%dis%record_array(this%flowcorrect, this%iout, iprint, -ibinun, &
303  budtxt(2), cdatafmp, nvaluesp, &
304  nwidthp, editdesc, dinact)
305  end if
306  end subroutine fmi_ot_flow
307 
308  !> @brief Deallocate variables
309  !!
310  !! Deallocate memory associated with FMI object
311  !<
312  subroutine gwtfmi_da(this)
313  ! -- modules
315  ! -- dummy
316  class(tspfmitype) :: this
317  ! -- todo: finalize hfr and bfr either here or in a finalize routine
318  !
319  ! -- deallocate any memory stored with gwfpackages
320  call this%deallocate_gwfpackages()
321  !
322  ! -- deallocate fmi arrays
323  if (associated(this%datp)) then
324  deallocate (this%datp)
325  deallocate (this%gwfpackages)
326  deallocate (this%flowpacknamearray)
327  call mem_deallocate(this%iatp)
328  call mem_deallocate(this%igwfmvrterm)
329  end if
330 
331  deallocate (this%aptbudobj)
332  call mem_deallocate(this%flowcorrect)
333  call mem_deallocate(this%ibdgwfsat0)
334  if (this%flows_from_file) then
335  call mem_deallocate(this%gwfstrgss)
336  call mem_deallocate(this%gwfstrgsy)
337  end if
338  !
339  ! -- special treatment, these could be from mem_checkin
340  call mem_deallocate(this%gwfhead, 'GWFHEAD', this%memoryPath)
341  call mem_deallocate(this%gwfsat, 'GWFSAT', this%memoryPath)
342  call mem_deallocate(this%gwfspdis, 'GWFSPDIS', this%memoryPath)
343  call mem_deallocate(this%gwfflowja, 'GWFFLOWJA', this%memoryPath)
344  !
345  ! -- deallocate scalars
346  call mem_deallocate(this%flows_from_file)
347  call mem_deallocate(this%iflowsupdated)
348  call mem_deallocate(this%iflowerr)
349  call mem_deallocate(this%igwfstrgss)
350  call mem_deallocate(this%igwfstrgsy)
351  call mem_deallocate(this%iubud)
352  call mem_deallocate(this%iuhds)
353  call mem_deallocate(this%iumvr)
354  call mem_deallocate(this%nflowpack)
355  call mem_deallocate(this%idryinactive)
356  !
357  ! -- deallocate parent
358  call this%NumericalPackageType%da()
359  end subroutine gwtfmi_da
360 
361  !> @ brief Allocate scalars
362  !!
363  !! Allocate scalar variables for an FMI object
364  !<
365  subroutine gwtfmi_allocate_scalars(this)
366  ! -- modules
368  ! -- dummy
369  class(tspfmitype) :: this
370  ! -- local
371  !
372  ! -- allocate scalars in parent
373  call this%FlowModelInterfaceType%allocate_scalars()
374  !
375  ! -- Allocate
376  call mem_allocate(this%iflowerr, 'IFLOWERR', this%memoryPath)
377  !
378  ! -- Although not a scalar, allocate the advanced package transport
379  ! budget object to zero so that it can be dynamically resized later
380  allocate (this%aptbudobj(0))
381  !
382  ! -- Initialize
383  this%iflowerr = 0
384  end subroutine gwtfmi_allocate_scalars
385 
386  !> @ brief Allocate arrays for FMI object
387  !!
388  !! Method to allocate arrays for the FMI package.
389  !<
390  subroutine gwtfmi_allocate_arrays(this, nodes)
392  ! -- modules
393  use constantsmodule, only: dzero
394  ! -- dummy
395  class(tspfmitype) :: this
396  integer(I4B), intent(in) :: nodes
397  ! -- local
398  integer(I4B) :: n
399  !
400  ! -- allocate parent arrays
401  call this%FlowModelInterfaceType%allocate_arrays(nodes)
402  !
403  ! -- Allocate variables needed for all cases
404  if (this%iflowerr == 0) then
405  call mem_allocate(this%flowcorrect, 1, 'FLOWCORRECT', this%memoryPath)
406  else
407  call mem_allocate(this%flowcorrect, nodes, 'FLOWCORRECT', this%memoryPath)
408  end if
409  do n = 1, size(this%flowcorrect)
410  this%flowcorrect(n) = dzero
411  end do
412  end subroutine gwtfmi_allocate_arrays
413 
414  !> @brief Set gwt transport cell status
415  !!
416  !! Dry GWF cells are treated differently by GWT and GWE. Transport does not
417  !! occur in deactivated GWF cells; however, GWE still simulates conduction
418  !! through dry cells.
419  !<
420  subroutine set_active_status(this, cnew)
421  ! -- modules
422  use constantsmodule, only: dhdry
423  ! -- dummy
424  class(tspfmitype) :: this
425  real(DP), intent(inout), dimension(:) :: cnew
426  ! -- local
427  integer(I4B) :: n
428  integer(I4B) :: m
429  integer(I4B) :: ipos
430  real(DP) :: crewet, tflow, flownm
431  character(len=15) :: nodestr
432  ! -- formats
433  character(len=*), parameter :: fmtoutmsg1 = &
434  "(1x,'WARNING: DRY CELL ENCOUNTERED AT ', a,'; RESET AS INACTIVE WITH &
435  &DRY ', a, '=', G13.5)"
436  character(len=*), parameter :: fmtoutmsg2 = &
437  &"(1x,'DRY CELL REACTIVATED AT', a, 'WITH STARTING', a, '=', G13.5)"
438  !
439  do n = 1, this%dis%nodes
440  ! -- Calculate the ibound-like array that has 0 if saturation
441  ! is zero and 1 otherwise
442  if (this%gwfsat(n) > dzero) then
443  this%ibdgwfsat0(n) = 1
444  else
445  this%ibdgwfsat0(n) = 0
446  end if
447  !
448  ! -- Check if active transport cell is inactive for flow
449  if (this%ibound(n) > 0) then
450  if (this%gwfhead(n) == dhdry) then
451  ! -- transport cell should be made inactive
452  this%ibound(n) = 0
453  cnew(n) = dhdry
454  call this%dis%noder_to_string(n, nodestr)
455  write (this%iout, fmtoutmsg1) &
456  trim(nodestr), trim(adjustl(this%depvartype)), dhdry
457  end if
458  end if
459  end do
460  !
461  ! -- if flow cell is dry, then set gwt%ibound = 0 and conc to dry
462  do n = 1, this%dis%nodes
463  !
464  ! -- Convert dry transport cell to active if flow has rewet
465  if (cnew(n) == dhdry) then
466  if (this%gwfhead(n) /= dhdry) then
467  !
468  ! -- obtain weighted concentration/temperature
469  crewet = dzero
470  tflow = dzero
471  do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1
472  m = this%dis%con%ja(ipos)
473  flownm = this%gwfflowja(ipos)
474  if (flownm > 0) then
475  if (this%ibound(m) /= 0) then
476  crewet = crewet + cnew(m) * flownm ! kluge note: apparently no need to multiply flows by eqnsclfac
477  tflow = tflow + this%gwfflowja(ipos) ! since it will divide out below anyway
478  end if
479  end if
480  end do
481  if (tflow > dzero) then
482  crewet = crewet / tflow
483  else
484  crewet = dzero
485  end if
486  !
487  ! -- cell is now wet
488  this%ibound(n) = 1
489  cnew(n) = crewet
490  call this%dis%noder_to_string(n, nodestr)
491  write (this%iout, fmtoutmsg2) &
492  trim(nodestr), trim(adjustl(this%depvartype)), crewet
493  end if
494  end if
495  end do
496  end subroutine set_active_status
497 
498  !> @brief Calculate the previous saturation level
499  !!
500  !! Calculate the groundwater cell head saturation for the end of
501  !! the last time step
502  !<
503  function gwfsatold(this, n, delt) result(satold)
504  ! -- modules
505  ! -- dummy
506  class(tspfmitype) :: this
507  integer(I4B), intent(in) :: n
508  real(dp), intent(in) :: delt
509  ! -- result
510  real(dp) :: satold
511  ! -- local
512  real(dp) :: vcell
513  real(dp) :: vnew
514  real(dp) :: vold
515  !
516  ! -- calculate the value
517  vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n))
518  vnew = vcell * this%gwfsat(n)
519  vold = vnew
520  if (this%igwfstrgss /= 0) vold = vold + this%gwfstrgss(n) * delt
521  if (this%igwfstrgsy /= 0) vold = vold + this%gwfstrgsy(n) * delt
522  satold = vold / vcell
523  end function gwfsatold
524 
525  !> @brief Read options from input file
526  !<
527  subroutine gwtfmi_read_options(this)
528  ! -- modules
529  use constantsmodule, only: linelength, dem6
532  ! -- dummy
533  class(tspfmitype) :: this
534  ! -- local
535  character(len=LINELENGTH) :: keyword
536  integer(I4B) :: ierr
537  logical :: isfound, endOfBlock
538  character(len=*), parameter :: fmtisvflow = &
539  "(4x,'CELL-BY-CELL FLOW INFORMATION WILL BE SAVED TO BINARY FILE &
540  &WHENEVER ICBCFL IS NOT ZERO AND FLOW IMBALANCE CORRECTION ACTIVE.')"
541  character(len=*), parameter :: fmtifc = &
542  &"(4x,'MASS WILL BE ADDED OR REMOVED TO COMPENSATE FOR FLOW IMBALANCE.')"
543  !
544  ! -- get options block
545  call this%parser%GetBlock('OPTIONS', isfound, ierr, blockrequired=.false., &
546  supportopenclose=.true.)
547  !
548  ! -- parse options block if detected
549  if (isfound) then
550  write (this%iout, '(1x,a)') 'PROCESSING FMI OPTIONS'
551  do
552  call this%parser%GetNextLine(endofblock)
553  if (endofblock) exit
554  call this%parser%GetStringCaps(keyword)
555  select case (keyword)
556  case ('SAVE_FLOWS')
557  this%ipakcb = -1
558  write (this%iout, fmtisvflow)
559  case ('FLOW_IMBALANCE_CORRECTION')
560  write (this%iout, fmtifc)
561  this%iflowerr = 1
562  case default
563  write (errmsg, '(a,a)') 'Unknown FMI option: ', &
564  trim(keyword)
565  call store_error(errmsg)
566  call this%parser%StoreErrorUnit()
567  end select
568  end do
569  write (this%iout, '(1x,a)') 'END OF FMI OPTIONS'
570  end if
571  end subroutine gwtfmi_read_options
572 
573  !> @brief Read PACKAGEDATA block
574  !!
575  !! Read packagedata block from input file
576  !<
577  subroutine gwtfmi_read_packagedata(this)
578  ! -- modules
579  use openspecmodule, only: access, form
583  ! -- dummy
584  class(tspfmitype) :: this
585  ! -- local
586  type(budgetobjecttype), pointer :: budobjptr
587  character(len=LINELENGTH) :: keyword, fname
588  character(len=LENPACKAGENAME) :: pname
589  integer(I4B) :: i
590  integer(I4B) :: ierr
591  integer(I4B) :: inunit
592  integer(I4B) :: iapt
593  logical :: isfound, endOfBlock
594  logical :: blockrequired
595  logical :: exist
596  type(budobjptrarray), dimension(:), allocatable :: tmpbudobj
597  !
598  ! -- initialize
599  iapt = 0
600  blockrequired = .true.
601  !
602  ! -- get options block
603  call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
604  blockrequired=blockrequired, &
605  supportopenclose=.true.)
606  !
607  ! -- parse options block if detected
608  if (isfound) then
609  write (this%iout, '(1x,a)') 'PROCESSING FMI PACKAGEDATA'
610  do
611  call this%parser%GetNextLine(endofblock)
612  if (endofblock) exit
613  call this%parser%GetStringCaps(keyword)
614  select case (keyword)
615  case ('GWFBUDGET')
616  call this%parser%GetStringCaps(keyword)
617  if (keyword /= 'FILEIN') then
618  call store_error('GWFBUDGET keyword must be followed by '// &
619  '"FILEIN" then by filename.')
620  call this%parser%StoreErrorUnit()
621  end if
622  call this%parser%GetString(fname)
623  inunit = getunit()
624  inquire (file=trim(fname), exist=exist)
625  if (.not. exist) then
626  call store_error('Could not find file '//trim(fname))
627  call this%parser%StoreErrorUnit()
628  end if
629  call openfile(inunit, this%iout, fname, 'DATA(BINARY)', form, &
630  access, 'UNKNOWN')
631  this%iubud = inunit
632  call this%initialize_bfr()
633  case ('GWFHEAD')
634  call this%parser%GetStringCaps(keyword)
635  if (keyword /= 'FILEIN') then
636  call store_error('GWFHEAD keyword must be followed by '// &
637  '"FILEIN" then by filename.')
638  call this%parser%StoreErrorUnit()
639  end if
640  call this%parser%GetString(fname)
641  inquire (file=trim(fname), exist=exist)
642  if (.not. exist) then
643  call store_error('Could not find file '//trim(fname))
644  call this%parser%StoreErrorUnit()
645  end if
646  inunit = getunit()
647  call openfile(inunit, this%iout, fname, 'DATA(BINARY)', form, &
648  access, 'UNKNOWN')
649  this%iuhds = inunit
650  call this%initialize_hfr()
651  case ('GWFMOVER')
652  call this%parser%GetStringCaps(keyword)
653  if (keyword /= 'FILEIN') then
654  call store_error('GWFMOVER keyword must be followed by '// &
655  '"FILEIN" then by filename.')
656  call this%parser%StoreErrorUnit()
657  end if
658  call this%parser%GetString(fname)
659  inunit = getunit()
660  call openfile(inunit, this%iout, fname, 'DATA(BINARY)', form, &
661  access, 'UNKNOWN')
662  this%iumvr = inunit
663  call budgetobject_cr_bfr(this%mvrbudobj, 'MVT', this%iumvr, &
664  this%iout)
665  call this%mvrbudobj%fill_from_bfr(this%dis, this%iout)
666  case default
667  !
668  ! --expand the size of aptbudobj, which stores a pointer to the budobj
669  allocate (tmpbudobj(iapt))
670  do i = 1, size(this%aptbudobj)
671  tmpbudobj(i)%ptr => this%aptbudobj(i)%ptr
672  end do
673  deallocate (this%aptbudobj)
674  allocate (this%aptbudobj(iapt + 1))
675  do i = 1, size(tmpbudobj)
676  this%aptbudobj(i)%ptr => tmpbudobj(i)%ptr
677  end do
678  deallocate (tmpbudobj)
679  !
680  ! -- Open the budget file and start filling it
681  iapt = iapt + 1
682  pname = keyword(1:lenpackagename)
683  call this%parser%GetStringCaps(keyword)
684  if (keyword /= 'FILEIN') then
685  call store_error('Package name must be followed by '// &
686  '"FILEIN" then by filename.')
687  call this%parser%StoreErrorUnit()
688  end if
689  call this%parser%GetString(fname)
690  inunit = getunit()
691  call openfile(inunit, this%iout, fname, 'DATA(BINARY)', form, &
692  access, 'UNKNOWN')
693  call budgetobject_cr_bfr(budobjptr, pname, inunit, &
694  this%iout, colconv2=['GWF '])
695  call budobjptr%fill_from_bfr(this%dis, this%iout)
696  this%aptbudobj(iapt)%ptr => budobjptr
697  end select
698  end do
699  write (this%iout, '(1x,a)') 'END OF FMI PACKAGEDATA'
700  end if
701  end subroutine gwtfmi_read_packagedata
702 
703  !> @brief Set the pointer to a budget object
704  !!
705  !! An advanced transport can pass in a name and a
706  !! pointer budget object, and this routine will look through the budget
707  !! objects managed by FMI and point to the one with the same name, such as
708  !! LAK-1, SFR-1, etc.
709  !<
710  subroutine set_aptbudobj_pointer(this, name, budobjptr)
711  ! -- modules
712  class(tspfmitype) :: this
713  ! -- dumm
714  character(len=*), intent(in) :: name
715  type(budgetobjecttype), pointer :: budobjptr
716  ! -- local
717  integer(I4B) :: i
718  !
719  ! -- find and set the pointer
720  do i = 1, size(this%aptbudobj)
721  if (this%aptbudobj(i)%ptr%name == name) then
722  budobjptr => this%aptbudobj(i)%ptr
723  exit
724  end if
725  end do
726  end subroutine set_aptbudobj_pointer
727 
728  !> @brief Initialize the groundwater flow terms based on the budget file
729  !! reader
730  !!
731  !! Initialize terms and figure out how many different terms and packages
732  !! are contained within the file
733  !<
735  ! -- modules
738  ! -- dummy
739  class(tspfmitype) :: this
740  ! -- local
741  integer(I4B) :: nflowpack
742  integer(I4B) :: i, ip
743  integer(I4B) :: naux
744  logical :: found_flowja
745  logical :: found_dataspdis
746  logical :: found_datasat
747  logical :: found_stoss
748  logical :: found_stosy
749  integer(I4B), dimension(:), allocatable :: imap
750  !
751  ! -- Calculate the number of gwf flow packages
752  allocate (imap(this%bfr%nbudterms))
753  imap(:) = 0
754  nflowpack = 0
755  found_flowja = .false.
756  found_dataspdis = .false.
757  found_datasat = .false.
758  found_stoss = .false.
759  found_stosy = .false.
760  do i = 1, this%bfr%nbudterms
761  select case (trim(adjustl(this%bfr%budtxtarray(i))))
762  case ('FLOW-JA-FACE')
763  found_flowja = .true.
764  case ('DATA-SPDIS')
765  found_dataspdis = .true.
766  case ('DATA-SAT')
767  found_datasat = .true.
768  case ('STO-SS')
769  found_stoss = .true.
770  this%igwfstrgss = 1
771  case ('STO-SY')
772  found_stosy = .true.
773  this%igwfstrgsy = 1
774  case default
775  nflowpack = nflowpack + 1
776  imap(i) = 1
777  end select
778  end do
779  !
780  ! -- allocate gwfpackage arrays (gwfpackages, iatp, datp, ...)
781  call this%allocate_gwfpackages(nflowpack)
782  !
783  ! -- Copy the package name and aux names from budget file reader
784  ! to the gwfpackages derived-type variable
785  ip = 1
786  do i = 1, this%bfr%nbudterms
787  if (imap(i) == 0) cycle
788  call this%gwfpackages(ip)%set_name(this%bfr%dstpackagenamearray(i), &
789  this%bfr%budtxtarray(i))
790  naux = this%bfr%nauxarray(i)
791  call this%gwfpackages(ip)%set_auxname(naux, &
792  this%bfr%auxtxtarray(1:naux, i))
793  ip = ip + 1
794  end do
795  !
796  ! -- Copy just the package names for the boundary packages into
797  ! the flowpacknamearray
798  ip = 1
799  do i = 1, size(imap)
800  if (imap(i) == 1) then
801  this%flowpacknamearray(ip) = this%bfr%dstpackagenamearray(i)
802  ip = ip + 1
803  end if
804  end do
805  !
806  ! -- Error if specific discharge, saturation or flowja not found
807  if (.not. found_dataspdis) then
808  write (errmsg, '(a)') 'Specific discharge not found in &
809  &budget file. SAVE_SPECIFIC_DISCHARGE and &
810  &SAVE_FLOWS must be activated in the NPF package.'
811  call store_error(errmsg)
812  end if
813  if (.not. found_datasat) then
814  write (errmsg, '(a)') 'Saturation not found in &
815  &budget file. SAVE_SATURATION and &
816  &SAVE_FLOWS must be activated in the NPF package.'
817  call store_error(errmsg)
818  end if
819  if (.not. found_flowja) then
820  write (errmsg, '(a)') 'FLOWJA not found in &
821  &budget file. SAVE_FLOWS must &
822  &be activated in the NPF package.'
823  call store_error(errmsg)
824  end if
825  if (count_errors() > 0) then
826  call this%parser%StoreErrorUnit()
827  end if
828  end subroutine initialize_gwfterms_from_bfr
829 
830  !> @brief Initialize groundwater flow terms from the groundwater budget
831  !!
832  !! Flows are coming from a gwf-gwt exchange object
833  !<
835  ! -- modules
836  use bndmodule, only: bndtype, getbndfromlist
837  ! -- dummy
838  class(tspfmitype) :: this
839  ! -- local
840  integer(I4B) :: ngwfpack
841  integer(I4B) :: ngwfterms
842  integer(I4B) :: ip
843  integer(I4B) :: imover
844  integer(I4B) :: ntomvr
845  integer(I4B) :: iterm
846  character(len=LENPACKAGENAME) :: budtxt
847  class(bndtype), pointer :: packobj => null()
848  !
849  ! -- determine size of gwf terms
850  ngwfpack = this%gwfbndlist%Count()
851  !
852  ! -- Count number of to-mvr terms, but do not include advanced packages
853  ! as those mover terms are not losses from the cell, but rather flows
854  ! within the advanced package
855  ntomvr = 0
856  do ip = 1, ngwfpack
857  packobj => getbndfromlist(this%gwfbndlist, ip)
858  imover = packobj%imover
859  if (packobj%isadvpak /= 0) imover = 0
860  if (imover /= 0) then
861  ntomvr = ntomvr + 1
862  end if
863  end do
864  !
865  ! -- Allocate arrays in fmi of size ngwfterms, which is the number of
866  ! packages plus the number of packages with mover terms.
867  ngwfterms = ngwfpack + ntomvr
868  call this%allocate_gwfpackages(ngwfterms)
869  !
870  ! -- Assign values in the fmi package
871  iterm = 1
872  do ip = 1, ngwfpack
873  !
874  ! -- set and store names
875  packobj => getbndfromlist(this%gwfbndlist, ip)
876  budtxt = adjustl(packobj%text)
877  call this%gwfpackages(iterm)%set_name(packobj%packName, budtxt)
878  this%flowpacknamearray(iterm) = packobj%packName
879  call this%gwfpackages(iterm)%set_auxname(packobj%naux, &
880  packobj%auxname)
881  iterm = iterm + 1
882  !
883  ! -- if this package has a mover associated with it, then add another
884  ! term that corresponds to the mover flows
885  imover = packobj%imover
886  if (packobj%isadvpak /= 0) imover = 0
887  if (imover /= 0) then
888  budtxt = trim(adjustl(packobj%text))//'-TO-MVR'
889  call this%gwfpackages(iterm)%set_name(packobj%packName, budtxt)
890  this%flowpacknamearray(iterm) = packobj%packName
891  call this%gwfpackages(iterm)%set_auxname(packobj%naux, &
892  packobj%auxname)
893  this%igwfmvrterm(iterm) = 1
894  iterm = iterm + 1
895  end if
896  end do
898 
899  !> @brief Initialize an array for storing PackageBudget objects.
900  !!
901  !! This routine allocates gwfpackages (an array of PackageBudget
902  !! objects) to the proper size and initializes member variables.
903  !<
904  subroutine gwtfmi_allocate_gwfpackages(this, ngwfterms)
905  ! -- modules
906  use constantsmodule, only: lenmempath
908  ! -- dummy
909  class(tspfmitype) :: this
910  integer(I4B), intent(in) :: ngwfterms
911  ! -- local
912  integer(I4B) :: n
913  character(len=LENMEMPATH) :: memPath
914  !
915  ! -- direct allocate
916  allocate (this%gwfpackages(ngwfterms))
917  allocate (this%flowpacknamearray(ngwfterms))
918  allocate (this%datp(ngwfterms))
919  !
920  ! -- mem_allocate
921  call mem_allocate(this%iatp, ngwfterms, 'IATP', this%memoryPath)
922  call mem_allocate(this%igwfmvrterm, ngwfterms, 'IGWFMVRTERM', this%memoryPath)
923  !
924  ! -- initialize
925  this%nflowpack = ngwfterms
926  do n = 1, this%nflowpack
927  this%iatp(n) = 0
928  this%igwfmvrterm(n) = 0
929  this%flowpacknamearray(n) = ''
930  !
931  ! -- Create a mempath for each individual flow package data set
932  ! of the form, MODELNAME/FMI-FTn
933  write (mempath, '(a, i0)') trim(this%memoryPath)//'-FT', n
934  call this%gwfpackages(n)%initialize(mempath)
935  end do
936  end subroutine gwtfmi_allocate_gwfpackages
937 
938  !> @brief Deallocate memory
939  !!
940  !! Deallocate memory that stores the gwfpackages array
941  !<
943  ! -- modules
944  ! -- dummy
945  class(tspfmitype) :: this
946  ! -- local
947  integer(I4B) :: n
948  !
949  ! -- initialize
950  do n = 1, this%nflowpack
951  call this%gwfpackages(n)%da()
952  end do
953  end subroutine gwtfmi_deallocate_gwfpackages
954 
955 end module tspfmimodule
This module contains the base boundary package.
class(bndtype) function, pointer, public getbndfromlist(list, idx)
Get boundary from package list.
This module contains the BudgetModule.
Definition: Budget.f90:20
subroutine, public rate_accumulator(flow, rin, rout)
@ brief Rate accumulator subroutine
Definition: Budget.f90:632
subroutine, public budgetobject_cr_bfr(this, name, ibinun, iout, colconv1, colconv2)
Create a new budget object from a binary flow file.
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter linelength
maximum length of a standard line
Definition: Constants.f90:45
real(dp), parameter dhdry
real dry cell constant
Definition: Constants.f90:94
integer(i4b), parameter lenpackagename
maximum length of the package name
Definition: Constants.f90:23
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dem6
real constant 1e-6
Definition: Constants.f90:109
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter lenbudtxt
maximum length of a budget component names
Definition: Constants.f90:37
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
subroutine, public urdaux(naux, inunit, iout, lloc, istart, istop, auxname, line, text)
Read auxiliary variables from an input line.
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
This module defines variable data types.
Definition: kind.f90:8
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
This module contains the PackageBudgetModule Module.
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
Definition: Sim.f90:168
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29
real(dp) function gwfsatold(this, n, delt)
Calculate the previous saturation level.
Definition: tsp-fmi.f90:504
integer(i4b), parameter nbditems
Definition: tsp-fmi.f90:24
subroutine fmi_bd(this, isuppress_output, model_budget)
Calculate budget terms associated with FMI object.
Definition: tsp-fmi.f90:252
subroutine gwtfmi_deallocate_gwfpackages(this)
Deallocate memory.
Definition: tsp-fmi.f90:943
character(len=lenbudtxt), dimension(nbditems) budtxt
Definition: tsp-fmi.f90:25
subroutine gwtfmi_allocate_scalars(this)
@ brief Allocate scalars
Definition: tsp-fmi.f90:366
subroutine gwtfmi_allocate_arrays(this, nodes)
@ brief Allocate arrays for FMI object
Definition: tsp-fmi.f90:391
subroutine initialize_gwfterms_from_gwfbndlist(this)
Initialize groundwater flow terms from the groundwater budget.
Definition: tsp-fmi.f90:835
subroutine gwtfmi_allocate_gwfpackages(this, ngwfterms)
Initialize an array for storing PackageBudget objects.
Definition: tsp-fmi.f90:905
subroutine fmi_fc(this, nodes, cold, nja, matrix_sln, idxglo, rhs)
Calculate coefficients and fill matrix and rhs terms associated with FMI object.
Definition: tsp-fmi.f90:190
subroutine initialize_gwfterms_from_bfr(this)
Initialize the groundwater flow terms based on the budget file reader.
Definition: tsp-fmi.f90:735
subroutine fmi_rp(this, inmvr)
Read and prepare.
Definition: tsp-fmi.f90:110
subroutine gwtfmi_read_options(this)
Read options from input file.
Definition: tsp-fmi.f90:528
subroutine set_aptbudobj_pointer(this, name, budobjptr)
Set the pointer to a budget object.
Definition: tsp-fmi.f90:711
subroutine set_active_status(this, cnew)
Set gwt transport cell status.
Definition: tsp-fmi.f90:421
subroutine fmi_ot_flow(this, icbcfl, icbcun)
Save budget terms associated with FMI object.
Definition: tsp-fmi.f90:273
character(len=lenpackagename) text
Definition: tsp-fmi.f90:22
subroutine fmi_ad(this, cnew)
Advance routine for FMI object.
Definition: tsp-fmi.f90:139
subroutine fmi_cq(this, cnew, flowja)
Calculate flow correction.
Definition: tsp-fmi.f90:223
subroutine gwtfmi_da(this)
Deallocate variables.
Definition: tsp-fmi.f90:313
subroutine gwtfmi_read_packagedata(this)
Read PACKAGEDATA block.
Definition: tsp-fmi.f90:578
subroutine, public fmi_cr(fmiobj, name_model, inunit, iout, eqnsclfac, depvartype)
Create a new FMI object.
Definition: tsp-fmi.f90:75
@ brief BndType
Derived type for the Budget object.
Definition: Budget.f90:39
A generic heterogeneous doubly-linked list.
Definition: List.f90:14
Derived type for storing flows.