/**************************************************************************/
// Hydrologic functions used to distribute surface water horizontally
//@ Alexey Voinov -- voinov@cbl.umces.edu
/**************************************************************************/

#include "SWater.h"

#define MAXCELLS 200 //number of cells that can be considered in open water

float n_out_sum, n_out_conc;  // global variables for RALF


//*************************************************************************
// Combines fluxing over steep areas and fluxing by equilibrating in open water
// Also moves 'Stuff' together with water. 

// Called from SURFACE_WATER variable. AvailWater is SURFACE_WATER
// HydMap - map of streams (usually HYDRO). It has the mouth of the river
//   marked by MOUTH (read from SWdata file) and stations where to monitor 
//   marked by STAT (read from SWdata file).
// Stuff - is the stuff to be moved with water


void SWTransport4( CVariable& AvailWater, CVariable& HydMap,  CVariable& HabMap, 
	CVariable& Elev, CVariable& Stuff, CVariable& outFlux )
	
{	
	static FILE *file, *fileN, *fileIn;
	static int date = 0;
	float TotalOut = 0.;
	float TotNout = 0.;   // the amount of stuff that is fluxed out from the last cell
	float Nconc, Ntot;
	
	int niter;
	float flow_rate = 1.;
	float df, w_sum, ftmp, move;
	Pix p, rp;
	Pix rPix[MAXCELLS];

	static float MAXDRUN, HSHEAD, WATLEV, NDIF, transport_st, flux_parm;	
	static int  Nmax, OPENW, GSTAT, MOUTH, ROWOUT, COLOUT, ONOFF;

	int il, rr, rc, NN;
	char ccc[1000];


	DistributedGrid& grid = AvailWater.Grid();
	grid.SetPointOrdering(0);	                 // Sets grid ordering to default

// Read control Data from file SWData
	if (fileIn == NULL) 
	{
		CPathString pathName(Env::DataPath()); 
    		pathName.Add("SWData");
		if ( (fileIn = fopen (pathName, "r")) == NULL)  gFatal( "Can't open SWData file! " );
		
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &MAXDRUN );	//number of iterations across the area over 1 time step
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &HSHEAD );	//half-saturation for routing
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &WATLEV );	//water level at outlet		
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &NDIF );

		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &OPENW );	//category for open water
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &GSTAT );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &MOUTH );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &Nmax );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &transport_st );

		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &flux_parm );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &ROWOUT );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &COLOUT );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &ONOFF );	//switches function on/off
	}

	if (!ONOFF ) return;

	if (fileN == NULL) 
	{ // open output file for stuff
		CPathString pathName(Env::ArchivePath()); 
		pathName.Add("NOut.txt");
		if ( (fileN = fopen (pathName, "w")) == NULL) gFatal( "Can't open NOut file! " ); 
	}

	if (file == NULL) 
	{ // open output file for water
		CPathString pathName(Env::ArchivePath()); 
		pathName.Add("WSOut.txt");
		if ( (file = fopen (pathName, "w")) == NULL) gFatal( "Can't open WSOut file! " ); 
	  
	  // prepare output file header
		fprintf ( file, "transport_st=%f flux_parm=%f MAXDRUN=%f HSHEAD=%f WATLEV=%f NDIF=%f\n", 
			      transport_st, flux_parm, MAXDRUN, HSHEAD, WATLEV, NDIF );

		for( p = grid.first(); p; grid.next(p) )
		{
			const OrderedPoint& pt = grid.GetPoint(p);
			                                      
			if( HydMap(pt) == GSTAT ) 
			{
				int ir0 = pt.row(), ic0 = pt.col();
				fprintf ( file, "\t%9dx%d       \t", ir0, ic0 );
				fprintf ( fileN, "\t%9dx%d\t            \t            ", ir0, ic0 );
			}
	        
			if ( pt.row()==ROWOUT && pt.col()==COLOUT)
               			fprintf ( fileN, "\t%9dx%d\t", ROWOUT, COLOUT );
			
		}
        
		fprintf ( file, "\tTotal from area" );
		fprintf ( fileN, "\t From area\t Conc at out\t Total in area" );
		
	}

	fprintf ( file, "\n%d", date );
	fprintf ( fileN, "\n%d", date++ );
	
	static CVariable* swFlux = NULL;	
	if(swFlux == NULL ) { swFlux = AvailWater.GetSimilarVariable("SWFlux"); }
	    // intermediate increment to stage

	static CVariable* NFlux = NULL;	
	if(NFlux == NULL ) { NFlux = AvailWater.GetSimilarVariable("NFlux"); }
	    // intermediate increment to stuff
	                      
	static CVariable* NTot = NULL;	
	if(NTot == NULL ) { NTot = AvailWater.GetSimilarVariable("NTot"); }
	    // array that is used to accumulate stuff fluxed through gaging points

	static CVariable* swTot = NULL;	
	if(swTot == NULL ) { swTot = AvailWater.GetSimilarVariable("SWTot"); }
	    // array that is used to accumulate water fluxed through gaging points
	
	swTot -> Set(0.0);
	NTot -> Set(0.0);
	                      
	AvailWater.LinkEdges();
	Stuff.LinkEdges();

	swFlux->Set(0.0);
	NFlux->Set(0.0);


