GRASS Programmer's Manual  6.5.svn(2014)-r66266
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
N_gwflow.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: groundwater flow in porous media
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2000 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18 
19 #include "grass/N_gwflow.h"
20 
21 /* *************************************************************** */
22 /* ***************** N_gwflow_data3d ***************************** */
23 /* *************************************************************** */
38 N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
39  int river, int drain)
40 {
42 
43  data = (N_gwflow_data3d *) G_calloc(1, sizeof(N_gwflow_data3d));
44 
45  data->phead = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
46  data->phead_start = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
47  data->status = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
48  data->hc_x = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
49  data->hc_y = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
50  data->hc_z = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
51  data->q = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
52  data->s = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
53  data->nf = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
54  data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
55 
56  if (river) {
57  data->river_head =
58  N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
59  data->river_leak =
60  N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
61  data->river_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
62  }
63  else {
64  data->river_head = NULL;
65  data->river_leak = NULL;
66  data->river_bed = NULL;
67  }
68 
69  if (drain) {
70  data->drain_leak =
71  N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
72  data->drain_bed = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
73  }
74  else {
75  data->drain_leak = NULL;
76  data->drain_bed = NULL;
77  }
78 
79  return data;
80 }
81 
82 /* *************************************************************** */
83 /* ********************* N_free_gwflow_data3d ******************** */
84 /* *************************************************************** */
93 {
94  if (data->phead)
95  N_free_array_3d(data->phead);
96  if (data->phead_start)
98  if (data->status)
99  N_free_array_3d(data->status);
100  if (data->hc_x)
101  N_free_array_3d(data->hc_x);
102  if (data->hc_y)
103  N_free_array_3d(data->hc_y);
104  if (data->hc_z)
105  N_free_array_3d(data->hc_z);
106  if (data->q)
107  N_free_array_3d(data->q);
108  if (data->s)
109  N_free_array_3d(data->s);
110  if (data->nf)
111  N_free_array_3d(data->nf);
112  if (data->r)
113  N_free_array_2d(data->r);
114  if (data->river_head)
116  if (data->river_leak)
118  if (data->river_bed)
119  N_free_array_3d(data->river_bed);
120  if (data->drain_leak)
122  if (data->drain_bed)
123  N_free_array_3d(data->drain_bed);
124 
125  G_free(data);
126 
127  data = NULL;
128 
129  return;
130 }
131 
132 /* *************************************************************** */
133 /* ******************** N_alloc_gwflow_data2d ******************** */
134 /* *************************************************************** */
148 N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river,
149  int drain)
150 {
152 
153  data = (N_gwflow_data2d *) G_calloc(1, sizeof(N_gwflow_data2d));
154 
155  data->phead = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
156  data->phead_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
157  data->status = N_alloc_array_2d(cols, rows, 1, CELL_TYPE);
158  data->hc_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
159  data->hc_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
160  data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
161  data->s = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
162  data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
163  data->r = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
164  data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
165  data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
166 
167  if (river) {
168  data->river_head = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
169  data->river_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
170  data->river_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
171  }
172  else {
173  data->river_head = NULL;
174  data->river_leak = NULL;
175  data->river_bed = NULL;
176  }
177 
178  if (drain) {
179  data->drain_leak = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
180  data->drain_bed = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
181  }
182  else {
183  data->drain_leak = NULL;
184  data->drain_bed = NULL;
185  }
186 
187 
188  return data;
189 }
190 
191 /* *************************************************************** */
192 /* ****************** N_free_gwflow_data2d *********************** */
193 /* *************************************************************** */
201 {
202  if (data->phead)
203  N_free_array_2d(data->phead);
204  if (data->phead_start)
206  if (data->status)
207  N_free_array_2d(data->status);
208  if (data->hc_x)
209  N_free_array_2d(data->hc_x);
210  if (data->hc_y)
211  N_free_array_2d(data->hc_y);
212  if (data->q)
213  N_free_array_2d(data->q);
214  if (data->s)
215  N_free_array_2d(data->s);
216  if (data->nf)
217  N_free_array_2d(data->nf);
218  if (data->r)
219  N_free_array_2d(data->r);
220  if (data->top)
221  N_free_array_2d(data->top);
222  if (data->bottom)
223  N_free_array_2d(data->bottom);
224  if (data->river_head)
226  if (data->river_leak)
228  if (data->river_bed)
229  N_free_array_2d(data->river_bed);
230  if (data->drain_leak)
232  if (data->drain_bed)
233  N_free_array_2d(data->drain_bed);
234 
235  G_free(data);
236 
237  data = NULL;;
238 
239  return;
240 }
241 
242 /* *************************************************************** */
243 /* ***************** N_callback_gwflow_3d ************************ */
244 /* *************************************************************** */
263 N_data_star *N_callback_gwflow_3d(void *gwdata, N_geom_data * geom, int col,
264  int row, int depth)
265 {
266  double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
267 
268  double dx, dy, dz, Ax, Ay, Az;
269 
270  double hc_x, hc_y, hc_z;
271 
272  double hc_xw, hc_yn, hc_zt;
273 
274  double hc_xe, hc_ys, hc_zb;
275 
276  double hc_start;
277 
278  double Ss, r, nf, q;
279 
280  double C, W, E, N, S, T, B, V;
281 
282  N_data_star *mat_pos;
283 
285 
286  /*cast the void pointer to the right data structure */
287  data = (N_gwflow_data3d *) gwdata;
288 
289  dx = geom->dx;
290  dy = geom->dy;
291  dz = geom->dz;
292  Az = N_get_geom_data_area_of_cell(geom, row);
293  Ay = geom->dx * geom->dz;
294  Ax = geom->dz * geom->dy;
295 
296  /*read the data from the arrays */
297  hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
298 
299  hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
300  hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
301  hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
302 
303  hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
304  hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
305  hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
306  hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
307  hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
308  hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
309 
310  hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
311  hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
312  hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
313  hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
314  hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
315  hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
316 
317  /*inner sources */
318  q = N_get_array_3d_d_value(data->q, col, row, depth);
319  /*specific yield */
320  Ss = N_get_array_3d_d_value(data->s, col, row, depth);
321  /*porosity */
322  nf = N_get_array_3d_d_value(data->nf, col, row, depth);
323 
324  /*mass balance center cell to western cell */
325  W = -1 * Ax * hc_w / dx;
326  /*mass balance center cell to eastern cell */
327  E = -1 * Ax * hc_e / dx;
328  /*mass balance center cell to northern cell */
329  N = -1 * Ay * hc_n / dy;
330  /*mass balance center cell to southern cell */
331  S = -1 * Ay * hc_s / dy;
332  /*mass balance center cell to top cell */
333  T = -1 * Az * hc_t / dz;
334  /*mass balance center cell to bottom cell */
335  B = -1 * Az * hc_b / dz;
336 
337  /*specific yield */
338  Ss = Az * dz * Ss;
339 
340  /*the diagonal entry of the matrix */
341  C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
342 
343  /*the entry in the right side b of Ax = b */
344  V = (q + hc_start * Ss / data->dt * Az);
345 
346  /*only the top cells will have recharge */
347  if (depth == geom->depths - 2) {
348  r = N_get_array_2d_d_value(data->r, col, row);
349  V += r * Az;
350  }
351 
352  G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
353 
354  /*create the 7 point star entries */
355  mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
356 
357  return mat_pos;
358 }
359 
360 /* *************************************************************** */
361 /* ****************** N_callback_gwflow_2d *********************** */
362 /* *************************************************************** */
379 N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
380  int row)
381 {
382  double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
383 
384  double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
385 
386  double dx, dy, Az;
387 
388  double hc_x, hc_y;
389 
390  double z, top;
391 
392  double hc_xw, hc_yn;
393 
394  double z_xw, z_yn;
395 
396  double hc_xe, hc_ys;
397 
398  double z_xe, z_ys;
399 
400  double hc, hc_start;
401 
402  double Ss, r, nf, q;
403 
404  double C, W, E, N, S, V;
405 
407 
408  N_data_star *mat_pos;
409 
410  double river_vect = 0; /*entry in vector */
411 
412  double river_mat = 0; /*entry in matrix */
413 
414  double drain_vect = 0; /*entry in vector */
415 
416  double drain_mat = 0; /*entry in matrix */
417 
418  /*cast the void pointer to the right data structure */
419  data = (N_gwflow_data2d *) gwdata;
420 
421  dx = geom->dx;
422  dy = geom->dy;
423  Az = N_get_geom_data_area_of_cell(geom, row);
424 
425  /*read the data from the arrays */
426  hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
427  hc = N_get_array_2d_d_value(data->phead, col, row);
428  top = N_get_array_2d_d_value(data->top, col, row);
429 
430 
431  if (hc > top) { /*If the aquifer is confined */
432  z = N_get_array_2d_d_value(data->top, col,
433  row) -
434  N_get_array_2d_d_value(data->bottom, col, row);
435  z_xw =
436  N_get_array_2d_d_value(data->top, col - 1,
437  row) -
438  N_get_array_2d_d_value(data->bottom, col - 1, row);
439  z_xe =
440  N_get_array_2d_d_value(data->top, col + 1,
441  row) -
442  N_get_array_2d_d_value(data->bottom, col + 1, row);
443  z_yn =
444  N_get_array_2d_d_value(data->top, col,
445  row - 1) -
446  N_get_array_2d_d_value(data->bottom, col, row - 1);
447  z_ys =
448  N_get_array_2d_d_value(data->top, col,
449  row + 1) -
450  N_get_array_2d_d_value(data->bottom, col, row + 1);
451  }
452  else { /* the aquifer is unconfined */
453 
454  /* If the aquifer is unconfied use an explicite scheme to solve
455  * the nonlinear equation. We use the phead from the first iteration */
456  z = N_get_array_2d_d_value(data->phead, col, row) -
457  N_get_array_2d_d_value(data->bottom, col, row);
458  z_xw = N_get_array_2d_d_value(data->phead, col - 1, row) -
459  N_get_array_2d_d_value(data->bottom, col - 1, row);
460  z_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
461  N_get_array_2d_d_value(data->bottom, col + 1, row);
462  z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
463  N_get_array_2d_d_value(data->bottom, col, row - 1);
464  z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
465  N_get_array_2d_d_value(data->bottom, col, row + 1);
466  }
467 
468  /*geometrical mean of cell height */
469  if (z_w > 0 || z_w < 0 || z_w == 0)
470  z_w = N_calc_arith_mean(z_xw, z);
471  else
472  z_w = z;
473  if (z_e > 0 || z_e < 0 || z_e == 0)
474  z_e = N_calc_arith_mean(z_xe, z);
475  else
476  z_e = z;
477  if (z_n > 0 || z_n < 0 || z_n == 0)
478  z_n = N_calc_arith_mean(z_yn, z);
479  else
480  z_n = z;
481  if (z_s > 0 || z_s < 0 || z_s == 0)
482  z_s = N_calc_arith_mean(z_ys, z);
483  else
484  z_s = z;
485 
486  /* Inner sources */
487  q = N_get_array_2d_d_value(data->q, col, row);
488  nf = N_get_array_2d_d_value(data->nf, col, row);
489 
490  /* specific yield of current cell face */
491  Ss = N_get_array_2d_d_value(data->s, col, row) * Az;
492  /* recharge */
493  r = N_get_array_2d_d_value(data->r, col, row) * Az;
494 
495  /*get the surrounding permeabilities */
496  hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
497  hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
498  hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
499  hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
500  hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
501  hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
502 
503  /* calculate the transmissivities */
504  T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
505  T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
506  T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
507  T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
508 
509  /*compute the river leakage, this is an explicit method */
510  if (data->river_leak &&
511  (N_get_array_2d_d_value(data->river_leak, col, row) != 0)) {
512  if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
513  river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
514  N_get_array_2d_d_value(data->river_leak, col, row);
515  river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
516  }
517  else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
518  river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
519  N_get_array_2d_d_value(data->river_bed, col, row))
520  * N_get_array_2d_d_value(data->river_leak, col, row);
521  river_mat = 0;
522  }
523  }
524 
525  /*compute the drainage, this is an explicit method */
526  if (data->drain_leak &&
527  (N_get_array_2d_d_value(data->drain_leak, col, row) != 0)) {
528  if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
529  drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
530  N_get_array_2d_d_value(data->drain_leak, col, row);
531  drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
532  }
533  else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
534  drain_vect = 0;
535  drain_mat = 0;
536  }
537  }
538 
539  /*mass balance center cell to western cell */
540  W = -1 * T_w * dy / dx;
541  /*mass balance center cell to eastern cell */
542  E = -1 * T_e * dy / dx;
543  /*mass balance center cell to northern cell */
544  N = -1 * T_n * dx / dy;
545  /*mass balance center cell to southern cell */
546  S = -1 * T_s * dx / dy;
547 
548  /*the diagonal entry of the matrix */
549  C = -1 * (W + E + N + S - Ss / data->dt - river_mat * Az -
550  drain_mat * Az);
551 
552  /*the entry in the right side b of Ax = b */
553  V = (q + hc_start * Ss / data->dt) + r + river_vect * Az +
554  drain_vect * Az;
555 
556  G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
557 
558  /*create the 5 point star entries */
559  mat_pos = N_create_5star(C, W, E, N, S, V);
560 
561  return mat_pos;
562 }
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
Definition: N_arrays.c:127
N_array_2d * status
Definition: N_gwflow.h:90
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
N_array_2d * top
Definition: N_gwflow.h:87
Matrix entries for a mass balance 5/7/9 star system.
Definition: N_pde.h:342
N_array_2d * s
Definition: N_gwflow.h:74
tuple q
Definition: forms.py:2019
N_array_2d * hc_y
Definition: N_gwflow.h:71
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
Definition: N_arrays.c:777
float r
Definition: named_colr.c:8
N_array_2d * river_bed
Definition: N_gwflow.h:80
N_array_3d * river_head
Definition: N_gwflow.h:48
This data structure contains all data needed to compute the groundwater mass balance in three dimensi...
Definition: N_gwflow.h:34
N_array_2d * phead_start
Definition: N_gwflow.h:69
N_array_3d * river_leak
Definition: N_gwflow.h:47
N_array_2d * phead
Definition: N_gwflow.h:68
N_array_2d * q
Definition: N_gwflow.h:72
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
Definition: N_arrays.c:69
N_array_3d * s
Definition: N_gwflow.h:43
N_array_3d * nf
Definition: N_gwflow.h:44
int depths
Definition: N_pde.h:139
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
Definition: N_arrays.c:723
N_array_3d * river_bed
Definition: N_gwflow.h:49
double dy
Definition: N_pde.h:134
N_array_3d * hc_z
Definition: N_gwflow.h:40
N_array_2d * r
Definition: N_gwflow.h:73
N_array_2d * r
Definition: N_gwflow.h:42
#define C
Definition: intr_char.c:17
Geometric information about the structured grid.
Definition: N_pde.h:127
N_gwflow_data3d * N_alloc_gwflow_data3d(int cols, int rows, int depths, int river, int drain)
Allocate memory for the groundwater calculation data structure in 3 dimensions.
Definition: N_gwflow.c:38
double dz
Definition: N_pde.h:135
N_array_3d * status
Definition: N_gwflow.h:55
tuple data
N_data_star * N_callback_gwflow_3d(void *gwdata, N_geom_data *geom, int col, int row, int depth)
This callback function creates the mass balance of a 7 point star.
Definition: N_gwflow.c:263
N_data_star * N_create_7star(double C, double W, double E, double N, double S, double T, double B, double V)
allocate and initialize a 7 point star data structure
N_array_3d * drain_leak
Definition: N_gwflow.h:52
void N_free_gwflow_data3d(N_gwflow_data3d *data)
Release the memory of the groundwater flow data structure in three dimensions.
Definition: N_gwflow.c:92
N_array_3d * phead_start
Definition: N_gwflow.h:37
N_array_3d * hc_y
Definition: N_gwflow.h:39
N_data_star * N_callback_gwflow_2d(void *gwdata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
Definition: N_gwflow.c:379
N_array_2d * hc_x
Definition: N_gwflow.h:70
N_data_star * N_create_5star(double C, double W, double E, double N, double S, double V)
allocate and initialize a 5 point star data structure
N_array_2d * river_head
Definition: N_gwflow.h:79
N_array_2d * nf
Definition: N_gwflow.h:75
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
Get the areay size in square meter of one cell (x*y) at row.
Definition: N_geom.c:199
double N_calc_arith_mean(double a, double b)
Calculate the arithmetic mean of values a and b.
Definition: N_tools.c:33
N_array_2d * drain_bed
Definition: N_gwflow.h:84
N_array_3d * hc_x
Definition: N_gwflow.h:38
N_array_2d * drain_leak
Definition: N_gwflow.h:83
N_gwflow_data2d * N_alloc_gwflow_data2d(int cols, int rows, int river, int drain)
Alllocate memory for the groundwater calculation data structure in 2 dimensions.
Definition: N_gwflow.c:148
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
Definition: N_arrays.c:987
return NULL
Definition: dbfopen.c:1394
This data structure contains all data needed to compute the groundwater mass balance in two dimension...
Definition: N_gwflow.h:66
tuple cols
void N_free_gwflow_data2d(N_gwflow_data2d *data)
Release the memory of the groundwater flow data structure in two dimensions.
Definition: N_gwflow.c:200
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
N_array_3d * drain_bed
Definition: N_gwflow.h:53
N_array_3d * phead
Definition: N_gwflow.h:36
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
Definition: N_arrays.c:375
#define N
Definition: inverse.c:8
N_array_2d * bottom
Definition: N_gwflow.h:88
N_array_3d * q
Definition: N_gwflow.h:41
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
Definition: N_tools.c:121
double dx
Definition: N_pde.h:133
N_array_2d * river_leak
Definition: N_gwflow.h:78