← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to raster/simwe/simlib/hydro.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/simwe/simlib/hydro.c
1
 
 
2
 
/****************************************************************************
3
 
 *
4
 
 * MODULE:       simwe library
5
 
 * AUTHOR(S):    Helena Mitasova, Jaro Hofierka, Lubos Mitas:
6
 
 * PURPOSE:      Hydrologic and sediment transport simulation (SIMWE)
7
 
 *
8
 
 * COPYRIGHT:    (C) 2002 by the GRASS Development Team
9
 
 *
10
 
 *               This program is free software under the GNU General Public
11
 
 *               License (>=v2). Read the file COPYING that comes with GRASS
12
 
 *               for details.
13
 
 *
14
 
 *****************************************************************************/
15
 
 
16
 
/* hydro.c (simlib), 20.nov.2002, JH */
17
 
 
18
 
#include <stdio.h>
19
 
#include <stdlib.h>
20
 
#include <math.h>
21
 
#include <grass/gis.h>
22
 
/* #include <grass/site.h> */
23
 
#include <grass/bitmap.h>
24
 
#include <grass/linkm.h>
25
 
#include <grass/glocale.h>
26
 
 
27
 
#define MAIN
28
 
#include <grass/waterglobs.h>
29
 
 
30
 
/* **************************************************** */
31
 
/*       create walker representation of si */
32
 
/* ******************************************************** */
33
 
/*                       .......... iblock loop */
34
 
 
35
 
void main_loop(void)
36
 
