MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
prt-fmi.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b
5  use simmodule, only: store_error
6  use simvariablesmodule, only: errmsg
8  use basedismodule, only: disbasetype
10 
11  implicit none
12  private
13  public :: prtfmitype
14  public :: fmi_cr
15 
16  character(len=LENPACKAGENAME) :: text = ' PRTFMI'
17 
19 
20  integer(I4B) :: max_faces !< maximum number of faces for grid cell polygons
21  double precision, allocatable, public :: sourceflows(:) ! cell source flows array
22  double precision, allocatable, public :: sinkflows(:) ! cell sink flows array
23  double precision, allocatable, public :: storageflows(:) ! cell storage flows array
24  double precision, allocatable, public :: boundaryflows(:) ! cell boundary flows array
25 
26  contains
27 
28  procedure :: fmi_ad
29  procedure :: fmi_df => prtfmi_df
30  procedure, private :: accumulate_flows
31 
32  end type prtfmitype
33 
34 contains
35 
36  !> @brief Create a new PrtFmi object
37  subroutine fmi_cr(fmiobj, name_model, input_mempath, inunit, iout)
38  ! dummy
39  type(prtfmitype), pointer :: fmiobj
40  character(len=*), intent(in) :: name_model
41  character(len=*), intent(in) :: input_mempath
42  integer(I4B), intent(inout) :: inunit
43  integer(I4B), intent(in) :: iout
44  !
45  ! Create the object
46  allocate (fmiobj)
47  !
48  ! create name and memory path
49  call fmiobj%set_names(1, name_model, 'FMI', 'FMI', input_mempath)
50  fmiobj%text = text
51  !
52  ! Allocate scalars
53  call fmiobj%allocate_scalars()
54  !
55  ! Set variables
56  fmiobj%inunit = inunit
57  fmiobj%iout = iout
58  !
59  ! Assign dependent variable label
60  fmiobj%depvartype = 'TRACKS '
61 
62  end subroutine fmi_cr
63 
64  !> @brief Time step advance
65  subroutine fmi_ad(this)
66  ! modules
67  use constantsmodule, only: dhdry
68  ! dummy
69  class(prtfmitype) :: this
70  ! local
71  integer(I4B) :: n
72  character(len=15) :: nodestr
73  character(len=*), parameter :: fmtdry = &
74  &"(/1X,'WARNING: DRY CELL ENCOUNTERED AT ',a,'; RESET AS INACTIVE')"
75  character(len=*), parameter :: fmtrewet = &
76  &"(/1X,'DRY CELL REACTIVATED AT ', a)"
77  !
78  ! Set flag to indicated that flows are being updated. For the case where
79  ! flows may be reused (only when flows are read from a file) then set
80  ! the flag to zero to indicated that flows were not updated
81  this%iflowsupdated = 1
82  !
83  ! If reading flows from a budget file, read the next set of records
84  if (this%iubud /= 0) then
85  call this%advance_bfr()
86  end if
87  !
88  ! If reading heads from a head file, read the next set of records
89  if (this%iuhds /= 0) then
90  call this%advance_hfr()
91  end if
92  !
93  ! If mover flows are being read from file, read the next set of records
94  if (this%iumvr /= 0) then
95  call this%mvrbudobj%bfr_advance(this%dis, this%iout)
96  end if
97  !
98  ! Accumulate flows
99  call this%accumulate_flows()
100  !
101  ! if flow cell is dry, then set this%ibound = 0
102  do n = 1, this%dis%nodes
103  !
104  ! Calculate the ibound-like array that has 0 if saturation
105  ! is zero and 1 otherwise
106  if (this%gwfsat(n) > dzero) then
107  this%ibdgwfsat0(n) = 1
108  else
109  this%ibdgwfsat0(n) = 0
110  end if
111  !
112  ! Check if active model cell is inactive for flow
113  if (this%ibound(n) > 0) then
114  if (this%gwfhead(n) == dhdry) then
115  ! cell should be made inactive
116  this%ibound(n) = 0
117  call this%dis%noder_to_string(n, nodestr)
118  write (this%iout, fmtdry) trim(nodestr)
119  end if
120  end if
121  !
122  ! Convert dry model cell to active if flow has rewet
123  if (this%ibound(n) == 0) then
124  if (this%gwfhead(n) /= dhdry) then
125  ! cell is now wet
126  this%ibound(n) = 1
127  call this%dis%noder_to_string(n, nodestr)
128  write (this%iout, fmtrewet) trim(nodestr)
129  end if
130  end if
131  end do
132 
133  end subroutine fmi_ad
134 
135  !> @brief Define the flow model interface
136  subroutine prtfmi_df(this, dis, idryinactive)
137  ! modules
138  use simmodule, only: store_error
139  ! dummy
140  class(prtfmitype) :: this
141  class(disbasetype), pointer, intent(in) :: dis
142  integer(I4B), intent(in) :: idryinactive
143  !
144  ! Call parent class define
145  call this%FlowModelInterfaceType%fmi_df(dis, idryinactive)
146  !
147  ! Allocate arrays
148  this%max_faces = this%dis%get_max_npolyverts() + 2
149  allocate (this%StorageFlows(this%dis%nodes))
150  allocate (this%SourceFlows(this%dis%nodes))
151  allocate (this%SinkFlows(this%dis%nodes))
152  allocate (this%BoundaryFlows(this%dis%nodes * this%max_faces))
153 
154  end subroutine prtfmi_df
155 
156  !> @brief Accumulate flows
157  subroutine accumulate_flows(this)
158  implicit none
159  ! dummy
160  class(prtfmitype) :: this
161  ! local
162  integer(I4B) :: j, i, ip, ib
163  integer(I4B) :: ioffset, iflowface, iauxiflowface
164  real(DP) :: qbnd
165  character(len=LENAUXNAME) :: auxname
166  integer(I4B) :: naux
167 
168  this%StorageFlows = dzero
169  if (this%igwfstrgss /= 0) &
170  this%StorageFlows = this%StorageFlows + &
171  this%gwfstrgss
172  if (this%igwfstrgsy /= 0) &
173  this%StorageFlows = this%StorageFlows + &
174  this%gwfstrgsy
175 
176  this%SourceFlows = dzero
177  this%SinkFlows = dzero
178  this%BoundaryFlows = dzero
179  do ip = 1, this%nflowpack
180  iauxiflowface = 0
181  naux = this%gwfpackages(ip)%naux
182  if (naux > 0) then
183  do j = 1, naux
184  auxname = this%gwfpackages(ip)%auxname(j)
185  if (trim(adjustl(auxname)) == "IFLOWFACE") then
186  iauxiflowface = j
187  exit
188  end if
189  end do
190  end if
191  do ib = 1, this%gwfpackages(ip)%nbound
192  i = this%gwfpackages(ip)%nodelist(ib)
193  if (i <= 0) cycle
194  if (this%ibound(i) <= 0) cycle
195  qbnd = this%gwfpackages(ip)%get_flow(ib)
196  ! todo, after initial release: default iflowface values for different packages
197  iflowface = 0
198  if (iauxiflowface > 0) then
199  iflowface = nint(this%gwfpackages(ip)%auxvar(iauxiflowface, ib))
200  ! maps bot -2 -> max_faces - 1, top -1 -> max_faces
201  if (iflowface < 0) iflowface = iflowface + this%max_faces + 1
202  end if
203  if (iflowface .gt. 0) then
204  ioffset = (i - 1) * this%max_faces
205  this%BoundaryFlows(ioffset + iflowface) = &
206  this%BoundaryFlows(ioffset + iflowface) + qbnd
207  else if (qbnd .gt. dzero) then
208  this%SourceFlows(i) = this%SourceFlows(i) + qbnd
209  else if (qbnd .lt. dzero) then
210  this%SinkFlows(i) = this%SinkFlows(i) + qbnd
211  end if
212  end do
213  end do
214 
215  end subroutine accumulate_flows
216 
217 end module prtfmimodule
This module contains simulation constants.
Definition: Constants.f90:9
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 lenauxname
maximum length of a aux variable
Definition: Constants.f90:35
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
This module defines variable data types.
Definition: kind.f90:8
subroutine accumulate_flows(this)
Accumulate flows.
Definition: prt-fmi.f90:158
subroutine fmi_ad(this)
Time step advance.
Definition: prt-fmi.f90:66
subroutine prtfmi_df(this, dis, idryinactive)
Define the flow model interface.
Definition: prt-fmi.f90:137
character(len=lenpackagename) text
Definition: prt-fmi.f90:16
subroutine, public fmi_cr(fmiobj, name_model, input_mempath, inunit, iout)
Create a new PrtFmi object.
Definition: prt-fmi.f90:38
This module contains simulation methods.
Definition: Sim.f90:10
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
This module contains simulation variables.
Definition: SimVariables.f90:9
character(len=maxcharlen) errmsg
error message string