faster/cleaner (all messages appear when verbose level above std)

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@73520 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.eb.hsebal01/main.c

@@ -43,6 +43,8 @@ double **G_alloc_matrix(int rows, int cols)
     return m;
 }
 
+
+
 int main(int argc, char *argv[])
 {
     struct Cell_head cellhd;
@@ -85,8 +87,6 @@ int main(int argc, char *argv[])
     int rowDry = 0, colDry = 0, rowWet = 0, colWet = 0;
 
     /********************************/
-
-    /********************************/
     xp = yp;
     yp = xp;
     xmin = ymin;
@@ -135,6 +135,11 @@ int main(int argc, char *argv[])
     input_t0dem->description =
 	_("Name of altitude corrected surface temperature raster map [K]");
 
+    input_eact = G_define_standard_option(G_OPT_R_INPUT);
+    input_eact->key = "vapourpressureactual";
+    input_eact->description =
+	_("Name of the actual vapour pressure (e_act) map [KPa]");
+
     input_ustar = G_define_option();
     input_ustar->key = "frictionvelocitystar";
     input_ustar->type = TYPE_DOUBLE;
@@ -145,12 +150,6 @@ int main(int argc, char *argv[])
 	_("Value of the height independent friction velocity (u*) [m/s]");
     input_ustar->guisection = _("Parameters");
 
-    input_eact = G_define_standard_option(G_OPT_R_INPUT);
-    input_eact->key = "vapourpressureactual";
-    input_eact->required = YES;
-    input_eact->description =
-	_("Name of the actual vapour pressure (e_act) map [KPa]");
-
     input_row_wet = G_define_option();
     input_row_wet->key = "row_wet_pixel";
     input_row_wet->type = TYPE_DOUBLE;
@@ -209,7 +208,7 @@ int main(int argc, char *argv[])
 
     /*If automatic flag, just forget the rest of options */
     if (flag2->answer)
-	G_message(_("Automatic mode selected"));
+	G_verbose_message(_("Automatic mode selected"));
     /*If not automatic & all pixels locations in col/row given */
     else if (!flag2->answer &&
 	     input_row_wet->answer &&
@@ -222,9 +221,9 @@ int main(int argc, char *argv[])
 	m_col_dry = atof(input_col_dry->answer);
 	/*If pixels locations are in projected coordinates */
 	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);
+	    G_verbose_message(_("Manual wet/dry pixels in image coordinates"));
+	G_verbose_message(_("Wet Pixel=> x:%f y:%f"), m_col_wet, m_row_wet);
+	G_verbose_message(_("Dry Pixel=> x:%f y:%f"), m_col_dry, m_row_dry);
     }
     /*If not automatic & missing any of the pixel location, Fatal Error */
     else {
@@ -278,17 +277,29 @@ int main(int argc, char *argv[])
 
     /***************************************************/
     /* Allocate memory for temporary images            */
+
+    /***************************************************/
     double **d_Roh, **d_Rah;
+    double **d_z0m, **d_t0dem, **d_eact;
 
     if ((d_Roh = G_alloc_matrix(nrows, ncols)) == NULL)
-	G_message("Unable to allocate memory for temporary d_Roh image");
+	G_verbose_message("Unable to allocate memory for temporary d_Roh image");
     if ((d_Rah = G_alloc_matrix(nrows, ncols)) == NULL)
-	G_message("Unable to allocate memory for temporary d_Rah image");
+	G_verbose_message("Unable to allocate memory for temporary d_Rah image");
+    if ((d_z0m = G_alloc_matrix(nrows, ncols)) == NULL)
+	G_verbose_message("Unable to allocate memory for temporary d_z0m image");
+    if ((d_t0dem = G_alloc_matrix(nrows, ncols)) == NULL)
+	G_verbose_message("Unable to allocate memory for temporary d_t0dem image");
+    if ((d_eact = G_alloc_matrix(nrows, ncols)) == NULL)
+	G_verbose_message("Unable to allocate memory for temporary d_eact image");
 
     /***************************************************/
-
     /* MANUAL T0DEM WET/DRY PIXELS */
+    DCELL d_u5 = 0.0;
+    DCELL d_roh1 = 0.0;
+    DCELL d_rah1 = 0.0;
     DCELL d_Rn_dry = 0.0, d_g0_dry = 0.0;
+    DCELL d_Rah_dry = 0.0, d_Roh_dry = 0.0;
     DCELL d_t0dem_dry = 0.0, d_t0dem_wet = 0.0;
 
     if (flag2->answer) {
@@ -300,176 +311,187 @@ int main(int argc, char *argv[])
 
 	/*********************/
 	for (row = 0; row < nrows; row++) {
-	    DCELL d_t0dem;
-	    DCELL d_z0m;
-	    DCELL d_eact;
-
 	    G_percent(row, nrows, 2);
 	    Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 	    Rast_get_d_row(infd_z0m, inrast_z0m, row);
-	    Rast_get_d_row(infd_eact, inrast_eact, row);
 	    Rast_get_d_row(infd_Rn, inrast_Rn, row);
 	    Rast_get_d_row(infd_g0, inrast_g0, row);
 	    /*process the data */
 	    for (col = 0; col < ncols; col++) {
-		d_t0dem = ((DCELL *) inrast_t0dem)[col];
-		d_z0m = ((DCELL *) inrast_z0m)[col];
-		d_eact = ((DCELL *) inrast_eact)[col];
+		d_t0dem[row][col] = (double) ((DCELL *) inrast_t0dem)[col];
+		d_z0m[row][col] = (double) ((DCELL *) inrast_z0m)[col];
 		d_Rn = ((DCELL *) inrast_Rn)[col];
 		d_g0 = ((DCELL *) inrast_g0)[col];
-		if (Rast_is_d_null_value(&d_t0dem) ||
-		    Rast_is_d_null_value(&d_eact) ||
-		    Rast_is_d_null_value(&d_z0m) ||
-		    Rast_is_d_null_value(&d_Rn) ||
+		if (Rast_is_d_null_value(&d_Rn) ||
 		    Rast_is_d_null_value(&d_g0)) {
+			d_Roh[row][col] = -999.9;
+			d_Rah[row][col] = -999.9;
 		    /* do nothing */
 		}
 		else {
-		    if (d_t0dem <= 250.0) {
+		    if (d_t0dem[row][col] <= 250.0 || d_z0m[row][col] < 0.01) {
+			d_Roh[row][col] = -999.9;
+			d_Rah[row][col] = -999.9;
 			/* do nothing */
 		    }
 		    else {
 			d_h0 = d_Rn - d_g0;
-			if (d_t0dem < t0dem_min &&
+			d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+		        d_rah1 =(1/(d_u5*pow(0.41,2)))*log(5/d_z0m[row][col])*log(5/(d_z0m[row][col]*0.1));
+		        d_roh1 =((998-d_eact[row][col])/(d_t0dem[row][col]*2.87))+(d_eact[row][col]/(d_t0dem[row][col]*4.61));
+		        if (d_roh1 > 5)
+		            d_roh1 = 1.0;
+		        else
+		        d_roh1 =((1000 - 4.65) / (d_t0dem[row][col] * 2.87))+(4.65 / (d_t0dem[row][col] * 4.61));
+
+			d_Roh[row][col] = d_roh1;
+			d_Rah[row][col] = d_rah1;
+
+			if (d_t0dem[row][col] < t0dem_min &&
 			    d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 0.0 &&
-			    d_h0 < 100.0) {
-			    t0dem_min = d_t0dem;
-			    d_t0dem_wet = d_t0dem;
+			    d_h0 < 100.0 && d_roh1 > 0.001 && d_rah1 > 0.001) {
+			    t0dem_min = d_t0dem[row][col];
+			    d_t0dem_wet = d_t0dem[row][col];
 			    d_Rn_wet = d_Rn;
 			    d_g0_wet = d_g0;
-			    m_col_wet = col;
-			    m_row_wet = row;
+			    colWet = col;
+			    rowWet = row;
 			}
-			if (d_t0dem > t0dem_max &&
+			if (d_t0dem[row][col] > t0dem_max &&
 			    d_Rn > 0.0 && d_g0 > 0.0 && d_h0 > 100.0 &&
-			    d_h0 < 500.0) {
-			    t0dem_max = d_t0dem;
-			    d_t0dem_dry = d_t0dem;
+			    d_h0 < 500.0 && d_roh1 > 0.001 && d_rah1 > 0.001) {
+			    t0dem_max = d_t0dem[row][col];
+			    d_t0dem_dry = d_t0dem[row][col];
 			    d_Rn_dry = d_Rn;
 			    d_g0_dry = d_g0;
-			    m_col_dry = col;
-			    m_row_dry = row;
+			    colDry = col;
+			    rowDry = row;
+			    d_Roh_dry = d_roh1;
+			    d_Rah_dry = d_rah1;
 			}
 		    }
 		}
 	    }
 	}
-	G_message("row_wet=%d\tcol_wet=%d", row_wet, col_wet);
-	G_message("row_dry=%d\tcol_dry=%d", row_dry, col_dry);
-	G_message("g0_wet=%f", d_g0_wet);
-	G_message("Rn_wet=%f", d_Rn_wet);
-	G_message("LE_wet=%f", d_Rn_wet - d_g0_wet);
-	G_message("t0dem_dry=%f", d_t0dem_dry);
-	G_message("rnet_dry=%f", d_Rn_dry);
-	G_message("g0_dry=%f", d_g0_dry);
-	G_message("h0_dry=%f", d_Rn_dry - d_g0_dry);
-	G_message("auto config completed");
+	G_verbose_message("row_wet=%d\tcol_wet=%d", rowWet, colWet);
+	G_verbose_message("row_dry=%d\tcol_dry=%d", rowDry, colDry);
+	G_verbose_message("g0_wet=%f", d_g0_wet);
+	G_verbose_message("Rn_wet=%f", d_Rn_wet);
+	G_verbose_message("LE_wet=%f", d_Rn_wet - d_g0_wet);
+	G_verbose_message("t0dem_wet=%f", d_t0dem_wet);
+	G_verbose_message("t0dem_dry=%f", d_t0dem_dry);
+	G_verbose_message("rnet_dry=%f", d_Rn_dry);
+	G_verbose_message("g0_dry=%f", d_g0_dry);
+	G_verbose_message("h0_dry=%f", d_Rn_dry - d_g0_dry);
+	G_verbose_message("Rah_dry=%f", d_Rah_dry);
+	G_verbose_message("Roh_dry=%f", d_Roh_dry);
+	G_verbose_message("auto config completed");
     }	/* END OF FLAG2 */
 
-    G_message("Passed here");
-
     /* MANUAL T0DEM WET/DRY PIXELS */
     /*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);
+	G_verbose_message("Manual: Dry Pixel | row:%i col:%i", row, col);
+    	rowDry = row;
+    	colDry = col;
     }
     else {
-	row = (int)m_row_dry;
-	col = (int)m_col_dry;
-	G_message("Dry Pixel | row:%i col:%i", row, col);
+	row = rowDry;
+	col = colDry;
+	G_verbose_message("Dry Pixel | row:%i col:%i", row, col);
     }
-    rowDry = row;
-    colDry = col;
-    Rast_get_d_row(infd_Rn, inrast_Rn, row);
-    Rast_get_d_row(infd_g0, inrast_g0, row);
-    Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
-    d_Rn_dry = ((DCELL *) inrast_Rn)[col];
-    d_g0_dry = ((DCELL *) inrast_g0)[col];
-    d_t0dem_dry = ((DCELL *) inrast_t0dem)[col];
     /*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);
+	G_verbose_message("Manual: Wet Pixel | row:%i col:%i", row, col);
+    	rowWet = row;
+    	colWet = col;
     }
     else {
-	row = m_row_wet;
-	col = m_col_wet;
-	G_message("Wet Pixel | row:%i col:%i", row, col);
+	row = rowWet;
+	col = colWet;
+	G_verbose_message("Wet Pixel | row:%i col:%i", row, col);
     }
-    rowWet = row;
-    colWet = col;
-    Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
-    d_t0dem_wet = ((DCELL *) inrast_t0dem)[col];
     /* END OF MANUAL WET/DRY PIXELS */
-    double h_dry;
 
+    /* Extract end-members */
+    Rast_get_d_row(infd_Rn, inrast_Rn, rowDry);
+    Rast_get_d_row(infd_g0, inrast_g0, rowDry);
+    Rast_get_d_row(infd_t0dem, inrast_t0dem, rowDry);
+    d_Rn_dry = ((DCELL *) inrast_Rn)[colDry];
+    d_g0_dry = ((DCELL *) inrast_g0)[colDry];
+    d_t0dem_dry = ((DCELL *) inrast_t0dem)[colDry];
+
+    Rast_get_d_row(infd_t0dem, inrast_t0dem, rowWet);
+    d_t0dem_wet = ((DCELL *) inrast_t0dem)[colWet];
+
+    double h_dry;
     h_dry = d_Rn_dry - d_g0_dry;
-    G_message("h_dry = %f", h_dry);
-    G_message("t0dem_dry = %f", d_t0dem_dry);
-    G_message("t0dem_wet = %f", d_t0dem_wet);
+    G_verbose_message("h_dry = %f", h_dry);
+    G_verbose_message("t0dem_dry = %f", d_t0dem_dry);
+    G_verbose_message("t0dem_wet = %f", d_t0dem_wet);
 
-    DCELL d_rah_dry = 0.0;
-    DCELL d_roh_dry = 0.0;
+    DCELL d_rah_dry = d_Rah_dry;
+    DCELL d_roh_dry = d_Roh_dry;
 
-    DCELL d_t0dem, d_z0m, d_eact;
-    DCELL d_u5;
-    DCELL d_roh1;
     DCELL d_h1, d_h2, d_h3;
-    DCELL d_rah1, d_rah2, d_rah3;
+    DCELL d_rah2, d_rah3;
     DCELL d_L, d_x, d_psih, d_psim;
 
     /* INITIALIZATION */
+    /******************************************************/
+    /*If d_rah_dry and d_roh_dry are not auto found, then:*/
+    if(d_Rah_dry==0.0 && d_Roh_dry==0.0){
+    /***************************************************/
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 2);
-	/* read a line input maps into buffers */
-	Rast_get_d_row(infd_z0m, inrast_z0m, row);
-	Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
-	Rast_get_d_row(infd_eact, inrast_eact, row);
 	/* read every cell in the line buffers */
 	for (col = 0; col < ncols; col++) {
-	    d_z0m = ((DCELL *) inrast_z0m)[col];
-	    d_t0dem = ((DCELL *) inrast_t0dem)[col];
-	    d_eact = ((DCELL *) inrast_eact)[col];
-	    if (Rast_is_d_null_value(&d_t0dem) ||
-		Rast_is_d_null_value(&d_eact) ||
-		Rast_is_d_null_value(&d_z0m)) {
+	    d_eact[row][col] = (double) ((DCELL *) inrast_eact)[col];
+	    d_z0m[row][col] = (double) ((DCELL *) inrast_z0m)[col];
+	    if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+		Rast_is_d_null_value(&d_eact[row][col]) ||
+		Rast_is_d_null_value(&d_z0m[row][col])) {
 		Rast_set_d_null_value(&outrast[col], 1);
 		d_Roh[row][col] = -999.9;
 		d_Rah[row][col] = -999.9;
 		if (row == rowDry && col == colDry) {	/*collect dry pix info */
-		    d_rah_dry = d_rah1;
-		    d_roh_dry = d_roh1;
-		    G_message("d_rah_dry=%f d_roh_dry=%f",d_rah_dry,d_roh_dry);
+		    d_rah_dry = d_Rah[row][col];
+		    d_roh_dry = d_Roh[row][col];
+		    G_verbose_message("Init: d_rah_dry=%f d_roh_dry=%f",d_rah_dry,d_roh_dry);
 		}
 	    }
 	    else {
-		d_u5 = (ustar / 0.41) * log(5 / d_z0m);
-		d_rah1 =
-		    (1 / (d_u5 * pow(0.41, 2))) * log(5 / d_z0m) * log(5/(d_z0m*0.1));
-		d_roh1 =
-		    ((998 - d_eact) / (d_t0dem * 2.87)) + (d_eact / (d_t0dem * 4.61));
+		d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+		d_rah1 =(1/(d_u5*pow(0.41,2)))*log(5/d_z0m[row][col])*log(5/(d_z0m[row][col]*0.1));
+		d_roh1 =((998-d_eact[row][col])/(d_t0dem[row][col]*2.87))+(d_eact[row][col]/(d_t0dem[row][col]*4.61));
 		if (d_roh1 > 5)
 		    d_roh1 = 1.0;
 		else
-		    d_roh1 =
-			((1000 - 4.65) / (d_t0dem * 2.87)) +
-			(4.65 / (d_t0dem * 4.61));
-		if (row == rowDry && col == colDry) {	/*collect dry pix info */
+		    d_roh1 =((1000 - 4.65) / (d_t0dem[row][col] * 2.87))+(4.65 / (d_t0dem[row][col] * 4.61));
+		if (row == rowDry && col == colDry) {	
+		    /*collect dry pix info */
 		    d_rah_dry = d_rah1;
 		    d_roh_dry = d_roh1;
-		    G_message("d_rah_dry=%f d_roh_dry=%f", d_rah_dry,
-			      d_roh_dry);
+		    G_verbose_message("row=%d col=%d", row, col);
+		    G_verbose_message("ustar=%f", ustar);
+		    G_verbose_message("d_u5=%f", d_u5);
+		    G_verbose_message("d_t0dem_dry=%f", d_t0dem[row][col]);
+		    G_verbose_message("d_z0m_dry=%f", d_z0m[row][col]);
+		    G_verbose_message("d_rah_dry=%f", d_rah_dry);
+		    G_verbose_message("d_roh_dry=%f", d_roh_dry);
 		}
 		d_Roh[row][col] = d_roh1;
 		d_Rah[row][col] = d_rah1;
 	    }
 	}
     }
+    } /*END OF if !d_Rah_dry and !d_Roh_dry*/
     DCELL d_dT_dry;
 
     /*Calculate dT_dry */
@@ -486,78 +508,91 @@ int main(int argc, char *argv[])
 
     a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
     b = (sumy - (a * sumx)) / 2.0;
-    G_message("d_dT_dry=%f", d_dT_dry);
-    G_message("dT1=%f * t0dem + (%f)", a, b);
+    G_verbose_message("d_dT_dry=%f", d_dT_dry);
+    G_verbose_message("dT1=%f * t0dem + (%f)", a, b);
     DCELL d_h_dry = 0.0;
     if(isnan(a) || isnan(b)){
-	    G_fatal_error("Delta T Convergence failed, exiting prematurily, please check output");
+    	G_free(outrast);
+    	Rast_close(outfd);
+	G_fatal_error("Delta T Convergence failed, exiting prematurily, please check output");
     }
 
     /* ITERATION 1 */
+    /***************************************************/
+    /***************************************************/
+    //outfd = Rast_open_old(h0, "");
+    //Rast_get_cellhd(h0, "", &cellhd);
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 2);
-	/* read a line input maps into buffers */
-	Rast_get_d_row(infd_z0m, inrast_z0m, row);
-	Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 	/* read every cell in the line buffers */
 	for (col = 0; col < ncols; col++) {
-	    d_z0m = ((DCELL *) inrast_z0m)[col];
-	    d_t0dem = ((DCELL *) inrast_t0dem)[col];
 	    d_rah1 = d_Rah[row][col];
 	    d_roh1 = d_Roh[row][col];
-	    if (Rast_is_d_null_value(&d_t0dem) ||
-		Rast_is_d_null_value(&d_z0m)) {
+	    if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+		Rast_is_d_null_value(&d_z0m[row][col])) {
 		Rast_set_d_null_value(&outrast[col], 1);
 	    }
 	    else {
 		if (d_rah1 < 1.0)
 		    d_h1 = 0.0;
 		else
-		    d_h1 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah1;
+		    d_h1 = (1004 * d_roh1) * (a * d_t0dem[row][col] + b) / d_rah1;
 		if (d_h1 < 0 && d_h1 > -50) {
 		    d_h1 = 0.0;
 		}
 		if (d_h1 < -50 || d_h1 > 1000) {
 		    Rast_set_d_null_value(&outrast[col], 1);
+		}else{
+			outrast[col] = d_h1;
+			d_L =-1004*d_roh1*pow(ustar,3)*d_t0dem[row][col]/(d_h1*9.81*0.41);
+			d_x = pow((1 - 16 * (5 / d_L)), 0.25);
+			d_psim =2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
+		    		2 * atan(d_x) + 0.5 * M_PI;
+			d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
+			d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+			d_rah2 =(1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m[row][col]) - d_psim)
+		    	* log((5 / (d_z0m[row][col] * 0.1)) - d_psih);
+			if (row == rowDry && col == colDry) {	/*collect dry pix info */
+		    		d_h1 = (1004 * d_roh1) * (a * d_t0dem[row][col] + b) / d_rah_dry;
+				d_L =-1004*d_roh1*pow(ustar,3)*d_t0dem[row][col]/(d_h1*9.81*0.41);
+				d_x = pow((1 - 16 * (5 / d_L)), 0.25);
+				d_psim =2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
+		    			2 * atan(d_x) + 0.5 * M_PI;
+				d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
+				d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
+				d_rah2 =(1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m[row][col]) - d_psim)
+		    			* log((5 / (d_z0m[row][col] * 0.1)) - d_psih);
+		    		d_rah_dry = d_rah2;
+		    		d_h_dry = d_h1;
+				G_verbose_message("d_z0m (dry)=%f",d_z0m[row][col]);
+				G_verbose_message("d_rah1 (dry)=%f",d_rah1);
+				G_verbose_message("d_rah2 (dry)=%f",d_rah2);
+				G_verbose_message("d_h1 (dry)=%f",d_h1);
+			}
+			d_Rah[row][col] = d_rah1;
 		}
-		outrast[col] = d_h1;
-		d_L =
-		    -1004 * d_roh1 * pow(ustar,
-		    3) * d_t0dem / (d_h1 * 9.81 * 0.41);
-		d_x = pow((1 - 16 * (5 / d_L)), 0.25);
-		d_psim =
-		    2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
-		    2 * atan(d_x) + 0.5 * M_PI;
-		d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
-		d_u5 = (ustar / 0.41) * log(5 / d_z0m);
-		d_rah2 =
-		    (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m) - d_psim)
-		    * log((5 / (d_z0m * 0.1)) - d_psih);
-		if (row == rowDry && col == colDry) {	/*collect dry pix info */
-		    d_rah_dry = d_rah2;
-		    d_h_dry = d_h1;
-		}
-		d_Rah[row][col] = d_rah1;
 	    }
 	}
