← Back to branch summary

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

« back to all changes in this revision

Viewing changes to imagery/i.eb.hsebal01/main.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:
imagery/i.eb.hsebal01/main.c
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       i.eb.hSEBAL01
 
5
 * AUTHOR(S):    Yann Chemin - yann.chemin@gmail.com
 
6
 * PURPOSE:      Calculates sensible heat flux by SEBAL iteration
 
7
 *               Delta T will be reassessed in the iterations !
 
8
 *               This has been seen in Bastiaanssen (1995), 
 
9
 *               later modified by Chemin and Alexandridis (2001).
 
10
 *               This code is implemented in Alexandridis et al. (2009).
 
11
 *
 
12
 * COPYRIGHT:    (C) 2002-2009 by the GRASS Development Team
 
13
 *
 
14
 *               This program is free software under the GNU General Public
 
15
 *               License (>=v2). Read the file COPYING that comes with GRASS
 
16
 *               for details.
 
17
 *
 
18
 * CHANGELOG:   
 
19
 *
 
20
 *****************************************************************************/
 
21
 
 
22
#include <stdio.h>
 
23
#include <string.h>
 
24
#include <stdlib.h>
 
25
#include <math.h>
 
26
#include <grass/gis.h>
 
27
#include <grass/raster.h>
 
28
#include <grass/glocale.h>
 
29
 
 
30
double **G_alloc_matrix(int rows, int cols)
 
31
{
 
32
    double **m;
 
33
    int i;
 
34
 
 
35
    m = (double **)G_calloc(rows, sizeof(double *));
 
36
    m[0] = (double *)G_calloc(rows * cols, sizeof(double));
 
37
    for (i = 1; i < rows; i++)
 
38
        m[i] = m[i - 1] + cols;
 
39
 
 
40
    return m;
 
41
}
 
42
 
 
43
int main(int argc, char *argv[])
 
