From 320e246aa0535f9cb375fdc17942d462a6584ada Mon Sep 17 00:00:00 2001
From: Zac Flamig <zac.flamig@noaa.gov>
Date: Fri, 26 May 2017 15:53:39 +0000
Subject: [PATCH 1/3] grid averaging code

---
 src/BasicGrids.cpp     |   8 +++
 src/GridNode.h         |   2 +
 src/KinematicRoute.cpp |  17 ++++--
 src/Simulator.cpp      | 118 ++++++++++++++++++++++++++++++++++++++++-
 src/Simulator.h        |   3 +-
 src/TimeSeries.cpp     |   2 +-
 6 files changed, 143 insertions(+), 7 deletions(-)

diff --git a/src/BasicGrids.cpp b/src/BasicGrids.cpp
index cc9d190..a4e4695 100755
--- a/src/BasicGrids.cpp
+++ b/src/BasicGrids.cpp
@@ -844,6 +844,8 @@ void CarveBasin(BasinConfigSection *basin, std::vector<GridNode> *nodes, std::ma
 		}
 		currentN->area = g_Projection->GetArea(currentN->refLoc.x, currentN->refLoc.y);
 		currentN->contribArea = currentN->area;
+    currentN->relief = g_DEM->data[currentN->y][currentN->x];
+		currentN->riverLen = currentN->horLen;
 		currentN->fac = g_FAM->data[currentN->y][currentN->x];
 		walkNodes.push(currentN);
 		
@@ -934,6 +936,8 @@ void CarveBasin(BasinConfigSection *basin, std::vector<GridNode> *nodes, std::ma
 					nextN->slope = ((DEMDiff < 1.0) ? 1.0 : DEMDiff) / nextN->horLen;
 					nextN->area = g_Projection->GetArea(nextN->refLoc.x, nextN->refLoc.y);
 					nextN->contribArea = nextN->area;
+					nextN->relief = g_DEM->data[nextN->y][nextN->x];
+					nextN->riverLen = nextN->horLen;
 					//printf("Pushing node %i %i (%i, %i) from %i %i (%i, %i) %i %i\n", nextN->x, nextN->y, g_DDM->data[nextN->y][nextN->x], g_FAM->data[nextN->y][nextN->x], currentN->x, currentN->y, g_DDM->data[currentN->y][currentN->x], g_FAM->data[currentN->y][currentN->x], currentNode, nodes->size());
 					walkNodes.push(nextN);
 				}
@@ -951,6 +955,10 @@ void CarveBasin(BasinConfigSection *basin, std::vector<GridNode> *nodes, std::ma
 		GridNode *node = &nodes->at(i);
 		if (node->downStreamNode != INVALID_DOWNSTREAM_NODE) {	
 			nodes->at(node->downStreamNode).contribArea += node->contribArea;
+			nodes->at(node->downStreamNode).riverLen += node->riverLen;
+      if (nodes->at(node->downStreamNode).relief < node->relief) {
+				nodes->at(node->downStreamNode).relief = node->relief;
+			}
 		}
 	}
 