-	Rast_put_d_row(outfd, outrast);
     }
 
     /*Calculate dT_dry */
+    G_verbose_message("d_h_dry=%f",d_h_dry);
+    G_verbose_message("d_rah_dry=%f",d_rah_dry);
+    G_verbose_message("d_roh_dry=%f",d_roh_dry);
     d_dT_dry = (d_h_dry * d_rah_dry) / (1004 * d_roh_dry);
     /*Calculate coefficients for next dT equation */
     /*      a = (d_dT_dry-0)/(d_t0dem_dry-d_t0dem_wet); */
     /*      b = (-1.0) * ( a * d_t0dem_wet ); */
-    /*      G_message("d_dT_dry=%f",d_dT_dry); */
-    /*      G_message("dT2=%f * t0dem + (%f)", a, b); */
+    /*      G_verbose_message("d_dT_dry=%f",d_dT_dry); */
+    /*      G_verbose_message("dT2=%f * t0dem + (%f)", a, b); */
     sumx = d_t0dem_wet + d_t0dem_dry;
     sumy = d_dT_dry + 0.0;
     sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
     sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
     a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
     b = (sumy - (a * sumx)) / 2.0;
-    G_message("d_dT_dry=%f", d_dT_dry);
-    G_message("dT1=%f * t0dem + (%f)", a, b);
+    G_verbose_message("d_dT_dry=%f", d_dT_dry);
+    G_verbose_message("dT2=%f * t0dem + (%f)", a, b);
     if(isnan(a) || isnan(b)){
     	G_free(outrast);
     	Rast_close(outfd);
@@ -565,23 +600,16 @@ int main(int argc, char *argv[])
     }
 
     /* ITERATION 2 */
