[GRASS-SVN] r31643 - grass-addons/gipe/i.eb.h_iter
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 5 23:36:32 EDT 2008
Author: ychemin
Date: 2008-06-05 23:36:32 -0400 (Thu, 05 Jun 2008)
New Revision: 31643
Modified:
grass-addons/gipe/i.eb.h_iter/fixed_deltat.c
grass-addons/gipe/i.eb.h_iter/main.c
Log:
Updated i.eb.h_iter to new programming standard, added iteration number input option
Modified: grass-addons/gipe/i.eb.h_iter/fixed_deltat.c
===================================================================
--- grass-addons/gipe/i.eb.h_iter/fixed_deltat.c 2008-06-05 19:57:03 UTC (rev 31642)
+++ grass-addons/gipe/i.eb.h_iter/fixed_deltat.c 2008-06-06 03:36:32 UTC (rev 31643)
@@ -4,7 +4,7 @@
#define PI 3.14159265358979323846
-double fixed_deltat(double u2m, double roh_air,double cp,double dt,double disp,double z0m,double z0h,double tempk){
+double fixed_deltat(double u2m, double roh_air,double cp,double dt,double disp,double z0m,double z0h,double tempk,int iteration){
int i;
double ublend;
double length;
@@ -14,20 +14,20 @@
double rah;
double h_in;
- ublend=u2m*(log(100-disp)-log(z0m))/(log(2-disp)-log(z0m));
+ ublend=u2m*(log10(100-disp)-log10(z0m))/(log10(2-disp)-log10(z0m));
psim=0.0;
psih=0.0;
- for(i=0;i<10;i++){
- ustar = 0.41*ublend/(log((100-disp)/z0m)-psim);
- rah = (log((2-disp)/z0h)-psih)/(0.41*ustar);
+ for(i=0;i<iteration;i++){
+ ustar = 0.41*ublend/(log10((100-disp)/z0m)-psim);
+ rah = (log10((2-disp)/z0h)-psih)/(0.41*ustar);
h_in = roh_air * cp * dt / rah;
length= -roh_air*cp*pow(ustar,3)*tempk/(0.41*9.81*h_in);
xm = pow(1.0-16.0*((100-disp)/length),0.25);
xh = pow(1.0-16.0*((2-disp)/length),0.25);
- psim = 2.0*log((1.0+xm)/2.0)+log((1+xm*xm)-2*atan(xm)+0.5*PI);
- psih = 2.0*log((1.0+xh*xh)/2.0);
+ psim = 2.0*log10((1.0+xm)/2.0)+log10((1+xm*xm)-2*atan(xm)+0.5*PI);
+ psih = 2.0*log10((1.0+xh*xh)/2.0);
}
return rah;
Modified: grass-addons/gipe/i.eb.h_iter/main.c
===================================================================
--- grass-addons/gipe/i.eb.h_iter/main.c 2008-06-05 19:57:03 UTC (rev 31642)
+++ grass-addons/gipe/i.eb.h_iter/main.c 2008-06-06 03:36:32 UTC (rev 31643)
@@ -3,7 +3,8 @@
* MODULE: i.eb.h_iter
* AUTHOR(S): Yann Chemin - yann.chemin at gmail.com
* PURPOSE: Calculates sensible heat flux by SEBAL iteration
- * if your DeltaT (or eq. to make it from surf.temp) is validated.
+ * if your DeltaT (or eq. to make it from surf.temp)
+ * is validated.
* ! Delta T will not be reassessed in the iterations !
* A flag allows the Bastiaanssen (1995) affine transform
* of surface temperature as used in his SEBAL model.
@@ -15,8 +16,10 @@
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
- * CHANGELOG: 28 October 06: Added u2m input (wind speed at 2 meters height)
- * 30 June 07: Debugging a Seg Fault
+ * CHANGELOG: 06 June 08: Added iteration number input
+ * 28 October 06: Added u2m input
+ * (wind speed at 2 meters height)
+ * 30 June 07: Debugging a Seg Fault
*
*****************************************************************************/
@@ -27,7 +30,7 @@
#include <grass/glocale.h>
-double fixed_deltat(double u2m, double roh_air,double cp,double dt,double disp,double z0m,double z0h,double tempk);
+double fixed_deltat(double u2m, double roh_air,double cp,double dt,double disp,double z0m,double z0h,double tempk, int iteration);
double h0(double roh_air, double cp, double rah, double dtair);
int main(int argc, char *argv[])
@@ -41,7 +44,8 @@
int sebal=0;//SEBAL Flag for affine transform of surf. temp.
struct GModule *module;
struct Option *input1, *input2, *input3, *input4, *input5;
- struct Option *input6, *input7, *input8, *input9, *input10, *output1;
+ struct Option *input6, *input7, *input8, *input9, *input10;
+ struct Option *input11, *output1;
struct Flag *flag1, *flag2;
struct History history; //metadata
@@ -60,10 +64,11 @@
double cp; //air specific heat
int i=0,j=0;
double a,b; //SEBAL slope and intercepts of surf. temp.
+ int iteration;
void *inrast_rohair, *inrast_tempk, *inrast_dtair;
void *inrast_disp, *inrast_z0m, *inrast_z0h, *inrast_u2m;
- unsigned char *outrast;
+ DCELL *outrast;
RASTER_MAP_TYPE data_type_output=DCELL_TYPE;
RASTER_MAP_TYPE data_type_rohair;
RASTER_MAP_TYPE data_type_dtair;
@@ -136,6 +141,14 @@
input10->key =_("u2m");
input10->description=_("Name of the wind speed at 2m height input layer (m/s)");
input10->answer =_("u2m");
+
+ input11 = G_define_option() ;
+ input11->key =_("iteration");
+ input11->type = TYPE_INTEGER;
+ input11->required = NO;
+ input11->gisprompt =_("parameter, integer number");
+ input11->description=_("number of iteration of rah stabilization");
+ input11->answer =_("3");
output1 = G_define_standard_option(G_OPT_R_OUTPUT) ;
output1->key =_("h0");
@@ -168,10 +181,14 @@
z0m = input8->answer;
z0h = input9->answer;
u2m = input10->answer;
+ if(input11->answer){
+ iteration = atoi(input11->answer);
+ }else{
+ iteration = 3;
+ }
- result = output1->answer;
- sebal = flag1->answer;
-
+ result = output1->answer;
+ sebal = flag1->answer;
/***************************************************/
/* TEST FOR -s FLAG COEFS
* Return error if not proper flag coefs
@@ -383,36 +400,30 @@
d_u2m = ((DCELL *) inrast_u2m)[col];
break;
}
- if(G_is_d_null_value(&d_rohair)){
- ((DCELL *) outrast)[col] = -999.99;
- }else if((!sebal)&&G_is_d_null_value(&d_dtair)){
- ((DCELL *) outrast)[col] = -999.99;
- }else if(G_is_d_null_value(&d_tempk)){
- ((DCELL *) outrast)[col] = -999.99;
- }else if(G_is_d_null_value(&d_disp)){
- ((DCELL *) outrast)[col] = -999.99;
- }else if(G_is_d_null_value(&d_z0m)){
- ((DCELL *) outrast)[col] = -999.99;
- }else if(G_is_d_null_value(&d_z0h)){
- ((DCELL *) outrast)[col] = -999.99;
- }else if(G_is_d_null_value(&d_u2m)){
- ((DCELL *) outrast)[col] = -999.99;
+ if(G_is_d_null_value(&d_rohair)||
+ ((!sebal)&&G_is_d_null_value(&d_dtair))||
+ G_is_d_null_value(&d_tempk)||
+ G_is_d_null_value(&d_disp)||
+ G_is_d_null_value(&d_z0m)||
+ G_is_d_null_value(&d_z0h)||
+ G_is_d_null_value(&d_u2m)){
+ G_set_d_null_value(&outrast[col],1);
}else {
/************************************/
/* calculate sensible heat flux */
if(sebal){/*if -s flag then calculate Delta T from Coefs*/
d_affine=a*d_tempk+b;
/* Run iterations to find Rah*/
- d_rah=fixed_deltat(d_u2m,d_rohair,cp,d_affine,d_disp,d_z0m,d_z0h,d_tempk);
+ d_rah=fixed_deltat(d_u2m,d_rohair,cp,d_affine,d_disp,d_z0m,d_z0h,d_tempk,iteration);
/*Process h*/
d = h0(d_rohair,cp,d_rah,d_affine);
}else{/* not -s flag then take delta T map input directly*/
/* Run iterations to find Rah*/
- d_rah=fixed_deltat(d_u2m,d_rohair,cp,d_dtair,d_disp,d_z0m,d_z0h,d_tempk);
+ d_rah=fixed_deltat(d_u2m,d_rohair,cp,d_dtair,d_disp,d_z0m,d_z0h,d_tempk,iteration);
/*Process h*/
d = h0(d_rohair,cp,d_rah,d_dtair);
}
- ((DCELL *) outrast)[col] = d;
+ outrast[col] = d;
}
}
if (G_put_raster_row (outfd, outrast, data_type_output) < 0)
More information about the grass-commit
mailing list