// *SC*, Type, class, and routine to support lakes for branching flows
	// *SC*, 2/29/08
	
	// flooded local minima --> non-flooded outlets
	typedef struct
	{
	tLNode *outlet;
	tLNode **drains;
	int nDrains;
	} lakeRouteType;

	// sort outlets by elevation
	int compareLakeRouteTypes(const void *p1, const void *p2)
	{
	lakeRouteType *lp1 = (lakeRouteType *) p1;
	lakeRouteType *lp2 = (lakeRouteType *) p2;
	double z1 = lp1->outlet->getZ();
	double z2 = lp2->outlet->getZ();

	if (z1 > z2) return -1;
	else if (z1 < z2) return 1;
	else return 0;
	}

	// sort flooded minima by elevation
	int compareDrains(const void *p1, const void *p2)
	{
	tLNode *np1 = *((tLNode **) p1);
	tLNode *np2 = *((tLNode **) p2);
	double z1 = np1->getZ();
	double z2 = np2->getZ();

	if (z1 > z2) return -1;
	else if (z1 < z2) return 1;
	else return 0;
	}

	// list of lake bottoms and corresponding outlets
	class LakeRoutes
	{
	public:
		tStreamNet *net;
		lakeRouteType *routes;
		int nRoutes;

		LakeRoutes(tStreamNet *net) { this->net = net; routes = 0; nRoutes = 0; };
		
		~LakeRoutes()
		{
		for (int i = 0; i < nRoutes; i ++) free(routes[i].drains);
		free(routes);
		}
		
		// add a (drain, outlet) pair
		void add(tLNode *drain, tLNode *outlet) 
		{
		if (nRoutes == 0) 
			{
			routes = (lakeRouteType *) malloc(sizeof(lakeRouteType));
			routes[nRoutes].outlet = outlet;
			routes[nRoutes].drains = (tLNode **) malloc(sizeof(tLNode *));
			routes[nRoutes].drains[0] = drain;
			routes[nRoutes].nDrains = 1;
			nRoutes ++;
			}
		else 
			{
			int i;
			for (i = 0; i < nRoutes; i ++) if (routes[i].outlet == outlet) break;
			if (i < nRoutes)
				{
				routes[i].drains = (tLNode **) realloc(routes[i].drains, (routes[i].nDrains + 1) * sizeof(tLNode *));
				routes[i].drains[routes[i].nDrains] = drain;
				routes[i].nDrains ++;
				}
			else
				{
				routes = (lakeRouteType *) realloc(routes, (nRoutes + 1) * sizeof(lakeRouteType));
				routes[nRoutes].outlet = outlet;
				routes[nRoutes].drains = (tLNode **) malloc(sizeof(tLNode *));
				routes[nRoutes].drains[0] = drain;
				routes[nRoutes].nDrains = 1;
				nRoutes ++;
				}
			}
		};
		
		// sort outlets and drains by elevation
		void sort(void) 
		{ 
		qsort(routes, nRoutes, sizeof(lakeRouteType), &compareLakeRouteTypes); 
		for (int i = 0; i < nRoutes; i ++)
			qsort(routes[i].drains, routes[i].nDrains, sizeof(tLNode *), &compareDrains);
		};
		
		// make sure outlets above drains
		void check(void) 
		{
		bool isOK = true;
		for (int i = 0; i < nRoutes; i ++)
			{
			for (int j = 0; j < routes[i].nDrains; j ++)
				{
				if (routes[i].drains[j]->getZ() >= routes[i].outlet->getZ()) 
					{
					printf("LakeRoutes:check(): !!! drain %12.4f %12.4f %12.4f >= outlet %12.4f %12.4f %12.4f !!!\n", 
						routes[i].drains[j]->getX(), routes[i].drains[j]->getY(), routes[i].drains[j]->getZ(),
						routes[i].outlet->getX(), routes[i].outlet->getY(), routes[i].outlet->getZ());	
					isOK = false;
					}
				}
			}
		if (isOK) printf("   check() OK\n");
		};

		// print outlets and drains
		void print(void) 
		{
		printf("   %d lake routes:\n", nRoutes);
		for (int i = 0; i < nRoutes; i ++)
			{
			printf("   %12.4f %12.4f %12.4f   -- %d drains:\n", 
				routes[i].outlet->getX(), routes[i].outlet->getY(), routes[i].outlet->getZ(), routes[i].nDrains);		
			for (int j = 0; j < routes[i].nDrains; j ++)
				{
				printf("      %12.4f %12.4f %12.4f\n", 
					routes[i].drains[j]->getX(), routes[i].drains[j]->getY(), routes[i].drains[j]->getZ());		
				}
			}
		};
		
		// from top to bottom (in elevation):
		//   * calculate flow at drains (into bottoms of lakes)
		//   * add to corresponding outlets
		//   * remove lake from further calculations (to avoid loops)
		void eval(void)
		{
		for (int i = 0; i < nRoutes; i ++)
			{
			for (int j = 0; j < routes[i].nDrains; j ++)
				net->RouteFlowArea(routes[i].drains[j], 0.0);
			for (int j = 0; j < routes[i].nDrains; j ++)
				routes[i].outlet->AddDrArea(routes[i].drains[j]->getDrArea());
			routes[i].outlet->setIgnoreFloodedEdgesFlag(true);
			}
		};
	};

	// *SC*, Drain lake bottoms to branching outlet
	// *SC*, 2/28/08
	// * build list of lake drains and corresponding outlets
	// * from highest to lowest, evaluate flow into lakes and add to outlets
	// * remove lakes from further calculations
	void tStreamNet::DrainLakes(void)
	{
	tLNode * curnode;
	tMesh< tLNode >::nodeListIter_t nodIter( meshPtr->getNodeList() );

	int nodes = 0;
	int notFlooded = 0;
	int flooded = 0;
	int minimum = 0;
	LakeRoutes lakeRoutes(this);

	printf("DrainLakes()\n");
	for (curnode = nodIter.FirstP(); nodIter.IsActive(); curnode = nodIter.NextP() )
		{
		nodes ++;
		
		if (curnode->getFloodStatus() == tLNode::kNotFlooded)
			{
			notFlooded ++;
			if (0) printf("DrainLakes(): %12.4f %12.4f %12.4f   -- not flooded --\n", 
				curnode->getX(), curnode->getY(), curnode->getZ());
			}

		// look for flooded nodes
		if (curnode->getFloodStatus() == tLNode::kFlooded)
			{
			flooded ++;
			bool wasMinimum = false;
			tLNode *startNode = curnode;
			
			// look for bottoms
			if (curnode->isLocalMinimum()) 
				{
				minimum ++;
				wasMinimum = true;
				printf("DrainLakes(): %12.4f %12.4f %12.4f   -- flooded and minimum --\n", 
					curnode->getX(), curnode->getY(), curnode->getZ());
				}
			else if (0) printf("DrainLakes(): %12.4f %12.4f %12.4f   -- flooded --\n", 
				curnode->getX(), curnode->getY(), curnode->getZ());
				
			// find outlet (flow path was constructed in FillLakes())
			while (curnode && curnode->getFloodStatus() != tLNode::kNotFlooded) 
				curnode = curnode->getDownstrmNbr();
				
			// error: no outlet
			if (!curnode) printf("DrainLakes(): !!! flooded node does not lead to outlet !!!\n");
			
			// find outlet of bottoms
			if (curnode && wasMinimum) 
				{
				if (curnode->getBoundaryFlag() != kNonBoundary) printf("   outlet at: %12.4f %12.4f %12.4f   -- boundary --\n", 
					curnode->getX(), curnode->getY(), curnode->getZ());
				else 
					{
					printf("   outlet at: %12.4f %12.4f %12.4f\n", 
						curnode->getX(), curnode->getY(), curnode->getZ());
					// add non-boundary (drain, outlet) pair
					lakeRoutes.add(startNode, curnode);
					}
				}
			}
		}
		
	printf("DrainLakes(): nodes = %d, notFlooded = %d, flooded = %d, minimum = %d\n", nodes, notFlooded, flooded, minimum);
	// sort top to bottom
	lakeRoutes.sort();
	lakeRoutes.check();
	lakeRoutes.print();
	// evaluate flows, leaving at outlets for further access in RouteFlowArea() below
	lakeRoutes.eval();
	}

