MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
prt-fmi.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b, lgp
4  use errorutilmodule, only: pstop
6  use simmodule, only: store_error
7  use simvariablesmodule, only: errmsg
9  use basedismodule, only: disbasetype
11 
12  implicit none
13  private
14  public :: prtfmitype
15  public :: fmi_cr
16 
17  character(len=LENPACKAGENAME) :: text = ' PRTFMI'
18 
20  private
21  integer(I4B), public :: max_faces !< maximum number of faces for grid cell polygons
22  real(dp), allocatable, public :: sourceflows(:) ! cell source flows array
23  real(dp), allocatable, public :: sinkflows(:) ! cell sink flows array
24  real(dp), allocatable, public :: storageflows(:) ! cell storage flows array
25  real(dp), allocatable, public :: boundaryflows(:) ! cell boundary flows array
26  integer(I4B), allocatable, public :: boundaryfaces(:) ! bitmask of assigned boundary faces
27 
28  contains
29 
30  procedure :: fmi_ad
31  procedure :: fmi_df => prtfmi_df
32  procedure, private :: accumulate_flows
33  procedure :: mark_boundary_face
34  procedure :: is_boundary_face
36  procedure, private :: iflowface_to_icellface
37 
38  end type prtfmitype
39 
40 contains
41 
42  !> @brief Create a new PrtFmi object
43  subroutine fmi_cr(fmiobj, name_model, input_mempath, inunit, iout)
44  ! dummy
45  type(prtfmitype), pointer :: fmiobj
46  character(len=*), intent(in) :: name_model
47  character(len=*), intent(in) :: input_mempath
48  integer(I4B), intent(inout) :: inunit
49  integer(I4B), intent(in) :: iout
50 
51  ! Create the object
52  allocate (fmiobj)
53 
54  ! create name and memory path
55  call fmiobj%set_names(1, name_model, 'FMI', 'FMI', input_mempath)
56  fmiobj%text = text
57 
58  ! Allocate scalars
59  call fmiobj%allocate_scalars()
60 
61  ! Set variables
62  fmiobj%inunit = inunit
63  fmiobj%iout = iout
64 
65  ! Assign dependent variable label
66  fmiobj%depvartype = 'TRACKS '
67 
68  end subroutine fmi_cr
69 
70  !> @brief Time step advance
71  subroutine fmi_ad(this)
72  ! modules
73  use constantsmodule, only: dhdry
74  ! dummy
75  class(prtfmitype) :: this
76  ! local
77  integer(I4B) :: n
78  character(len=15) :: nodestr
79  character(len=*), parameter :: fmtdry = &
80  &"(/1X,'WARNING: DRY CELL ENCOUNTERED AT ',a,'; RESET AS INACTIVE')"
81  character(len=*), parameter :: fmtrewet = &
82  &"(/1X,'DRY CELL REACTIVATED AT ', a)"
83 
84  ! Set flag to indicated that flows are being updated. For the case where
85  ! flows may be reused (only when flows are read from a file) then set
86  ! the flag to zero to indicated that flows were not updated
87  this%iflowsupdated = 1
88 
89  ! If reading flows from a budget file, read the next set of records
90  if (this%iubud /= 0) call this%advance_bfr()
91 
92  ! If reading heads from a head file, read the next set of records
93  if (this%iuhds /= 0) call this%advance_hfr()
94 
95  ! If mover flows are being read from file, read the next set of records
96  if (this%iumvr /= 0) &
97  call this%mvrbudobj%bfr_advance(this%dis, this%iout)
98 
99  ! Accumulate flows
100  call this%accumulate_flows()
101 
102  ! if flow cell is dry, then set this%ibound = 0
103  do n = 1, this%dis%nodes
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  class(prtfmitype) :: this
138  class(disbasetype), pointer, intent(in) :: dis
139  integer(I4B), intent(in) :: idryinactive
140 
141  call this%FlowModelInterfaceType%fmi_df(dis, idryinactive)
142 
143  this%max_faces = this%dis%get_max_npolyverts() + 2
144  allocate (this%StorageFlows(this%dis%nodes))
145  allocate (this%SourceFlows(this%dis%nodes))
146  allocate (this%SinkFlows(this%dis%nodes))
147  allocate (this%BoundaryFlows(this%dis%nodes * this%max_faces))
148  allocate (this%BoundaryFaces(this%dis%nodes))
149 
150  end subroutine prtfmi_df
151 
152  !> @brief Accumulate flows
153  subroutine accumulate_flows(this)
154  ! dummy
155  class(prtfmitype) :: this
156  ! local
157  integer(I4B) :: j, i, ip, ib
158  integer(I4B) :: ioffset, iflowface, iauxiflowface, icellface
159  real(DP) :: qbnd
160  character(len=LENAUXNAME) :: auxname
161  integer(I4B) :: naux
162 
163  this%StorageFlows = dzero
164  if (this%igwfstrgss /= 0) &
165  this%StorageFlows = this%StorageFlows + this%gwfstrgss
166  if (this%igwfstrgsy /= 0) &
167  this%StorageFlows = this%StorageFlows + this%gwfstrgsy
168 
169  this%SourceFlows = dzero
170  this%SinkFlows = dzero
171  this%BoundaryFlows = dzero
172  this%BoundaryFaces = 0
173  do ip = 1, this%nflowpack
174  iauxiflowface = 0
175  naux = this%gwfpackages(ip)%naux
176  if (naux > 0) then
177  do j = 1, naux
178  auxname = this%gwfpackages(ip)%auxname(j)
179  if (trim(adjustl(auxname)) == "IFLOWFACE") then
180  iauxiflowface = j
181  exit
182  end if
183  end do
184  end if
185  do ib = 1, this%gwfpackages(ip)%nbound
186  i = this%gwfpackages(ip)%nodelist(ib)
187  if (i <= 0) cycle
188  if (this%ibound(i) <= 0) cycle
189  qbnd = this%gwfpackages(ip)%get_flow(ib)
190  ! todo, after initial release: default iflowface values for different packages
191  iflowface = 0
192  icellface = 0
193  if (iauxiflowface > 0) then
194  iflowface = nint(this%gwfpackages(ip)%auxvar(iauxiflowface, ib))
195  icellface = this%iflowface_to_icellface(iflowface)
196  end if
197  if (icellface > 0) then
198  call this%mark_boundary_face(i, icellface)
199  ioffset = (i - 1) * this%max_faces
200  this%BoundaryFlows(ioffset + icellface) = &
201  this%BoundaryFlows(ioffset + icellface) + qbnd
202  else if (qbnd .gt. dzero) then
203  this%SourceFlows(i) = this%SourceFlows(i) + qbnd
204  else if (qbnd .lt. dzero) then
205  this%SinkFlows(i) = this%SinkFlows(i) + qbnd
206  end if
207  end do
208  end do
209 
210  end subroutine accumulate_flows
211 
212  !> @brief Mark a face as a boundary face.
213  subroutine mark_boundary_face(this, ic, iface)
214  class(prtfmitype) :: this
215  integer(I4B), intent(in) :: ic !< node number (reduced)
216  integer(I4B), intent(in) :: iface !< cell face number
217  ! local
218  integer(I4B) :: bit_pos
219 
220  if (ic <= 0 .or. ic > this%dis%nodes) then
221  print *, 'Invalid cell number: ', ic
222  print *, 'Expected a value in range [1, ', this%dis%nodes, ']'
223  call pstop(1)
224  end if
225  if (iface <= 0) then
226  print *, 'Invalid face number: ', iface
227  print *, 'Expected a value in range [1, ', this%max_faces, ']'
228  call pstop(1)
229  end if
230  bit_pos = iface - 1 ! bit position 0-based
231  if (bit_pos < 0 .or. bit_pos > 31) then
232  print *, 'Invalid bitmask position: ', iface
233  print *, 'Expected a value in range [0, 31]'
234  call pstop(1)
235  end if
236  this%BoundaryFaces(ic) = ibset(this%BoundaryFaces(ic), bit_pos)
237  end subroutine mark_boundary_face
238 
239  !> @brief Check if a face is assigned to a boundary package.
240  function is_boundary_face(this, ic, iface) result(is_boundary)
241  class(prtfmitype) :: this
242  integer(I4B), intent(in) :: ic !< node number (reduced)
243  integer(I4B), intent(in) :: iface !< cell face number
244  logical(LGP) :: is_boundary
245  ! local
246  integer(I4B) :: bit_pos
247 
248  is_boundary = .false.
249  if (ic <= 0 .or. ic > this%dis%nodes) then
250  print *, 'Invalid cell number: ', ic
251  print *, 'Expected a value in range [1, ', this%dis%nodes, ']'
252  call pstop(1)
253  end if
254  if (iface <= 0) then
255  print *, 'Invalid face number: ', iface
256  print *, 'Expected a value in range [1, ', this%max_faces, ']'
257  call pstop(1)
258  end if
259  bit_pos = iface - 1 ! bit position 0-based
260  if (bit_pos < 0 .or. bit_pos > 31) then
261  print *, 'Invalid bitmask position: ', iface
262  print *, 'Expected a value in range [0, 31]'
263  call pstop(1)
264  end if
265  is_boundary = btest(this%BoundaryFaces(ic), bit_pos)
266  end function is_boundary_face
267 
268  !> @brief Check if a face is an assigned boundary with net outflow.
269  function is_net_out_boundary_face(this, ic, iface) result(is_net_out_boundary)
270  class(prtfmitype) :: this
271  integer(I4B), intent(in) :: ic !< node number (reduced)
272  integer(I4B), intent(in) :: iface !< cell face number
273  logical(LGP) :: is_net_out_boundary
274  ! local
275  integer(I4B) :: ioffset
276 
277  is_net_out_boundary = .false.
278  if (.not. this%is_boundary_face(ic, iface)) return
279  ioffset = (ic - 1) * this%max_faces
280  if (this%BoundaryFlows(ioffset + iface) < dzero) is_net_out_boundary = .true.
281  end function is_net_out_boundary_face
282 
283  !> @brief Convert an iflowface number to a cell face number.
284  !! Maps bottom (-2) -> max_faces - 1, top (-1) -> max_faces.
285  function iflowface_to_icellface(this, iflowface) result(iface)
286  class(prtfmitype), intent(inout) :: this
287  integer(I4B), intent(in) :: iflowface
288  integer(I4B) :: iface
289 
290  iface = iflowface
291  if (iface < 0) iface = iface + this%max_faces + 1
292  end function iflowface_to_icellface
293 
294 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 lenvarname
maximum length of a variable name
Definition: Constants.f90:17
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
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
Definition: ErrorUtil.f90:24
This module defines variable data types.
Definition: kind.f90:8
subroutine accumulate_flows(this)
Accumulate flows.
Definition: prt-fmi.f90:154
integer(i4b) function iflowface_to_icellface(this, iflowface)
Convert an iflowface number to a cell face number. Maps bottom (-2) -> max_faces - 1,...
Definition: prt-fmi.f90:286
subroutine fmi_ad(this)
Time step advance.
Definition: prt-fmi.f90:72
subroutine mark_boundary_face(this, ic, iface)
Mark a face as a boundary face.
Definition: prt-fmi.f90:214
subroutine prtfmi_df(this, dis, idryinactive)
Define the flow model interface.
Definition: prt-fmi.f90:137
logical(lgp) function is_net_out_boundary_face(this, ic, iface)
Check if a face is an assigned boundary with net outflow.
Definition: prt-fmi.f90:270
character(len=lenpackagename) text
Definition: prt-fmi.f90:17
subroutine, public fmi_cr(fmiobj, name_model, input_mempath, inunit, iout)
Create a new PrtFmi object.
Definition: prt-fmi.f90:44
logical(lgp) function is_boundary_face(this, ic, iface)
Check if a face is assigned to a boundary package.
Definition: prt-fmi.f90:241
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