{
37
 
 
38
 
    int i, ii, l, k;
39
 
/*    int icoub, lwout, nmult; */
40
 
    int icoub, nmult;
41
 
    int iw, iblock, lw;
42
 
    int itime, iter1;
43
 
    int nfiterh, nfiterw;
44
 
    int mgen, mgen2, mgen3;
45
 
    int nblock;
46
 
    int icfl;
47
 
    int mitfac;
48
 
/*  int mitfac, p; */
49
 
    double x, y;
50
 
    double velx, vely, stxm, stym;
51
 
    double factor, conn, gaux, gauy;
52
 
    double d1, addac, decr;
53
 
    double barea, sarea, walkwe;
54
 
    double gen, gen2, wei2, wei3, wei, weifac;
55
 
    float eff;
56
 
 
57
 
    nblock = 1;
58
 
    icoub = 0;
59
 
    icfl = 0;
60
 
/*    nstack = 0; */
61
 
 
62
 
    if (maxwa > (MAXW - mx * my)) {
63
 
        mitfac = maxwa / (MAXW - mx * my);
64
 
        nblock = mitfac + 1;
65
 
        maxwa = maxwa / nblock;
66
 
    }
67
 
 
68
 
    G_debug(2, " maxwa, nblock %d %d", maxwa, nblock);
69
 
 
70
 
    for (iblock = 1; iblock <= nblock; iblock++) {
71
 
        ++icoub;
72
 
 
73
 
        lw = 0;
74
 
        walkwe = 0.;
75
 
        barea = stepx * stepy;
76
 
        sarea = bresx * bresy;
77
 
        G_debug(2, " barea,sarea,rwalk,sisum: %f %f %f %f", barea, sarea,
78
 
                rwalk, sisum);
79
 
        lwwfin = 0;
80
 
        /* write hh.walkers0 */
81
 
 
82
 
        for (k = 0; k < my; k++) {
83
 
            for (l = 0; l < mx; l++) {  /* run thru the whole area */
84
 
                if (zz[k][l] != UNDEF) {
85
 
 
86
 
                    x = xp0 + stepx * (double)(l);
87
 
                    y = yp0 + stepy * (double)(k);
88
 
 
89
 
                    gen = rwalk * si[k][l] / sisum;
90
 
                    mgen = (int)gen;
91
 
                    wei = gen / (double)(mgen + 1);
92
 
 
93
 
                    /*if (si[k][l] != 0.) { */
94
 
                    /* this stuff later for multiscale */
95
 
 
96
 
                    gen2 =
97
 
                        (double)maxwab *si[k][l] / (si0 *
98
 
                                                    (double)(mx2o * my2o));
99
 
                    gen2 = gen2 * (barea / sarea);
100
 
                    mgen2 = (int)gen2;
101
 
                    wei2 = gen2 / (double)(mgen2 + 1);
102
 
                    mgen3 =
103
 
                        (int)((double)mgen2 * wei2 / ((double)mgen * wei));
104
 
                    nmult = mgen3 + 1;
105
 
                    wei3 = gen2 / (double)((mgen + 1) * (mgen2 + 1));
106
 
                    weifac = wei3 / wei;
107
 
                    /*              } else {
108
 
                       nmult = 1;
109
 
                       weifac = 1.;
110
 
                       fprintf(stderr, "\n zero rainfall excess in cell"); 
111
 
                       } */
112
 
 
113
 
                    /*G_debug(2, " gen,gen2,wei,wei2,mgen3,nmult: %f %f %f %f %d %d",gen,gen2,wei,wei2,mgen3,nmult);
114
 
                     */
115
 
                    for (iw = 1; iw <= mgen + 1; iw++) {        /* assign walkers */
116
 
                        ++lw;
117
 
 
118
 
                        if (lw > MAXW)
119
 
                            G_fatal_error(_("nwalk (%d) > maxw (%d)!"), lw, MAXW);
120
 
 
121
 
                        w[lw][1] = x + stepx * (ulec() - 0.5);
122
 
                        w[lw][2] = y + stepy * (ulec() - 0.5);
123
 
                        w[lw][3] = wei;
124
 
 
125
 
                        walkwe += w[lw][3];
126
 
                        vavg[lw][1] = v1[k][l];
127
 
                        vavg[lw][2] = v2[k][l];
128
 
                        if (w[lw][1] >= xmin && w[lw][2] >= ymin &&
129
 
                            w[lw][1] <= xmax && w[lw][2] <= ymax) {
130
 
                            iflag[lw] = 0;
131
 
                        }
132
 
                        else {
133
 
                            iflag[lw] = 1;
134
 
                        }
135
 
 
136
 
/*
137
 
                        lwout = lw / ldemo;
138
 
                        lwout *= ldemo;
139
 
 
140
 
                        if (lwout == lw) {
141
 
                            ++lwwfin;
142
 
                        }
143
 
*/
144
 
                    }
145
 
                }               /*DEFined area */
146
 
            }
147
 
        }
148
 
        nwalk = lw;
149
 
        G_debug(2, " number of written walkers: %d", lwwfin);
150
 
        G_debug(2, " nwalk, maxw %d %d", nwalk, MAXW);
151
 
        G_debug(2, " walkwe (walk weight),frac %f %f", walkwe, frac);
152
 
 
153
 
        stxm = stepx * (double)(mx + 1) - xmin;
154
 
        stym = stepy * (double)(my + 1) - ymin;
155
 
        nwalka = 0;
156
 
        deldif = sqrt(deltap) * frac;   /* diffuse factor */
157
 
 
158
 
 
159
 
        factor = deltap * sisum / (rwalk * (double)nblock);
160
 
 
161
 
        G_debug(2, " deldif,factor %f %e", deldif, factor);
162
 
 
163
 
        /* ********************************************************** */
164
 
        /*       main loop over the projection time */
165
 
        /* *********************************************************** */
166
 
 
167
 
 
168
 
        G_debug(2, "main loop over the projection time... ");
169
 
 
170
 
        for (i = 1; i <= miter; i++) {  /* iteration loop depending on simulation time and deltap */
171
 
            G_percent(i, miter, 1);
172
 
            iter1 = i / iterout;
173
 
            iter1 *= iterout;
174
 
            if (iter1 == i) {
175
 
                nfiterw = i / iterout + 10;
176
 
                nfiterh = i / iterout + 40;
177
 
                G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d",
178
 
                        iblock, i, miter, nwalk, nwalka);
179
 
            }
180
 
 
181
 
            if (nwalka == 0 && i > 1)
182
 
                goto L_800;
183
 
 
184
 
            /* ************************************************************ */
185
 
            /*                               .... propagate one step */
186
 
            /* ************************************************************ */
187
 
 
188
 
            addac = factor;
189
 
            conn = (double)nblock / (double)iblock;
190
 
            if (i == 1) {
191
 
                addac = factor * .5;
192
 
            }
193
 
            nwalka = 0;
194
 
 
195
 
            for (lw = 1; lw <= nwalk; lw++) {
196
 
                if (w[lw][3] > EPS) {   /* check the walker weight */
197
 
                    ++nwalka;
198
 
                    l = (int)((w[lw][1] + stxm) / stepx) - mx - 1;
199
 
                    k = (int)((w[lw][2] + stym) / stepy) - my - 1;
200
 
 
201
 
                    if (l > mx - 1 || k > my - 1 || k < 0 || l < 0) {
202
 
 
203
 
                        G_debug(2, " k,l=%d,%d", k, l);
204
 
                        printf("    lw,w=%d %f %f", lw, w[lw][1], w[lw][2]);
205
 
                        G_debug(2, "    stxym=%f %f", stxm, stym);
206
 
                        printf("    step=%f %f", stepx, stepy);
207
 
                        G_debug(2, "    m=%d %d", my, mx);
208
 
                        printf("    nwalka,nwalk=%d %d", nwalka, nwalk);
209
 
                        G_debug(2, "  ");
210
 
                    }
211
 
 
212
 
                    if (zz[k][l] != UNDEF) {
213
 
 
214
 
                        if (infil != NULL) {    /* infiltration part */
215
 
 
216
 
                            if (inf[k][l] - si[k][l] > 0.) {
217
 
 
218
 
                                decr = pow(addac * w[lw][3], 3. / 5.);  /* decreasing factor in m */
219
 
                                if (inf[k][l] > decr) {
220
 
                                    inf[k][l] -= decr;  /* decrease infilt. in cell and eliminate the walker */
221
 
                                    w[lw][3] = 0.;
222
 
                                }
223
 
                                else {
224
 
                                    w[lw][3] -= pow(inf[k][l], 5. / 3.) / addac;        /* use just proportional part of the walker weight */
225
 
                                    inf[k][l] = 0.;
226
 
 
227
 
                                }
228
 
 
229
 
                            }
230
 
 
231
 
                        }
232
 
 
233
 
 
234
 
                        gama[k][l] += (addac * w[lw][3]);       /* add walker weigh to water depth or conc. */
235
 
 
236
 
                        d1 = gama[k][l] * conn;
237
 
                        hhc = pow(d1, 3. / 5.);
238
 
 
239
 
                        if (hhc > hhmax && wdepth == NULL) {    /* increased diffusion if w.depth > hhmax */
240
 
                            dif[k][l] = (halpha + 1) * deldif;
241
 
                            velx = vavg[lw][1];
242
 
                            vely = vavg[lw][2];
243
 
                        }
244
 
                        else {
245
 
                            dif[k][l] = deldif;
246
 
                            velx = v1[k][l];
247
 
                            vely = v2[k][l];
248
 
                        }
249
 
 
250
 
 
251
 
                        if (traps != NULL && trap[k][l] != 0.) {        /* traps */
252
 
 
253
 
                            eff = ulec();       /* random generator */
254
 
 
255
 
                            if (eff <= trap[k][l]) {
256
 
                                velx = -0.1 * v1[k][l]; /* move it slightly back */
257
 
                                vely = -0.1 * v2[k][l];
258
 
                            }
259
 
 
260
 
                        }
261
 
 
262
 
 
263
 
                        gaux = gasdev();
264
 
                        gauy = gasdev();
265
 
 
266
 
                        w[lw][1] += (velx + dif[k][l] * gaux);  /* move the walker */
267
 
                        w[lw][2] += (vely + dif[k][l] * gauy);
268
 
 
269
 
                        if (hhc > hhmax && wdepth == NULL) {
270
 
                            vavg[lw][1] = hbeta * (vavg[lw][1] + v1[k][l]);
271
 
                            vavg[lw][2] = hbeta * (vavg[lw][2] + v2[k][l]);
272
 
                        }
273
 
 
274
 
                        if (w[lw][1] <= xmin || w[lw][2] <= ymin || w[lw][1]
275
 
                            >= xmax || w[lw][2] >= ymax) {
276
 
                            w[lw][3] = 1e-10;   /* eliminate walker if it is out of area */
277
 
                        }
278
 
                        else {
279
 
                            if (wdepth != NULL) {
280
 
                                l = (int)((w[lw][1] + stxm) / stepx) - mx - 1;
281
 
                                k = (int)((w[lw][2] + stym) / stepy) - my - 1;
282
 
                                w[lw][3] *= sigma[k][l];
283
 
                            }
284
 
 
285
 
                        }       /* else */
286
 
                    }           /*DEFined area */
287
 
                    else {
288
 
                        w[lw][3] = 1e-10;       /* eliminate walker if it is out of area */
289
 
                    }
290
 
                }
291
 
 
292
 
/* output walkers commented out */
293
 
                /*        write the walkers which reach the box */
294
 
/*
295
 
                if (w[lw][1] >= xmin && w[lw][2] >= ymin && w[lw][1]
296
 
                    <= xmax && w[lw][2] <= ymax && iflag[lw] == 0) {
297
 
 
298
 
                    lwout = lw / ldemo;
299
 
                    lwout *= ldemo;
300
 
                    if (lwout == lw && (i == miter || i == iter1)) {    
301
 
                        stack[nstack][1] = mixx / conv + w[lw][1] / conv;
302
 
                        stack[nstack][2] = miyy / conv + w[lw][2] / conv;
303
 
                        stack[nstack][3] = w[lw][3];
304
 
 
305
 
                        nstack++;
306
 
*/
307
 
                        /*                  iflag[lw] = 1; */
308
 
/*
309
 
                    }
310
 
                }
311
 
*/
312
 
 
313
 
                /* walker is leaving the box */
314
 
 
315
 
                /*              if (w[lw][1] <= mixx && w[lw][2] <= bymi && w[lw][1] 
316
 
                   >= bxma && w[lw][2] >= byma && iflag[lw] == 1) {
317
 
                   iflag[lw] = 0;
318
 
                   } */
319
 
 
320
 
                /* every iterout iteration write out the selected ldemo walkers */
321
 
                /* and water depth and time */
322
 
 
323
 
 
324
 
            }                   /* lw loop */
325
 
 
326
 
            if (i == iter1) {
327
 
 
328
 
                /* call output for iteration output */
329
 
 
330
 
                if (ts == 1) {
331
 
                    if (erdep != NULL)
332
 
                        erod(gama);     /* divergence of gama field */
333
 
 
334
 
                    conn = (double)nblock / (double)iblock;
335
 
                    itime = (int)(i * deltap * timec);
336
 
                    ii = output_data(itime, conn);
337
 
/*                  nstack = 0; */
338
 
                    if (ii != 1)
339
 
                        G_fatal_error(_("Unable to write raster maps"));
340
 
                }
341
 
 
342
 
            }
343
 
 
344
 
/* ascii data site file output for gamma  - hydrograph or sediment*/
345
 
/* cchez incl. sqrt(sinsl) */
346
 
/* sediment */
347
 
/*defined area */
348
 
/*
349
 
            if (sfile != NULL) {        
350
 
 
351
 
                for (p = 0; p < npoints; p++) {
352
 
 
353
 
                    l = (int)((points[p].east - mixx + stxm) / stepx) - mx -
354
 
                        1;
355
 
                    k = (int)((points[p].north - miyy + stym) / stepy) - my -
356
 
                        1;
357
 
 
358
 
                    if (zz[k][l] != UNDEF) {
359
 
 
360
 
                        if (wdepth == NULL) {
361
 
                            points[p].z1 = step * gama[k][l] * cchez[k][l];     
362
 
                        }
363
 
                        else
364
 
                            points[p].z1 = gama[k][l] * slope[k][l];    
365
 
 
366
 
                        G_debug(2, " k,l,z1 %d %d %f", k, l, points[p].z1);
367
 
 
368
 
                        fprintf(fw, "%f %f %f\n", points[p].east / conv,
369
 
                                points[p].north / conv, points[p].z1);
370
 
                    }           
371
 
 
372
 
                }
373
 
 
374
 
            }
375
 
*/
376
 
 
377
 
        }                       /* miter */
378
 
 
379
 
      L_800:
380
 
        /*        if (iwrib != nblock) {
381
 
           icount = icoub / iwrib;
382
 
 
383
 
           if (icoub == (icount * iwrib)) {
384
 
           ++icfl;
385
 
           nflw = icfl + 50;
386
 
           conn = (double) nblock / (double) iblock;
387
 
 
388
 
           }
389
 
           } */
390
 
 
391
 
 
392
 
        if (err != NULL) {
393
 
            for (k = 0; k < my; k++) {
394
 
                for (l = 0; l < mx; l++) {
395
 
                    if (zz[k][l] != UNDEF) {
396
 
                        d1 = gama[k][l] * (double)conn;
397
 
                        gammas[k][l] += pow(d1, 3. / 5.);
398
 
                    }           /* DEFined area */
399
 
                }
400
 
            }
401
 
        }
402
 
        if (erdep != NULL)
403
 
            erod(gama);
404
 
    }
405
 
    /*                       ........ end of iblock loop */
406
 
 
407
 
}