[GRASS-SVN] r32297 - grass-addons/gipe/i.eb.h_SEBAL95

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Jul 25 10:20:56 EDT 2008


Author: ychemin
Date: 2008-07-25 10:20:56 -0400 (Fri, 25 Jul 2008)
New Revision: 32297

Added:
   grass-addons/gipe/i.eb.h_SEBAL95/psi_m.c
Modified:
   grass-addons/gipe/i.eb.h_SEBAL95/psi_h.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/sensi_h_noz0m.c
   grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_z0m.c
Log:
reworking of psichrometric dependencies of aerodynamic resitance to momentum heat

Modified: grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c	2008-07-25 13:51:00 UTC (rev 32296)
+++ grass-addons/gipe/i.eb.h_SEBAL95/psi_h.c	2008-07-25 14:20:56 UTC (rev 32297)
@@ -4,29 +4,23 @@
 
 #define PI 3.1415927
 
-double psi_h(double t0_dem, double h, double U_0, double roh_air)
+double psi_h(double t0_dem, double h, double ustar, double roh_air, double hu)
 {
-	double result;
-	double n5_temp, n11_mem, n12_mem;
-
+	double result, n5, n11, n12;
 //	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);
+		n5 = (-1004* roh_air*pow(ustar,3)* t0_dem)/(0.41*9.81* h);
 	} else {
-		n5_temp = -1000.0;
+		n5 = -1000.0;
 	}
-
-	if(n5_temp < 0.0){
-		n12_mem = pow((1-16*(2/n5_temp)),0.25);
-		n11_mem = (2*log((1+pow(n12_mem,2))/2));
+	if(n5 < 0.0){
+		n12 = pow((1-16*(hu/n5)),0.25);
+		n11 = (2*log((1+pow(n12,2))/2));
 	} else {
-		n12_mem = 1.0;
-		n11_mem = -5*2/n5_temp;
+		n12 = 1.0;
+		n11 = -5*hu/n5;
 	}
-
-	result = n11_mem;
-
-	return (result);
+	result = n11;
+	return result;
 }
 

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-07-25 14:20:56 UTC (rev 32297)
@@ -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
   + *

Modified: grass-addons/gipe/i.eb.h_SEBAL95/rah1.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/rah1.c	2008-07-25 13:51:00 UTC (rev 32296)
+++ grass-addons/gipe/i.eb.h_SEBAL95/rah1.c	2008-07-25 14:20:56 UTC (rev 32297)
@@ -2,14 +2,11 @@
 #include<math.h>
 #include"functions.h"
 
-#define PI 3.1415927
-
-double rah1(double psih, double ustar, double hu)
+double rah1(double zom_0, double psih, double psim, double ustar)
 {
-	double result;
-
-	result = (log(hu/0.01)-psih)/(ustar*0.41);
-
-	return (result);
+	double  u5, result;
+	u5 = (ustar/0.41)*log(5/zom_0);
+	result = (1/(pow(u5*0.41,2))) * (log(5/zom_0)-psim) * (log(5/(zom_0*0.1))-psih);
+	return result;
 }
 

Modified: grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c
===================================================================
--- grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c	2008-07-25 13:51:00 UTC (rev 32296)
+++ grass-addons/gipe/i.eb.h_SEBAL95/rah_0.c	2008-07-25 14:20:56 UTC (rev 32297)
@@ -2,12 +2,12 @@
 #include<math.h>
 #include"functions.h"
 
-double rah_0(double zom_0, double u_0, double hu)
+double rah_0(double zom_0, double u_hu, double hu)
 {
-	double result;
-
-	result = log(hu/0.01)/(u_0*0.41);
-
-	return (result);
+	double ustar, u5, result;
+	ustar = (u_hu*0.41)/(log(hu/zom_0));
+	u5 = (ustar/0.41)*log(5/zom_0);
+	result = (1/(pow(u5*0.41,2))) * log(5/zom_0) * log(5/(zom_0*0.1));
+	return result;
 }
 

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-25 13:51:00 UTC (rev 32296)
+++ grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_noz0m.c	2008-07-25 14:20:56 UTC (rev 32297)
@@ -18,8 +18,8 @@
 	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 psih, psim, psih_desert,psim_desert;
+	double ustar_desert,ustar_desertold,zom_desert;
 	double result;
 
 	/* Fat-free junk food */
@@ -48,7 +48,7 @@
 // 	printf("*****************************zom = %5.3f\n",zom0);
 	u_0 		= U_0(zom0, u_hu, hu);
 // 	printf("*****************************u0\n");
-	rah[0] 		= rah_0(zom0, u_0, hu);
+	rah[0] 		= rah_0(zom0, u_hu, hu);
 // 	printf("*****************************rah = %5.3f\n",rah[0]);
 	h[0] 		= h_0(roh_air[0], rah[0], dtair[0]);
 // 	printf("*****************************h\n");
@@ -70,23 +70,25 @@
 			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],u_hu,hu);
-		rah[ic] = rah1(psih, ustar[ic], 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]);	
 		/* 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,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;
-			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,u_hu,hu);
+			psih_desert 	= psi_h(t0_dem_desert,h_desert,ustar_desertold,roh_air_desert,hu);
 		}
-		rah_desert	= rah1(psih_desert,ustar_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...*/
 		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-25 13:51:00 UTC (rev 32296)
+++ grass-addons/gipe/i.eb.h_SEBAL95/sensi_h_z0m.c	2008-07-25 14:20:56 UTC (rev 32297)
@@ -18,8 +18,8 @@
 	double u_0;
 	double h_desert, rah_desert, roh_air_desert;
 	double dtair_desert;
-	double psih_desert,ustar_desert,ustar_desertold,zom_desert;
-	double psih;
+	double psih,psim,psih_desert,psim_desert;
+	double ustar_desert,ustar_desertold,zom_desert;
 	double result;
 
 	/* Fat-free junk food */
@@ -46,7 +46,7 @@
 // 	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_0, hu);
+	rah[0] 		= rah_0(zom0, u_hu, hu);
 // 	printf("*****************************rah = %5.3f\n",rah[0]);
 	h[0] 		= h_0(roh_air[0], rah[0], dtair[0]);
 // 	printf("*****************************h\n");
@@ -68,23 +68,26 @@
 			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],u_hu,hu);
-		rah[ic] = rah1(psih, ustar[ic],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]);	
 		/* 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,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;
-			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,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(psih_desert,ustar_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...*/
 		dtair[ic] 	= dt_air(t0_dem, tempk_water, tempk_desert, dtair_desert);



More information about the grass-commit mailing list