diff --git a/src/GridNode.h b/src/GridNode.h
index ba592bb..12c6079 100755
--- a/src/GridNode.h
+++ b/src/GridNode.h
@@ -17,6 +17,8 @@ struct GridNode {
   float area;
 	float contribArea;
   float horLen;
+	float riverLen;
+	float relief;
   bool channelGridCell;
   GaugeConfigSection *gauge;
   unsigned long downStreamNode;
diff --git a/src/KinematicRoute.cpp b/src/KinematicRoute.cpp
index fda8860..3d2eca9 100755
--- a/src/KinematicRoute.cpp
+++ b/src/KinematicRoute.cpp
@@ -387,6 +387,15 @@ void KWRoute::InitializeParameters(std::map<GaugeConfigSection *, float *> *para
   //This pass distributes parameters
   size_t numNodes = nodes->size();
 	size_t unused = 0;
+
+	if (paramGrids->at(PARAM_KINEMATIC_ALPHA) && g_DEM->IsSpatialMatch(paramGrids->at(PARAM_KINEMATIC_ALPHA))) {
+      printf("Alpha grid is match!\n");
+    }
+
+  if (paramGrids->at(PARAM_KINEMATIC_BETA) && g_DEM->IsSpatialMatch(paramGrids->at(PARAM_KINEMATIC_BETA))) {
+      printf("Beta grid is match!\n");
+    }
+
 	for (size_t i = 0; i < numNodes; i++) {
 		GridNode *node = &nodes->at(i);
 		KWGridNode *cNode = &(kwNodes[i]);
@@ -402,8 +411,8 @@ void KWRoute::InitializeParameters(std::map<GaugeConfigSection *, float *> *para
 		}
 		cNode->incomingWater[KW_LAYER_INTERFLOW] = 0.0;
     cNode->incomingWater[KW_LAYER_FASTFLOW] = 0.0;
-    
-		// Deal with the distributed parameters here
+		
+    // Deal with the distributed parameters here
 		GridLoc pt;
 		for (size_t paramI = 0; paramI < PARAM_KINEMATIC_QTY; paramI++) {
 			FloatGrid *grid = paramGrids->at(paramI);
@@ -430,7 +439,7 @@ void KWRoute::InitializeParameters(std::map<GaugeConfigSection *, float *> *para
     }
     
     if (cNode->params[PARAM_KINEMATIC_ALPHA] < 0.0) {
-      //printf("Node Alpha(%f) is less than 0, setting to 1.\n", cNode->params[PARAM_KINEMATIC_ALPHA]);
+      printf("Node Alpha(%f) is less than 0, setting to 1.\n", cNode->params[PARAM_KINEMATIC_ALPHA]);
       cNode->params[PARAM_KINEMATIC_ALPHA] = 1.0;
     }
 
@@ -446,7 +455,7 @@ void KWRoute::InitializeParameters(std::map<GaugeConfigSection *, float *> *para
 
     
     
-    if (node->fac > cNode->params[PARAM_KINEMATIC_TH]) {
+    if (node->contribArea > cNode->params[PARAM_KINEMATIC_TH]) {
 			node->channelGridCell = true;
       cNode->channelGridCell = true;
     } else {
diff --git a/src/Simulator.cpp b/src/Simulator.cpp
index 2a36d43..e968ce7 100755
--- a/src/Simulator.cpp
+++ b/src/Simulator.cpp
@@ -1,6 +1,7 @@
 #include <cstdio>
 #include <cstring>
 #include <sys/stat.h>
+#include <dirent.h>
 #if _OPENMP
 #include <omp.h>
 #endif
@@ -599,6 +600,121 @@ void Simulator::CleanUp() {
 }
 
 void Simulator::BasinAvg() {
+	char buffer[CONFIG_MAX_LEN*2];
+  std::vector<float> avgVals, areaVals, fileVals;
+  long numNodes = nodes.size();
+  avgVals.resize(numNodes);
+	areaVals.resize(numNodes);
+	fileVals.resize(numNodes);
+  gridWriter.Initialize();
+
+	DIR *dir;
+	struct dirent *ent;
+
+	INFO_LOGF("Running basin averaging over files in output folder %s", outputPath);
+
+
+	// Compute & Output these first so they get included in the averaging process
+  for (long i = numNodes - 1; i >= 0; i--) {
+      avgVals[i] = nodes[i].relief - g_DEM->data[nodes[i].y][nodes[i].x];
+  }
+  sprintf(buffer, "%s/relief.tif", outputPath);
+  gridWriter.WriteGrid(&nodes, &avgVals, buffer, false);
+
+  for (long i = numNodes - 1; i >= 0; i--) {
+      avgVals[i] = avgVals[i] / nodes[i].riverLen ;
+  }
+  sprintf(buffer, "%s/relief.ratio.tif", outputPath);
+  gridWriter.WriteGrid(&nodes, &avgVals, buffer, false);
+
+
+	if ((dir = opendir(outputPath)) == NULL) {
+		ERROR_LOGF("%s", "Failed to open output directory for reading files to average!");
+		return;
+	}
+
+	while ((ent = readdir(dir)) != NULL) {
+		if (ent->d_name[0] == '.') {
+			continue;
+		}
+		sprintf(buffer, "%s/%s", outputPath, ent->d_name);
+		INFO_LOGF("Averaging file %s", buffer);
+		FloatGrid *fileGrid = ReadFloatTifGrid(buffer);
+		if (!fileGrid) {
+			continue;
+		}
+
+		if (g_DEM->IsSpatialMatch(fileGrid)) {
+    // The grids are the same! Our life is easy!
+    #pragma omp parallel for
+    for (size_t i = 0; i < nodes.size(); i++) {
+      GridNode *node = &(nodes[i]);
+      if (fileGrid->data[node->y][node->x] != fileGrid->noData) {
+        fileVals[i] = fileGrid->data[node->y][node->x];
+      } else {
+        fileVals[i] = 0;
+      }
+    }
+
+  } else {
+    // The grids are different, we must do some resampling fun.
+    #pragma omp parallel for
+    for (size_t i = 0; i < nodes.size(); i++) {
+      GridLoc pt;
+                        GridNode *node = &(nodes[i]);
+      if (fileGrid->GetGridLoc(node->refLoc.x, node->refLoc.y, &pt) && fileGrid->data[pt.y][pt.x] != fileGrid->noData) {
+        fileVals[i] = fileGrid->data[pt.y][pt.x];
+      } else {
+        fileVals[i] = 0;
+      }
+                }
+
+  }
+
+  delete fileGrid;
+	
+	
+		for (long i = numNodes - 1; i >= 0; i--) {
+      GridNode *node = &(nodes[i]);
+      float addVal = avgVals[i] + (fileVals[i] * nodes[i].area);
+      float areaAdd = areaVals[i] + nodes[i].area;
+      avgVals[i] = addVal;
+			areaVals[i] = areaAdd;
+      if (node->downStreamNode != INVALID_DOWNSTREAM_NODE) {
+        avgVals[node->downStreamNode] += addVal;
+				areaVals[node->downStreamNode] += areaAdd;
+      }
+    }	
+	
+		for (long i = numNodes - 1; i >= 0; i--) {
+      avgVals[i] = avgVals[i] / areaVals[i];
+			areaVals[i] = 0.0;
+    }
+	
+		sprintf(buffer, "%s/%s.avg.tif", outputPath, ent->d_name);
+    gridWriter.WriteGrid(&nodes, &avgVals, buffer, false);	
+		for (long i = numNodes - 1; i >= 0; i--) {
+			avgVals[i] = 0.0;
+		}
+	}
+
+	closedir(dir);
+
+	for (long i = numNodes - 1; i >= 0; i--) {
+      avgVals[i] = nodes[i].contribArea;
+	}
+  sprintf(buffer, "%s/basin.area.tif", outputPath);
+  gridWriter.WriteGrid(&nodes, &avgVals, buffer, false);
+
+	for (long i = numNodes - 1; i >= 0; i--) {
+      avgVals[i] = nodes[i].riverLen;
+  }
+  sprintf(buffer, "%s/river.length.tif", outputPath);
+  gridWriter.WriteGrid(&nodes, &avgVals, buffer, false);
+
+}
+
+void Simulator::BasinAvgPrecip() {
   PrecipReader precipReader;
   char buffer[CONFIG_MAX_LEN*2];
 #if _OPENMP
@@ -1708,7 +1824,7 @@ void Simulator::PreloadForcings(char *file, bool cali) {
       // therefore only care about averages!
       //
       readVec.resize(nodes.size());
-      std::vector<float> *vec, *vecPrev;
+      std::vector<float> *vec, *vecPrev = NULL;
       
       sprintf(buffer, "%s/%s", precipSec->GetLoc(), precipFile->GetName());
       vec = &(currentPrecipCali[tsIndex]);
diff --git a/src/Simulator.h b/src/Simulator.h
index 7951e08..d5e343a 100755
--- a/src/Simulator.h
+++ b/src/Simulator.h
@@ -24,6 +24,7 @@ class Simulator {
   
   void CleanUp();
 	void BasinAvg();
+  void BasinAvgPrecip();
   void Simulate(bool trackPeaks = false);
   float SimulateForCali(float *testParams);
   float *SimulateForCaliTS(float *testParams);
@@ -73,7 +74,7 @@ class Simulator {
   TimeUnit *timeStep, *timeStepSR, *timeStepLR, *timeStepPrecip, *timeStepQPF, *timeStepPET, *timeStepTemp, *timeStepTempF;
   float precipConvert, qpfConvert, petConvert, timeStepHours, timeStepHoursLR;
   TimeVar currentTime, currentTimePrecip, currentTimeQPF, currentTimePET, currentTimeTemp, currentTimeTempF, beginTime, endTime, warmEndTime, beginLRTime;
-  DatedName *precipFile, *qpfFile, *petFile, *tempFile, *tempFFile, currentTimeText, currentTimeTextOutput, currentDayTextOutput;
+  DatedName *precipFile, *qpfFile, *petFile, *tempFile, *tempFFile, currentTimeText, currentTimeTextOutput;
   std::vector<float> currentFF, currentSF, currentQ, avgPrecip, avgPET, avgSWE, currentSWE, avgT, avgSM, avgFF, avgSF, currentDepth;
   std::vector<FloatGrid *> paramGrids, paramGridsRoute, paramGridsSnow, paramGridsInundation;
   bool hasQPF, hasTempF, wantsDA;
diff --git a/src/TimeSeries.cpp b/src/TimeSeries.cpp
index 9ebb129..690ba41 100755
--- a/src/TimeSeries.cpp
+++ b/src/TimeSeries.cpp
@@ -26,7 +26,7 @@ void TimeSeries::LoadTimeSeries(char *file) {
 	while (!feof(tsFile)) {
 		char buffer[CONFIG_MAX_LEN];
 		float dataValue;
-		if (fscanf(tsFile, "%[^,],%f%*c", &(buffer[0]), &dataValue) == 2) {
+		if (fscanf(tsFile, "%[^,],%f ", &(buffer[0]), &dataValue) == 2) {
 			TSDataPoint *pt = new TSDataPoint;
 			pt->time.LoadTimeExcel(buffer);
 			pt->value = dataValue;

From fa17ad5450463a4ecf6cc9676fb8d9e8d1135814 Mon Sep 17 00:00:00 2001
From: Zac Flamig <zac.flamig@noaa.gov>
Date: Mon, 13 Nov 2017 06:02:43 +0000
Subject: [PATCH 2/3] Max inundation

---
 src/BasicGrids.cpp     |  3 +++
 src/GriddedOutput.cpp  |  2 ++
 src/GriddedOutput.h    |  3 ++-
 src/KinematicRoute.cpp | 28 --------------------
 src/LAEAProjection.cpp |  6 ++---
 src/Simulator.cpp      | 60 ++++++++++++++++++++++++++----------------
 src/Simulator.h        |  2 +-
 7 files changed, 49 insertions(+), 55 deletions(-)

diff --git a/src/BasicGrids.cpp b/src/BasicGrids.cpp
index a4e4695..1bf3224 100755
--- a/src/BasicGrids.cpp
+++ b/src/BasicGrids.cpp
@@ -1267,6 +1267,9 @@ void ReclassifyDDM() {
 				case 32:
           g_DDM->data[row][col] = FLOW_NORTHWEST;
           break;
+				default:
+					g_DDM->data[row][col] = g_DDM->noData;
+					break;
 			}
     }
   }
diff --git a/src/GriddedOutput.cpp b/src/GriddedOutput.cpp
index ee4aaeb..e1ebecf 100755
--- a/src/GriddedOutput.cpp
+++ b/src/GriddedOutput.cpp
@@ -23,6 +23,7 @@ const char *GriddedOutputText[] = {
 	"maxthresholdexceedance",
 	"maxthresholdexceedancep",
 	"precipaccum",
+  "maxinundation",
 };
 
 const int GriddedOutputFlags[] = {
@@ -45,4 +46,5 @@ const int GriddedOutputFlags[] = {
 	OG_MAXTHRES,
 	OG_MAXTHRESP,
 	OG_PRECIPACCUM,
+  OG_MAXDEPTH,
 };
diff --git a/src/GriddedOutput.h b/src/GriddedOutput.h
index a3cb780..fc23a77 100755
--- a/src/GriddedOutput.h
+++ b/src/GriddedOutput.h
@@ -21,9 +21,10 @@ enum SUPPORTED_OUTPUT_GRIDS {
 	OG_MAXTHRES = 32768,
 	OG_MAXTHRESP = 65536,
 	OG_PRECIPACCUM = 131072,
+  OG_MAXDEPTH = 262144,
 };
 
-#define OG_QTY 19
+#define OG_QTY 20
 
 extern const char *GriddedOutputText[];
 extern const int GriddedOutputFlags[];
diff --git a/src/KinematicRoute.cpp b/src/KinematicRoute.cpp
index 3d2eca9..d0a85b3 100755
--- a/src/KinematicRoute.cpp
+++ b/src/KinematicRoute.cpp
@@ -29,34 +29,6 @@ float KWRoute::SetObsInflow(long index, float inflow) {
 		cNode->incomingWaterOverland = inflow / node->horLen;
 	} else {
 		prev = cNode->states[STATE_KW_PQ];
-		float diff = 0.0;
-		if (inflow > 1.0 && prev > 1.0) {
-			diff = inflow/prev - 1.0;
-		}
-		GaugeConfigSection *thisGauge = node->gauge;
-		float multiplier = 1.0;
-		if (diff > 0.2) {
-			multiplier = 0.5;
-		} else if (diff < -0.2) {
-			multiplier = 2.0;
-		}
-		printf(" Multiplier %f,%f,%f,%f ", multiplier, diff, inflow, prev);
-		if (multiplier != 1.0) {
-			size_t numNodes = nodes->size();
-			for (size_t i = 0; i < numNodes; i++) {
-				node = &nodes->at(i);
-				if (node->gauge == thisGauge) {
-    			KWGridNode *cNode = &(kwNodes[i]);
-					cNode->params[PARAM_KINEMATIC_ALPHA] *= multiplier;
-					if (cNode->params[PARAM_KINEMATIC_ALPHA] < 0.01) {
-						cNode->params[PARAM_KINEMATIC_ALPHA] = 0.01;
-					} else if (cNode->params[PARAM_KINEMATIC_ALPHA] > 200.0) {
-						cNode->params[PARAM_KINEMATIC_ALPHA] = 200.0;
-					}
-				}
-			}
-
-		}
 		cNode->states[STATE_KW_PQ] = inflow;
 		cNode->incomingWaterChannel = inflow;
 	}
diff --git a/src/LAEAProjection.cpp b/src/LAEAProjection.cpp
index d05e1db..cf4a3d9 100755
--- a/src/LAEAProjection.cpp
+++ b/src/LAEAProjection.cpp
@@ -5,9 +5,9 @@
 #include "LAEAProjection.h"
 
 LAEAProjection::LAEAProjection() {
-	stan_par = TORADIANS(45.0);
-	cent_lon = TORADIANS(-100.0);
-	a = 6370997.0; 
+	stan_par = TORADIANS(48.0);  // TORADIANS(45.0);
+	cent_lon = TORADIANS(9.0); //TORADIANS(-100.0);
+	a = 6378388.0; //6370997.0; 
 }
 
 LAEAProjection::~LAEAProjection() {
diff --git a/src/Simulator.cpp b/src/Simulator.cpp
index e968ce7..b05bf51 100755
--- a/src/Simulator.cpp
+++ b/src/Simulator.cpp
@@ -330,6 +330,7 @@ bool Simulator::InitializeSimu(TaskConfigSection *task) {
     currentQ.resize(nodes.size());
     currentSWE.resize(nodes.size());
     currentDepth.resize(nodes.size());
+    computeVec.resize(nodes.size());
   } else {
     outputRP = false;
     currentFF.resize(lumpedNodes.size());
@@ -1318,13 +1319,14 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
   }
   
   // Hard coded RP counting
-  std::vector<float> count2, rpGrid, rpMaxGrid, maxGrid;
+  std::vector<float> count2, rpGrid, rpMaxGrid, maxGrid, maxDepthGrid;
   std::vector<float> SM;
   std::vector<float> dailyMaxQ, dailyMinSM, dailyMaxQHour;
   count2.resize(currentFF.size());
   rpGrid.resize(currentFF.size());
   rpMaxGrid.resize(currentFF.size());
   maxGrid.resize(currentFF.size());
+  maxDepthGrid.resize(currentFF.size());
   SM.resize(currentFF.size());
   
   for (size_t i = 0; i < currentFF.size(); i++) {
@@ -1332,6 +1334,7 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
     rpGrid[i] = 0.0;
     rpMaxGrid[i] = 0.0;
     maxGrid[i] = 0.0;
+		maxDepthGrid[i] = 0.0;
     currentFF[i] = 0.0;
     currentSF[i] = 0.0;
     currentQ[i] = 0.0;
@@ -1462,6 +1465,10 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
         SaveTSOutput();
       }
       
+			if (iModel) {
+        iModel->Inundation(&currentQ, &currentDepth);
+			}
+
       if (trackPeaks) {
         for (size_t i = 0; i < currentFF.size(); i++) {
           if (currentQ[i] > peakVals[indexYear][i]) {
@@ -1482,6 +1489,11 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
             rpMaxGrid[i] = rpGrid[i];
           }
         }
+				if (iModel) {
+					if (currentDepth[i] > maxDepthGrid[i]) {
+						maxDepthGrid[i] = currentDepth[i];
+					}
+				}
       }
       
       /*int month = currentTime.GetTM()->tm_mon;
@@ -1499,9 +1511,9 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
         sprintf(buffer, "%s/q.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), wbModel->GetName());
 	for (size_t i = 0; i < currentQ.size(); i++) {
         	float val = floorf(currentQ[i] * 10.0f + 0.5f) / 10.0f;
-        	currentDepth[i] = val;
+        	computeVec[i] = val;
     	}
-        gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+        gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
       }
       if ((griddedOutputs & OG_SM) == OG_SM) {
         sprintf(buffer, "%s/sm.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), wbModel->GetName());
@@ -1528,25 +1540,24 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
         gridWriter.WriteGrid(&nodes, &currentTempSimu, buffer, false);
       }
       if (iModel && (griddedOutputs & OG_DEPTH) == OG_DEPTH) {
-        iModel->Inundation(&currentQ, &currentDepth);
         sprintf(buffer, "%s/depth.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), iModel->GetName());
         gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
       }
       if ((griddedOutputs & OG_UNITQ) == OG_UNITQ) {
         for (size_t i = 0; i < currentQ.size(); i++) {
-          currentDepth[i] = currentQ[i] / nodes[i].contribArea;
-	  float val = floorf(currentDepth[i] * 10.0f + 0.5f) / 10.0f;
-      	  currentDepth[i] = val;
+          computeVec[i] = currentQ[i] / nodes[i].contribArea;
+	  float val = floorf(computeVec[i] * 10.0f + 0.5f) / 10.0f;
+      	  computeVec[i] = val;
         }
         sprintf(buffer, "%s/unitq.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), wbModel->GetName());
-        gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+        gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
       }
       if (outputThres && (griddedOutputs & OG_THRES) == OG_THRES) {
         for (size_t i = 0; i < currentQ.size(); i++) {
-          currentDepth[i] = ComputeThresValue(currentQ[i], actionVals[i], minorVals[i], moderateVals[i], majorVals[i]);
+          computeVec[i] = ComputeThresValue(currentQ[i], actionVals[i], minorVals[i], moderateVals[i], majorVals[i]);
         }
         sprintf(buffer, "%s/thres.%s.%s.tif", outputPath, currentTimeTextOutput.GetName(), wbModel->GetName());
-        gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+        gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
         
       }
       
@@ -1614,9 +1625,9 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
     sprintf(buffer, "%s/maxq.%04i%02i%02i.%02i%02i%02i.tif", outputPath, ctWE->tm_year + 1900, ctWE->tm_mon + 1, ctWE->tm_mday, ctWE->tm_hour, ctWE->tm_min, ctWE->tm_sec);
     for (size_t i = 0; i < currentQ.size(); i++) {
     	float val = floorf(maxGrid[i] * 10.0f + 0.5f) / 10.0f;
-	currentDepth[i] = val;
+	computeVec[i] = val;
     }
-    gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+    gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
   }
   
   if (sModel && (griddedOutputs & OG_MAXSWE) == OG_MAXSWE) {
@@ -1626,32 +1637,37 @@ void Simulator::SimulateDistributed(bool trackPeaks) {
   
   if ((griddedOutputs & OG_MAXUNITQ) == OG_MAXUNITQ) {
     for (size_t i = 0; i < currentQ.size(); i++) {
-      currentDepth[i] = maxGrid[i] / nodes[i].contribArea;
-      float val = floorf(currentDepth[i] * 10.0f + 0.5f) / 10.0f;
-      currentDepth[i] = val;
+      computeVec[i] = maxGrid[i] / nodes[i].contribArea;
+      float val = floorf(computeVec[i] * 10.0f + 0.5f) / 10.0f;
+      computeVec[i] = val;
     }
     sprintf(buffer, "%s/maxunitq.%04i%02i%02i.%02i%02i%02i.tif", outputPath, ctWE->tm_year + 1900, ctWE->tm_mon + 1, ctWE->tm_mday, ctWE->tm_hour, ctWE->tm_min, ctWE->tm_sec);
-    gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+    gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
   }
+
+  if ((griddedOutputs & OG_MAXDEPTH) == OG_MAXDEPTH) {
+    sprintf(buffer, "%s/maxdepth.%04i%02i%02i.%02i%02i%02i.tif", outputPath, ctWE->tm_year + 1900, ctWE->tm_mon + 1, ctWE->tm_mday, ctWE->tm_hour, ctWE->tm_min, ctWE->tm_sec);
+    gridWriter.WriteGrid(&nodes, &maxDepthGrid, buffer, false);
+  }	
   
   if (outputThres && (griddedOutputs & OG_MAXTHRES) == OG_MAXTHRES) {
     for (size_t i = 0; i < currentQ.size(); i++) {
-      currentDepth[i] = floorf(ComputeThresValue(maxGrid[i], actionVals[i], minorVals[i], moderateVals[i], majorVals[i]) * 10.0f + 0.5f) / 10.0f;
-      if (currentDepth[i] < 1.0) {
-	currentDepth[i] = 0.0;
+      computeVec[i] = floorf(ComputeThresValue(maxGrid[i], actionVals[i], minorVals[i], moderateVals[i], majorVals[i]) * 10.0f + 0.5f) / 10.0f;
+      if (computeVec[i] < 1.0) {
+	computeVec[i] = 0.0;
 	}
     }
     sprintf(buffer, "%s/maxthres.%04i%02i%02i.%02i%02i%02i.tif", outputPath, ctWE->tm_year + 1900, ctWE->tm_mon + 1, ctWE->tm_mday, ctWE->tm_hour, ctWE->tm_min, ctWE->tm_sec);
-    gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+    gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
     
   }
   
   if (outputThresP && (griddedOutputs & OG_MAXTHRESP) == OG_MAXTHRESP) {
     for (size_t i = 0; i < currentQ.size(); i++) {
-      currentDepth[i] = floorf(ComputeThresValueP(maxGrid[i], actionVals[i], actionSDVals[i], minorVals[i], minorSDVals[i], moderateVals[i], moderateSDVals[i], majorVals[i], majorSDVals[i]) * 10.0f + 0.5f) / 10.0f;
+      computeVec[i] = floorf(ComputeThresValueP(maxGrid[i], actionVals[i], actionSDVals[i], minorVals[i], minorSDVals[i], moderateVals[i], moderateSDVals[i], majorVals[i], majorSDVals[i]) * 10.0f + 0.5f) / 10.0f;
     }
     sprintf(buffer, "%s/maxthresp.%04i%02i%02i.%02i%02i%02i.tif", outputPath, ctWE->tm_year + 1900, ctWE->tm_mon + 1, ctWE->tm_mday, ctWE->tm_hour, ctWE->tm_min, ctWE->tm_sec);
-    gridWriter.WriteGrid(&nodes, &currentDepth, buffer, false);
+    gridWriter.WriteGrid(&nodes, &computeVec, buffer, false);
     
   }
   
diff --git a/src/Simulator.h b/src/Simulator.h
index d5e343a..598eb45 100755
--- a/src/Simulator.h
+++ b/src/Simulator.h
@@ -75,7 +75,7 @@ class Simulator {
   float precipConvert, qpfConvert, petConvert, timeStepHours, timeStepHoursLR;
   TimeVar currentTime, currentTimePrecip, currentTimeQPF, currentTimePET, currentTimeTemp, currentTimeTempF, beginTime, endTime, warmEndTime, beginLRTime;
   DatedName *precipFile, *qpfFile, *petFile, *tempFile, *tempFFile, currentTimeText, currentTimeTextOutput;
-  std::vector<float> currentFF, currentSF, currentQ, avgPrecip, avgPET, avgSWE, currentSWE, avgT, avgSM, avgFF, avgSF, currentDepth;
+  std::vector<float> currentFF, currentSF, currentQ, avgPrecip, avgPET, avgSWE, currentSWE, avgT, avgSM, avgFF, avgSF, currentDepth, computeVec;
   std::vector<FloatGrid *> paramGrids, paramGridsRoute, paramGridsSnow, paramGridsInundation;
   bool hasQPF, hasTempF, wantsDA;
 	bool inLR;

From 73711cd695531f8c1c219cb91cee9b13a39400ee Mon Sep 17 00:00:00 2001
From: Zac Flamig <zac.flamig@noaa.gov>
Date: Tue, 10 Apr 2018 18:47:41 +0000
Subject: [PATCH 3/3] average grids

---
 src/BasicGrids.cpp | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/BasicGrids.cpp b/src/BasicGrids.cpp
index 1bf3224..8a52baa 100755
--- a/src/BasicGrids.cpp
+++ b/src/BasicGrids.cpp
@@ -317,7 +317,7 @@ void ClipBasicGrids(long x, long y, long search, const char *output) {
 	maxY = loc.y + search;
 	minY = loc.y - search;
   
-	FindIndBasins(-124.73, -66.95, 50.00, 24.5);
+	FindIndBasins(-180.0, 180.0, 90.0, -90.0);
 	//FindIndBasins(-84.84, -67.26, 44.17, 35.76);
   
 	printf("Search box is %ld, %ld, %ld, %ld, %ld, %ld\n", minX, maxX, minY, maxY, x, y);
@@ -726,6 +726,10 @@ void CarveBasin(BasinConfigSection *basin, std::vector<GridNode> *nodes, std::ma
 			gauge->SetLat(ref.y);
 			gauge->SetLon(ref.x);
 		}
+    if (g_DEM->data[loc->y][loc->x] == g_DEM->noData || g_FAM->data[loc->y][loc->x] == g_FAM->noData ||  g_FAM->data[loc->y][loc->x] < 0) {
+      ERROR_LOGF("Gauge \"%s\" is located in a no data grid cell!", gauge->GetName());
+      return;
+    }
 		gauge->SetFlowAccum(g_FAM->data[loc->y][loc->x]);
 		gauge->SetUsed(false);
 		INFO_LOGF("Gauge %s (%f, %f; %ld, %ld): FAM %ld", gauge->GetName(), gauge->GetLat(), gauge->GetLon(), loc->y, loc->x, gauge->GetFlowAccum());
@@ -918,7 +922,8 @@ void CarveBasin(BasinConfigSection *basin, std::vector<GridNode> *nodes, std::ma
 			//We compute slope here too!
 			if (keepGoing) {
 			for (int i = 1; i < FLOW_QTY; i++) {
-				if (TestUpstream(currentN, (FLOW_DIR)i, &nextNode)) {
+				if (TestUpstream(currentN, (FLOW_DIR)i, &nextNode) && currentNode < totalAccum) {
+          //printf("currentNode is %li\n", currentNode);
 					GridNode *nextN = &(*nodes)[currentNode];
 					nextN->index = currentNode;
 					currentNode++;
@@ -1071,7 +1076,7 @@ bool TestUpstream(long nextX, long nextY, FLOW_DIR dir, GridLoc *loc) {
 			return false;
 	}
   
-	if (nextX >= 0 && nextY >= 0 && nextX < g_DDM->numCols && nextY < g_DDM->numRows && g_DDM->data[nextY][nextX] == wantDir) { // && g_FAM->data[nextY][nextX] <= currentFAC) {
+	if (nextX >= 0 && nextY >= 0 && nextX < g_DDM->numCols && nextY < g_DDM->numRows && g_DDM->data[nextY][nextX] == wantDir && g_DEM->data[nextY][nextX] != g_DEM->noData && g_FAM->data[nextY][nextX] != g_FAM->noData) { // && g_FAM->data[nextY][nextX] <= currentFAC) {
 		loc->x = nextX;
 		loc->y = nextY;
 		return true;