-
     /***************************************************/
-
     /***************************************************/
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 2);
-	/* read a line input maps into buffers */
-	Rast_get_d_row(infd_z0m, inrast_z0m, row);
-	Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 	/* read every cell in the line buffers */
 	for (col = 0; col < ncols; col++) {
-	    d_z0m = ((DCELL *) inrast_z0m)[col];
-	    d_t0dem = ((DCELL *) inrast_t0dem)[col];
 	    d_rah2 = d_Rah[row][col];
 	    d_roh1 = d_Roh[row][col];
-	    if (Rast_is_d_null_value(&d_t0dem) ||
-		Rast_is_d_null_value(&d_z0m)) {
+	    if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+		Rast_is_d_null_value(&d_z0m[row][col])) {
 		Rast_set_d_null_value(&outrast[col], 1);
 	    }
 	    else {
@@ -589,7 +617,7 @@ int main(int argc, char *argv[])
 		    d_h2 = 0.0;
 		}
 		else {
-		    d_h2 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah2;
+		    d_h2 = (1004 * d_roh1) * (a * d_t0dem[row][col] + b) / d_rah2;
 		}
 		if (d_h2 < 0 && d_h2 > -50) {
 		    d_h2 = 0.0;
@@ -600,15 +628,15 @@ int main(int argc, char *argv[])
 		outrast[col] = d_h2;
 		d_L =
 		    -1004 * d_roh1 * pow(ustar,
-		    3) * d_t0dem / (d_h2 * 9.81 * 0.41);
+		    3) * d_t0dem[row][col] / (d_h2 * 9.81 * 0.41);
 		d_x = pow((1 - 16 * (5 / d_L)), 0.25);
 		d_psim = 2 * log((1 + d_x) / 2) + log((1 + pow(d_x, 2)) / 2) -
 		    2 * atan(d_x) + 0.5 * M_PI;
 		d_psih = 2 * log((1 + pow(d_x, 2)) / 2);
-		d_u5 = (ustar / 0.41) * log(5 / d_z0m);
+		d_u5 = (ustar / 0.41) * log(5 / d_z0m[row][col]);
 		d_rah3 =
-		    (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m) -
-		    d_psim) * log((5 / (d_z0m * 0.1)) - d_psih);
+		    (1 / (d_u5 * pow(0.41, 2))) * log((5 / d_z0m[row][col]) -
+		    d_psim) * log((5 / (d_z0m[row][col] * 0.1)) - d_psih);
 		if (row == rowDry && col == colDry) {	/*collect dry pix info */
 		    d_rah_dry = d_rah3;
 		    d_h_dry = d_h2;
@@ -616,7 +644,6 @@ int main(int argc, char *argv[])
 		d_Rah[row][col] = d_rah2;
 	    }
 	}