// *SC*, Added define to allow branching flows
// *SC*, 2/23/08
// *SC*, Now superseded by global var 'gBranchingFlows' (2/28/08)
#define MULTIPLE_DOWNSTREAM_FLOWS 2
// *SC*, Added define to inhibit routing flows other than from inlet
// *SC*, 2/27/08
// *SC*, Should only be set to 0 if MULTIPLE_DOWNSTREAM_FLOWS == 1
// *SC*, Already obsolete (2/28/08)
#define ROUTE_ALL_FLOWS 1

void tStreamNet::DrainAreaVoronoi()
{
   if (0) //DEBUG
     std::cout << "DrainAreaVoronoi()..." << std::endl;

   tLNode * curnode;
   tMesh< tLNode >::nodeListIter_t nodIter( meshPtr->getNodeList() );

   // Reset drainage areas to zero
   for( curnode = nodIter.FirstP(); nodIter.IsActive();
        curnode = nodIter.NextP() )
        {
       	curnode->setDrArea( 0. );
			// *SC*, Clear calculated flow flag for doing upstream branching calcs
			// *SC*, 2/27/08
       	//if (MULTIPLE_DOWNSTREAM_FLOWS == 2) curnode->setFlowCalculated(false);
			// *SC*, Added global var
			// *SC*, 2/28/08
       	if (gBranchingFlows) curnode->setFlowCalculatedFlag(false);
       	}

	// process lakes for branching flows
	if (gBranchingFlows) DrainLakes();
	
   // send voronoi area for each node to the node at the other end of the
   // flowedge and downstream
	// *SC*, Added define to inhibit routing flows other than from inlet
	// *SC*, 2/27/08
   if (ROUTE_ALL_FLOWS) for( curnode = nodIter.FirstP(); nodIter.IsActive();
        curnode = nodIter.NextP() )
   {

      // Debug Quintijn
      if(curnode->getVArea() < 0.0){
        std::cout<< "Voronoi area <  0.0 at \n";
        std::cout<< curnode->getX() <<' '<<curnode->getY()<<std::endl;
        std::cout<< "Area= "<<curnode->getVArea()<<std::endl;
      	exit(1);
      }

      RouteFlowArea( curnode, curnode->getVArea() );
   }
   if( inlet.innode != 0 )
   {
#if 1
      //inlet.FindNewInlet();            QC commented this out -> avoid channel hopping around
#endif
      RouteFlowArea( inlet.innode, inlet.inDrArea );
   }
	
	// add unbranched flows (since lake bottoms 'jump' to outlets otherwise)
	if (gBranchingFlows)
		{
		gBranchingFlows = false;
		
		tMesh< tLNode >::nodeListIter_t nodIter( meshPtr->getNodeList() );
		for( curnode = nodIter.FirstP(); nodIter.IsActive(); curnode = nodIter.NextP() )
			RouteFlowArea( curnode, curnode->getVArea() );
		
		if( inlet.innode != 0 ) RouteFlowArea( inlet.innode, inlet.inDrArea );
		
		gBranchingFlows = true;
		}
		
   if (0) //DEBUG
     std::cout << "DrainAreaVoronoi() finished" << std::endl;
}