44
{
 
45
    struct Cell_head cellhd;
 
46
 
 
47
    /* buffer for in, tmp and out raster */
 
48
    void *inrast_Rn, *inrast_g0;
 
49
    void *inrast_z0m, *inrast_t0dem;
 
50
    DCELL *outrast;
 
51
    int nrows = 0, ncols = 0;
 
52
    int row = 0, col = 0;
 
53
    int row_wet = 0, col_wet = 0;
 
54
    int row_dry = 0, col_dry = 0;
 
55
    double m_row_wet = 0.0, m_col_wet = 0.0;
 
56
    double m_row_dry = 0.0, m_col_dry = 0.0;
 
57
    int infd_Rn, infd_g0;
 
58
    int infd_z0m, infd_t0dem;
 
59
    int outfd;
 
60
    char *Rn, *g0;
 
61
    char *z0m, *t0dem;
 
62
    char *h0;
 
63
 
 
64
    double ustar, ea;
 
65
    struct History history;
 
66
    struct GModule *module;
 
67
    struct Option *input_Rn, *input_g0;
 
68
    struct Option *input_z0m, *input_t0dem, *input_ustar;
 
69
    struct Option *input_ea, *output;
 
70
    struct Option *input_row_wet, *input_col_wet;
 
71
    struct Option *input_row_dry, *input_col_dry;
 
72
    struct Flag *flag2, *flag3;
 
73
 
 
74
    /********************************/
 
75
    double xp = 0.0, yp = 0.0;
 
76
    double xmin = 0.0, ymin = 0.0;
 
77
    double xmax = 0.0, ymax = 0.0;
 
78
    double stepx = 0.0, stepy = 0.0;
 
79
    double latitude = 0.0, longitude = 0.0;
 
80
    int rowDry = 0, colDry = 0, rowWet = 0, colWet = 0;
 
81
 
 
82
    /********************************/
 
83
 
 
84
    /********************************/
 
85
    xp = yp;
 
86
    yp = xp;
 
87
    xmin = ymin;
 
88
    ymin = xmin;
 
89
    xmax = ymax;
 
90
    ymax = xmax;
 
91
    stepx = stepy;
 
92
    stepy = stepx;
 
93
    latitude = longitude;
 
94
    longitude = latitude;
 
95
    rowDry = colDry;
 
96
    colDry = rowDry;
 
97
    rowWet = colWet;
 
98
    colWet = rowWet;
 
99
 
 
100
    /********************************/
 
101
    G_gisinit(argv[0]);
 
102
 
 
103
    module = G_define_module();
 
104
    G_add_keyword(_("imagery"));
 
105
    G_add_keyword(_("energy balance"));
 
106
    G_add_keyword(_("soil moisture"));
 
107
    G_add_keyword(_("evaporative fraction"));
 
108
    G_add_keyword(_("SEBAL"));
 
109
    module->description =
 
110
        _("Computes sensible heat flux iteration SEBAL 01.");
 
111
 
 
112
    /* Define different options */
 
113
    input_Rn = G_define_standard_option(G_OPT_R_INPUT);
 
114
    input_Rn->key = "netradiation";
 
115
    input_Rn->description =
 
116
        _("Name of instantaneous net radiation raster map [W/m2]");
 
117
 
 
118
    input_g0 = G_define_standard_option(G_OPT_R_INPUT);
 
119
    input_g0->key = "soilheatflux";
 
120
    input_g0->description =
 
121
        _("Name of instantaneous soil heat flux raster map [W/m2]");
 
122
 
 
123
    input_z0m = G_define_standard_option(G_OPT_R_INPUT);
 
124
    input_z0m->key = "aerodynresistance";
 
125
    input_z0m->description =
 
126
        _("Name of aerodynamic resistance to heat momentum raster map [s/m]");
 
127
 
 
128
    input_t0dem = G_define_standard_option(G_OPT_R_INPUT);
 
129
    input_t0dem->key = "temperaturemeansealevel";
 
130
    input_t0dem->description =
 
131
        _("Name of altitude corrected surface temperature raster map [K]");
 
132
 
 
133
    input_ustar = G_define_option();
 
134
    input_ustar->key = "frictionvelocitystar";
 
135
    input_ustar->type = TYPE_DOUBLE;
 
136
    input_ustar->required = YES;
 
137
    input_ustar->gisprompt = "old,value";
 
138
    input_ustar->answer = "0.32407";
 
139
    input_ustar->description =
 
140
        _("Value of the height independent friction velocity (u*) [m/s]");
 
141
    input_ustar->guisection = _("Parameters");
 
142
 
 
143
    input_ea = G_define_option();
 
144
    input_ea->key = "vapourpressureactual";
 
145
    input_ea->type = TYPE_DOUBLE;
 
146
    input_ea->required = YES;
 
147
    input_ea->answer = "1.511";
 
148
    input_ea->description =
 
149
        _("Value of the actual vapour pressure (e_act) [KPa]");
 
150
    input_ea->guisection = _("Parameters");
 
151
 
 
152
    input_row_wet = G_define_option();
 
153
    input_row_wet->key = "row_wet_pixel";
 
154
    input_row_wet->type = TYPE_DOUBLE;
 
155
    input_row_wet->required = NO;
 
156
    input_row_wet->description = _("Row value of the wet pixel");
 
157
    input_row_wet->guisection = _("Parameters");
 
158
 
 
159
    input_col_wet = G_define_option();
 
160
    input_col_wet->key = "column_wet_pixel";
 
161
    input_col_wet->type = TYPE_DOUBLE;
 
162
    input_col_wet->required = NO;
 
163
    input_col_wet->description = _("Column value of the wet pixel");
 
164
    input_col_wet->guisection = _("Parameters");
 
165
 
 
166
    input_row_dry = G_define_option();
 
167
    input_row_dry->key = "row_dry_pixel";
 
168
    input_row_dry->type = TYPE_DOUBLE;
 
169
    input_row_dry->required = NO;
 
170
    input_row_dry->description = _("Row value of the dry pixel");
 
171
    input_row_dry->guisection = _("Parameters");
 
172
 
 
173
    input_col_dry = G_define_option();
 
174
    input_col_dry->key = "column_dry_pixel";
 
175
    input_col_dry->type = TYPE_DOUBLE;
 
176
    input_col_dry->required = NO;
 
177
    input_col_dry->description = _("Column value of the dry pixel");
 
178
    input_col_dry->guisection = _("Parameters");
 
179
 
 
180
    output = G_define_standard_option(G_OPT_R_OUTPUT);
 
181
    output->description =
 
182
        _("Name for output sensible heat flux raster map [W/m2]");
 
183
 
 
184
    /* Define the different flags */
 
185
    flag2 = G_define_flag();
 
186
    flag2->key = 'a';
 
187
    flag2->description = _("Automatic wet/dry pixel (careful!)");
 
188
 
 
189
    flag3 = G_define_flag();
 
190
    flag3->key = 'c';
 
191
    flag3->description =
 
192
        _("Dry/Wet pixels coordinates are in image projection, not row/col");
 
193
 
 
194
    if (G_parser(argc, argv))
 
195
        exit(EXIT_FAILURE);
 
196
 
 
197
    /* get entered parameters */
 
198
    Rn = input_Rn->answer;
 
199
    g0 = input_g0->answer;
 
200
    z0m = input_z0m->answer;
 
201
    t0dem = input_t0dem->answer;
 
202
 
 
203
    h0 = output->answer;
 
204
 
 
205
    ustar = atof(input_ustar->answer);
 
206
    ea = atof(input_ea->answer);
 
207
 
 
208
    /*If automatic flag, just forget the rest of options */
 
209
    if (flag2->answer)
 
210
        G_message(_("Automatic mode selected"));
 
211
    /*If not automatic & all pixels locations in col/row given */
 
212
    else if (!flag2->answer &&
 
213
             input_row_wet->answer &&
 
214
             input_col_wet->answer &&
 
215
             input_row_dry->answer && input_col_dry->answer) {
 
216
        m_row_wet = atof(input_row_wet->answer);
 
217
        m_col_wet = atof(input_col_wet->answer);
 
218
        m_row_dry = atof(input_row_dry->answer);
 
219
        m_col_dry = atof(input_col_dry->answer);
 
220
        /*If pixels locations are in projected coordinates */
 
221
        if (flag3->answer)
 
222
            G_message(_("Manual wet/dry pixels in image coordinates"));
 
223
        G_message(_("Wet Pixel=> x:%f y:%f"), m_col_wet, m_row_wet);
 
224
        G_message(_("Dry Pixel=> x:%f y:%f"), m_col_dry, m_row_dry);
 
225
    }
 
226
    /*If not automatic & missing any of the pixel location, Fatal Error */
 
227
    else {
 
228
        G_fatal_error(_("Either auto-mode either wet/dry pixels coordinates should be provided!"));
 
229
    }
 
230
 
 
231
    /* check legal output name */
 
232
    if (G_legal_filename(h0) < 0)
 
233
        G_fatal_error(_("<%s> is an illegal name"), h0);
 
234
 
 
235
    infd_Rn = Rast_open_old(Rn, "");
 
236
    infd_g0 = Rast_open_old(g0, "");
 
237
    infd_z0m = Rast_open_old(z0m, "");
 
238
    infd_t0dem = Rast_open_old(t0dem, "");
 
239
 
 
240
    Rast_get_cellhd(Rn, "", &cellhd);
 
241
    Rast_get_cellhd(g0, "", &cellhd);
 
242
    Rast_get_cellhd(z0m, "", &cellhd);
 
243
    Rast_get_cellhd(t0dem, "", &cellhd);
 
244
 
 
245
    /* Allocate input buffer */
 
246
    inrast_Rn = Rast_allocate_d_buf();
 
247
    inrast_g0 = Rast_allocate_d_buf();
 
248
    inrast_z0m = Rast_allocate_d_buf();
 
249
    inrast_t0dem = Rast_allocate_d_buf();
 
250
 
 
251
    /***************************************************/
 
252
    /* Setup pixel location variables */
 
253
 
 
254
    /***************************************************/
 
255
    stepx = cellhd.ew_res;
 
256
    stepy = cellhd.ns_res;
 
257
 
 
258
    xmin = cellhd.west;
 
259
    xmax = cellhd.east;
 
260
    ymin = cellhd.south;
 
261
    ymax = cellhd.north;
 
262
 
 
263
    nrows = Rast_window_rows();
 
264
    ncols = Rast_window_cols();
 
265
 
 
266
    /***************************************************/
 
267
    /* Allocate output buffer */
 
268
 
 
269
    /***************************************************/
 
270
    outrast = Rast_allocate_d_buf();
 
271
    outfd = Rast_open_new(h0, DCELL_TYPE);
 
272
 
 
273
    /***************************************************/
 
274
    /* Allocate memory for temporary images            */
 
275
    double **d_Roh, **d_Rah;
 
276
 
 
277
    if ((d_Roh = G_alloc_matrix(nrows, ncols)) == NULL)
 
278
        G_message("Unable to allocate memory for temporary d_Roh image");
 
279
    if ((d_Rah = G_alloc_matrix(nrows, ncols)) == NULL)
 
280
        G_message("Unable to allocate memory for temporary d_Rah image");
 
281
 
 
282
    /***************************************************/
 
283
 
 
284
    /* MANUAL T0DEM WET/DRY PIXELS */
 
285
    DCELL d_Rn_dry = 0.0, d_g0_dry = 0.0;
 
286
    DCELL d_t0dem_dry = 0.0, d_t0dem_wet = 0.0;
 
287
 
 
288
    if (flag2->answer) {
 
289
        /* Process tempk min / max pixels */
 
290
        /* Internal use only */
 
291
        DCELL d_Rn_wet = 0.0, d_g0_wet = 0.0;
 
292
        DCELL d_Rn, d_g0, d_h0;
 
293
        DCELL t0dem_min = 1000.0, t0dem_max = 0.0;
 
294
 
 
295
        /*********************/
 
296
        for (row = 0; row < nrows; row++) {
 
297
            DCELL d_t0dem;
 
298
 
 
299
            G_percent(row, nrows, 2);
 
300
            Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
301
            Rast_get_d_row(infd_Rn, inrast_Rn, row);
 
302
            Rast_get_d_row(infd_g0, inrast_g0, row);
 
303
            /*process the data */
 
304
            for (col = 0; col < ncols; col++) {
 
305
                d_t0dem = ((DCELL *) inrast_t0dem)[col];
 
306
                d_Rn = ((DCELL *) inrast_Rn)[col];
 
307
                d_g0 = ((DCELL *) inrast_g0)[col];
 
308
                if (Rast_is_d_null_value(&d_t0dem) ||
 
309
                    Rast_is_d_null_value(&d_Rn) ||
 
310
                    Rast_is_d_null_value(&d_g0)) {
 
311
                    /* do nothing */
 
312
                }
 
313
                else {
 
314
                    if (d_t0dem <= 250.0) {
 
315
                        /* do nothing */
 
316
                    }
 
317
                    else {
 
318
                        d_h0 = d_Rn - d_g0;
 
319
                        if (d_t0dem < t0dem_min &&
 
320
                            d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 0.0 &&
 
321
                            d_h0 < 100.0) {
 
322
                            t0dem_min = d_t0dem;
 
323
                            d_t0dem_wet = d_t0dem;
 
324
                            d_Rn_wet = d_Rn;
 
325
                            d_g0_wet = d_g0;
 
326
                            m_col_wet = col;
 
327
                            m_row_wet = row;
 
328
                        }
 
329
                        if (d_t0dem > t0dem_max &&
 
330
                            d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 100.0 &&
 
331
                            d_h0 < 500.0) {
 
332
                            t0dem_max = d_t0dem;
 
333
                            d_t0dem_dry = d_t0dem;
 
334
                            d_Rn_dry = d_Rn;
 
335
                            d_g0_dry = d_g0;
 
336
                            m_col_dry = col;
 
337
                            m_row_dry = row;
 
338
                        }
 
339
                    }
 
340
                }
 
341
            }
 
342
        }
 
343
        G_message("row_wet=%d\tcol_wet=%d", row_wet, col_wet);
 
344
        G_message("row_dry=%d\tcol_dry=%d", row_dry, col_dry);
 
345
        G_message("g0_wet=%f", d_g0_wet);
 
346
        G_message("Rn_wet=%f", d_Rn_wet);
 
347
        G_message("LE_wet=%f", d_Rn_wet - d_g0_wet);
 
348
        G_message("t0dem_dry=%f", d_t0dem_dry);
 
349
        G_message("rnet_dry=%f", d_Rn_dry);
 
350
        G_message("g0_dry=%f", d_g0_dry);
 
351
        G_message("h0_dry=%f", d_Rn_dry - d_g0_dry);
 
352
    }                           /* END OF FLAG2 */
 
353
 
 
354
    G_message("Passed here");
 
355
 
 
356
    /* MANUAL T0DEM WET/DRY PIXELS */
 
357
    /*DRY PIXEL */
 
358
    if (flag3->answer) {
 
359
        /*Calculate coordinates of row/col from projected ones */
 
360
        row = (int)((ymax - m_row_dry) / (double)stepy);
 
361
        col = (int)((m_col_dry - xmin) / (double)stepx);
 
362
        G_message("Dry Pixel | row:%i col:%i", row, col);
 
363
    }
 
364
    else {
 
365
        row = (int)m_row_dry;
 
366
        col = (int)m_col_dry;
 
367
        G_message("Dry Pixel | row:%i col:%i", row, col);
 
368
    }
 
369
    rowDry = row;
 
370
    colDry = col;
 
371
    Rast_get_d_row(infd_Rn, inrast_Rn, row);
 
372
    Rast_get_d_row(infd_g0, inrast_g0, row);
 
373
    Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
374
    d_Rn_dry = ((DCELL *) inrast_Rn)[col];
 
375
    d_g0_dry = ((DCELL *) inrast_g0)[col];
 
376
    d_t0dem_dry = ((DCELL *) inrast_t0dem)[col];
 
377
    /*WET PIXEL */
 
378
    if (flag3->answer) {
 
379
        /*Calculate coordinates of row/col from projected ones */
 
380
        row = (int)((ymax - m_row_wet) / (double)stepy);
 
381
        col = (int)((m_col_wet - xmin) / (double)stepx);
 
382
        G_message("Wet Pixel | row:%i col:%i", row, col);
 
383
    }
 
384
    else {
 
385
        row = m_row_wet;
 
386
        col = m_col_wet;
 
387
        G_message("Wet Pixel | row:%i col:%i", row, col);
 
388
    }
 
389
    rowWet = row;
 
390
    colWet = col;
 
391
    Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
392
    d_t0dem_wet = ((DCELL *) inrast_t0dem)[col];
 
393
    /* END OF MANUAL WET/DRY PIXELS */
 
394
    double h_dry;
 
395
 
 
396
    h_dry = d_Rn_dry - d_g0_dry;
 
397
    G_message("h_dry = %f", h_dry);
 
398
    G_message("t0dem_dry = %f", d_t0dem_dry);
 
399
    G_message("t0dem_wet = %f", d_t0dem_wet);
 
400
 
 
401
    DCELL d_rah_dry = 0.0;
 
402
    DCELL d_roh_dry = 0.0;
 
403
 
 
404
    DCELL d_t0dem, d_z0m;
 
405
    DCELL d_u5;
 
406
    DCELL d_roh1;
 
407
    DCELL d_h1, d_h2, d_h3;
 
408
    DCELL d_rah1, d_rah2, d_rah3;
 
409
    DCELL d_L, d_x, d_psih, d_psim;
 
410
 
 
411
    /* INITIALIZATION */
 
412
    for (row = 0; row < nrows; row++) {
 
413
        G_percent(row, nrows, 2);
 
414
        /* read a line input maps into buffers */
 
415
        Rast_get_d_row(infd_z0m, inrast_z0m, row);
 
416
        Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
417
        /* read every cell in the line buffers */
 
418
        for (col = 0; col < ncols; col++) {
 
419
            d_z0m = ((DCELL *) inrast_z0m)[col];
 
420
            d_t0dem = ((DCELL *) inrast_t0dem)[col];
 
421
            if (Rast_is_d_null_value(&d_t0dem) ||
 
422
                Rast_is_d_null_value(&d_z0m)) {
 
423
                /* do nothing */
 
424
                d_Roh[row][col] = -999.9;
 
425
                d_Rah[row][col] = -999.9;
 
426
            }
 
427
            else {
 
428
                d_u5 = (ustar / 0.41) * log(5 / d_z0m);
 
429
                d_rah1 =
 
430
                    (1 / (d_u5 * pow(0.41, 2))) * log(5 / d_z0m) * log(5 /
 
431
                                                                       (d_z0m
 
432
                                                                        *
 
433
                                                                        0.1));
 
434
                d_roh1 =
 
435
                    ((998 - ea) / (d_t0dem * 2.87)) + (ea / (d_t0dem * 4.61));
 
436
                if (d_roh1 > 5)
 
437
                    d_roh1 = 1.0;
 
438
                else
 
439
                    d_roh1 =
 
440
                        ((1000 - 4.65) / (d_t0dem * 2.87)) +
 
441
                        (4.65 / (d_t0dem * 4.61));
 
442
                if (row == rowDry && col == colDry) {   /*collect dry pix info */
 
443
                    d_rah_dry = d_rah1;
 
444
                    d_roh_dry = d_roh1;
 
445
                    G_message("d_rah_dry=%f d_roh_dry=%f", d_rah_dry,
 
446
                              d_roh_dry);
 
447
                }
 
448
                d_Roh[row][col] = d_roh1;
 
449
                d_Rah[row][col] = d_rah1;
 
450
            }
 
451
        }
 
452
    }
 
453
    DCELL d_dT_dry;
 
454
 
 
455
    /*Calculate dT_dry */
 
456
    d_dT_dry = (h_dry * d_rah_dry) / (1004 * d_roh_dry);
 
457
    double a, b;
 
458
 
 
459
    /*Calculate coefficients for next dT equation */
 
460
    /*a = 1.0/ ((d_dT_dry-0.0) / (d_t0dem_dry-d_t0dem_wet)); */
 
461
    /*b = ( a * d_t0dem_wet ) * (-1.0); */
 
462
    double sumx = d_t0dem_wet + d_t0dem_dry;
 
463
    double sumy = d_dT_dry + 0.0;
 
464
    double sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
 
465
    double sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
 
466
 
 
467
    a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
 
468
    b = (sumy - (a * sumx)) / 2.0;
 
469
    G_message("d_dT_dry=%f", d_dT_dry);
 
470
    G_message("dT1=%f * t0dem + (%f)", a, b);
 
471
    DCELL d_h_dry = 0.0;
 
472
 
 
473
    /* ITERATION 1 */
 
474
    for (row = 0; row < nrows; row++) {
 
475
        G_percent(row, nrows, 2);
 
476
        /* read a line input maps into buffers */
 
477
        Rast_get_d_row(infd_z0m, inrast_z0m, row);
 
478
        Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
479
        /* read every cell in the line buffers */
 
480
        for (col = 0; col < ncols; col++) {
 
481
            d_z0m = ((DCELL *) inrast_z0m)[col];
 
482
            d_t0dem = ((DCELL *) inrast_t0dem)[col];
 
483
            d_rah1 = d_Rah[row][col];
 
484
            d_roh1 = d_Roh[row][col];
 
485
            if (Rast_is_d_null_value(&d_t0dem) ||
 
486
                Rast_is_d_null_value(&d_z0m)) {
 
487
                /* do nothing */
 
488
            }
 
489
            else {
 
490
                if (d_rah1 < 1.0)
 
491
                    d_h1 = 0.0;
 
492
                else
 
493
                    d_h1 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah1;
 
494
                d_L =
 
495
                    -1004 * d_roh1 * pow(ustar,
 
496
                                         3) * d_t0dem / (d_h1 * 9.81 * 0.41);
 
497
                d_x = pow((1 - 16 * (5 / d_L)), 0.25);
 
498
                d_psim =
 
499
                    2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
 
500
                    2 * atan(d_x) + 0.5 * M_PI;
 
501
                d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
 
502
                d_u5 = (ustar / 0.41) * log(5 / d_z0m);
 
503
                d_rah2 =
 
504
                    (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m) - d_psim)
 
505
                    * log((5 / (d_z0m * 0.1)) - d_psih);
 
506
                if (row == rowDry && col == colDry) {   /*collect dry pix info */
 
507
                    d_rah_dry = d_rah2;
 
508
                    d_h_dry = d_h1;
 
509
                }
 
510
                d_Rah[row][col] = d_rah1;
 
511
            }
 
512
        }
 
513
    }
 
