[GRASS-SVN] r32480 - in grass-addons/gipe: . i.eb.h_SEBAL95
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Aug 3 05:46:40 EDT 2008
Author: ychemin
Date: 2008-08-03 05:46:40 -0400 (Sun, 03 Aug 2008)
New Revision: 32480
Added:
grass-addons/gipe/i.eb.h_SEBAL95/
grass-addons/gipe/i.eb.h_SEBAL95/Makefile
grass-addons/gipe/i.eb.h_SEBAL95/U_0.c
grass-addons/gipe/i.eb.h_SEBAL95/description.html
grass-addons/gipe/i.eb.h_SEBAL95/dtair.c
grass-addons/gipe/i.eb.h_SEBAL95/dtair_desert.c
grass-addons/gipe/i.eb.h_SEBAL95/functions.h
grass-addons/gipe/i.eb.h_SEBAL95/h1.c
grass-addons/gipe/i.eb.h_SEBAL95/h_0.c
grass-addons/gipe/i.eb.h_SEBAL95/main.c
grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c
grass-addons/gipe/i.eb.h_SEBAL95/psi_m.c
grass-addons/gipe/i.eb.h_SEBAL95/rah1.c
grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c
grass-addons/gipe/i.eb.h_SEBAL95/roh_air.c
grass-addons/gipe/i.eb.h_SEBAL95/roh_air_0.c
grass-addons/gipe/i.eb.h_SEBAL95/sensi_h.c
grass-addons/gipe/i.eb.h_SEBAL95/u_star.c
grass-addons/gipe/i.eb.h_SEBAL95/zom_0.c
Log:
Restored clean-up version of i.eb.h_SEBAL95
Added: grass-addons/gipe/i.eb.h_SEBAL95/Makefile
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/Makefile (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/Makefile 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.eb.h_SEBAL95
+
+LIBES = $(GISLIB)
+DEPENDENCIES = $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/gipe/i.eb.h_SEBAL95/U_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/U_0.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/U_0.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,15 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double U_0(double zom_0, double u2m)
+{
+ double u_0;
+
+ u_0 = u2m*0.41*log10(200/(0.15/7))/(log10(2/(0.15/7))*log10(200/zom_0));
+
+// printf("u_0 = %5.3f\n", u_0);
+
+ return (u_0);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/U_0.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/description.html
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/description.html (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/description.html 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,161 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>r.evapo.PM</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass.smlogo.gif" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<H2>NAME</H2> <B><I>i.eb.h_SEBAL95 </I></B>- computation of <i>sensible heat flux</i> [W/m2] after Bastiaanssen, 1995 in [1].
+
+<P><I>(GRASS Raster Program)</I>
+
+<H2>SYNOPSIS</H2>
+<B>i.eb.h_SEBAL95</B>
+<BR>
+<B>i.eb.h_SEBAL95</B> help</br>
+<BR>
+
+<B>i.eb.h_SEBAL95</b> <b>[ -qzn ]</b>
+
+<b>DEM</b>=name
+<b>T</b>=name
+<b>RH</b>=name
+<b>WS</b>=name
+<b>NSR</b>=name
+<b>Vh</b>=name
+<b>ETP</b>=name
+
+<H2>DESCRIPTION</H2>
+
+<p><EM>i.eb.h_SEBAL95</EM> given the vegetation height (hc), humidity (RU),
+wind speed at two meters height (WS), temperature (T), digital terrain model (DEM),
+and net radiation (NSR) raster input maps,
+calculates the sensible heat flux map (h0).
+
+<p>Optionally the user can activate a flag (-z)
+that allows him setting to zero all of the negative evapotranspiration cells;
+in fact these negative values motivated by the condensation of the air water
+vapour content, are sometime undesired because they can produce computational
+problems. The usage of the flag -n detect that the module is run in night hours
+and the appropriate soil heat flux is calculated.
+
+<p>The algorithm implements well known approaches: the hourly Penman-Monteith method as presented in Allen et al. (1998) for land surfaces and the Penman method (Penman, 1948) for water surfaces.<br>
+
+<p>Land and water surfaces are idenfyied by Vh:<br>
+- where Vh>0 vegetation is present and evapotranspiration is calculated;<br>
+- where Vh=0 bare ground is present and evapotranspiration is calculated;<br>
+- where Vh<0 water surface is present and evaporation is calculated;<br>
+
+<p>For more details on the algorithms see [1].
+
+
+<H2>OPTIONS</H2>
+
+The program will run non-interactively if the user specifies program
+arguments and flag settings on the command line using the following
+form:
+
+<P><B>i.eb.h_SEBAL95</b> <b>[ -qzd ]</b>
+<b>DEM</b>=name
+<b>T</b>=name
+<b>RH</b>=name
+<b>WS</b>=name
+<b>NSR</b>=name
+<b>Vh</b>=name
+<b>ETP</b>=name
+
+
+
+<P>Alternatively, the user can simply type <EM>i.eb.h_SEBAL95</EM> on the
+command line and the program will ask for parameter values and flag
+settings interactively, using the standard GRASS parser interface.
+
+
+<H3>Flags:</H3>
+<dl>
+ <dt><B>-q</B>
+ <dd>Run quietly (do not display status messages). By default
+ <EM>r.evapo.PM</EM> is run verbosely.
+ <dt><B>-z</B>
+ <dd>Set negative calculated evapotranspiration values to zero.
+ <dt><B>-n</B>
+ <dd>Calculate soil heat flux for night time. By default
+ <EM>r.evapo.PM</EM> calculate it for day time.
+</dl>
+
+
+<H3>Parameters:</H3>
+<dl>
+ <dt><B>DEM</B>=<I>name</I>
+ <dd>Input elevation raster [m a.s.l.]. Required.</dd>
+
+ <dt><b>T</b>=<I>name</I>
+ <dd>Input temperature raster [�C]. Required.</dd>
+
+ <dt><b>RH</b> =<I>name</I>
+ <dd>Input relative humidity raster [%]. Required.</dd>
+
+ <dt><b>WS</b> =<I>name</I>
+ <dd>Input wind speed at two meters raster [m/s]. Required.</dd>
+
+ <dt><b>NSR</b> =<I>name</I>
+ <dd>Input net solar radiation raster [MJ/(m2*h)]. Required.</dd>
+
+ <dt><b>Vh</b> =<I>name</I>
+ <dd>Input vegetation heigth raster [m]. Required.</dd>
+
+ <dt><b>ETP</b> =<I>name</I>
+ <dd>Output evapotranspiration raster [mm/h]. Required.</dd>
+
+</dl>
+
+
+<H2>NOTES</H2>
+
+<p>Net solar radiation map in MJ/(m2*h) can be computed from the combination of the r.sun ,
+run in mode 1, and the r.mapcalc commands.
+
+<p>The sum of the three radiation components outputted by r.sun (beam, diffuse, and reflected)
+multiplied by the Wh->Mj conversion factor (0.0036) and optionally by a
+clear sky factor [0-1] allows the generation of a map to be used as
+an NSR input for the <EM>r.evapo.PM</EM> command.
+<dt>example:
+<br><dd>r.sun -s elevin=dem aspin=aspect slopein=slope lin=2 albedo=alb_Mar incidout=out beam_rad=beam diff_rad=diffuse refl_rad=reflected day=73 time=13:00 dist=100;
+<br><dd>r.mapcalc 'NSR=0.0036*(beam+diffuse+reflected)';
+
+
+<H2>SEE ALSO</H2>
+<ul>
+ <li><a href=i.eb.h_iter.html>i.eb.h_iter</a>,
+ <a href=i.eb.h0.html>i.eb.h0</a>
+</ul>
+
+
+
+<H2>AUTHORS</H2>
+ <p>
+ <i>
+ <br>Yann Chemin, International Rice Research Institute, Los Banos, The Philippines.
+ </i>
+ <p>Contact: <a href="mailto:y.chemin at cgiar.org"> Yann chemin</a>
+
+
+<H2>REFERENCES</H2>
+
+ <p>[1] Bastiaanssen, W.G.M., 1995.
+ Estimation of Land surface paramters by remote sensing under clear-sky conditions. PhD thesis, Wageningen University, Wageningen, The Netherlands.
+
+ <p>[2] Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998.
+ Crop Evapotranspiration: Guidelines for computing crop water requirements.
+ Irrigation and Drainage Paper 56, Food and Agriculture Organization of the
+ United Nations, Rome, pp. 300
+
+ <p>[3] Penman, H. L. 1948. Natural evaporation from open water,
+ bare soil and grass. Proc. Roy. Soc. London, A193, pp. 120-146.
+<p><i>Last changed: $Date: 2007/07/29 19:30:00 $</i>
+</body>
+</html>
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/description.html
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/dtair.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/dtair.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/dtair.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,20 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+/* Pixel-based input required are: tempk water & desert
+ * additionally, dtair in Desert point should be given
+ */
+
+
+double dt_air(double t0_dem, double tempk_water, double tempk_desert, double dtair_desert)
+{
+ double a, b, result;
+ a = (dtair_desert-0.0)/(tempk_desert-tempk_water);
+ b = dtair_desert - a * tempk_desert;
+ result = t0_dem * a + b;
+// printf("__________________________________");
+// printf("dtair=%5.3f * Tempk + (%5.3f)\n",a,b);
+ return result;
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/dtair.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/dtair_desert.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/dtair_desert.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/dtair_desert.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,14 @@
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include"functions.h"
+
+double dt_air_desert(double h_desert, double roh_air_desert, double rah_desert)
+{
+ double result;
+
+ result = (h_desert * rah_desert)/(roh_air_desert * 1004);
+
+ return result;
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/dtair_desert.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/functions.h
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/functions.h (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/functions.h 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,51 @@
+/* These are the headers for the RS functions used in SEBAL */
+/* ychemin at yahoo.com - yann.chemin at ait.ac.th */
+/* LGPL - 2004 */
+
+/* Initial functions */
+double DN2Rad_landsat7( double Lmin, double LMax, double QCalMax, double QCalmin, int DN );
+double Rad2Ref_landsat7( double radiance, double doy,double sun_elevation, double k_exo );
+double tempk_landsat7( double chan6 );
+double bb_alb_aster( double greenchan, double redchan, double nirchan, double swirchan1, double swirchan2, double swirchan3, double swirchan4, double swirchan5, double swirchan6 );
+double bb_alb_landsat( double bluechan, double greenchan, double redchan, double nirchan, double chan5, double chan7 );
+double bb_alb_noaa( double redchan, double nirchan );
+double bb_alb_modis( double redchan, double nirchan, double chan3, double chan4, double chan5, double chan6, double chan7 );
+double nd_vi( double redchan, double nirchan );
+double emissivity_generic( double ndvi );
+double emissivity_modis( double e31, double e32 );
+double t0_dem( double dem, double tempk );
+double t_air_approx( double temp_k );
+
+/* Instantaneous rnet and g0 */
+double r_net( double bbalb, double ndvi, double tempk, double tair, double e0, double tsw, double doy, double utc, double sunzangle );
+double g_0(double bbalb, double ndvi, double tempk, double rnet, double utc);
+
+/* Diurnal r_net and etpot */
+double solar_day(double lat, double doy, double tsw );
+double r_net_day( double bbalb, double solar, double tsw );
+double et_pot_day( double bbalb, double solar, double tempk, double tsw );
+
+/* Sensible heat flux functions */
+double sensi_h( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u2m, double dem_desert);
+double roh_air_0(double tempk);
+double zom_0(double ndvi, double ndvi_max);
+double U_0(double zom_0, double u2m);
+double rah_0(double zom_0, double u_0);
+double h_0(double roh_air, double rah, double dtair);
+double dt_air_approx( double temp_k );
+double dt_air_0(double t0_dem, double tempk_water, double tempk_desert);
+double dt_air_desert(double h_desert, double roh_air_desert, double rah_desert);
+double dt_air(double t0_dem, double tempk_water, double tempk_desert, double dtair_desert);
+double rohair(double dem, double tempk, double dtair);
+double h1(double roh_air, double rah, double dtair);
+double u_star(double t0_dem, double h, double ustar, double roh_air, double zom, double u2m);
+double psi_h(double t0_dem, double h, double U_0, double roh_air);
+double rah1(double psih, double u_star);
+
+/* Final outputs */
+double evap_fr(double r_net, double g0, double h);
+double et_a(double r_net_day, double evap_f, double tempk);
+double biomass( double ndvi, double solar_day, double evap_fr, double light_use_ef );
+double soilmoisture( double ndvi, double evap_fr );
+
+
Added: grass-addons/gipe/i.eb.h_SEBAL95/h1.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/h1.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/h1.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,13 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double h1(double roh_air, double rah, double dtair)
+{
+ double result;
+
+ result = roh_air*1004*dtair / rah;
+
+ return result;
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/h1.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/h_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/h_0.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/h_0.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,13 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double h_0(double roh_air, double rah, double dtair)
+{
+ double result;
+
+ result = roh_air*1004*(dtair) / rah;
+
+ return (result);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/h_0.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/main.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/main.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/main.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,991 @@
+/****************************************************************************
+ *
+ * MODULE: i.eb.h_SEBAL95
+ * AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
+ * PURPOSE: Calculates sensible heat flux by SEBAL iteration
+ * Delta T will be reassessed in the iterations !
+ * This has been seen in Bastiaanssen (1995).
+ *
+ * COPYRIGHT: (C) 2002-2009 by the GRASS Development Team
+ *
+ * This program is free software under the GNU General Public
+ * License (>=v2). Read the file COPYING that comes with GRASS
+ * for details.
+ *
+ * CHANGELOG:
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <math.h>
+#include <grass/glocale.h>
+/*#include <omp.h>*/
+
+
+double sensi_h( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u2m, double dem_desert);
+
+int main(int argc, char *argv[])
+{
+ struct Cell_head cellhd;
+ char *mapset; // mapset name
+
+ /* buffer for in out raster */
+ void *inrast_T,*inrast_ndvi,*inrast_dem, *inrast_u2m;
+ void *inrast_Rn,*inrast_g0,*inrast_albedo;
+ DCELL *outrast;
+
+ int nrows, ncols;
+ int row, col;
+ int row_wet, col_wet;
+ int row_dry, col_dry;
+ double m_row_wet, m_col_wet;
+ double m_row_dry, m_col_dry;
+ int infd_T,infd_ndvi,infd_dem,infd_u2m;
+ int infd_Rn,infd_g0,infd_albedo;
+ int outfd;
+
+ char *mapset_T,*mapset_ndvi,*mapset_dem,*mapset_u2m;
+ char *mapset_Rn,*mapset_g0,*mapset_albedo;
+ char *T, *ndvi, *dem, *u2m, *Rn, *g0, *albedo;
+ char *h0;
+
+ struct History history;
+ struct GModule *module;
+ struct Option *input_T, *input_ndvi, *input_dem, *input_u2m;
+ struct Option *input_Rn, *input_g0, *input_albedo;
+ struct Option *output;
+ struct Option *input_row_wet, *input_col_wet;
+ struct Option *input_row_dry, *input_col_dry;
+ struct Option *input_iter;
+ struct Flag *flag1, *flag2, *flag3, *day, *zero;
+ /*******************************/
+ RASTER_MAP_TYPE data_type_T;
+ RASTER_MAP_TYPE data_type_ndvi;
+ RASTER_MAP_TYPE data_type_dem;
+ RASTER_MAP_TYPE data_type_u2m;
+ RASTER_MAP_TYPE data_type_Rn;
+ RASTER_MAP_TYPE data_type_g0;
+ RASTER_MAP_TYPE data_type_albedo;
+ RASTER_MAP_TYPE data_type_output=DCELL_TYPE;
+ /*******************************/
+ int iteration = 10; /*SEBAL95 loop number*/
+ /********************************/
+ /* Stats for dry/wet pixels */
+ double t0dem_min=400.0,t0dem_max=200.0;
+ double tempk_min=400.0,tempk_max=200.0;
+ /********************************/
+ double xp, yp;
+ double xmin, ymin;
+ double xmax, ymax;
+ double stepx,stepy;
+ double latitude, longitude;
+ /********************************/
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ module->description = _("Sensible Heat Flux iteration SEBAL 95");
+
+ /* Define different options */
+ input_T = G_define_standard_option(G_OPT_R_INPUT);
+ input_T->key = "tempk";
+ input_T->description = _("Name of Surface Skin Temperature input map [K]");
+
+ input_dem = G_define_standard_option(G_OPT_R_INPUT);
+ input_dem->key = "dem";
+ input_dem->description = _("Name of dem input map [m a.s.l.]");
+
+ input_u2m = G_define_standard_option(G_OPT_R_INPUT);
+ input_u2m->key = "u2m";
+ input_u2m->description = _("Name of Wind speed at 2m height input map [m/s]");
+
+ input_ndvi = G_define_standard_option(G_OPT_R_INPUT);
+ input_ndvi->key = "ndvi";
+ input_ndvi->description = _("Name of NDVI input map [-]");
+
+ input_albedo = G_define_standard_option(G_OPT_R_INPUT);
+ input_albedo->key= "albedo";
+ input_albedo->description = _("Name of Albedo input map [-]");
+
+ input_Rn = G_define_standard_option(G_OPT_R_INPUT);
+ input_Rn->key = "rnet";
+ input_Rn->description = _("Name of instantaneous Net Solar Radiation input map [W/m2]");
+
+ input_g0 = G_define_standard_option(G_OPT_R_INPUT);
+ input_g0->key = "g0";
+ input_g0->description = _("Name of instantaneous Soil Heat Flux input map [W/m2]");
+
+ input_iter = G_define_option();
+ input_iter->key = "iteration";
+ input_iter->type = TYPE_INTEGER;
+ input_iter->required = NO;
+ input_iter->gisprompt = "old,value";
+ input_iter->description = _("Value of the number of SEBAL95 loops (default is 10)");
+ input_iter->guisection = _("Optional");
+
+ input_row_wet = G_define_option();
+ input_row_wet->key = "row_wet";
+ input_row_wet->type = TYPE_INTEGER;
+ input_row_wet->required = NO;
+ input_row_wet->gisprompt = "old,value";
+ input_row_wet->description = _("Row value of the wet pixel");
+ input_row_wet->guisection = _("Optional");
+
+ input_col_wet = G_define_option();
+ input_col_wet->key = "col_wet";
+ input_col_wet->type = TYPE_INTEGER;
+ input_col_wet->required = NO;
+ input_col_wet->gisprompt = "old,value";
+ input_col_wet->description = _("Column value of the wet pixel");
+ input_col_wet->guisection = _("Optional");
+
+ input_row_dry = G_define_option();
+ input_row_dry->key = "row_dry";
+ input_row_dry->type = TYPE_INTEGER;
+ input_row_dry->required = NO;
+ input_row_dry->gisprompt = "old,value";
+ input_row_dry->description = _("Row value of the dry pixel");
+ input_row_dry->guisection = _("Optional");
+
+ input_col_dry = G_define_option();
+ input_col_dry->key = "col_dry";
+ input_col_dry->type = TYPE_INTEGER;
+ input_col_dry->required = NO;
+ input_col_dry->gisprompt = "old,value";
+ input_col_dry->description = _("Column value of the dry pixel");
+ input_col_dry->guisection = _("Optional");
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT) ;
+ output->key = "h0";
+ output->description= _("Name of output sensible heat flux layer [W/m2]");
+
+ /* Define the different flags */
+ flag1 = G_define_flag() ;
+ flag1->key = 't' ;
+ flag1->description = _("Temperature histogram check (careful!)") ;
+
+ flag2 = G_define_flag() ;
+ flag2->key = 'a' ;
+ flag2->description = _("Automatic wet/dry pixel (careful!)") ;
+
+ flag3 = G_define_flag() ;
+ flag3->key = 'c' ;
+ flag3->description = _("Coordinates of manual dry/wet pixels are in image projection and not row/col") ;
+
+ zero = G_define_flag() ;
+ zero->key = 'z' ;
+ zero->description = _("set negative evapo to zero");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* get entered parameters */
+ T = input_T->answer;
+ dem = input_dem->answer;
+ u2m = input_u2m->answer;
+ ndvi = input_ndvi->answer;
+ Rn = input_Rn->answer;
+ g0 = input_g0->answer;
+ albedo = input_albedo->answer;
+
+ h0 = output->answer;
+
+ if(input_iter->answer){
+ iteration=atoi(input_iter->answer);
+ }
+ if(input_row_wet->answer&&input_row_dry&&
+ input_col_wet->answer&&input_col_dry){
+ m_row_wet = atof(input_row_wet->answer);
+ m_col_wet = atof(input_col_wet->answer);
+ m_row_dry = atof(input_row_dry->answer);
+ m_col_dry = atof(input_col_dry->answer);
+ if(flag3->answer){
+ G_message("Manual wet/dry pixels in image coordinates");
+ G_message("Wet Pixel=> x:%f y:%f",m_col_wet,m_row_wet);
+ G_message("Dry Pixel=> x:%f y:%f",m_col_dry,m_row_dry);
+ } else {
+ G_message("Wet Pixel=> row:%.0f col:%.0f",m_row_wet,m_col_wet);
+ G_message("Dry Pixel=> row:%.0f col:%.0f",m_row_dry,m_col_dry);
+ }
+ }
+ /* find maps in mapset */
+ mapset_T = G_find_cell2 (T, "");
+ if (mapset_T == NULL)
+ G_fatal_error (_("cell file [%s] not found"), T);
+ mapset_dem = G_find_cell2 (dem, "");
+ if (mapset_dem == NULL)
+ G_fatal_error (_("cell file [%s] not found"), dem);
+ mapset_u2m = G_find_cell2 (u2m, "");
+ if (mapset_u2m == NULL)
+ G_fatal_error (_("cell file [%s] not found"), u2m);
+ mapset_ndvi = G_find_cell2 (ndvi, "");
+ if (mapset_ndvi == NULL)
+ G_fatal_error (_("cell file [%s] not found"), ndvi);
+ mapset_Rn = G_find_cell2 (Rn, "");
+ if (mapset_Rn == NULL)
+ G_fatal_error (_("cell file [%s] not found"), Rn);
+ mapset_g0 = G_find_cell2 (g0, "");
+ if (mapset_g0 == NULL)
+ G_fatal_error (_("cell file [%s] not found"), g0);
+ mapset_albedo = G_find_cell2 (albedo, "");
+ if (mapset_albedo == NULL)
+ G_fatal_error(_("cell file [%s] not found"),albedo);
+
+ /* check legal output name */
+ if (G_legal_filename (h0) < 0)
+ G_fatal_error (_("[%s] is an illegal name"), h0);
+
+ /* determine the input map type (CELL/FCELL/DCELL) */
+ data_type_T = G_raster_map_type(T, mapset_T);
+ data_type_dem = G_raster_map_type(dem, mapset_dem);
+ data_type_u2m = G_raster_map_type(u2m, mapset_u2m);
+ data_type_ndvi = G_raster_map_type(ndvi, mapset_ndvi);
+ data_type_Rn = G_raster_map_type(Rn, mapset_Rn);
+ data_type_g0 = G_raster_map_type(g0, mapset_g0);
+ data_type_albedo = G_raster_map_type(albedo, mapset_albedo);
+
+ if ( (infd_T = G_open_cell_old (T, mapset_T)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"), T);
+ if ( (infd_dem = G_open_cell_old (dem, mapset_dem)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"),dem);
+ if ( (infd_u2m = G_open_cell_old (u2m, mapset_u2m)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"),u2m);
+ if ( (infd_ndvi = G_open_cell_old (ndvi, mapset_ndvi)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"),ndvi);
+ if ( (infd_Rn = G_open_cell_old (Rn, mapset_Rn)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"),Rn);
+ if ( (infd_g0 = G_open_cell_old (g0, mapset_g0)) < 0)
+ G_fatal_error (_("Cannot open cell file [%s]"),g0);
+ if((infd_albedo=G_open_cell_old (albedo,mapset_albedo)) < 0)
+ G_fatal_error(_("Cannot open cell file [%s]"),albedo);
+
+ if (G_get_cellhd (T, mapset_T, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), T);
+ if (G_get_cellhd (dem, mapset_dem, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), dem);
+ if (G_get_cellhd (u2m, mapset_u2m, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), u2m);
+ if (G_get_cellhd (ndvi, mapset_ndvi, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), ndvi);
+ if (G_get_cellhd (Rn, mapset_Rn, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), Rn);
+ if (G_get_cellhd (g0, mapset_g0, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), g0);
+ if (G_get_cellhd (albedo, mapset_albedo, &cellhd) < 0)
+ G_fatal_error (_("Cannot read file header of [%s]"), albedo);
+
+ /* Allocate input buffer */
+ inrast_T = G_allocate_raster_buf(data_type_T);
+ inrast_dem = G_allocate_raster_buf(data_type_dem);
+ inrast_u2m = G_allocate_raster_buf(data_type_u2m);
+ inrast_ndvi = G_allocate_raster_buf(data_type_ndvi);
+ inrast_Rn = G_allocate_raster_buf(data_type_Rn);
+ inrast_g0 = G_allocate_raster_buf(data_type_g0);
+ inrast_albedo = G_allocate_raster_buf(data_type_albedo);
+
+ /***************************************************/
+ G_debug(3, "number of rows %d",cellhd.rows);
+ /***************************************************/
+ /* Setup pixel location variables */
+ /***************************************************/
+ stepx=cellhd.ew_res;
+ stepy=cellhd.ns_res;
+
+ xmin=cellhd.west;
+ xmax=cellhd.east;
+ ymin=cellhd.south;
+ ymax=cellhd.north;
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+ /***************************************************/
+ /* Allocate output buffer */
+ /***************************************************/
+ outrast = G_allocate_d_raster_buf();
+
+ if((outfd = G_open_raster_new (h0,DCELL_TYPE)) < 0)
+ G_fatal_error (_("Could not open <%s>"),h0);
+
+ DCELL d_ndvi; /* Input raster */
+ DCELL d_ndvi_max = 0.0; /* Generated here */
+ /* THREAD 1 */
+ /* NDVI Max */
+ for (row = 0; row < nrows; row++)
+ {
+ G_percent (row, nrows, 2);
+ if (G_get_d_raster_row (infd_ndvi, inrast_ndvi, row) < 0)
+ G_fatal_error (_("Could not read from <%s>"),ndvi);
+ /* #pragma omp for parallel default (shared) \
+ shared(ncols,nrows) \
+ private (col,row,d_ndvi)*/
+ for (col=0; col < ncols; col++)
+ {
+ switch(data_type_ndvi){
+ case CELL_TYPE:
+ d_ndvi = (double) ((CELL *) inrast_ndvi)[col];
+ break;
+ case FCELL_TYPE:
+ d_ndvi = (double) ((FCELL *) inrast_ndvi)[col];
+ break;
+ case DCELL_TYPE:
+ d_ndvi = (double) ((DCELL *) inrast_ndvi)[col];
+ break;
+ }
+ //d_ndvi = ((DCELL *) inrast_ndvi)[col];
+ if(G_is_d_null_value(&d_ndvi)){
+ /* do nothing */
+ } else if ((d_ndvi)>d_ndvi_max&&(d_ndvi)<0.999){
+ d_ndvi_max = d_ndvi;
+ }
+ }
+ }
+ G_message("ndvi_max=%f\n",d_ndvi_max);
+ /* Pick up wet and dry pixel values */
+ DCELL d_Rn; /* Input raster */
+ DCELL d_g0; /* Input raster */
+ DCELL d_tempk_wet;
+ DCELL d_tempk_dry;
+ DCELL d_Rn_dry;
+ DCELL d_g0_dry;
+ DCELL d_t0dem_dry;
+ DCELL d_dem_dry;
+ DCELL d_dT_dry;
+ DCELL d_ndvi_dry;
+ DCELL d_dT;
+ /*START Temperature minimum search */
+ /* THREAD 1 */
+ /*This is correcting for un-Earthly temperatures*/
+ /*It finds when histogram is actually starting to pull up...*/
+ int peak1, peak2, peak3;
+ int i_peak1, i_peak2, i_peak3;
+ int bottom1a, bottom1b;
+ int bottom2a, bottom2b;
+ int bottom3a, bottom3b;
+ int i_bottom1a, i_bottom1b;
+ int i_bottom2a, i_bottom2b;
+ int i_bottom3a, i_bottom3b;
+ if(flag1->answer){
+ int i=0;
+ int histogram[400];
+ for (i=0;i<400;i++){
+ histogram[i]=0;
+ }
+ DCELL d_T;
+ /****************************/
+ /* Process pixels histogram */
+ for (row = 0; row < nrows; row++)
+ {
+ G_percent (row, nrows, 2);
+ /* read input map */
+ if( (G_get_raster_row(infd_T,inrast_T,row,data_type_T)) < 0){
+ G_fatal_error (_("Could not read from <%s>"),T);
+ }
+ /*process the data */
+ for (col=0; col < ncols; col++)
+ {
+ switch(data_type_T)
+ {
+ case CELL_TYPE:
+ d_T = (double) ((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_T = (double) ((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_T = (double) ((DCELL *) inrast_T)[col];
+ break;
+ }
+ if(G_is_d_null_value(&d_T)){
+ /*Do nothing*/
+ } else {
+ int temp;
+ temp = (int) d_T;
+ if(temp>0){
+ histogram[temp]=histogram[temp]+1.0;
+ }
+ }
+ }
+ }
+ int sum=0;
+ double average=0.0;
+ for(i=0;i<400;i++){
+ sum += histogram[i];
+ }
+ average = (double) sum / 400.0;
+ G_message("Histogram of Temperature map (if it has rogue values to clean)");
+ //int peak1, peak2, peak3;
+ //int i_peak1, i_peak2, i_peak3;
+ peak1=0;
+ peak2=0;
+ peak3=0;
+ i_peak1=0;
+ i_peak2=0;
+ i_peak3=0;
+ //int bottom1a, bottom1b;
+ //int bottom2a, bottom2b;
+ //int bottom3a, bottom3b;
+ //int i_bottom1a, i_bottom1b;
+ //int i_bottom2a, i_bottom2b;
+ //int i_bottom3a, i_bottom3b;
+ bottom1a=100000;
+ bottom1b=100000;
+ bottom2a=100000;
+ bottom2b=100000;
+ bottom3a=100000;
+ bottom3b=100000;
+ i_bottom1a=1000;
+ i_bottom1b=1000;
+ i_bottom2a=1000;
+ i_bottom2b=1000;
+ i_bottom3a=1000;
+ i_bottom3b=1000;
+ for(i=0;i<400;i++){
+ /* Search for highest peak of dataset (2) */
+ /* Highest Peak */
+ if(histogram[i]>peak2){
+ peak2 = histogram[i];
+ i_peak2=i;
+ }
+ }
+ int stop=0;
+ for(i=i_peak2;i>5;i--){
+ if(((histogram[i]+histogram[i-1]+histogram[i-2]+histogram[i-3]+histogram[i-4])/5)<histogram[i]&&stop==0){
+ bottom2a = histogram[i];
+ i_bottom2a = i;
+ } else if(((histogram[i]+histogram[i-1]+histogram[i-2]+histogram[i-3]+histogram[i-4])/5)>histogram[i]&&stop==0){
+ /*Search for peaks of datasets (1)*/
+ peak1 = histogram[i];
+ i_peak1=i;
+ stop=1;
+ }
+ }
+ stop=0;
+ for(i=i_peak2;i<395;i++){
+ if(((histogram[i]+histogram[i+1]+histogram[i+2]+histogram[i+3]+histogram[i+4])/5)<histogram[i]&&stop==0){
+ bottom2b = histogram[i];
+ i_bottom2b = i;
+ } else if(((histogram[i]+histogram[i+1]+histogram[i+2]+histogram[i+3]+histogram[i+4])/5)>histogram[i]&&stop==0){
+ /*Search for peaks of datasets (3)*/
+ peak3 = histogram[i];
+ i_peak3=i;
+ stop=1;
+ }
+ }
+ /* First histogram lower bound */
+ for(i=250;i<i_peak1;i++){
+ if(histogram[i]<bottom1a){
+ bottom1a = histogram[i];
+ i_bottom1a = i;
+ }
+ }
+ /* First histogram higher bound */
+ for(i=i_peak2;i>i_peak1;i--){
+ if(histogram[i]<=bottom1b){
+ bottom1b = histogram[i];
+ i_bottom1b = i;
+ }
+ }
+ /* Third histogram lower bound */
+ for(i=i_peak2;i<i_peak3;i++){
+ if(histogram[i]<bottom3a){
+ bottom3a = histogram[i];
+ i_bottom3a = i;
+ }
+ }
+ /* Third histogram higher bound */
+ for(i=399;i>i_peak3;i--){
+ if(histogram[i]<bottom3b){
+ bottom3b = histogram[i];
+ i_bottom3b = i;
+ }
+ }
+ G_message("bottom1a: [%i]=>%i\n",i_bottom1a, bottom1a);
+ G_message("peak1: [%i]=>%i\n",i_peak1, peak1);
+ G_message("bottom1b: [%i]=>%i\n",i_bottom1b, bottom1b);
+ G_message("bottom2a: [%i]=>%i\n",i_bottom2a, bottom2a);
+ G_message("peak2: [%i]=>%i\n",i_peak2, peak2);
+ G_message("bottom2b: [%i]=>%i\n",i_bottom2b, bottom2b);
+ G_message("bottom3a: [%i]=>%i\n",i_bottom3a, bottom3a);
+ G_message("peak3: [%i]=>%i\n",i_peak3, peak3);
+ G_message("bottom3b: [%i]=>%i\n",i_bottom3b, bottom3b);
+ }//END OF FLAG1
+
+ /* End of processing histogram*/
+ /*******************/
+
+ /*Process wet pixel values*/
+ /* FLAG2 */
+ if(flag2->answer){
+ /* THREAD 3 */
+ /* Process tempk min / max pixels */
+ /* #pragma omp for parallel default (shared) \
+ shared(ncols,nrows) \
+ private (col,row,d_albedo,d_tempk,d_dem,d_t0dem,d_Rn,d_g0)*/
+ /* Internal use only */
+ DCELL d_Rn_wet;
+ DCELL d_g0_wet;
+ /*********************/
+ for (row = 0; row < nrows; row++){
+ DCELL d_albedo;
+ DCELL d_tempk;
+ DCELL d_dem;
+ DCELL d_t0dem;
+ DCELL d_h0=100.0;//for flag 1
+ DCELL d_h0_max=100.0;//for flag 1
+ G_percent(row,nrows,2);
+ if(G_get_raster_row(infd_albedo,inrast_albedo,row,data_type_albedo)<0)
+ G_fatal_error(_("Could not read from <%s>"),albedo);
+ if(G_get_raster_row(infd_T,inrast_T,row,data_type_T)<0)
+ G_fatal_error(_("Could not read from <%s>"),T);
+ if(G_get_raster_row(infd_dem,inrast_dem,row,data_type_dem)<0)
+ G_fatal_error(_("Could not read from <%s>"),dem);
+ if(G_get_raster_row(infd_Rn, inrast_Rn, row,data_type_Rn) < 0)
+ G_fatal_error(_("Could not read from <%s>"),Rn);
+ if(G_get_raster_row(infd_g0, inrast_g0, row,data_type_g0) < 0)
+ G_fatal_error(_("Could not read from <%s>"),g0);
+ if (G_get_d_raster_row (infd_ndvi, inrast_ndvi, row) < 0)
+ G_fatal_error (_("Could not read from <%s>"),ndvi);
+ /*process the data */
+ for (col=0; col < ncols; col++)
+ {
+ switch(data_type_albedo){
+ case CELL_TYPE:
+ d_albedo = (double) ((CELL *) inrast_albedo)[col];
+ break;
+ case FCELL_TYPE:
+ d_albedo = (double) ((FCELL *) inrast_albedo)[col];
+ break;
+ case DCELL_TYPE:
+ d_albedo = (double) ((DCELL *) inrast_albedo)[col];
+ break;
+ }
+ switch(data_type_T){
+ case CELL_TYPE:
+ d_tempk = (double) ((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double) ((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double) ((DCELL *) inrast_T)[col];
+ break;
+ }
+ switch(data_type_dem){
+ case CELL_TYPE:
+ d_dem = (double) ((CELL *) inrast_dem)[col];
+ break;
+ case FCELL_TYPE:
+ d_dem = (double) ((FCELL *) inrast_dem)[col];
+ break;
+ case DCELL_TYPE:
+ d_dem = (double) ((DCELL *) inrast_dem)[col];
+ break;
+ }
+ switch(data_type_Rn){
+ case CELL_TYPE:
+ d_Rn = (double) ((CELL *) inrast_Rn)[col];
+ break;
+ case FCELL_TYPE:
+ d_Rn = (double) ((FCELL *) inrast_Rn)[col];
+ break;
+ case DCELL_TYPE:
+ d_Rn = (double) ((DCELL *) inrast_Rn)[col];
+ break;
+ }
+ switch(data_type_g0){
+ case CELL_TYPE:
+ d_g0 = (double) ((CELL *) inrast_g0)[col];
+ break;
+ case FCELL_TYPE:
+ d_g0 = (double) ((FCELL *) inrast_g0)[col];
+ break;
+ case DCELL_TYPE:
+ d_g0 = (double) ((DCELL *) inrast_g0)[col];
+ break;
+ }
+ switch(data_type_ndvi){
+ case CELL_TYPE:
+ d_ndvi = (double) ((CELL *) inrast_ndvi)[col];
+ break;
+ case FCELL_TYPE:
+ d_ndvi = (double) ((FCELL *) inrast_ndvi)[col];
+ break;
+ case DCELL_TYPE:
+ d_ndvi = (double) ((DCELL *) inrast_ndvi)[col];
+ break;
+ }
+ if(G_is_d_null_value(&d_albedo)||
+ G_is_d_null_value(&d_tempk)||
+ G_is_d_null_value(&d_dem)||
+ G_is_d_null_value(&d_Rn)||
+ G_is_d_null_value(&d_g0)||
+ G_is_d_null_value(&d_ndvi)){
+ /* do nothing */
+ }else{
+ d_t0dem = d_tempk + 0.00627*d_dem;
+ if(d_t0dem<=250.0||d_tempk<=250.0){
+ /* do nothing */
+ } else {
+ d_h0=d_Rn-d_g0;
+ if(d_t0dem<t0dem_min&&d_albedo<0.15){
+ t0dem_min=d_t0dem;
+ tempk_min=d_tempk;
+ d_tempk_wet=d_tempk;
+ d_Rn_wet=d_Rn;
+ d_g0_wet=d_g0;
+ col_wet=col;
+ row_wet=row;
+ }
+ if(flag1->answer&&
+ d_tempk>=(double)i_peak1-5.0&&
+ d_tempk<(double)i_peak1+1.0){
+ tempk_min=d_tempk;
+ d_tempk_wet=d_tempk;
+ d_Rn_wet=d_Rn;
+ d_g0_wet=d_g0;
+ col_wet=col;
+ row_wet=row;
+ }
+ if(d_t0dem>t0dem_max){
+ t0dem_max=d_t0dem;
+ d_t0dem_dry=d_t0dem;
+ tempk_max=d_tempk;
+ d_tempk_dry=d_tempk;
+ d_Rn_dry=d_Rn;
+ d_g0_dry=d_g0;
+ d_dem_dry=d_dem;
+ col_dry=col;
+ row_dry=row;
+ d_ndvi_dry=d_ndvi;
+ }
+ if(flag1->answer&&
+ d_tempk>=(double)i_peak3-0.0&&
+ d_tempk<(double)i_peak3+7.0&&
+ d_h0>100.0&&d_h0>d_h0_max&&
+ d_g0>10.0&&d_Rn>100.0&&
+ d_albedo>0.3){
+ tempk_max=d_tempk;
+ d_tempk_dry=d_tempk;
+ d_Rn_dry=d_Rn;
+ d_g0_dry=d_g0;
+ d_h0_max=d_h0;
+ d_dem_dry=d_dem;
+ col_dry=col;
+ row_dry=row;
+ d_ndvi_dry=d_ndvi;
+ }
+ }
+ }
+ }
+ }
+ G_message("tempk_min=%f\ntempk_max=%f\n",tempk_min,tempk_max);
+ G_message("row_wet=%d\tcol_wet=%d\n",row_wet,col_wet);
+ G_message("row_dry=%d\tcol_dry=%d\n",row_dry,col_dry);
+ G_message("tempk_wet=%f\n",d_tempk_wet);
+ G_message("g0_wet=%f\n",d_g0_wet);
+ G_message("Rn_wet=%f\n",d_Rn_wet);
+ G_message("LE_wet=%f\n",d_Rn_wet-d_g0_wet);
+ G_message("tempk_dry=%f\n",d_tempk_dry);
+ G_message("dem_dry=%f\n",d_dem_dry);
+ G_message("t0dem_dry=%f\n",d_t0dem_dry);
+ G_message("rnet_dry=%f\n",d_Rn_dry);
+ G_message("g0_dry=%f\n",d_g0_dry);
+ G_message("h0_dry=%f\n",d_Rn_dry-d_g0_dry);
+ } /* END OF FLAG2 */
+ /*MPI_BARRIER*/
+
+ /* MANUAL WET/DRY PIXELS */
+ if(input_row_wet->answer&&input_row_dry->answer&&
+ input_col_wet->answer&&input_col_dry->answer){
+ /*DRY PIXEL*/
+ if(flag3->answer){
+ /*Calculate coordinates of row/col from projected ones*/
+ row = (int) (( ymax - m_row_dry ) / (double) stepy) ;
+ col = (int) (( m_col_dry - xmin ) / (double) stepx) ;
+ G_message("Dry Pixel | row:%i col:%i",row,col);
+ m_row_dry = row;
+ m_col_dry = col;
+ }
+ row = (int) m_row_dry;
+ col = (int) m_col_dry;
+ G_message("Dry Pixel | row:%i col:%i",row,col);
+ DCELL d_tempk;
+ DCELL d_dem;
+ DCELL d_t0dem;
+ if(G_get_raster_row(infd_T,inrast_T,row,data_type_T)<0)
+ G_fatal_error(_("Could not read from <%s>"),T);
+ if(G_get_raster_row(infd_dem,inrast_dem,row,data_type_dem)<0)
+ G_fatal_error(_("Could not read from <%s>"),dem);
+ if(G_get_raster_row(infd_Rn, inrast_Rn, row,data_type_Rn) < 0)
+ G_fatal_error(_("Could not read from <%s>"),Rn);
+ if(G_get_raster_row(infd_g0, inrast_g0, row,data_type_g0) < 0)
+ G_fatal_error(_("Could not read from <%s>"),g0);
+ switch(data_type_T){
+ case CELL_TYPE:
+ d_tempk = (double) ((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double) ((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double) ((DCELL *) inrast_T)[col];
+ break;
+ }
+ switch(data_type_dem){
+ case CELL_TYPE:
+ d_dem = (double) ((CELL *) inrast_dem)[col];
+ break;
+ case FCELL_TYPE:
+ d_dem = (double) ((FCELL *) inrast_dem)[col];
+ break;
+ case DCELL_TYPE:
+ d_dem = (double) ((DCELL *) inrast_dem)[col];
+ break;
+ }
+ switch(data_type_Rn){
+ case CELL_TYPE:
+ d_Rn = (double) ((CELL *) inrast_Rn)[col];
+ break;
+ case FCELL_TYPE:
+ d_Rn = (double) ((FCELL *) inrast_Rn)[col];
+ break;
+ case DCELL_TYPE:
+ d_Rn = (double) ((DCELL *) inrast_Rn)[col];
+ break;
+ }
+ switch(data_type_g0){
+ case CELL_TYPE:
+ d_g0 = (double) ((CELL *) inrast_g0)[col];
+ break;
+ case FCELL_TYPE:
+ d_g0 = (double) ((FCELL *) inrast_g0)[col];
+ break;
+ case DCELL_TYPE:
+ d_g0 = (double) ((DCELL *) inrast_g0)[col];
+ break;
+ }
+ d_t0dem = d_tempk + 0.00627 * d_dem;
+ d_t0dem_dry = d_t0dem;
+ d_tempk_dry = d_tempk;
+ d_Rn_dry = d_Rn;
+ d_g0_dry = d_g0;
+ d_dem_dry = d_dem;
+ /*WET PIXEL*/
+ if(flag3->answer){
+ /*Calculate coordinates of row/col from projected ones*/
+ row = (int) (( ymax - m_row_wet ) / (double) stepy) ;
+ col = (int) (( m_col_wet - xmin ) / (double) stepx) ;
+ G_message("Wet Pixel | row:%i col:%i",row,col);
+ m_row_wet = row;
+ m_col_wet = col;
+ }
+ row = m_row_wet;
+ col = m_col_wet;
+ G_message("Wet Pixel | row:%i col:%i",row,col);
+ if(G_get_raster_row(infd_T,inrast_T,row,data_type_T)<0)
+ G_fatal_error(_("Could not read from <%s>"),T);
+ switch(data_type_T){
+ case CELL_TYPE:
+ d_tempk = (double) ((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double) ((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double) ((DCELL *) inrast_T)[col];
+ break;
+ }
+ d_tempk_wet = d_tempk;
+ G_message("Manual Pixels\n");
+ G_message("***************************\n");
+ G_message("row_wet=%d\tcol_wet=%d\n",(int) m_row_wet,(int) m_col_wet);
+ G_message("row_dry=%d\tcol_dry=%d\n",(int) m_row_dry,(int) m_col_dry);
+ G_message("tempk_wet=%f\n",d_tempk_wet);
+ G_message("tempk_dry=%f\n",d_tempk_dry);
+ G_message("dem_dry=%f\n",d_dem_dry);
+ G_message("t0dem_dry=%f\n",d_t0dem_dry);
+ G_message("rnet_dry=%f\n",d_Rn_dry);
+ G_message("g0_dry=%f\n",d_g0_dry);
+ G_message("h0_dry=%f\n",d_Rn_dry-d_g0_dry);
+ }
+ /* END OF MANUAL WET/DRY PIXELS */
+
+ for (row = 0; row < nrows; row++)
+ {
+ DCELL d_albedo;
+ DCELL d_tempk;
+ DCELL d_dem;
+ DCELL d_u2m;
+ DCELL d_t0dem;
+ DCELL d; /* Output pixel */
+ G_percent(row,nrows,2);
+ /* read a line input maps into buffers*/
+ if(G_get_raster_row(infd_albedo,inrast_albedo,row,data_type_albedo)<0)
+ G_fatal_error(_("Could not read from <%s>"),albedo);
+ if(G_get_raster_row(infd_T,inrast_T,row,data_type_T)<0)
+ G_fatal_error(_("Could not read from <%s>"),T);
+ if(G_get_raster_row(infd_dem,inrast_dem,row,data_type_dem)<0)
+ G_fatal_error(_("Could not read from <%s>"),dem);
+ if(G_get_raster_row(infd_u2m,inrast_u2m,row,data_type_u2m)<0)
+ G_fatal_error(_("Could not read from <%s>"),u2m);
+ if(G_get_raster_row(infd_ndvi,inrast_ndvi,row,data_type_ndvi) < 0)
+ G_fatal_error(_("Could not read from <%s>"),ndvi);
+ if(G_get_raster_row(infd_Rn,inrast_Rn,row,data_type_Rn) < 0)
+ G_fatal_error(_("Could not read from <%s>"),Rn);
+ if(G_get_raster_row(infd_g0,inrast_g0,row,data_type_g0) < 0)
+ G_fatal_error(_("Could not read from <%s>"),g0);
+ /* read every cell in the line buffers */
+ for (col=0; col < ncols; col++){
+ switch(data_type_albedo){
+ case CELL_TYPE:
+ d_albedo = (double) ((CELL *) inrast_albedo)[col];
+ break;
+ case FCELL_TYPE:
+ d_albedo = (double) ((FCELL *) inrast_albedo)[col];
+ break;
+ case DCELL_TYPE:
+ d_albedo = (double) ((DCELL *) inrast_albedo)[col];
+ break;
+ }
+ switch(data_type_T){
+ case CELL_TYPE:
+ d_tempk = (double) ((CELL *) inrast_T)[col];
+ break;
+ case FCELL_TYPE:
+ d_tempk = (double) ((FCELL *) inrast_T)[col];
+ break;
+ case DCELL_TYPE:
+ d_tempk = (double) ((DCELL *) inrast_T)[col];
+ break;
+ }
+ switch(data_type_dem){
+ case CELL_TYPE:
+ d_dem = (double) ((CELL *) inrast_dem)[col];
+ break;
+ case FCELL_TYPE:
+ d_dem = (double) ((FCELL *) inrast_dem)[col];
+ break;
+ case DCELL_TYPE:
+ d_dem = (double) ((DCELL *) inrast_dem)[col];
+ break;
+ }
+ switch(data_type_u2m){
+ case CELL_TYPE:
+ d_u2m = (double) ((CELL *) inrast_u2m)[col];
+ break;
+ case FCELL_TYPE:
+ d_u2m = (double) ((FCELL *) inrast_u2m)[col];
+ break;
+ case DCELL_TYPE:
+ d_u2m = (double) ((DCELL *) inrast_u2m)[col];
+ break;
+ }
+ switch(data_type_ndvi){
+ case CELL_TYPE:
+ d_ndvi = (double) ((CELL *) inrast_ndvi)[col];
+ break;
+ case FCELL_TYPE:
+ d_ndvi = (double) ((FCELL *) inrast_ndvi)[col];
+ break;
+ case DCELL_TYPE:
+ d_ndvi = (double) ((DCELL *) inrast_ndvi)[col];
+ break;
+ }
+ switch(data_type_Rn){
+ case CELL_TYPE:
+ d_Rn = (double) ((CELL *) inrast_Rn)[col];
+ break;
+ case FCELL_TYPE:
+ d_Rn = (double) ((FCELL *) inrast_Rn)[col];
+ break;
+ case DCELL_TYPE:
+ d_Rn = (double) ((DCELL *) inrast_Rn)[col];
+ break;
+ }
+ switch(data_type_g0){
+ case CELL_TYPE:
+ d_g0 = (double) ((CELL *) inrast_g0)[col];
+ break;
+ case FCELL_TYPE:
+ d_g0 = (double) ((FCELL *) inrast_g0)[col];
+ break;
+ case DCELL_TYPE:
+ d_g0 = (double) ((DCELL *) inrast_g0)[col];
+ break;
+ }
+ if(G_is_d_null_value(&d_tempk) ||
+ G_is_d_null_value(&d_dem) ||
+ G_is_d_null_value(&d_u2m) ||
+ G_is_d_null_value(&d_ndvi) ||
+ G_is_d_null_value(&d_Rn) ||
+ G_is_d_null_value(&d_g0) ||
+ d_g0<0.0 || d_Rn<0.0 ||
+ d_dem<=-100.0 || d_dem>9000.0 ||
+ d_tempk<200.0){
+ G_set_d_null_value(&outrast[col],1);
+ } else {
+ /* Albedo < 0*/
+ if(d_albedo<0.001){
+ d_albedo=0.001;
+ }
+ /* Calculate T0dem */
+ d_t0dem = (double)d_tempk + 0.00627 * (double) d_dem;
+ /* G_message("**InLoop d_t0dem=%5.3f",d_t0dem);
+ G_message(" d_dem=%5.3f",d_dem);
+ G_message(" d_tempk=%5.3f",d_tempk);
+ G_message(" d_albedo=%5.3f",d_albedo);
+ G_message(" d_Rn=%5.3f",d_Rn);
+ G_message(" d_g0=%5.3f",d_g0);
+ G_message(" d_ndvi=%5.3f",d_ndvi);
+ G_message(" d_u_hu=%5.3f",d_u_hu);
+ */ /* Calculate sensible heat flux */
+ d = sensi_h(iteration,
+ d_tempk_wet,
+ d_tempk_dry,
+ d_t0dem,
+ d_tempk,
+ d_ndvi,
+ d_ndvi_max,
+ d_dem,
+ d_Rn_dry,
+ d_g0_dry,
+ d_t0dem_dry,
+ d_u2m,
+ d_dem_dry);
+ // G_message(" d_h0=%5.3f",d);
+ if (zero->answer && d<0.0){
+ d=0.0;
+ }
+ outrast[col] = d;
+ }
+ }
+ if (G_put_d_raster_row (outfd, outrast) < 0)
+ G_fatal_error (_("Cannot write to <%s>"),h0);
+ }
+ G_free(inrast_T);
+ G_close_cell (infd_T);
+ G_free(inrast_dem);
+ G_close_cell (infd_dem);
+ G_free(inrast_u2m);
+ G_close_cell (infd_u2m);
+ G_free(inrast_ndvi);
+ G_close_cell (infd_ndvi);
+ G_free(inrast_Rn);
+ G_close_cell (infd_Rn);
+ G_free(inrast_g0);
+ G_close_cell (infd_g0);
+ G_free(inrast_albedo);
+ G_close_cell (infd_albedo);
+
+ G_free(outrast);
+ G_close_cell (outfd);
+ /* add command line incantation to history file */
+ G_short_history(h0, "raster", &history);
+ G_command_history(&history);
+ G_write_history(h0, &history);
+
+ exit(EXIT_SUCCESS);
+}
Added: grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,32 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double psi_h(double t0_dem, double h, double U_0, double roh_air)
+{
+ double result;
+ double n5_temp, n11_mem, n12_mem;
+
+// printf("h input to psih() is %5.3f\n",h);
+
+ if(h != 0.0){
+ n5_temp = (-1004* roh_air*pow(U_0,3)* t0_dem)/(0.41*9.81* h);
+ } else {
+ n5_temp = -1000.0;
+ }
+
+ if(n5_temp < 0.0){
+ n12_mem = pow((1-16*(2/n5_temp)),0.25);
+ n11_mem = (2*log10((1+pow(n12_mem,2))/2));
+ } else {
+ n12_mem = 1.0;
+ n11_mem = -5*2/n5_temp;
+ }
+
+ result = n11_mem;
+
+ return (result);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/psi_m.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/psi_m.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/psi_m.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,27 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double psi_m(double t0_dem, double h, double ustar, double roh_air, double hu)
+{
+ double result;
+ double n5, n10, n12;
+// printf("h input to psim() is %5.3f\n",h);
+ if(h != 0.0){
+ n5 = (-1004* roh_air*pow(ustar,3)* t0_dem)/(0.41*9.81* h);
+ } else {
+ n5 = -1000.0;
+ }
+ if(n5 < 0.0){
+ n12 = pow((1-16*(hu/n5)),0.25);
+ n10 = (2*log((1+n12)/2))+log((1+pow(n12,2))/2)-2*atan(n12)+0.5*PI;
+ } else {
+ n12 = 1.0;
+ n10 = -5*hu/n5;
+ }
+ result = n10;
+ return (result);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/psi_m.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/rah1.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/rah1.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/rah1.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,15 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double rah1(double psih, double ustar)
+{
+ double result;
+
+ result = (log10(2/0.01)-psih)/(ustar*0.41);
+
+ return (result);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/rah1.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,13 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double rah_0(double zom_0, double u_0)
+{
+ double result;
+
+ result = log10(2/0.01)/(u_0*0.41);
+
+ return (result);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/roh_air.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/roh_air.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/roh_air.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,15 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double rohair(double dem, double tempk, double dtair)
+{
+ double a, b, result;
+
+ a = tempk - dtair;
+ b = (( a - 0.00627*dem)/a);
+ result = 349.467 * pow( b , 5.26)/ a ;
+
+ return result;
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/roh_air.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/roh_air_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/roh_air_0.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/roh_air_0.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,14 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double roh_air_0(double tempk)
+{
+ double A, B, result;
+
+ A = (18*(6.11*exp(17.27*34.8/(34.8+237.3)))/100.0);
+ B = (18*(6.11*exp(17.27*34.8/(34.8+237.3)))/100.0);
+ result = (1000.0 - A)/(tempk*2.87)+ B/(tempk*4.61);
+ return result;
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/roh_air_0.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/sensi_h.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/sensi_h.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/sensi_h.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,112 @@
+/* This is the main loop used in SEBAL */
+/* ychemin at yahoo.com - yann.chemin at ait.ac.th */
+/* GPL >= 2 - April 2004 */
+
+#include<stdio.h>
+#include<math.h>
+#include<stdlib.h>
+#include "functions.h"
+
+/* Arrays Declarations */
+#define ITER_MAX 10
+
+double sensi_h( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u2m, double dem_desert)
+{
+ /* Arrays Declarations */
+ double dtair[ITER_MAX], roh_air[ITER_MAX], rah[ITER_MAX];
+ double h[ITER_MAX];
+ double ustar[ITER_MAX], zom[ITER_MAX];
+
+ /* Declarations */
+ int i, j, ic, debug=0;
+ double u_0, zom0;
+ double h_desert, rah_desert, roh_air_desert;
+ double dtair_desert;
+ double psih_desert,ustar_desert,ustar_desertold,zom_desert;
+ double psih;
+ double result;
+
+ /* Fat-free junk food */
+ if (iteration>ITER_MAX){
+ iteration=ITER_MAX;
+ }
+
+ if(debug==1){
+ printf("*****************************\n");
+ printf("t0_dem = %5.3f\n",t0_dem);
+ printf("ndvi = %5.3f ndvimax = %5.3f\n",ndvi,ndvi_max);
+ printf("*****************************\n");
+ }
+// dtair[0] = dt_air_0(t0_dem, tempk_water, tempk_desert);
+ dtair[0] = 5.0;
+// printf("*****************************dtair = %5.3f\n",dtair[0]);
+ roh_air[0] = roh_air_0(tempk);
+// printf("*****************************rohair=%5.3f\n",roh_air[0]);
+ roh_air_desert = roh_air_0(tempk_desert);
+// printf("**rohairdesert = %5.3f\n",roh_air_desert);
+ zom0 = zom_0(ndvi, ndvi_max);
+// printf("*****************************zom = %5.3f\n",zom0);
+ u_0 = U_0(zom0, u2m);
+// printf("*****************************u0\n");
+ rah[0] = rah_0(zom0, u_0);
+// printf("*****************************rah = %5.3f\n",rah[0]);
+ h[0] = h_0(roh_air[0], rah[0], dtair[0]);
+// printf("*****************************h\n");
+ if(debug==1){
+ printf("dtair[0] = %5.3f K\n", dtair[0]);
+ printf("roh_air[0] = %5.3f kg/m3\n", roh_air[0]);
+ printf("roh_air_desert0 = %5.3f kg/m3\n", roh_air_desert);
+ printf("zom_0 = %5.3f\n", zom0);
+ printf("u_0 = %5.3f\n", u_0);
+ printf("rah[0] = %5.3f s/m\n", rah[0]);
+ printf("h[0] = %5.3f W/m2\n", h[0]);
+ }
+
+/*----------------------------------------------------------------*/
+/*Main iteration loop of SEBAL*/
+ zom[0] = zom0;
+ for(ic=1;ic<iteration+1;ic++){
+ if(debug==1){
+ printf("\n ******** ITERATION %i *********\n",ic);
+ }
+ /* Where is roh_air[i]? */
+ psih = psi_h(t0_dem,h[ic-1],u_0,roh_air[ic-1]);
+ ustar[ic] = u_star(t0_dem,h[ic-1],u_0,roh_air[ic-1],zom[0],u2m);
+ rah[ic] = rah1(psih, ustar[ic]);
+ /* get desert point values from maps */
+ if(ic==1){
+ h_desert = rnet_desert - g0_desert;
+ zom_desert = 0.002;
+ psih_desert = psi_h(t0_dem_desert,h_desert,u_0,roh_air_desert);
+ ustar_desert = u_star(t0_dem_desert,h_desert,u_0,roh_air_desert,zom_desert,u2m);
+ } else {
+ roh_air_desert = rohair(dem_desert,tempk_desert,dtair_desert);
+ h_desert = h1(roh_air_desert,rah_desert,dtair_desert);
+ ustar_desertold = ustar_desert;
+ psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert);
+ ustar_desert = u_star(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,zom_desert,u2m);
+ }
+ rah_desert = rah1(psih_desert,ustar_desert);
+ dtair_desert = dt_air_desert(h_desert, roh_air_desert, rah_desert);
+ /* This should find the new dtair from inversed h equation...*/
+ dtair[ic] = dt_air(t0_dem, tempk_water, tempk_desert, dtair_desert);
+ /* This produces h[ic] and roh_air[ic+1] */
+ roh_air[ic] = rohair(dem, tempk, dtair[ic]);
+ h[ic] = h1(roh_air[ic], rah[ic], dtair[ic]);
+ /* Output values of the iteration parameters */
+ if(debug==1){
+ printf("psih[%i] = %5.3f\n", ic, psih);
+ printf("ustar[%i] = %5.3f\n", ic, ustar[ic]);
+ printf("rah[%i] = %5.3f s/m\n",ic, rah[ic]);
+ printf("h_desert = %5.3f\n", h_desert);
+ printf("rohair_desert = %5.3f\n", roh_air_desert);
+ printf("psih_desert = %5.3f\tustar_desert = %5.3f\trah_desert = %5.3f\n", psih_desert, ustar_desert, rah_desert);
+ printf("dtair_desert = %8.5f\n", dtair_desert);
+ printf("dtair[%i] = %5.3f K\n", ic, dtair[ic]);
+ printf("roh_air[%i] = %5.3f kg/m3\n", ic, roh_air[ic]);
+ printf("h[%i] = %5.3f W/m2\n",ic, h[ic]);
+ }
+ }
+ return h[iteration];
+}
+
Added: grass-addons/gipe/i.eb.h_SEBAL95/u_star.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/u_star.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/u_star.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,45 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+#define PI 3.1415927
+
+double u_star(double t0_dem, double h, double ustar, double roh_air, double zom, double u2m)
+{
+ double result;
+ double n5_temp; /* Monin-Obukov Length */
+ double n10_mem; /* psi m */
+ double n31_mem; /* x for U200 (that is bh...) */
+ double hv=0.15; /* crop height (m) */
+ double bh=200; /* blending height (m) */
+
+// printf("t0dem = %5.3f\n", t0_dem);
+// printf("h = %5.3f\n", h);
+// printf("U_0 = %5.3f\n", U_0);
+// printf("roh_air = %5.3f\n", roh_air);
+// printf("zom = %5.3f\n", zom);
+// printf("u2m = %5.3f\n", u2m);
+
+ if(h != 0.0){
+ n5_temp = (-1004* roh_air*pow(ustar,3)* t0_dem)/(0.41*9.81* h);
+ } else {
+ n5_temp = -1000.0;
+ }
+
+ if(n5_temp < 0.0){
+
+ n31_mem = pow((1-16*(200/n5_temp)),0.25);
+ n10_mem = (2*log10((1+n31_mem)/2)+log10((1+pow(n31_mem,2))/2)-2*atan(n31_mem)+0.5*PI);
+
+ } else {
+
+// n31_mem = 1.0;
+ n10_mem = -5*2/n5_temp;
+
+ }
+
+ result = ((u2m*0.41/log10(2/(hv/7)))/0.41*log10(bh /(hv/7)*0.41))/(log10(bh / zom)-n10_mem);
+
+ return (result);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/u_star.c
___________________________________________________________________
Name: svn:executable
+ *
Added: grass-addons/gipe/i.eb.h_SEBAL95/zom_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/zom_0.c (rev 0)
+++ grass-addons/gipe/i.eb.h_SEBAL95/zom_0.c 2008-08-03 09:46:40 UTC (rev 32480)
@@ -0,0 +1,19 @@
+#include<stdio.h>
+#include<math.h>
+#include"functions.h"
+
+double zom_0(double ndvi, double ndvi_max)
+{
+ double a, b, zom;
+// double hv_ndvimax=1.5; /* crop vegetation height (m) */
+// double hv_desert=0.002; /* desert base vegetation height (m) */
+
+// a = (log(hv_desert)-((log(hv_ndvimax/7)-log(hv_desert))/(ndvi_max-0.02)*0.02));
+// b = ((log(hv_ndvimax/7)-log(hv_desert))/(ndvi_max-0.02)) * ndvi;
+// zom = exp(a+b);
+
+// printf("****************zom = %5.3f\n", zom);
+ zom = exp(3.3219*ndvi-3.9939);
+ return (zom);
+}
+
Property changes on: grass-addons/gipe/i.eb.h_SEBAL95/zom_0.c
___________________________________________________________________
Name: svn:executable
+ *
More information about the grass-commit
mailing list