// etc..........................................................................


/*****************************************************************************\
**
**  tStreamNet::RouteFlowArea
**
**  Starting with the current node, this routine increments
**  the drainage area of the node and each node downstream by _addedArea_.
**
\*****************************************************************************/

// *SC*, No need for fn to be inline, hopefully compiler would just ignore this
// *SC*, 2/26/08
// inline
void tStreamNet::RouteFlowArea( tLNode *curnode, double addedArea )
{
// *SC*, Added experimental support for branching flows
// *SC*, 2/23/08
// *SC*, Branching flows: 'right' way to do it: downstream-->upstream
// *SC*, 2/27/08
// *SC*, Added global var
// *SC*, 2/28/08
//if (MULTIPLE_DOWNSTREAM_FLOWS == 2)
if (gBranchingFlows)
	{
	// return if already done
	if (curnode->getFlowCalculatedFlag()) 
		{
		if (0 && !gSilentMode) printf("RouteFlowArea(2): %12.4f %12.4f %12.4f   -- done --\n", 
			curnode->getX(), curnode->getY(), curnode->getZ());
		return;
		}
	
	// return if boundary or sink(?)
	if (curnode->getBoundaryFlag() != kNonBoundary 
		/*|| curnode->getFloodStatus() == tLNode::kSink*/) 
		{
		if (0 && !gSilentMode) printf("RouteFlowArea(2): %12.4f %12.4f %12.4f   -- boundary/sink --\n", 
			curnode->getX(), curnode->getY(), curnode->getZ());
		return;
		}
	
	if (!gSilentMode) printf("RouteFlowArea(2): %12.4f %12.4f %12.4f\n", 
		curnode->getX(), curnode->getY(), curnode->getZ());
	
	// initialize flow to Voronoi area or inlet area
	double sumFlow;
	if (curnode == inlet.innode) sumFlow = inlet.inDrArea;
	else sumFlow = curnode->getVArea();
	tEdge *e0 = curnode->getEdg();
	tEdge *e = e0;
	
	// loop through upstream edges
	do
		{
		// note that slopes are first calculated with CalcSlope() and 
		// then can be accessed with getSlope() without recalculating,
		// not that there's very much math involved.  slopes are also
		// signed the wrong way: up is negative, down is positive
		
		// upstream edge has negative slope(!)
		if (e->CalcSlope() >= 0) continue;
		
		// need support for lakes: upstream node is lower in elevation, 
		// but has flow edge to here (and from every other flooded node).
		// issue is avoiding infinite loop when evaluating lower node flow
		// by calling back to this node, etc...
		// 2/29/08: actually more complicated than this; see DrainLakes() above
		
		// disregard upstream boundaries and sinks(?)
		tLNode *upstreamNode = (tLNode *) e->getDestinationPtrNC();
		if (upstreamNode->getBoundaryFlag() != kNonBoundary 
			/*|| upstreamNode->getFloodStatus() == tLNode::kSink*/) continue;
			
		// make sure flow OK from upstream node
		tEdge *ec = e->getComplementEdge();
		if (!ec->FlowAllowed() || ec->CalcSlope() <= 0) continue;
		
		// recursively calculate upstream flow (recursion depth 
		// is only max flow length, e.g. length of mesh)
		RouteFlowArea(upstreamNode, 0.0);  // second argument is ignored
		
		// sum slopes of all downstream edges from upstream node
		double sumSlope = 0.0;
		tEdge *e2 = ec;
		do
			{
			if (e2->FlowAllowed() && e2->CalcSlope() > 0)
				{
				if (upstreamNode->getIgnoreFloodedEdgesFlag() 
					&& ((tLNode *) e2->getDestinationPtrNC())->getFloodStatus() == tLNode::kFlooded) continue;
				sumSlope += e2->getSlope();
				}
			}
		while (e2 = e2->getCCWEdg(), e2 != ec);
			
		// added flow is weighted by slope to here
		if (sumSlope > 0) sumFlow += upstreamNode->getDrArea() * (ec->getSlope() / sumSlope);
		}
	while (e = e->getCCWEdg(), e != e0);
	
	// set flow area and flag
	curnode->setDrArea(sumFlow);
	curnode->setFlowCalculatedFlag(true);
	}  // More below ('else' clause)...

// *SC*, Branching flows: wrong way to do it: upstream-->downstream
// *SC*, 2/23/08
// *SC*, Already obsolete (2/28/08)
/*else if (MULTIPLE_DOWNSTREAM_FLOWS == 1)
	{
	if (!gSilentMode) std::cout << "RouteFlowArea(1): " << curnode->getX() << " " << curnode->getY() 
		<< " " << curnode->getZ() << " " << addedArea << std::endl;
	
	if (curnode->getBoundaryFlag() != kNonBoundary 
		|| curnode->getFloodStatus() == tLNode::kSink) return;
		
	curnode->AddDrArea( addedArea );
		
	// loop over all edges
	double sumSlope = 0.0;
	tEdge *e0 = curnode->getEdg();
	tEdge *e = e0;
	do
		{
		// look for allowed downstream (positive slope) flows, add slope to sum
		if (e->FlowAllowed() && e->CalcSlope() > 0)
			sumSlope += e->getSlope();
		e = e->getCCWEdg();
		}
	while (e != e0);

	// loop over edges again
	e = e0;
	do
		{
		// look for allowed downstream (positive slope) flows
		if (e->FlowAllowed() && e->getSlope() > 0)
			{
			// weight added area and send downstream
			double flowArea = addedArea * (e->getSlope() / sumSlope);
			RouteFlowArea((tLNode *) e->getDestinationPtrNC(), flowArea);
			}
		e = e->getCCWEdg();
		}
	while (e != e0);
	}*/
	
else
	{
   if (0) //DEBUG
     std::cout << "RouteFlowArea(0)..." << std::endl;
//#if DEBUG
   int niterations=0;  // Safety feature: prevents std::endless loops
//#endif

   // As long as the current node is neither a boundary nor a sink, add
   // _addedArea_ to its total drainage area and advance to the next node
   // downstream
   while( (curnode->getBoundaryFlag() == kNonBoundary) &&
	  (curnode->getFloodStatus()!=tLNode::kSink) )
   {
      curnode->AddDrArea( addedArea );
      curnode = curnode->getDownstrmNbr();
//#if DEBUG
      niterations++;
      if( unlikely(niterations>9990) )
	RouteError( curnode );
//#endif
   }
   if (0) //DEBUG
     std::cout << "RouteFlowArea() finished" << std::endl;
	}
}