514
 
 
515
    /*Calculate dT_dry */
 
516
    d_dT_dry = (d_h_dry * d_rah_dry) / (1004 * d_roh_dry);
 
517
    /*Calculate coefficients for next dT equation */
 
518
    /*      a = (d_dT_dry-0)/(d_t0dem_dry-d_t0dem_wet); */
 
519
    /*      b = (-1.0) * ( a * d_t0dem_wet ); */
 
520
    /*      G_message("d_dT_dry=%f",d_dT_dry); */
 
521
    /*      G_message("dT2=%f * t0dem + (%f)", a, b); */
 
522
    sumx = d_t0dem_wet + d_t0dem_dry;
 
523
    sumy = d_dT_dry + 0.0;
 
524
    sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
 
525
    sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
 
526
    a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
 
527
    b = (sumy - (a * sumx)) / 2.0;
 
528
    G_message("d_dT_dry=%f", d_dT_dry);
 
529
    G_message("dT1=%f * t0dem + (%f)", a, b);
 
530
 
 
531
    /* ITERATION 2 */
 
532
 
 
533
    /***************************************************/
 
534
 
 
535
    /***************************************************/
 
536
    for (row = 0; row < nrows; row++) {
 
537
        G_percent(row, nrows, 2);
 
538
        /* read a line input maps into buffers */
 
539
        Rast_get_d_row(infd_z0m, inrast_z0m, row);
 
540
        Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
541
        /* read every cell in the line buffers */
 
542
        for (col = 0; col < ncols; col++) {
 
543
            d_z0m = ((DCELL *) inrast_z0m)[col];
 
544
            d_t0dem = ((DCELL *) inrast_t0dem)[col];
 
545
            d_rah2 = d_Rah[row][col];
 
546
            d_roh1 = d_Roh[row][col];
 
547
            if (Rast_is_d_null_value(&d_t0dem) ||
 
548
                Rast_is_d_null_value(&d_z0m)) {
 
549
                /* do nothing */
 
550
            }
 
551
            else {
 
552
                if (d_rah2 < 1.0) {
 
553
                    d_h2 = 0.0;
 
554
                }
 
555
                else {
 
556
                    d_h2 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah2;
 
557
                }
 
558
                d_L =
 
559
                    -1004 * d_roh1 * pow(ustar,
 
560
                                         3) * d_t0dem / (d_h2 * 9.81 * 0.41);
 
561
                d_x = pow((1 - 16 * (5 / d_L)), 0.25);
 
562
                d_psim = 2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
 
563
                    2 * atan(d_x) + 0.5 * M_PI;
 
564
                d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
 
565
                d_u5 = (ustar / 0.41) * log(5 / d_z0m);
 
566
                d_rah3 =
 
567
                    (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m) -
 
568
                                                      d_psim) * log((5 /
 
569
                                                                     (d_z0m *
 
570
                                                                      0.1)) -
 
571
                                                                    d_psih);
 
572
                if (row == rowDry && col == colDry) {   /*collect dry pix info */
 
573
                    d_rah_dry = d_rah3;
 
574
                    d_h_dry = d_h2;
 
575
                }
 
576
                d_Rah[row][col] = d_rah2;
 
577
            }
 
578
        }
 
