[GRASS-SVN] r46539 - in grass-addons/imagery/gipe: . i.evapo.zk
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jun 3 06:16:36 EDT 2011
Author: ychemin
Date: 2011-06-03 03:16:36 -0700 (Fri, 03 Jun 2011)
New Revision: 46539
Added:
grass-addons/imagery/gipe/i.evapo.zk/
grass-addons/imagery/gipe/i.evapo.zk/Makefile
grass-addons/imagery/gipe/i.evapo.zk/i.evapo.zk.html
grass-addons/imagery/gipe/i.evapo.zk/main.c
grass-addons/imagery/gipe/i.evapo.zk/zk_daily_et.c
Log:
Added i.evapo.zk (html file not formatted)
Added: grass-addons/imagery/gipe/i.evapo.zk/Makefile
===================================================================
--- grass-addons/imagery/gipe/i.evapo.zk/Makefile (rev 0)
+++ grass-addons/imagery/gipe/i.evapo.zk/Makefile 2011-06-03 10:16:36 UTC (rev 46539)
@@ -0,0 +1,10 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.evapo.pt
+
+LIBES = $(RASTERLIB) $(GISLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Added: grass-addons/imagery/gipe/i.evapo.zk/i.evapo.zk.html
===================================================================
--- grass-addons/imagery/gipe/i.evapo.zk/i.evapo.zk.html (rev 0)
+++ grass-addons/imagery/gipe/i.evapo.zk/i.evapo.zk.html 2011-06-03 10:16:36 UTC (rev 46539)
@@ -0,0 +1,92 @@
+<H2>DESCRIPTION</H2>
+
+<EM>i.evapo.zk</EM> Calculates the global diurnal evapotranspiration after Zhang, Kimball, Nemani and Running (2010).
+
+Zhang, K., Kimball, J.S., Nemani, R.R., Running, S.W. (2010). A continuous satellite-derived global record of land surface evapotranspiration from 1983 to 2006. WRR 46, W09522
+
+
+
+<H2>NOTES</H2>
+Main function:
+ETa (biome_type, ndvi, tday, sh, patm, Rn, G)
+
+Biome_type\t as defined below
+ndvi\t\t NDVI value\t\t [-]
+tday\t\t day temperature\t [C]
+sh\t\t specific humidity\t [-]
+patm\t\t atmospheric pressure\t [Pa]
+Rn\t\t day net radiation\t [MJ/m2/d]
+G\t\t day soil heat flux\t [MJ/m2/d]
+
+IGBP Biome types with ID used in this model
+-----------------------------------
+Code\t ID\t Description
+BENF\t 0\t Boreal Evergreen Needleleaf Forest (less or eq 212 frost-free days)
+TENF\t 1\t Temperate Evergreen Needleleaf Forest (more than 212 frost-free days)
+EBF\t 2\t Evergreen Broadleaf Forest
+DBF\t 4\t Deciduous Broadleaf Forest
+CSH\t 6\t Closed Shrubland
+OSH\t 7\t Open Shrubland
+WSV\t 8\t Woody Savannah
+SV\t 9\t Savannah
+GRS\t 10\t Grassland
+CRP\t 12\t Cropland
+
+IGBP Classification
+1 Evergreen Needleleaf Forest
+2 Evergreen Broadleaf Forest
+3* Deciduous Needleleaf Forest
+4 Deciduous Broadleaf Forest
+5* Mixed Forest (stats of other classes areas used in article)
+6 Closed Shrublands
+7 Open Shrublands
+8 Woody Savannas
+9 Savannas
+10 Grasslands
+11* Permanent Wetlands
+12 Croplands
+13* Urban and Built-Up
+14* Cropland/Natural Vegetation Mosaic
+15* Snow and Ice
+16* Barren or Sparsely Vegetated
+17* Water Bodies (Evaporation Priestley-Taylor used in article)
+
+* Not used in this model
+
+IGBP Biome types and configuration of internal parameters of the model"
+----------------------------------------------------------------"
+#Code Description TcloseMinC TopenMaxC VPDClosePa VPDOpenPa ToptC BetaC kPa GaMs-1 GtotMs-1 GchMs-1 B1Sm-1 B2Sm-1 B3 b1 b2 b3 b4"
+#BENF Boreal Evergreen Needleleaf Forest -8 40 2800 500 12 25 150 0.03 0.002 0.08 208.3 8333.3 10"
+#TENF Temperate Evergreen Needleleaf Forest -8 40 2800 500 25 25 200 0.03 0.004 0.08 133.3 888.9 6"
+#EBF Evergreen Broadleaf Forest -8 50 4000 500 40 40 300 0.03 0.006 0.01 57.7 769.2 4.5"
+#DBF Deciduous Broadleaf Forest -6 45 2800 650 28 25 200 0.04 0.002 0.01 85.8 694.7 4"
+#MF Mixed Forest"
+#CSH Closed Shrubland -8 45 3300 500 19 20 400 0.01 0.001 0.04 202 4040.4 6.5"
+#OSH Open Shrubland -8 40 3700 500 10 30 50 0.005 0.012 0.04 178.6 178.6 8"
+#WSV Woody Savannah -8 50 3200 500 32 28 900 0.002 0.0018 0.04 0.2 24000 6.5 57.1 3333.3 8 -0.01035"
+#SV Savannah -8 40 5000 650 32 30 800 0.001 0.001 0.04 790.9 8181.8 10"
+#GRS Grassland -8 40 3800 650 20 30 500 0.001 0.001 0.04 175 2000 6"
+#CRP Cropland -8 45 3800 650 20 30 450 0.005 0.003 0.04 105 3000 3"
+ #For WSV when NDVI>0.64"
+
+
+<H2>TODO</H2>
+
+
+<H2>SEE ALSO</H2>
+
+<em>
+<A HREF="i.evapo.pm.html">i.evapo.PM</A><br>
+<A HREF="i.evapo.potrad.html">i.evapo.potrad</A><br>
+<A HREF="i.eb.netrad.html">i.eb.netrad</A><br>
+<A HREF="i.eb.g0.html">i.eb.g0</A><br>
+</em>
+
+
+<H2>AUTHORS</H2>
+
+Yann Chemin, GRASS Development team, 2011<BR>
+
+
+<p>
+<i>Last changed: $Date: 2011-05-12 21:29:43 +0530 (Wed, 12 Jan 2011) $</i>
Added: grass-addons/imagery/gipe/i.evapo.zk/main.c
===================================================================
--- grass-addons/imagery/gipe/i.evapo.zk/main.c (rev 0)
+++ grass-addons/imagery/gipe/i.evapo.zk/main.c 2011-06-03 10:16:36 UTC (rev 46539)
@@ -0,0 +1,232 @@
+/*****************************************************************************
+*
+* MODULE: i.evapo.pt
+* AUTHOR: Yann Chemin yann.chemin at gmail.com
+*
+* PURPOSE: To estimate the daily evapotranspiration by means
+* of Prestley and Taylor method (1972).
+*
+* COPYRIGHT: (C) 2007-2011 by the GRASS Development Team
+*
+* This program is free software under the GNU General Public
+* Licence (>=2). Read the file COPYING that comes with GRASS
+* for details.
+*
+***************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+
+double eta (double biome_type, double ndvi, double tday, double sh, double patm,double Rn, double G, double dem);
+
+int main(int argc, char *argv[])
+{
+ /* buffer for input-output rasters */
+ void *inrast_biomt, *inrast_ndvi, *inrast_tday, *inrast_sh;
+ void *inrast_patm, *inrast_rnetd, *inrast_g0, *inrast_dem;
+ DCELL *outrast;
+
+ /* pointers to input-output raster files */
+ int infd_biomt, infd_ndvi, infd_tday, infd_sh;
+ int infd_patm, infd_rnetd, infd_g0, infd_dem;
+ int outfd;
+
+ /* names of input-output raster files */
+ char *biomt, *ndvi, *tday, *sh;
+ char *patm, *rnetd, *g0, *dem;
+ char *eta;
+
+ /* input-output cell values */
+ DCELL d_biomt, d_ndvi, d_tday, d_sh;
+ DCELL d_patm, d_rnetd, d_g0, d_dem;
+ DCELL d_daily_et;
+
+ /* region informations and handler */
+ struct Cell_head cellhd;
+ int nrows, ncols;
+ int row, col;
+
+ /* parser stuctures definition */
+ struct GModule *module;
+ struct Option *input_biomt, *input_ndvi, *input_tday, *input_sh;
+ struct Option *input_patm, *input_rnetd, *input_g0, *input_dem;
+ struct Option *output;
+ struct Flag *zero;
+ struct Colors color;
+ struct History history;
+
+ /* Initialize the GIS calls */
+ G_gisinit(argv[0]);
+
+ module = G_define_module();
+ G_add_keyword(_("imagery"));
+ G_add_keyword(_("evapotranspiration"));
+ module->description =
+ _("Computes global evapotranspiration calculation "
+ "Zhang, Kimball, Nemani and Running formulation, 2010.");
+
+ /* Define different options */
+
+ input_biomt = G_define_standard_option(G_OPT_R_INPUT);
+ input_biomt->key = "biome_type";
+ input_biomt->description = _("Name of input IGBP biome type raster map [-]");
+
+ input_ndvi = G_define_standard_option(G_OPT_R_INPUT);
+ input_ndvi->key = "ndvi";
+ input_ndvi->description = _("Name of input Normalized Difference Vegetation Index raster map [-]");
+
+ input_tday = G_define_standard_option(G_OPT_R_INPUT);
+ input_tday->key = "airtemperature";
+ input_tday->description = _("Name of input air temperature raster map [C]");
+
+ input_sh = G_define_standard_option(G_OPT_R_INPUT);
+ input_sh->key = "specifichumidity";
+ input_sh->description = _("Name of input specific humidity raster map [-]");
+
+ input_patm = G_define_standard_option(G_OPT_R_INPUT);
+ input_patm->key = "atmosphericpressure";
+ input_patm->description = _("Name of input atmospheric pressure raster map [Pa]");
+
+ input_rnetd = G_define_standard_option(G_OPT_R_INPUT);
+ input_rnetd->key = "netradiation";
+ input_rnetd->description = _("Name of input net radiation raster map [MJ/m2/d]");
+
+ input_g0 = G_define_standard_option(G_OPT_R_INPUT);
+ input_g0->key = "soilheatflux";
+ input_g0->description = _("Name of input soil heat flux raster map [MJ/m2/d]");
+
+ input_dem = G_define_standard_option(G_OPT_R_INPUT);
+ input_dem->key = "elevation";
+ input_dem->description = _("Name of input elevation raster map [m]");
+
+ output = G_define_standard_option(G_OPT_R_OUTPUT);
+ output->description = _("Name of output evapotranspiration raster map [mm/d]");
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* get entered parameters */
+ biomt = input_biomt->answer;
+ ndvi = input_ndvi->answer;
+ tday = input_tday->answer;
+ sh = input_sh->answer;
+ patm = input_patm->answer;
+ rnetd = input_rnetd->answer;
+ g0 = input_g0->answer;
+ dem = input_dem->answer;
+
+ eta = output->answer;
+
+ /* open pointers to input raster files */
+ infd_biomt = Rast_open_old(biomt, "");
+ infd_ndvi = Rast_open_old(ndvi, "");
+ infd_tday = Rast_open_old(tday, "");
+ infd_sh = Rast_open_old(sh, "");
+ infd_patm = Rast_open_old(patm, "");
+ infd_rnetd = Rast_open_old(rnetd, "");
+ infd_g0 = Rast_open_old(g0, "");
+ infd_dem = Rast_open_old(dem, "");
+
+ /* read headers of raster files */
+ Rast_get_cellhd(biomt, "", &cellhd);
+ Rast_get_cellhd(ndvi, "", &cellhd);
+ Rast_get_cellhd(tday, "", &cellhd);
+ Rast_get_cellhd(sh, "", &cellhd);
+ Rast_get_cellhd(patm, "", &cellhd);
+ Rast_get_cellhd(rnetd, "", &cellhd);
+ Rast_get_cellhd(g0, "", &cellhd);
+ Rast_get_cellhd(dem, "", &cellhd);
+
+ /* Allocate input buffer */
+ inrast_biomt = Rast_allocate_d_buf();
+ inrast_ndvi = Rast_allocate_d_buf();
+ inrast_tday = Rast_allocate_d_buf();
+ inrast_sh = Rast_allocate_d_buf();
+ inrast_patm = Rast_allocate_d_buf();
+ inrast_rnetd = Rast_allocate_d_buf();
+ inrast_g0 = Rast_allocate_d_buf();
+ inrast_dem = Rast_allocate_d_buf();
+
+ /* get rows and columns number of the current region */
+ nrows = Rast_window_rows();
+ ncols = Rast_window_cols();
+
+ /* allocate output buffer */
+ outrast = Rast_allocate_d_buf();
+
+ /* open pointers to output raster files */
+ outfd = Rast_open_new(eta, DCELL_TYPE);
+
+ /* start the loop through cells */
+ for (row = 0; row < nrows; row++) {
+
+ G_percent(row, nrows, 2);
+ /* read input raster row into line buffer */
+ Rast_get_d_row(infd_biomt, inrast_biomt, row);
+ Rast_get_d_row(infd_ndvi, inrast_ndvi, row);
+ Rast_get_d_row(infd_tday, inrast_tday, row);
+ Rast_get_d_row(infd_sh, inrast_sh, row);
+ Rast_get_d_row(infd_patm, inrast_patm, row);
+ Rast_get_d_row(infd_rnetd, inrast_rnetd, row);
+ Rast_get_d_row(infd_g0, inrast_g0, row);
+ Rast_get_d_row(infd_dem, inrast_dem, row);
+
+ for (col = 0; col < ncols; col++) {
+ /* read current cell from line buffer */
+ d_biomt = ((DCELL *) inrast_biomt)[col];
+ d_ndvi = ((DCELL *) inrast_ndvi)[col];
+ d_tday = ((DCELL *) inrast_tday)[col];
+ d_sh = ((DCELL *) inrast_sh)[col];
+ d_patm = ((DCELL *) inrast_patm)[col];
+ d_rnetd = ((DCELL *) inrast_rnetd)[col];
+ d_g0 = ((DCELL *) inrast_g0)[col];
+ d_dem = ((DCELL *) inrast_dem)[col];
+
+ /*Calculate ET */
+ d_daily_et = zk_daily_et(d_biomt, d_ndvi, d_tday, d_sh, d_patm, d_rnetd, d_g0, d_dem);
+
+ /* write calculated ETP to output line buffer */
+ outrast[col] = d_daily_et;
+ }
+
+ /* write output line buffer to output raster file */
+ Rast_put_d_row(outfd, outrast);
+ }
+ /* free buffers and close input maps */
+
+ G_free(inrast_biomt);
+ G_free(inrast_ndvi);
+ G_free(inrast_tday);
+ G_free(inrast_sh);
+ G_free(inrast_patm);
+ G_free(inrast_rnetd);
+ G_free(inrast_g0);
+ G_free(inrast_dem);
+ Rast_close(infd_biomt);
+ Rast_close(infd_ndvi);
+ Rast_close(infd_tday);
+ Rast_close(infd_sh);
+ Rast_close(infd_patm);
+ Rast_close(infd_rnetd);
+ Rast_close(infd_g0);
+ Rast_close(infd_dem);
+
+ /* generate color table between -20 and 20 */
+ Rast_make_rainbow_colors(&color, -20, 20);
+ Rast_write_colors(eta, G_mapset(), &color);
+
+ Rast_short_history(eta, "raster", &history);
+ Rast_command_history(&history);
+ Rast_write_history(eta, &history);
+
+ /* free buffers and close output map */
+ G_free(outrast);
+ Rast_close(outfd);
+
+ return (EXIT_SUCCESS);
+}
Added: grass-addons/imagery/gipe/i.evapo.zk/zk_daily_et.c
===================================================================
--- grass-addons/imagery/gipe/i.evapo.zk/zk_daily_et.c (rev 0)
+++ grass-addons/imagery/gipe/i.evapo.zk/zk_daily_et.c 2011-06-03 10:16:36 UTC (rev 46539)
@@ -0,0 +1,316 @@
+#include<math.h>
+
+double ndvimin, ndvimax;
+/*Biome Type Configuration variables*/
+double tclosemin, topenmax;
+double vpdclose, vpdopen;
+double topt, beta, k;
+double ga, gtot, gch;
+double b1, b2, b3, b4;
+/*Biome Type Configuration is valid=1, is invalid=0*/
+int valid_conf;
+
+/*Variables*/
+double e_sat, b, t;
+double gs, g0, mtday;
+double rh, vpd, mvpd;
+double fracveg, latent, MaMw;
+double Cp, psi, gtotc, Delta;
+double rho, Esoil, Ecanopy;
+
+double fc( double ndvi ){
+ /*Fraction of vegetation cover*/
+ ndvimin = 0.05;
+ ndvimax = 0.95;
+ return ( ( ndvi - ndvimin ) / ( ndvimax - ndvimin ) );
+}
+
+double bdpc( double ndvi, double b1, double b2, double b3, double b4 ){
+ /*Biome dependent potential conductance (g0 in Zhang et al, 2010. WRR)*/
+ return (1.0 / (b1 + b2 * exp(-b3*ndvi))+b4);
+}
+
+double mTday( double tday, double tclosemin, double topenmax, double topt, double beta ){
+ /*Temperature dependent reduction function of the biome dependent potential conductance*/
+ if (tday <= tclosemin) return 0.01;
+ else if (tday >= topenmax) return 0.01;
+ else return ( exp( - ( (tday-topt)/beta )*( (tday-topt)/beta ) ) );
+}
+
+double mVPD( double vpd, double vpdclose, double vpdopen ){
+ /*VPD dependent reduction function of the biome dependent potential conductance*/
+ if (vpd <= vpdopen) return 1.0;
+ else if (vpd >= vpdclose) return 1.0;
+ else return ( ( vpdclose - vpd ) / ( vpdclose - vpdopen ) );
+}
+
+double esat( tday ){
+ /*From FAO56, output is in [Pa]*/
+ return( 1000 * 0.6108 * exp( 17.27 * tday / (tday + 237.3) ) );
+}
+
+double vpdeficit( double rh, double tday ){
+ /*From FAO56, vpd is esat - eact*/
+ e_sat = esat( tday );
+ return ( ( 1.0 - rh ) * e_sat );
+}
+
+double rhumidity( double sh, double tday, double patm ){
+ /*http://en.wikipedia.org/wiki/Humidity#Specific_humidity
+ Output [0.0-1.0]*/
+ e_sat = esat( tday );
+ /*printf("e_sat\t=%f\t[Pa]\n",e_sat);*/
+ b = ( 0.622 + sh ) * e_sat;
+ return ( sh * patm / b );
+}
+
+double slopesvpcurve( double tday ){
+ /*From FAO56, output is in [Pa/C]*/
+ e_sat = esat( tday );
+ return( 4098.0 * e_sat / pow( tday + 237.3, 2 ) );
+}
+
+double rhoair( double dem, double tday ){
+ /*Requires T in Kelvin*/
+ t = tday + 273.15;
+ b = ( ( t - ( 0.00627 * dem ) ) / t );
+ return( 349.467 * pow( b, 5.26 ) / t );
+}
+
+double zk_daily_et ( double biome_type, double ndvi, double tday, double sh, double patm, double rnetd, double soilHF, double dem){
+ /*ETa global model from Zhang et al 2010. A continuous satellite derived global record of land surface evapotranspiration from 1983 to 2006. WRR 46*/
+ b4 = 0.0; /*except for WSV with NDVI > 0.64*/
+ if (biome_type == 0){
+ /*BENF*/
+ tclosemin = -8.0;
+ topenmax = 40.0;
+ vpdclose = 2800.0;
+ vpdopen = 500.0;
+ topt = 12.0;
+ beta = 25.0;
+ k = 150.0;
+ ga = 0.03;
+ gtot = 0.002;
+ gch = 0.08;
+ b1 = 208.3;
+ b2 = 8333.3;
+ b3 = 10.0;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 1){
+ /*TENF*/
+ tclosemin = -8.0;
+ topenmax = 40.0;
+ vpdclose = 2800.0;
+ vpdopen = 500.0;
+ topt = 25.0;
+ beta = 25.0;
+ k = 200.0;
+ ga = 0.03;
+ gtot = 0.004;
+ gch = 0.08;
+ b1 = 133.3;
+ b2 = 888.9;
+ b3 = 6.0;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 2){
+ /*EBF*/
+ tclosemin = -8.0;
+ topenmax = 50.0;
+ vpdclose = 4000.0;
+ vpdopen = 500.0;
+ topt = 40.0;
+ beta = 40.0;
+ k = 300.0;
+ ga = 0.03;
+ gtot = 0.006;
+ gch = 0.01;
+ b1 = 57.7;
+ b2 = 769.2;
+ b3 = 4.5;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 4){
+ /*DBF*/
+ tclosemin = -6.0;
+ topenmax = 45.0;
+ vpdclose = 2800.0;
+ vpdopen = 650.0;
+ topt = 28.0;
+ beta = 25.0;
+ k = 200.0;
+ ga = 0.04;
+ gtot = 0.002;
+ gch = 0.01;
+ b1 = 85.8;
+ b2 = 694.7;
+ b3 = 4;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 6){
+ /*CSH*/
+ tclosemin = -8.0;
+ topenmax = 45.0;
+ vpdclose = 3300.0;
+ vpdopen = 500.0;
+ topt = 19.0;
+ beta = 20.0;
+ k = 400.0;
+ ga = 0.01;
+ gtot = 0.001;
+ gch = 0.04;
+ b1 = 202.0;
+ b2 = 4040.4;
+ b3 = 6.5;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 7){
+ /*OSH*/
+ tclosemin = -8.0;
+ topenmax = 40.0;
+ vpdclose = 3700.0;
+ vpdopen = 500.0;
+ topt = 10.0;
+ beta = 30.0;
+ k = 50.0;
+ ga = 0.005;
+ gtot = 0.012;
+ gch = 0.04;
+ b1 = 178.6;
+ b2 = 178.6;
+ b3 = 8;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 8 && ndvi < 0.64){
+ /*WSV && ndvi < 0.64*/
+ tclosemin = -8.0;
+ topenmax = 50.0;
+ vpdclose = 3200.0;
+ vpdopen = 500.0;
+ topt = 32.0;
+ beta = 28.0;
+ k = 900.0;
+ ga = 0.002;
+ gtot = 0.0018;
+ gch = 0.04;
+ b1 = 0.2;
+ b2 = 24000;
+ b3 = 6.5;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 8 && ndvi > 0.64){
+ /*WSV && ndvi > 0.64*/
+ tclosemin = -8.0;
+ topenmax = 50.0;
+ vpdclose = 3200.0;
+ vpdopen = 500.0;
+ topt = 32.0;
+ beta = 28.0;
+ k = 900.0;
+ ga = 0.002;
+ gtot = 0.0018;
+ gch = 0.04;
+ b1 = 57.1;
+ b2 = 3333.3;
+ b3 = 8.0;
+ b4 = -0.01035;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 9){
+ /*SV*/
+ tclosemin = -8.0;
+ topenmax = 40.0;
+ vpdclose = 5000.0;
+ vpdopen = 650.0;
+ topt = 32.0;
+ beta = 30.0;
+ k = 800.0;
+ ga = 0.001;
+ gtot = 0.001;
+ gch = 0.04;
+ b1 = 790.9;
+ b2 = 8181.8;
+ b3 = 10.0;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 10){
+ /*GRS*/
+ tclosemin = -8.0;
+ topenmax = 40.0;
+ vpdclose = 3800.0;
+ vpdopen = 650.0;
+ topt = 20.0;
+ beta = 30.0;
+ k = 500.0;
+ ga = 0.001;
+ gtot = 0.001;
+ gch = 0.04;
+ b1 = 175.0;
+ b2 = 2000;
+ b3 = 6.0;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else if (biome_type == 12){
+ /*CRP*/
+ tclosemin = -8.0;
+ topenmax = 45.0;
+ vpdclose = 3800.0;
+ vpdopen = 650.0;
+ topt = 20.0;
+ beta = 30.0;
+ k = 450.0;
+ ga = 0.005;
+ gtot = 0.003;
+ gch = 0.04;
+ b1 = 105.0;
+ b2 = 300.0;
+ b3 = 3.0;
+ /* Configuration is valid*/
+ valid_conf = 1;
+ } else {
+ /* Configuration is invalid*/
+ valid_conf = 0;
+ }
+ /*Compute potential conductance for this biome and this NDVI*/
+ g0 = bdpc(ndvi,b1,b2,b3,b4);
+ /*Preprocessing for Surface conductance (gs) in PM (FAO56), gc in this article*/
+ mtday = mTday( tday, tclosemin, topenmax, topt, beta );
+ /*relative humidity*/
+ rh = rhumidity( sh, tday, patm );
+ /*printf("rh\t=%f\t[-]\n",rh);*/
+ vpd = vpdeficit( rh, tday );
+ /*printf("vpd\t=%f\t\t[Pa]\n",vpd);*/
+ mvpd = mVPD( vpd, vpdclose, vpdopen );
+ /*Actually computing Surface conductance (gs) in PM (FAO56), gc in this article*/
+ gs = g0 * mtday * mvpd;
+ /*print "rs\t=%f\t[s/m]\n",1/gs);*/
+ /*Fraction of vegetation cover*/
+ fracveg = fc(ndvi);
+ /*printf("fc\t=%f\t[-]\n", fracveg);*/
+ /*preprocessing for soil Evaporation*/
+ latent = 2.45; /*MJ/Kg FAO56*/
+ MaMw = 0.622; /* - FAO56*/
+ Cp = 1.013 * 0.001; /* MJ/Kg/C FAO56 */
+ psi = patm * Cp / (MaMw * latent); /*psi = patm * 0.6647 / 1000*/
+ /*printf("psi\t=%f\t[Pa/C]\n",psi);*/
+ gtotc = gtot * ((273.15+tday) / 293.13) * (101300.0 / patm);
+ Delta = slopesvpcurve( tday ); /*slope in Pa/C*/
+ /*printf("Delta\t=%f\t[de/dt]\n",Delta);*/
+ rho = rhoair( dem, tday );
+ /*printf("rho\t=%f\t[kg/m3]\n",rho);*/
+ /*soil Evaporation*/
+ Esoil = pow(rh,vpd/k) * (Delta*(1-fracveg)*(rnetd-soilHF)+rho*Cp*vpd*ga) / (Delta+psi*ga/gtotc);
+ /*Canopy evapotranspiration*/
+ Ecanopy = (Delta*fracveg*(rnetd-soilHF)+rho*Cp*vpd*ga) / (Delta+psi*(1.0+ga/gs));
+ /*printf("------------------------------------------------------\n");*/
+ /*printf("Esoil\t=%f\t[mm/d]\n", Esoil);*/
+ /*printf("Ecanopy\t=%f\t[mm/d]\n", Ecanopy);*/
+ /*printf("------------------------------------------------------\n");*/
+ if (valid_conf){
+ return( (1-fracveg) * Esoil + fracveg * Ecanopy );
+ } else {
+ return( -28768 );
+ }
+}
+
More information about the grass-commit
mailing list