← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/gpde/N_gwflow.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/gpde/N_gwflow.c
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
 
/* *************************************************************** */
24
 
/*!
25
 
 * \brief Alllocate memory for the groundwater calculation data structure in 3 dimensions
26
 
 *
27
 
 * The groundwater calculation data structure will be allocated including
28
 
 * all appendant 3d and 2d arrays. The offset for the 3d arrays is one
29
 
 * to establish homogeneous Neumann boundary conditions at the calculation area border.
30
 
 * This data structure is used to create a linear equation system based on the computation of
31
 
 * groundwater flow in porous media with the finite volume method.
32
 
 *
33
 
 * \param cols   int
34
 
 * \param rows   int
35
 
 * \param depths int
36
 
 * \return N_gwflow_data3d *
37
 
 * */
38
 
N_gwflow_data3d *N_alloc_gwflow_data3d(int cols, int rows, int depths,
39
 
                                       int river, int drain)
40
 
{
41
 
    N_gwflow_data3d *data;
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
 
/* *************************************************************** */
85
 
/*!
86
 
 * \brief Release the memory of the groundwater flow data structure in three dimensions
87
 
 *
88
 
 * \param data N_gwflow_data3d *
89
 
 * \return void *
90
 
 * */
91
 
 
92
 
void N_free_gwflow_data3d(N_gwflow_data3d * data)
93
 
{
94
 
    if (data->phead)
95
 
        N_free_array_3d(data->phead);
96
 
    if (data->phead_start)
97
 
        N_free_array_3d(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)
115
 
        N_free_array_3d(data->river_head);
116
 
    if (data->river_leak)
117
 
        N_free_array_3d(data->river_leak);
118
 
    if (data->river_bed)
119
 
        N_free_array_3d(data->river_bed);
120
 
    if (data->drain_leak)
121
 
        N_free_array_3d(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
 
/* *************************************************************** */
135
 
/*!
136
 
 * \brief Alllocate memory for the groundwater calculation data structure in 2 dimensions
137
 
 * 
138
 
 * The groundwater calculation data structure will be allocated including
139
 
 * all appendant 2d arrays. The offset for the 3d arrays is one
140
 
 * to establish homogeneous Neumann boundary conditions at the calculation area border.
141
 
 * This data structure is used to create a linear equation system based on the computation of
142
 
 * groundwater flow in porous media with the finite volume method.
143
 
 *
144
 
 * \param cols int
145
 
 * \param rows int
146
 
 * \return N_gwflow_data2d *
147
 
 * */
148
 
N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river,
149
 
                                       int drain)
150
 
{
151
 
    N_gwflow_data2d *data;
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
 
/* *************************************************************** */
194
 
/*!
195
 
 * \brief Release the memory of the groundwater flow data structure in two dimensions
196
 
 *
197
 
 * \param data N_gwflow_data2d *
198
 
 * \return void
199
 
 * */
200
 
void N_free_gwflow_data2d(N_gwflow_data2d * data)
201
 
{
202
 
    if (data->phead)
203
 
        N_free_array_2d(data->phead);
204
 
    if (data->phead_start)
205
 
        N_free_array_2d(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)
225
 
        N_free_array_2d(data->river_head);
226
 
    if (data->river_leak)
227
 
        N_free_array_2d(data->river_leak);
228
 
    if (data->river_bed)
229
 
        N_free_array_2d(data->river_bed);
230
 
    if (data->drain_leak)
231
 
        N_free_array_2d(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
 
/* *************************************************************** */
245
 
/*!
246
 
 * \brief This callback function creates the mass balance of a 7 point star
247
 
 *
248
 
 * The mass balance is based on the common groundwater flow equation:
249
 
 *
250
 
 * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
251
 
 *
252
 
 * This equation is discretizised with the finite volume method in three dimensions.
253
 
 *
254
 
 *
255
 
 * \param gwdata N_gwflow_data3d *
256
 
 * \param geom N_geom_data *
257
 
 * \param col   int
258
 
 * \param row   int
259
 
 * \param depth int
260
 
 * \return N_data_star *
261
 
 *
262
 
 * */
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
 
    double dx, dy, dz, Ax, Ay, Az;
268
 
    double hc_x, hc_y, hc_z;
269
 
    double hc_xw, hc_yn, hc_zt;
270
 
    double hc_xe, hc_ys, hc_zb;
271
 
    double hc_start;
272
 
    double Ss, r, nf, q;
273
 
    double C, W, E, N, S, T, B, V;
274
 
    N_data_star *mat_pos;
275
 
    N_gwflow_data3d *data;
276
 
 
277
 
    /*cast the void pointer to the right data structure */
278
 
    data = (N_gwflow_data3d *) gwdata;
279
 
 
280
 
    dx = geom->dx;
281
 
    dy = geom->dy;
282
 
    dz = geom->dz;
283
 
    Az = N_get_geom_data_area_of_cell(geom, row);
284
 
    Ay = geom->dx * geom->dz;
285
 
    Ax = geom->dz * geom->dy;
286
 
 
287
 
    /*read the data from the arrays */
288
 
    hc_start = N_get_array_3d_d_value(data->phead_start, col, row, depth);
289
 
 
290
 
    hc_x = N_get_array_3d_d_value(data->hc_x, col, row, depth);
291
 
    hc_y = N_get_array_3d_d_value(data->hc_y, col, row, depth);
292
 
    hc_z = N_get_array_3d_d_value(data->hc_z, col, row, depth);
293
 
 
294
 
    hc_xw = N_get_array_3d_d_value(data->hc_x, col - 1, row, depth);
295
 
    hc_xe = N_get_array_3d_d_value(data->hc_x, col + 1, row, depth);
296
 
    hc_yn = N_get_array_3d_d_value(data->hc_y, col, row - 1, depth);
297
 
    hc_ys = N_get_array_3d_d_value(data->hc_y, col, row + 1, depth);
298
 
    hc_zt = N_get_array_3d_d_value(data->hc_z, col, row, depth + 1);
299
 
    hc_zb = N_get_array_3d_d_value(data->hc_z, col, row, depth - 1);
300
 
 
301
 
    hc_w = N_calc_harmonic_mean(hc_xw, hc_x);
302
 
    hc_e = N_calc_harmonic_mean(hc_xe, hc_x);
303
 
    hc_n = N_calc_harmonic_mean(hc_yn, hc_y);
304
 
    hc_s = N_calc_harmonic_mean(hc_ys, hc_y);
305
 
    hc_t = N_calc_harmonic_mean(hc_zt, hc_z);
306
 
    hc_b = N_calc_harmonic_mean(hc_zb, hc_z);
307
 
 
308
 
    /*inner sources */
309
 
    q = N_get_array_3d_d_value(data->q, col, row, depth);
310
 
    /*specific yield */
311
 
    Ss = N_get_array_3d_d_value(data->s, col, row, depth);
312
 
    /*porosity */
313
 
    nf = N_get_array_3d_d_value(data->nf, col, row, depth);
314
 
 
315
 
    /*mass balance center cell to western cell */
316
 
    W = -1 * Ax * hc_w / dx;
317
 
    /*mass balance center cell to eastern cell */
318
 
    E = -1 * Ax * hc_e / dx;
319
 
    /*mass balance center cell to northern cell */
320
 
    N = -1 * Ay * hc_n / dy;
321
 
    /*mass balance center cell to southern cell */
322
 
    S = -1 * Ay * hc_s / dy;
323
 
    /*mass balance center cell to top cell */
324
 
    T = -1 * Az * hc_t / dz;
325
 
    /*mass balance center cell to bottom cell */
326
 
    B = -1 * Az * hc_b / dz;
327
 
 
328
 
    /*specific yield */
329
 
    Ss = Az * dz * Ss;
330
 
 
331
 
    /*the diagonal entry of the matrix */
332
 
    C = -1 * (W + E + N + S + T + B - Ss / data->dt * Az);
333
 
 
334
 
    /*the entry in the right side b of Ax = b */
335
 
    V = (q + hc_start * Ss / data->dt * Az);
336
 
 
337
 
    /*only the top cells will have recharge */
338
 
    if (depth == geom->depths - 2) {
339
 
        r = N_get_array_2d_d_value(data->r, col, row);
340
 
        V += r * Az;
341
 
    }
342
 
 
343
 
    G_debug(5, "N_callback_gwflow_3d: called [%i][%i][%i]", depth, col, row);
344
 
 
345
 
    /*create the 7 point star entries */
346
 
    mat_pos = N_create_7star(C, W, E, N, S, T, B, V);
347
 
 
348
 
    return mat_pos;
349
 
}
350
 
 
351
 
/* *************************************************************** */
352
 
/* ****************** N_callback_gwflow_2d *********************** */
353
 
/* *************************************************************** */
354
 
/*!
355
 
 * \brief This callback function creates the mass balance of a 5 point star
356
 
 *
357
 
 * The mass balance is based on the common groundwater flow equation:
358
 
 *
359
 
 * \f[Ss \frac{\partial h}{\partial t} = \nabla {\bf K} \nabla h + q \f]
360
 
 *
361
 
 * This equation is discretizised with the finite volume method in two dimensions.
362
 
 *
363
 
 * \param gwdata N_gwflow_data2d *
364
 
 * \param geom N_geom_data *
365
 
 * \param col int
366
 
 * \param row int
367
 
 * \return N_data_star *
368
 
 *
369
 
 * */
370
 
N_data_star *N_callback_gwflow_2d(void *gwdata, N_geom_data * geom, int col,
371
 
                                  int row)
372
 
{
373
 
    double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
374
 
    double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
375
 
    double dx, dy, Az;
376
 
    double hc_x, hc_y;
377
 
    double z, top;
378
 
    double hc_xw, hc_yn;
379
 
    double z_xw, z_yn;
380
 
    double hc_xe, hc_ys;
381
 
    double z_xe, z_ys;
382
 
    double hc, hc_start;
383
 
    double Ss, r, nf, q;
384
 
    double C, W, E, N, S, V;
385
 
    N_gwflow_data2d *data;
386
 
    N_data_star *mat_pos;
387
 
    double river_vect = 0;      /*entry in vector */
388
 
    double river_mat = 0;       /*entry in matrix */
389
 
    double drain_vect = 0;      /*entry in vector */
390
 
    double drain_mat = 0;       /*entry in matrix */
391
 
 
392
 
    /*cast the void pointer to the right data structure */
393
 
    data = (N_gwflow_data2d *) gwdata;
394
 
 
395
 
    dx = geom->dx;
396
 
    dy = geom->dy;
397
 
    Az = N_get_geom_data_area_of_cell(geom, row);
398
 
 
399
 
    /*read the data from the arrays */
400
 
    hc_start = N_get_array_2d_d_value(data->phead_start, col, row);
401
 
    hc = N_get_array_2d_d_value(data->phead, col, row);
402
 
    top = N_get_array_2d_d_value(data->top, col, row);
403
 
 
404
 
 
405
 
    if (hc > top) {             /*If the aquifer is confined */
406
 
        z = N_get_array_2d_d_value(data->top, col,
407
 
                                   row) -
408
 
            N_get_array_2d_d_value(data->bottom, col, row);
409
 
        z_xw =
410
 
            N_get_array_2d_d_value(data->top, col - 1,
411
 
                                   row) -
412
 
            N_get_array_2d_d_value(data->bottom, col - 1, row);
413
 
        z_xe =
414
 
            N_get_array_2d_d_value(data->top, col + 1,
415
 
                                   row) -
416
 
            N_get_array_2d_d_value(data->bottom, col + 1, row);
417
 
        z_yn =
418
 
            N_get_array_2d_d_value(data->top, col,
419
 
                                   row - 1) -
420
 
            N_get_array_2d_d_value(data->bottom, col, row - 1);
421
 
        z_ys =
422
 
            N_get_array_2d_d_value(data->top, col,
423
 
                                   row + 1) -
424
 
            N_get_array_2d_d_value(data->bottom, col, row + 1);
425
 
    }
426
 
    else {                      /* the aquifer is unconfined */
427
 
 
428
 
        /* If the aquifer is unconfied use an explicite scheme to solve
429
 
         * the nonlinear equation. We use the phead from the first iteration */
430
 
        z = N_get_array_2d_d_value(data->phead, col, row) -
431
 
            N_get_array_2d_d_value(data->bottom, col, row);
432
 
        z_xw = N_get_array_2d_d_value(data->phead, col - 1, row) -
433
 
            N_get_array_2d_d_value(data->bottom, col - 1, row);
434
 
        z_xe = N_get_array_2d_d_value(data->phead, col + 1, row) -
435
 
            N_get_array_2d_d_value(data->bottom, col + 1, row);
436
 
        z_yn = N_get_array_2d_d_value(data->phead, col, row - 1) -
437
 
            N_get_array_2d_d_value(data->bottom, col, row - 1);
438
 
        z_ys = N_get_array_2d_d_value(data->phead, col, row + 1) -
439
 
            N_get_array_2d_d_value(data->bottom, col, row + 1);
440
 
    }
441
 
 
442
 
    /*geometrical mean of cell height */
443
 
    if (z_w > 0 || z_w < 0 || z_w == 0)
444
 
        z_w = N_calc_arith_mean(z_xw, z);
445
 
    else
446
 
        z_w = z;
447
 
    if (z_e > 0 || z_e < 0 || z_e == 0)
448
 
        z_e = N_calc_arith_mean(z_xe, z);
449
 
    else
450
 
        z_e = z;
451
 
    if (z_n > 0 || z_n < 0 || z_n == 0)
452
 
        z_n = N_calc_arith_mean(z_yn, z);
453
 
    else
454
 
        z_n = z;
455
 
    if (z_s > 0 || z_s < 0 || z_s == 0)
456
 
        z_s = N_calc_arith_mean(z_ys, z);
457
 
    else
458
 
        z_s = z;
459
 
 
460
 
    /* Inner sources */
461
 
    q = N_get_array_2d_d_value(data->q, col, row);
462
 
    nf = N_get_array_2d_d_value(data->nf, col, row);
463
 
 
464
 
    /* specific yield  of current cell face */
465
 
    Ss = N_get_array_2d_d_value(data->s, col, row) * Az;
466
 
    /* recharge */
467
 
    r = N_get_array_2d_d_value(data->r, col, row) * Az;
468
 
 
469
 
    /*get the surrounding permeabilities */
470
 
    hc_x = N_get_array_2d_d_value(data->hc_x, col, row);
471
 
    hc_y = N_get_array_2d_d_value(data->hc_y, col, row);
472
 
    hc_xw = N_get_array_2d_d_value(data->hc_x, col - 1, row);
473
 
    hc_xe = N_get_array_2d_d_value(data->hc_x, col + 1, row);
474
 
    hc_yn = N_get_array_2d_d_value(data->hc_y, col, row - 1);
475
 
    hc_ys = N_get_array_2d_d_value(data->hc_y, col, row + 1);
476
 
 
477
 
    /* calculate the transmissivities */
478
 
    T_w = N_calc_harmonic_mean(hc_xw, hc_x) * z_w;
479
 
    T_e = N_calc_harmonic_mean(hc_xe, hc_x) * z_e;
480
 
    T_n = N_calc_harmonic_mean(hc_yn, hc_y) * z_n;
481
 
    T_s = N_calc_harmonic_mean(hc_ys, hc_y) * z_s;
482
 
 
483
 
    /*compute the river leakage, this is an explicit method */
484
 
    if (data->river_leak &&
485
 
        (N_get_array_2d_d_value(data->river_leak, col, row) != 0)) {
486
 
        if (hc > N_get_array_2d_d_value(data->river_bed, col, row)) {
487
 
            river_vect = N_get_array_2d_d_value(data->river_head, col, row) *
488
 
                N_get_array_2d_d_value(data->river_leak, col, row);
489
 
            river_mat = N_get_array_2d_d_value(data->river_leak, col, row);
490
 
        }
491
 
        else if (hc < N_get_array_2d_d_value(data->river_bed, col, row)) {
492
 
            river_vect = (N_get_array_2d_d_value(data->river_head, col, row) -
493
 
                          N_get_array_2d_d_value(data->river_bed, col, row))
494
 
                * N_get_array_2d_d_value(data->river_leak, col, row);
495
 
            river_mat = 0;
496
 
        }
497
 
    }
498
 
 
499
 
    /*compute the drainage, this is an explicit method */
500
 
    if (data->drain_leak &&
501
 
        (N_get_array_2d_d_value(data->drain_leak, col, row) != 0)) {
502
 
        if (hc > N_get_array_2d_d_value(data->drain_bed, col, row)) {
503
 
            drain_vect = N_get_array_2d_d_value(data->drain_bed, col, row) *
504
 
                N_get_array_2d_d_value(data->drain_leak, col, row);
505
 
            drain_mat = N_get_array_2d_d_value(data->drain_leak, col, row);
506
 
        }
507
 
        else if (hc <= N_get_array_2d_d_value(data->drain_bed, col, row)) {
508
 
            drain_vect = 0;
509
 
            drain_mat = 0;
510
 
        }
511
 
    }
512
 
 
513
 
    /*mass balance center cell to western cell */
514
 
    W = -1 * T_w * dy / dx;
515
 
    /*mass balance center cell to eastern cell */
516
 
    E = -1 * T_e * dy / dx;
517
 
    /*mass balance center cell to northern cell */
518
 
    N = -1 * T_n * dx / dy;
519
 
    /*mass balance center cell to southern cell */
520
 
    S = -1 * T_s * dx / dy;
521
 
 
522
 
    /*the diagonal entry of the matrix */
523
 
    C = -1 * (W + E + N + S - Ss / data->dt - river_mat * Az -
524
 
              drain_mat * Az);
525
 
 
526
 
    /*the entry in the right side b of Ax = b */
527
 
    V = (q + hc_start * Ss / data->dt) + r + river_vect * Az +
528
 
        drain_vect * Az;
529
 
 
530
 
    G_debug(5, "N_callback_gwflow_2d: called [%i][%i]", row, col);
531
 
 
532
 
    /*create the 5 point star entries */
533
 
    mat_pos = N_create_5star(C, W, E, N, S, V);
534
 
 
535
 
    return mat_pos;
536
 
}