MODFLOW 6  version 6.7.0.dev0
USGS Modular Hydrologic Model
All Classes Namespaces Files Functions Variables Enumerator Pages
mf6bmiGrid.f90
Go to the documentation of this file.
1 !> @brief This module contains BMI routines to expose the MODFLOW 6 discretization
2 !!
3 !! NB: this module is experimental and still under development:
4 !! - add error handling
5 !! - add var address checks
6 !! - ...
7 !<
8 module mf6bmigrid
9  use mf6bmiutil
10  use mf6bmierror
11  use iso_c_binding, only: c_double, c_ptr, c_loc
13  use kindmodule, only: dp, i4b
16  implicit none
17 
18 contains
19 
20  ! Get the grid identifier for the given variable.
21  function get_var_grid(c_var_address, var_grid) result(bmi_status) &
22  bind(C, name="get_var_grid")
23  !DIR$ ATTRIBUTES DLLEXPORT :: get_var_grid
24  ! -- modules
25  use listsmodule, only: basemodellist
27  ! -- dummy variables
28  character(kind=c_char), intent(in) :: c_var_address(*)
29  integer(kind=c_int), intent(out) :: var_grid
30  integer(kind=c_int) :: bmi_status
31  ! -- local variables
32  character(len=LENMODELNAME) :: model_name
33  character(len=LENMEMPATH) :: var_address
34  integer(I4B) :: i
35  logical(LGP) :: success
36  class(basemodeltype), pointer :: basemodel
37 
38  var_grid = -1
39 
40  bmi_status = bmi_failure
41  var_address = char_array_to_string(c_var_address, &
42  strlen(c_var_address, lenmemaddress + 1))
43  model_name = extract_model_name(var_address, success)
44  if (.not. success) then
45  ! we failed
46  return
47  end if
48 
49  do i = 1, basemodellist%Count()
50  basemodel => getbasemodelfromlist(basemodellist, i)
51  if (basemodel%name == model_name) then
52  var_grid = basemodel%id
53  bmi_status = bmi_success
54  return
55  end if
56  end do
57  end function get_var_grid
58 
59  ! Get the grid type as a string.
60  function get_grid_type(grid_id, grid_type) result(bmi_status) &
61  bind(C, name="get_grid_type")
62  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_type
63  ! -- dummy variables
64  integer(kind=c_int), intent(in) :: grid_id
65  character(kind=c_char), intent(inout) :: grid_type(bmi_lengridtype)
66  integer(kind=c_int) :: bmi_status
67  ! -- local variables
68  character(len=LENGRIDTYPE) :: grid_type_f
69  character(len=LENMODELNAME) :: model_name
70 
71  bmi_status = bmi_failure
72  model_name = get_model_name(grid_id)
73  if (model_name == '') return
74 
75  call get_grid_type_model(model_name, grid_type_f)
76 
77  if (grid_type_f == "DIS") then
78  grid_type_f = "rectilinear"
79  else if ((grid_type_f == "DISV") .or. (grid_type_f == "DISU")) then
80  grid_type_f = "unstructured"
81  else
82  return
83  end if
84  grid_type = string_to_char_array(trim(grid_type_f), len_trim(grid_type_f))
85  bmi_status = bmi_success
86  end function get_grid_type
87 
88  ! Get number of dimensions of the computational grid.
89  function get_grid_rank(grid_id, grid_rank) result(bmi_status) &
90  bind(C, name="get_grid_rank")
91  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_rank
92  ! -- dummy variables
93  integer(kind=c_int), intent(in) :: grid_id
94  integer(kind=c_int), intent(out) :: grid_rank
95  integer(kind=c_int) :: bmi_status
96  ! -- local variables
97  character(len=LENMODELNAME) :: model_name
98  integer(I4B), dimension(:), pointer, contiguous :: grid_shape
99 
100  bmi_status = bmi_failure
101  ! TODO: It is currently only implemented for DIS grids
102  if (.not. confirm_grid_type(grid_id, "DIS")) return
103 
104  ! get shape array
105  model_name = get_model_name(grid_id)
106  call mem_setptr(grid_shape, "MSHAPE", create_mem_path(model_name, 'DIS'))
107 
108  if (grid_shape(1) == 1) then
109  grid_rank = 2
110  else
111  grid_rank = 3
112  end if
113  bmi_status = bmi_success
114  end function get_grid_rank
115 
116  ! Get the total number of elements in the computational grid.
117  function get_grid_size(grid_id, grid_size) result(bmi_status) &
118  bind(C, name="get_grid_size")
119  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_size
120  ! -- dummy variables
121  integer(kind=c_int), intent(in) :: grid_id
122  integer(kind=c_int), intent(out) :: grid_size
123  integer(kind=c_int) :: bmi_status
124  ! -- local variables
125  character(len=LENMODELNAME) :: model_name
126  integer(I4B), dimension(:), pointer, contiguous :: grid_shape
127  character(kind=c_char) :: grid_type(bmi_lengridtype)
128  character(len=LENGRIDTYPE) :: grid_type_f
129  integer(I4B) :: status
130 
131  bmi_status = bmi_failure
132 
133  if (get_grid_type(grid_id, grid_type) /= bmi_success) return
134  grid_type_f = char_array_to_string(grid_type, &
135  strlen(grid_type, lengridtype + 1))
136  model_name = get_model_name(grid_id)
137 
138  if (grid_type_f == "rectilinear") then
139  call mem_setptr(grid_shape, "MSHAPE", create_mem_path(model_name, 'DIS'))
140  grid_size = grid_shape(1) * grid_shape(2) * grid_shape(3)
141  bmi_status = bmi_success
142  return
143  else if (grid_type_f == "unstructured") then
144  status = get_grid_node_count(grid_id, grid_size)
145  bmi_status = bmi_success
146  return
147  end if
148  end function get_grid_size
149 
150  ! Get the dimensions of the computational grid.
151  function get_grid_shape(grid_id, grid_shape) result(bmi_status) &
152  bind(C, name="get_grid_shape")
153  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_shape
154  ! -- dummy variables
155  integer(kind=c_int), intent(in) :: grid_id
156  integer(kind=c_int), intent(out) :: grid_shape(*)
157  integer(kind=c_int) :: bmi_status
158  ! -- local variables
159  integer, dimension(:), pointer, contiguous :: grid_shape_ptr
160  character(len=LENMODELNAME) :: model_name
161  character(kind=c_char) :: grid_type(bmi_lengridtype)
162 
163  bmi_status = bmi_failure
164  ! make sure function is only used for implemented grid_types
165  if (get_grid_type(grid_id, grid_type) /= bmi_success) return
166 
167  ! get shape array
168  model_name = get_model_name(grid_id)
169  call mem_setptr(grid_shape_ptr, "MSHAPE", create_mem_path(model_name, 'DIS'))
170 
171  if (grid_shape_ptr(1) == 1) then
172  grid_shape(1:2) = grid_shape_ptr(2:3) ! 2D
173  else
174  grid_shape(1:3) = grid_shape_ptr ! 3D
175  end if
176  bmi_status = bmi_success
177  end function get_grid_shape
178 
179  ! Provides an array (whose length is the number of rows) that gives the x-coordinate for each row.
180  function get_grid_x(grid_id, grid_x) result(bmi_status) &
181  bind(C, name="get_grid_x")
182  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_x
183  ! -- dummy variables
184  integer(kind=c_int), intent(in) :: grid_id
185  real(kind=c_double), intent(out) :: grid_x(*)
186  integer(kind=c_int) :: bmi_status
187  ! -- local variables
188  integer(I4B) :: i
189  integer, dimension(:), pointer, contiguous :: grid_shape_ptr
190  character(len=LENMODELNAME) :: model_name
191  character(kind=c_char) :: grid_type(bmi_lengridtype)
192  real(dp), dimension(:, :), pointer, contiguous :: vertices_ptr
193  character(len=LENGRIDTYPE) :: grid_type_f
194  integer(I4B) :: x_size
195 
196  bmi_status = bmi_failure
197  ! make sure function is only used for implemented grid_types
198  if (get_grid_type(grid_id, grid_type) /= bmi_success) return
199  grid_type_f = char_array_to_string(grid_type, &
200  strlen(grid_type, lengridtype + 1))
201 
202  model_name = get_model_name(grid_id)
203  if (grid_type_f == "rectilinear") then
204  call mem_setptr(grid_shape_ptr, "MSHAPE", &
205  create_mem_path(model_name, 'DIS'))
206  ! The dimension of x is in the last element of the shape array.
207  ! + 1 because we count corners, not centers.
208  x_size = grid_shape_ptr(size(grid_shape_ptr)) + 1
209  grid_x(1:x_size) = [(i, i=0, x_size - 1)]
210  else if (grid_type_f == "unstructured") then
211  call mem_setptr(vertices_ptr, "VERTICES", &
212  create_mem_path(model_name, 'DIS'))
213  ! x-coordinates are in the 1st column
214  x_size = size(vertices_ptr(1, :))
215  grid_x(1:x_size) = vertices_ptr(1, :)
216  else
217  bmi_status = bmi_failure
218  return
219  end if
220  bmi_status = bmi_success
221  end function get_grid_x
222 
223  ! Provides an array (whose length is the number of rows) that gives the y-coordinate for each row.
224  function get_grid_y(grid_id, grid_y) result(bmi_status) &
225  bind(C, name="get_grid_y")
226  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_y
227  ! -- dummy variables
228  integer(kind=c_int), intent(in) :: grid_id
229  real(kind=c_double), intent(out) :: grid_y(*)
230  integer(kind=c_int) :: bmi_status
231  ! -- local variables
232  integer(I4B) :: i
233  integer, dimension(:), pointer, contiguous :: grid_shape_ptr
234  character(len=LENMODELNAME) :: model_name
235  character(kind=c_char) :: grid_type(bmi_lengridtype)
236  real(dp), dimension(:, :), pointer, contiguous :: vertices_ptr
237  character(len=LENGRIDTYPE) :: grid_type_f
238  integer(I4B) :: y_size
239 
240  bmi_status = bmi_failure
241  if (get_grid_type(grid_id, grid_type) /= bmi_success) return
242  grid_type_f = char_array_to_string(grid_type, &
243  strlen(grid_type, lengridtype + 1))
244 
245  model_name = get_model_name(grid_id)
246  if (grid_type_f == "rectilinear") then
247  call mem_setptr(grid_shape_ptr, "MSHAPE", &
248  create_mem_path(model_name, 'DIS'))
249  ! The dimension of y is in the second last element of the shape array.
250  ! + 1 because we count corners, not centers.
251  y_size = grid_shape_ptr(size(grid_shape_ptr - 1)) + 1
252  grid_y(1:y_size) = [(i, i=y_size - 1, 0, -1)]
253  else if (grid_type_f == "unstructured") then
254  call mem_setptr(vertices_ptr, "VERTICES", &
255  create_mem_path(model_name, 'DIS'))
256  ! y-coordinates are in the 2nd column
257  y_size = size(vertices_ptr(2, :))
258  grid_y(1:y_size) = vertices_ptr(2, :)
259  else
260  bmi_status = bmi_failure
261  return
262  end if
263  bmi_status = bmi_success
264  end function get_grid_y
265 
266  ! NOTE: node in BMI-terms is a vertex in Modflow terms
267  ! Get the number of nodes in an unstructured grid.
268  function get_grid_node_count(grid_id, count) result(bmi_status) &
269  bind(C, name="get_grid_node_count")
270  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_node_count
271  ! -- dummy variables
272  integer(kind=c_int), intent(in) :: grid_id
273  integer(kind=c_int), intent(out) :: count
274  integer(kind=c_int) :: bmi_status
275  ! -- local variables
276  character(len=LENMODELNAME) :: model_name
277  integer(I4B), pointer :: nvert_ptr
278 
279  ! make sure function is only used for DISU grids
280  bmi_status = bmi_failure
281  if (.not. confirm_grid_type(grid_id, "DISU")) return
282 
283  model_name = get_model_name(grid_id)
284  call mem_setptr(nvert_ptr, "NVERT", create_mem_path(model_name, 'DIS'))
285  count = nvert_ptr
286  bmi_status = bmi_success
287  end function get_grid_node_count
288 
289  ! TODO_JH: This currently only works for 2D DISU models
290  ! Get the number of faces in an unstructured grid.
291  function get_grid_face_count(grid_id, count) result(bmi_status) &
292  bind(C, name="get_grid_face_count")
293  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_face_count
294  ! -- modules
295  use listsmodule, only: basemodellist
297  ! -- dummy variables
298  integer(kind=c_int), intent(in) :: grid_id
299  integer(kind=c_int), intent(out) :: count
300  integer(kind=c_int) :: bmi_status
301  ! -- local variables
302  character(len=LENMODELNAME) :: model_name
303  integer(I4B) :: i
304  class(numericalmodeltype), pointer :: numericalmodel
305 
306  ! make sure function is only used for DISU grids
307  bmi_status = bmi_failure
308  if (.not. confirm_grid_type(grid_id, "DISU")) return
309 
310  model_name = get_model_name(grid_id)
311  do i = 1, basemodellist%Count()
312  numericalmodel => getnumericalmodelfromlist(basemodellist, i)
313  if (numericalmodel%name == model_name) then
314  count = numericalmodel%dis%nodes
315  end if
316  end do
317  bmi_status = bmi_success
318  end function get_grid_face_count
319 
320  ! Get the face-node connectivity.
321  function get_grid_face_nodes(grid_id, face_nodes) result(bmi_status) &
322  bind(C, name="get_grid_face_nodes")
323  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_face_nodes
324  ! -- dummy variables
325  integer(kind=c_int), intent(in) :: grid_id
326  integer(kind=c_int), intent(out) :: face_nodes(*)
327  integer(kind=c_int) :: bmi_status
328  ! -- local variables
329  character(len=LENMODELNAME) :: model_name
330  integer, dimension(:), pointer, contiguous :: javert_ptr
331  integer, dimension(:), allocatable :: nodes_per_face
332  integer :: face_count
333  integer :: face_nodes_count
334 
335  ! make sure function is only used for DISU grids
336  bmi_status = bmi_failure
337  if (.not. confirm_grid_type(grid_id, "DISU")) return
338 
339  model_name = get_model_name(grid_id)
340  call mem_setptr(javert_ptr, "JAVERT", create_mem_path(model_name, 'DIS'))
341 
342  bmi_status = get_grid_face_count(grid_id, face_count)
343  if (bmi_status == bmi_failure) return
344 
345  allocate (nodes_per_face(face_count))
346  bmi_status = get_grid_nodes_per_face(grid_id, nodes_per_face)
347  if (bmi_status == bmi_failure) return
348 
349  face_nodes_count = sum(nodes_per_face + 1)
350 
351  face_nodes(1:face_nodes_count) = javert_ptr(:)
352  bmi_status = bmi_success
353  end function get_grid_face_nodes
354 
355  ! Get the number of nodes for each face.
356  function get_grid_nodes_per_face(grid_id, nodes_per_face) result(bmi_status) &
357  bind(C, name="get_grid_nodes_per_face")
358  !DIR$ ATTRIBUTES DLLEXPORT :: get_grid_nodes_per_face
359  ! -- dummy variables
360  integer(kind=c_int), intent(in) :: grid_id
361  integer(kind=c_int), intent(out) :: nodes_per_face(*)
362  integer(kind=c_int) :: bmi_status
363  ! -- local variables
364  integer(I4B) :: i
365  character(len=LENMODELNAME) :: model_name
366  integer, dimension(:), pointer, contiguous :: iavert_ptr
367 
368  ! make sure function is only used for DISU grids
369  bmi_status = bmi_failure
370  if (.not. confirm_grid_type(grid_id, "DISU")) return
371 
372  model_name = get_model_name(grid_id)
373  call mem_setptr(iavert_ptr, "IAVERT", create_mem_path(model_name, 'DIS'))
374 
375  do i = 1, size(iavert_ptr) - 1
376  nodes_per_face(i) = iavert_ptr(i + 1) - iavert_ptr(i) - 1
377  end do
378  bmi_status = bmi_success
379  end function get_grid_nodes_per_face
380 end module mf6bmigrid
class(basemodeltype) function, pointer, public getbasemodelfromlist(list, idx)
Definition: BaseModel.f90:172
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter lenmodelname
maximum length of the model name
Definition: Constants.f90:22
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
This module defines variable data types.
Definition: kind.f90:8
type(listtype), public basemodellist
Definition: mf6lists.f90:16
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
Detailed error information for the BMI.
Definition: mf6bmiError.f90:6
integer, parameter bmi_failure
BMI status code for failure (taken from bmi.f90, CSDMS)
Definition: mf6bmiError.f90:12
integer, parameter bmi_success
BMI status code for success (taken from bmi.f90, CSDMS)
Definition: mf6bmiError.f90:13
This module contains BMI routines to expose the MODFLOW 6 discretization.
Definition: mf6bmiGrid.f90:8
integer(kind=c_int) function get_grid_rank(grid_id, grid_rank)
Definition: mf6bmiGrid.f90:91
integer(kind=c_int) function get_grid_face_count(grid_id, count)
Definition: mf6bmiGrid.f90:293
integer(kind=c_int) function get_grid_shape(grid_id, grid_shape)
Definition: mf6bmiGrid.f90:153
integer(kind=c_int) function get_grid_face_nodes(grid_id, face_nodes)
Definition: mf6bmiGrid.f90:323
integer(kind=c_int) function get_grid_type(grid_id, grid_type)
Definition: mf6bmiGrid.f90:62
integer(kind=c_int) function get_var_grid(c_var_address, var_grid)
Definition: mf6bmiGrid.f90:23
integer(kind=c_int) function get_grid_nodes_per_face(grid_id, nodes_per_face)
Definition: mf6bmiGrid.f90:358
integer(kind=c_int) function get_grid_size(grid_id, grid_size)
Definition: mf6bmiGrid.f90:119
integer(kind=c_int) function get_grid_node_count(grid_id, count)
Definition: mf6bmiGrid.f90:270
integer(kind=c_int) function get_grid_y(grid_id, grid_y)
Definition: mf6bmiGrid.f90:226
integer(kind=c_int) function get_grid_x(grid_id, grid_x)
Definition: mf6bmiGrid.f90:182
This module contains helper routines and parameters for the MODFLOW 6 BMI.
Definition: mf6bmiUtil.f90:4
character(len=lenmodelname) function get_model_name(grid_id)
Get the model name from the grid id.
Definition: mf6bmiUtil.f90:191
subroutine get_grid_type_model(model_name, grid_type_f)
Get the grid type for a named model as a fortran string.
Definition: mf6bmiUtil.f90:239
integer(c_int), bind(C, name="BMI_LENGRIDTYPE") bmi_lengridtype
max. length for grid type C-strings
Definition: mf6bmiUtil.f90:30
integer(i4b), parameter lengridtype
max length for Fortran grid type string
Definition: mf6bmiUtil.f90:24
logical function confirm_grid_type(grid_id, expected_type)
Confirm that grid is of an expected type.
Definition: mf6bmiUtil.f90:260
pure character(kind=c_char, len=1) function, dimension(length+1) string_to_char_array(string, length)
Convert Fortran string to C-style character string.
Definition: mf6bmiUtil.f90:149
pure integer(i4b) function strlen(char_array, max_len)
Returns the string length without the trailing null character.
Definition: mf6bmiUtil.f90:113
pure character(len=length) function char_array_to_string(char_array, length)
Convert C-style string to Fortran character string.
Definition: mf6bmiUtil.f90:133
character(len=lenmodelname) function extract_model_name(var_address, success)
Extract the model name from a memory address string.
Definition: mf6bmiUtil.f90:166
class(numericalmodeltype) function, pointer, public getnumericalmodelfromlist(list, idx)
Highest level model type. All models extend this parent type.
Definition: BaseModel.f90:13