-	Rast_put_d_row(outfd, outrast);
     }
 
     /*Calculate dT_dry */
@@ -624,16 +651,16 @@ int main(int argc, char *argv[])
     /*Calculate coefficients for next dT equation */
     /*      a = (d_dT_dry-0)/(d_t0dem_dry-d_t0dem_wet); */
     /*      b = (-1.0) * ( a * d_t0dem_wet ); */
-    /*      G_message("d_dT_dry=%f",d_dT_dry); */
-    /*      G_message("dT3=%f * t0dem + (%f)", a, b); */
+    /*      G_verbose_message("d_dT_dry=%f",d_dT_dry); */
+    /*      G_verbose_message("dT3=%f * t0dem + (%f)", a, b); */
     sumx = d_t0dem_wet + d_t0dem_dry;
     sumy = d_dT_dry + 0.0;
     sumx2 = pow(d_t0dem_wet, 2) + pow(d_t0dem_dry, 2);
     sumxy = (d_t0dem_wet * 0.0) + (d_t0dem_dry * d_dT_dry);
     a = (sumxy - ((sumx * sumy) / 2.0)) / (sumx2 - (pow(sumx, 2) / 2.0));
     b = (sumy - (a * sumx)) / 2.0;
-    G_message("d_dT_dry=%f", d_dT_dry);
-    G_message("dT1=%f * t0dem + (%f)", a, b);
+    G_verbose_message("d_dT_dry=%f", d_dT_dry);
+    G_verbose_message("dT3=%f * t0dem + (%f)", a, b);
     if(isnan(a) || isnan(b)){
     	G_free(outrast);
     	Rast_close(outfd);
@@ -641,24 +668,15 @@ int main(int argc, char *argv[])
     }
 
     /* ITERATION 3 */