579
    }
 
580
 
 
581
    /*Calculate dT_dry */
 
582
    d_dT_dry = (d_h_dry * d_rah_dry) / (1004 * d_roh_dry);
 
583
    /*Calculate coefficients for next dT equation */
 
584
    /*      a = (d_dT_dry-0)/(d_t0dem_dry-d_t0dem_wet); */
 
585
    /*      b = (-1.0) * ( a * d_t0dem_wet ); */
 
586
    /*      G_message("d_dT_dry=%f",d_dT_dry); */
 
587
    /*      G_message("dT3=%f * t0dem + (%f)", a, b); */
 
588
    sumx = d_t0dem_wet + d_t0dem_dry;
 
589
    sumy = d_dT_dry + 0.0;
 
590
    sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
 
591
    sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
 
592
    a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
 
593
    b = (sumy - (a * sumx)) / 2.0;
 
594
    G_message("d_dT_dry=%f", d_dT_dry);
 
595
    G_message("dT1=%f * t0dem + (%f)", a, b);
 
596
 
 
597
    /* ITERATION 3 */
 
598
 
 
599
    /***************************************************/
 
600
 
 
601
    /***************************************************/
 
602
 
 
603
    for (row = 0; row < nrows; row++) {
 
604
        G_percent(row, nrows, 2);
 
605
        /* read a line input maps into buffers */
 
606
        Rast_get_d_row(infd_z0m, inrast_z0m, row);
 
607
        Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 
608
        /* read every cell in the line buffers */
 
609
        for (col = 0; col < ncols; col++) {
 
610
            d_z0m = ((DCELL *) inrast_z0m)[col];
 
611
            d_t0dem = ((DCELL *) inrast_t0dem)[col];
 
612
            d_rah3 = d_Rah[row][col];
 
613
            d_roh1 = d_Roh[row][col];
 
614
            if (Rast_is_d_null_value(&d_t0dem) ||
 
615
                Rast_is_d_null_value(&d_z0m)) {
 
616
                Rast_set_d_null_value(&outrast[col], 1);
 
617
            }
 
618
            else {
 
619
                if (d_rah3 < 1.0) {
 
620
                    d_h3 = 0.0;
 
621
                }
 
622
                else {
 
623
                    d_h3 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah3;
 
624
                }
 
625
                if (d_h3 < 0 && d_h3 > -50) {
 
626
                    d_h3 = 0.0;
 
627
                }
 
628
                if (d_h3 < -50 || d_h3 > 1000) {
 
629
                    Rast_set_d_null_value(&outrast[col], 1);
 
630
                }
 
631
                outrast[col] = d_h3;
 
632
            }
 
633
        }
 
634
        Rast_put_d_row(outfd, outrast);
 
635
    }
 
636
 
 
637
 
 
638
    G_free(inrast_z0m);
 
639
    Rast_close(infd_z0m);
 
640
    G_free(inrast_t0dem);
 
641
    Rast_close(infd_t0dem);
 
642
 
 
643
    G_free(outrast);
 
644
    Rast_close(outfd);
 
645
 
 
646
    /* add command line incantation to history file */
 
647
    Rast_short_history(h0, "raster", &history);
 
648
    Rast_command_history(&history);
 
649
    Rast_write_history(h0, &history);
 
650
 
 
651
    exit(EXIT_SUCCESS);
 
652
}