// loop over the whole area, but skipping the open water cells	
	for( p = grid.first(); p; grid.next(p) )               
	{ 
		const OrderedPoint& pt = grid.GetPoint(p);     // pt same as p, but ordered
	     
		float& w0 = AvailWater(pt);
		float nf = 0.;

		if (date == 0) outFlux(pt) = 0;

		if (HabMap(pt) == OPENW && HydMap(pt) != 4) continue;
	//water flow in terrestrial areas and those that are outlets from reservoirs

		df = w0*flux_parm;

		(*swFlux)(pt) -= df;

		if (w0>transport_st)  nf = NDIF*Stuff(pt)*flux_parm;
			// it should be *df/w0, which is in fact = flux_parm
		
		(*NFlux)(pt) -= nf;
		outFlux(pt) += nf; //accumulate outflow of nutrient from cell

		Pix rp = p;
	    	
	// making the length of flow path dependent of the amount of water to move
	
		float stt = w0*w0; 
		float godo = 1 + MAXDRUN*stt/(stt+HSHEAD);
		niter = godo;
	    		    	
		for(int iter=0; iter<niter && rp; iter++) 
		{ // make niter iterations to run water further downhill
			rp = grid.LinkedPix( rp );           // gets downstream link rp

			if ( rp )
			{
				const OrderedPoint& rpt = grid.GetPoint(rp);
									
				if (HydMap(rpt) == GSTAT)
					// GSTAT mark the gaging stations on River map.
					// this is to sum up all the water and all the 
					// materialthat is moved
					// through the cells with gaging stations
				{
					(*swTot)(rpt) += df;
					(*NTot)(rpt) += nf;
				}
				if (HydMap(rpt) == MOUTH)  
					// MOUTH mark the mouth on River map.
         			{
					TotNout += nf;
					TotalOut += df;
          			}
				if (HabMap(rpt) == OPENW) break;  //stop if we hit open water cell
                        
			}
		} // end loop over cells along the path of flow
		   		
                // the following is when everything is dumped into one last recipient cell
		if ( rp ) 
		{  
		 	const OrderedPoint& rpt = grid.GetPoint(rp);
		 	(*swFlux)(rpt) += df;
		 	if (w0>transport_st) (*NFlux)(rpt) += nf;
		}
	  } //end loop over all area
			
	AvailWater.AddData(*swFlux);
	Stuff.AddData(*NFlux);