-
     /***************************************************/
-
     /***************************************************/
-
     for (row = 0; row < nrows; row++) {
 	G_percent(row, nrows, 2);
-	/* read a line input maps into buffers */
-	Rast_get_d_row(infd_z0m, inrast_z0m, row);
-	Rast_get_d_row(infd_t0dem, inrast_t0dem, row);
 	/* read every cell in the line buffers */
 	for (col = 0; col < ncols; col++) {
-	    d_z0m = ((DCELL *) inrast_z0m)[col];
-	    d_t0dem = ((DCELL *) inrast_t0dem)[col];
 	    d_rah3 = d_Rah[row][col];
-	    d_roh1 = d_Roh[row][col];
-	    if (Rast_is_d_null_value(&d_t0dem) ||
-		Rast_is_d_null_value(&d_z0m)) {
+	    if (Rast_is_d_null_value(&d_t0dem[row][col]) ||
+		Rast_is_d_null_value(&d_z0m[row][col])) {
 		Rast_set_d_null_value(&outrast[col], 1);
 	    }
 	    else {
@@ -666,21 +684,22 @@ int main(int argc, char *argv[])
 		    d_h3 = 0.0;
 		}
 		else {
-		    d_h3 = (1004 * d_roh1) * (a * d_t0dem + b) / d_rah3;
+		    d_h3 = (1004 * d_Roh[row][col]) * (a * d_t0dem[row][col] + b) / d_rah3;
 		}
 		if (d_h3 < 0 && d_h3 > -50) {
 		    d_h3 = 0.0;
 		}
 		if (d_h3 < -50 || d_h3 > 1000) {
 		    Rast_set_d_null_value(&outrast[col], 1);
-		}
+		} else {
 		outrast[col] = d_h3;
+		}
 	    }
 	}
 	Rast_put_d_row(outfd, outrast);
     }
-
-
+    G_free(inrast_eact);
+    Rast_close(infd_eact);
     G_free(inrast_z0m);
     Rast_close(infd_z0m);
     G_free(inrast_t0dem);