[GRASS-SVN] r32346 - grass-addons/gipe/i.eb.h_SEBAL95
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Jul 27 19:10:29 EDT 2008
Author: ychemin
Date: 2008-07-27 19:10:29 -0400 (Sun, 27 Jul 2008)
New Revision: 32346
Removed:
grass-addons/gipe/i.eb.h_SEBAL95/dtair_0.c
Modified:
grass-addons/gipe/i.eb.h_SEBAL95/main.c
grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_noz0m.c
grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_z0m.c
Log:
rearranged initialization parameters of cold and dry pixels, allowed external dT map to overrule internal initialization of cold and dry pixels
Deleted: grass-addons/gipe/i.eb.h_SEBAL95/dtair_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/dtair_0.c 2008-07-27 22:32:43 UTC (rev 32345)
+++ grass-addons/gipe/i.eb.h_SEBAL95/dtair_0.c 2008-07-27 23:10:29 UTC (rev 32346)
@@ -1,44 +0,0 @@
-#include<stdio.h>
-#include<math.h>
-#include"functions.h"
-
-/* Pixel-based input required are: tempk water & desert
- * additionally, dtair in Desert is vaguely initialized
- */
-#define ZERO 273.15
-
-double dt_air_0(double t0_dem, double tempk_water, double tempk_desert)
-{
- double a, b, result;
- double dtair_desert_0;
-
- if(tempk_desert > (ZERO+48.0)){
- dtair_desert_0 = 13.0;
- } else if(tempk_desert >= (ZERO+40.0) && tempk_desert < (ZERO+48.0)){
- dtair_desert_0 = 10.0;
- } else if(tempk_desert >= (ZERO+32.0) && tempk_desert < (ZERO+40.0)){
- dtair_desert_0 = 7.0;
- } else if(tempk_desert >= (ZERO+25.0) && tempk_desert < (ZERO+32.0)){
- dtair_desert_0 = 5.0;
- } else if(tempk_desert >= (ZERO+18.0) && tempk_desert < (ZERO+25.0)){
- dtair_desert_0 = 3.0;
- } else if(tempk_desert >= (ZERO+11.0) && tempk_desert < (ZERO+18.0)){
- dtair_desert_0 = 1.0;
- } else {
- dtair_desert_0 = 0.0;
-// printf("WARNING!!! dtair_desert_0 is NOT VALID!\n");
- }
-
-// printf("dtair0 = %.0f K\t",dtair_desert_0);
-
- a = (dtair_desert_0-0.0)/(tempk_desert-tempk_water);
- b = 0.0 - a * tempk_water;
-
-// printf("dt_air_0(a) = %5.3f Tempk(b) %5.3f\n",a,b);
-
- result = t0_dem * a + b;
-// printf("dt_air_0 = %5.3f\n",result);
-
- return result;
-}
-
Modified: grass-addons/gipe/i.eb.h_SEBAL95/main.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/main.c 2008-07-27 22:32:43 UTC (rev 32345)
+++ grass-addons/gipe/i.eb.h_SEBAL95/main.c 2008-07-27 23:10:29 UTC (rev 32346)
@@ -26,8 +26,8 @@
/*#include <omp.h>*/
-double sensi_h_z0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double zom0, double dtair0, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert);
-double sensi_h_noz0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double dtair0, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert);
+double sensi_h_z0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double zom0, double dtair0, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert, double dtair_desert);
+double sensi_h_noz0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double dtair0, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert, double dtair_desert);
int main(int argc, char *argv[])
{
@@ -427,6 +427,8 @@
DCELL d_g0_dry;
DCELL d_t0dem_dry;
DCELL d_dem_dry;
+ DCELL d_dT_dry;
+ DCELL d_dT;
/*START Temperature minimum search */
/* THREAD 1 */
/*This is correcting for un-Earthly temperatures*/
@@ -618,6 +620,10 @@
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(input_dT->answer){
+ if(G_get_raster_row(infd_dT,inrast_dT,row,data_type_dT)<0)
+ G_fatal_error(_("Could not read from <%s>"),T);
+ }
/*process the data */
for (col=0; col < ncols; col++)
{
@@ -676,11 +682,26 @@
d_g0 = (double) ((DCELL *) inrast_g0)[col];
break;
}
+ if(input_dT->answer){
+ switch(data_type_dT){
+ case CELL_TYPE:
+ d_dT=(double) ((CELL *) inrast_dT)[col];
+ break;
+ case FCELL_TYPE:
+ d_dT=(double) ((FCELL *)inrast_dT)[col];
+ break;
+ case DCELL_TYPE:
+ d_dT=(double) ((DCELL *)inrast_dT)[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_g0)||
+ ((input_dT->answer)&&
+ (G_is_d_null_value(&d_dT)))){
/* do nothing */
}else{
d_t0dem = d_tempk + 0.001649*d_dem;
@@ -717,6 +738,8 @@
d_dem_dry=d_dem;
col_dry=col;
row_dry=row;
+ if(input_dT->answer)
+ d_dT_dry=d_dT;
}
if(flag1->answer&&
d_tempk>=(double)i_peak3-0.0&&
@@ -732,6 +755,8 @@
d_dem_dry=d_dem;
col_dry=col;
row_dry=row;
+ if(input_dT->answer)
+ d_dT_dry=d_dT;
}
}
}
@@ -750,6 +775,8 @@
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);
+ if(input_dT->answer)
+ G_message("dT_dry=%f\n",d_dT_dry);
} /* END OF FLAG2 */
/*MPI_BARRIER*/
@@ -771,6 +798,7 @@
DCELL d_tempk;
DCELL d_dem;
DCELL d_t0dem;
+ DCELL d_dT;
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)
@@ -779,6 +807,10 @@
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(input_dT->answer){
+ if(G_get_raster_row(infd_dT,inrast_dT,row,data_type_dT)<0)
+ G_fatal_error(_("Could not read from <%s>"),T);
+ }
switch(data_type_T){
case CELL_TYPE:
d_tempk = (double) ((CELL *) inrast_T)[col];
@@ -823,13 +855,28 @@
d_g0 = (double) ((DCELL *) inrast_g0)[col];
break;
}
+ if(input_dT->answer){
+ switch(data_type_dT){
+ case CELL_TYPE:
+ d_dT=(double) ((CELL *) inrast_dT)[col];
+ break;
+ case FCELL_TYPE:
+ d_dT=(double) ((FCELL *)inrast_dT)[col];
+ break;
+ case DCELL_TYPE:
+ d_dT=(double) ((DCELL *)inrast_dT)[col];
+ break;
+ }
+ }
d_t0dem = d_tempk + 0.001649 * 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;
-
+ if(input_dT->answer){
+ d_dT_dry = d_dT;
+ }
/*WET PIXEL*/
if(flag3->answer){
/*Calculate coordinates of row/col from projected ones*/
@@ -867,6 +914,8 @@
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);
+ if(input_dT->answer)
+ G_message("dT_dry=%f\n",d_dT_dry);
}
/* END OF MANUAL WET/DRY PIXELS */
@@ -1038,8 +1087,8 @@
G_set_d_null_value(&outrast[col],1);
} else {
/* Albedo < 0*/
- if(d_albedo<0.01){
- d_albedo=0.01;
+ if(d_albedo<0.001){
+ d_albedo=0.001;
}
/* Calculate T0dem */
d_t0dem = (double)d_tempk + 0.001649 * (double) d_dem;
@@ -1056,18 +1105,22 @@
if(input_dT->answer){
/* do nothing */
} else {
- /* let the internal constant take over */
+ /* let the internal constant
+ * take over */
d_dT = -1.0;
+ d_dT_dry = -1.0;
}
- d = sensi_h_z0m(iteration,d_tempk_wet,d_tempk_dry,d_t0dem,d_tempk,d_z0m,d_dT,d_dem,d_Rn_dry,d_g0_dry,d_t0dem_dry,d_u_hu,d_hu,d_dem_dry);
+ d = sensi_h_z0m(iteration,d_tempk_wet,d_tempk_dry,d_t0dem,d_tempk,d_z0m,d_dT,d_dem,d_Rn_dry,d_g0_dry,d_t0dem_dry,d_u_hu,d_hu,d_dem_dry,d_dT_dry);
} else {
if(input_dT->answer){
/* do nothing */
} else {
- /* let the internal constant take over */
+ /* let the internal constant
+ * take over */
d_dT = -1.0;
+ d_dT_dry = -1.0;
}
- d = sensi_h_noz0m(iteration,d_tempk_wet,d_tempk_dry,d_t0dem,d_tempk,d_dT,d_ndvi,d_ndvi_max,d_dem,d_Rn_dry,d_g0_dry,d_t0dem_dry,d_u_hu,d_hu,d_dem_dry);
+ d = sensi_h_noz0m(iteration,d_tempk_wet,d_tempk_dry,d_t0dem,d_tempk,d_dT,d_ndvi,d_ndvi_max,d_dem,d_Rn_dry,d_g0_dry,d_t0dem_dry,d_u_hu,d_hu,d_dem_dry,d_dT_dry);
}
// G_message(" d_h0=%5.3f",d);
if (zero->answer && d<0.0){
Modified: grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_noz0m.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_noz0m.c 2008-07-27 22:32:43 UTC (rev 32345)
+++ grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_noz0m.c 2008-07-27 23:10:29 UTC (rev 32346)
@@ -6,7 +6,10 @@
/* Arrays Declarations */
#define ITER_MAX 10
-double sensi_h_noz0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double dtair0, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert)
+/* if dtair_desert < 0 then internal calculations */
+/* if dtair0 < 0 then internal calculations */
+
+double sensi_h_noz0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double dtair0, double ndvi, double ndvi_max, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert, double dtair_desert)
{
/* Arrays Declarations */
double dtair[ITER_MAX], roh_air[ITER_MAX], rah[ITER_MAX];
@@ -17,7 +20,6 @@
int i, j, ic, debug=0;
double u_0, zom0;
double h_desert, rah_desert, roh_air_desert;
- double dtair_desert;
double psih, psim, psih_desert,psim_desert;
double ustar_desert,ustar_desertold,zom_desert;
double result;
@@ -33,12 +35,8 @@
printf("ndvi = %5.3f ndvimax = %5.3f\n",ndvi,ndvi_max);
printf("*****************************\n");
}
-// dtair[0] = dt_air_0(t0_dem, tempk_water, tempk_desert);
- if(dtair0 < 0){
- dtair[0] = 5.0;
- } else {
- dtair[0] = dtair0;
- }
+ if(t0_dem<250.0)
+ t0_dem=250.0;
// printf("*****************************dtair = %5.3f\n",dtair[0]);
roh_air[0] = roh_air_0(tempk);
// printf("*****************************rohair=%5.3f\n",roh_air[0]);
@@ -50,6 +48,21 @@
// printf("*****************************u0\n");
rah[0] = rah_0(zom0, u_hu, hu);
// printf("*****************************rah = %5.3f\n",rah[0]);
+ if(dtair_desert < 0.0){
+ h_desert= rnet_desert - g0_desert;
+ zom_desert= 0.002;
+ ustar_desert= u_star(t0_dem_desert,h_desert,u_0,roh_air_desert,zom_desert,u_hu,hu);
+ psih_desert= psi_h(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
+ psim_desert= psi_m(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
+ rah_desert = rah1(zom_desert,psih_desert,psim_desert,ustar_desert);
+// printf("*****************************rah_desert = %5.3f\n",rah_desert);
+ dtair_desert = dt_air_desert(h_desert, roh_air_desert, rah_desert);
+ }
+ if(dtair0 < 0.0){
+ dtair[0] = dt_air(t0_dem,tempk_water,tempk_desert,dtair_desert);
+ } else {
+ dtair[0] = dtair0;
+ }
h[0] = h_0(roh_air[0], rah[0], dtair[0]);
// printf("*****************************h\n");
if(debug==1){
@@ -69,26 +82,20 @@
if(debug==1){
printf("\n ******** ITERATION %i *********\n",ic);
}
- /* Where is roh_air[i]? */
ustar[ic] = u_star(t0_dem,h[ic-1],u_0,roh_air[ic-1],zom[0],u_hu,hu);
psih = psi_h(t0_dem,h[ic-1],ustar[ic],roh_air[ic-1],hu);
psim = psi_m(t0_dem,h[ic-1],ustar[ic],roh_air[ic-1],hu);
rah[ic] = rah1(zom0, psih, psim, ustar[ic]);
+ if(rah[ic]<0.0)
+ rah[ic]=0.0;
/* get desert point values from maps */
- if(ic==1){
- h_desert = rnet_desert - g0_desert;
- zom_desert = 0.002;
- ustar_desert = u_star(t0_dem_desert,h_desert,u_0,roh_air_desert,zom_desert,u_hu,hu);
- psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
- psim_desert = psi_m(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
- } 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;
- ustar_desert = u_star(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,zom_desert,u_hu,hu);
- psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
- }
- rah_desert = rah1(zom0,psih_desert,psim_desert,ustar_desert);
+ roh_air_desert = rohair(dem_desert,tempk_desert,dtair_desert);
+ h_desert = h1(roh_air_desert,rah_desert,dtair_desert);
+ ustar_desertold = ustar_desert;
+ ustar_desert = u_star(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,zom_desert,u_hu,hu);
+ psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
+ psim_desert = psi_m(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
+ rah_desert = rah1(zom_desert,psih_desert,psim_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);
Modified: grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_z0m.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_z0m.c 2008-07-27 22:32:43 UTC (rev 32345)
+++ grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_z0m.c 2008-07-27 23:10:29 UTC (rev 32346)
@@ -6,7 +6,10 @@
/* Arrays Declarations */
#define ITER_MAX 10
-double sensi_h_z0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double zom0, double dtair0, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert)
+/* if dtair_desert < 0 then internal calculations */
+/* if dtair0 < 0 then internal calculations */
+
+double sensi_h_z0m( int iteration, double tempk_water, double tempk_desert, double t0_dem, double tempk, double zom0, double dtair0, double dem, double rnet_desert, double g0_desert, double t0_dem_desert, double u_hu, double hu, double dem_desert,double dtair_desert)
{
/* Arrays Declarations */
double dtair[ITER_MAX], roh_air[ITER_MAX], rah[ITER_MAX];
@@ -17,7 +20,6 @@
int i, j, ic, debug=0;
double u_0;
double h_desert, rah_desert, roh_air_desert;
- double dtair_desert;
double psih,psim,psih_desert,psim_desert;
double ustar_desert,ustar_desertold,zom_desert;
double result;
@@ -33,21 +35,31 @@
printf("z0m0 = %5.3f\n",zom0);
printf("*****************************\n");
}
-// dtair[0] = dt_air_0(t0_dem, tempk_water, tempk_desert);
- if(dtair0 < 0){
- dtair[0] = 5.0;
- } else {
- dtair[0] = dtair0;
- }
-// printf("*****************************dtair = %5.3f\n",dtair[0]);
+ if(t0_dem<250.0)
+ t0_dem=250.0;
+// printf("****************************dtair = %5.3f\n",dtair[0]);
roh_air[0] = roh_air_0(tempk);
-// printf("*****************************rohair=%5.3f\n",roh_air[0]);
+// printf("***************************rohair=%5.3f\n",roh_air[0]);
roh_air_desert = roh_air_0(tempk_desert);
// printf("**rohairdesert = %5.3f\n",roh_air_desert);
u_0 = U_0(zom0, u_hu, hu);
// printf("*****************************u0\n");
rah[0] = rah_0(zom0, u_hu, hu);
// printf("*****************************rah = %5.3f\n",rah[0]);
+ if(dtair_desert < 0.0){
+ h_desert= rnet_desert - g0_desert;
+ zom_desert= 0.002;
+ ustar_desert= u_star(t0_dem_desert,h_desert,u_0,roh_air_desert,zom_desert,u_hu,hu);
+ psih_desert= psi_h(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
+ psim_desert= psi_m(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
+ rah_desert= rah1(zom0,psih_desert,psim_desert,ustar_desert);
+ dtair_desert= dt_air_desert(h_desert, roh_air_desert, rah_desert);
+ }
+ if(dtair0 < 0.0){
+ dtair[0] = dt_air(t0_dem,tempk_water,tempk_desert,dtair_desert);
+ } else {
+ dtair[0] = dtair0;
+ }
h[0] = h_0(roh_air[0], rah[0], dtair[0]);
// printf("*****************************h\n");
if(debug==1){
@@ -67,26 +79,19 @@
if(debug==1){
printf("\n ******** ITERATION %i *********\n",ic);
}
- /* Where is roh_air[i]? */
ustar[ic] = u_star(t0_dem,h[ic-1],u_0,roh_air[ic-1],zom[0],u_hu,hu);
- psih = psi_h(t0_dem,h[ic-1],ustar[ic],roh_air[ic-1],hu);
- psim = psi_m(t0_dem,h[ic-1],ustar[ic],roh_air[ic-1],hu);
- rah[ic] = rah1(zom0,psih,psim,ustar[ic]);
+ psih= psi_h(t0_dem,h[ic-1],ustar[ic],roh_air[ic-1],hu);
+ psim= psi_m(t0_dem,h[ic-1],ustar[ic],roh_air[ic-1],hu);
+ rah[ic] = rah1(zom0,psih,psim,ustar[ic]);
+ if(rah[ic]<0.0)
+ rah[ic]=0.0;
/* get desert point values from maps */
- if(ic==1){
- h_desert = rnet_desert - g0_desert;
- zom_desert = 0.002;
- ustar_desert = u_star(t0_dem_desert,h_desert,u_0,roh_air_desert,zom_desert,u_hu,hu);
- psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
- psim_desert = psi_m(t0_dem_desert,h_desert,ustar_desert,roh_air_desert,hu);
- } 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;
- ustar_desert = u_star(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,zom_desert,u_hu,hu);
- psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
- psim_desert = psi_m(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
- }
+ roh_air_desert = rohair(dem_desert,tempk_desert,dtair_desert);
+ h_desert = h1(roh_air_desert,rah_desert,dtair_desert);
+ ustar_desertold = ustar_desert;
+ ustar_desert = u_star(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,zom_desert,u_hu,hu);
+ psih_desert = psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
+ psim_desert = psi_m(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
rah_desert = rah1(zom0,psih_desert,psim_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...*/
More information about the grass-commit
mailing list