// another loop over the whole area, but only for open water cells	
	for( p = grid.first(); p; grid.next(p) )               
	{ 
		const OrderedPoint& pt = grid.GetPoint(p);     // pt same as p, but ordered
		if (HabMap(pt) != OPENW) continue;

	// water flow in Open Water regions
		w_sum = Elev(pt)+AvailWater(pt);
		int count = 0, ow = 1;  // number of open water cells, presence of open water in the vicinity

			// calculate the average water level
			// loops until either the max allowed number of cells reached, 
			// or no adjacent water cell present in the vicinity
		for (NN = 1; count < Nmax-1 && ow; NN++)  //NN is number of concentric layer
		{ 
			ow = 0; 
			for (rr = -NN; rr <= NN; rr++)
			for (rc = -NN; rc <= NN; rc += 2*NN)

			{
				rp = grid.TranslateByGridCoord( p, rr, rc );

				if( rp )  
				{
					const OrderedPoint& rpt = grid.GetPoint(rp);
					if (HabMap(rpt) == OPENW) 
					{ 
						w_sum += Elev(rpt)+AvailWater(rpt);
						ow = 1;
						rPix[count++] = rp;
					}  
				} 
			}

			for (rc= -NN+1; rc<NN; rc++)
			for (rr = -NN; rr<=NN; rr+=2*NN)
			{  
				rp = grid.TranslateByGridCoord( p, rr, rc );

				if( rp )  
				{
					const OrderedPoint& rpt = grid.GetPoint(rp);
					if (HabMap(rpt) == OPENW) 
					{ 
						w_sum += Elev(rpt)+AvailWater(rpt);
						ow = 1;
						rPix[count++] = rp;
					}  
				} 
			}
		} // end loop over concentric neighborhood 

		float h0 = w_sum/(count+1);
	        
	        // make sure that the average water level does not expose ground
		if (h0 < Elev(pt)) h0 -= (Elev(pt)-h0)/count;
	        
		for( il = 0; il < count-1; il++ ) 
		{ 
			rp = rPix[il];
			const OrderedPoint& rpt = grid.GetPoint(rp);

			if (h0 < Elev(rpt)) h0 -= (Elev(rpt)-h0)/count;
		}
		// we can still potentially generate water in this case !!!	
				  
		// equilibrate  
		for( il = 0; il < count; il++ ) 
		{ 
			rp = rPix[il];
			const OrderedPoint& rpt = grid.GetPoint(rp);
		  
			ftmp = ( h0 - (Elev(rpt)+AvailWater(rpt)) ) * flow_rate;
			if (-ftmp > AvailWater(rpt)) ftmp = -AvailWater(rpt);
		      
		  // for exchange of stuff, I assume that the transport is occuring via the center
		  // point pt. So that at each step the stuff is moved from each of the 
		  // vicinity cells to/from the center point. This can potentially drive the
		  // amount of stuff in the center point below zero, but at the end of the loop
		  // over the vicinity it should be back to normal.
		     
			if (ftmp > AvailWater(pt)) 
			{ 
				ftmp = AvailWater(pt);
				move = Stuff(pt);
			}
			else 
			{ 
				if (ftmp > 0.)
				move = (AvailWater(pt) > 0.) ? Stuff(pt)*ftmp/AvailWater(pt) : 0.;
				else 
				move = (AvailWater(rpt) > 0.) ? Stuff(rpt)*ftmp/AvailWater(rpt) : 0.;
			}
		          
			Stuff(rpt) += move;
			AvailWater(rpt) += ftmp; 

			Stuff(pt) -= move;
			AvailWater(pt) -= ftmp;
		  
		} // end for equilibration loop
    
		//take care of the outflow from estuary
		if (HydMap(pt) == MOUTH && (ftmp = Elev(pt)+AvailWater(pt)-WATLEV) > 0) 
		{ 
			if (ftmp > AvailWater(pt)) ftmp = AvailWater(pt);
			  //if new level is below the water level, we remove all the water
			if (AvailWater(pt)>0)
				move = Stuff(pt)*ftmp/AvailWater(pt);
			else move = 0;  
			Stuff(pt) -= move;
			AvailWater(pt) -= ftmp;
			TotalOut += ftmp;
			TotNout += move;
		}  
		  	     
	} //end loop for OPEN Water areas
	

	
// printing results
	Ntot = 0;
		
	for( p = grid.first(); p; grid.next(p) )
	{
		const OrderedPoint& pt = grid.GetPoint(p);  
		Ntot += Stuff(pt);	// total amount of stuff in the area

		if( HydMap(pt) == GSTAT )
		{ 
			fprintf ( file, "\t%12.6f", (*swTot)(pt) );
			fprintf ( file, "\t%12.6f", AvailWater(pt) );
			fprintf ( fileN, "\t%12.6f", (*NTot)(pt) );
			fprintf ( fileN, "\t%12.6f", Stuff(pt) );
			if ( (*swTot)(pt) > 0 ) fprintf ( fileN, "\t%12.6f", (*NTot)(pt)/(*swTot)(pt) );
			else fprintf ( fileN, "\t    0.0     " );
		}
		
		if ( pt.row()==ROWOUT && pt.col()==COLOUT)
			fprintf ( fileN, "\t%12.6f", Stuff(pt) );
	}
	Nconc = (TotalOut > 0) ? TotNout/TotalOut : 0;

	fprintf ( fileN, "\t%12.6f\t%12.6f\t%12.6f", TotNout, Nconc, Ntot );
	fprintf ( file, "\t%12.4f", TotalOut );

// Handover Values of SWater to global variables (sorry) for printout for RALF
// and further analysis in UCode.cc(Goalfunction)

	n_out_sum = TotNout;
	n_out_conc = Nconc;
	
	fflush(file); 
	fflush(fileN); 
}


/***********************************************************************/