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

#include "SWater.h"

//*************************************************************************
// Fluxing procedure over steep areas, assuming that all the water 
// is moved downhill. Allows fluxing over a path of several cells. 
// 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).
// Fluxx - the precalculated rate of fluxing
// Stuff - is the stuff to be moved with water

void SWTransport3( CVariable& AvailWater, CVariable& HydMap, 
	CVariable& Slop, CVariable& Stuff )

{	
	static FILE *file, *fileN, *fileIn;
	static int date = 0;
	float TotalOut = 0.;
	float TotNout = 0.0;   // the amount of stuff that is fluxed out from the last cell
	float Nconc, Ntot;
	
	float df, w_sum, ftmp, move;
	Pix p, rp;

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

	char ccc[1000];

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

	if (fileIn == NULL)   // Read Data from file /data/SWData
	{
		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 );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &HSHEAD );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &WATLEV );
		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%f", &NDIF );

		fgets ( ccc, 1000, fileIn); 
		sscanf ( ccc, "%d", &OPENW );
		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 );
	}

	if (!ONOFF ) return;

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

	if (file == NULL) 
	{ // open output file
		CPathString pathName(Env::ArchivePath()); 
    		pathName.Add("WSOut");
		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%10dx%d\t", ir0, ic0 );
				fprintf ( fileN, "\t%10dx%d\t", ir0, ic0 );
			}
       
			if ( pt.row()==ROWOUT && pt.col()==COLOUT)
				fprintf ( fileN, "\t%10dx%d\t ", ROWOUT, COLOUT );
		
		}
        
		fprintf ( file, "\tTotal Outflow from area" );
		fprintf ( fileN, "\t Stuff flux from area\t Conc at outflow\t Total Stuff 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 nutrient
	                      
	static CVariable* NTot = NULL;	
	if(NTot == NULL ) { NTot = AvailWater.GetSimilarVariable("NTot"); }
	    // intermediate increment to nutrient

	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);
	                      
	Stuff.LinkEdges();
	AvailWater.LinkEdges();
	swFlux->Set(0.0);
	NFlux->Set(0.0);
	
// loop over the whole area
	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.;

		df = w0*flux_parm;

		(*swFlux)(pt) -= df;

		if (w0>transport_st)  nf = NDIF*Stuff(pt)*df/w0;
		
		(*NFlux)(pt) -= nf;

		Pix rp = p;
	    	
	// making the length of flow path dependent of the amount of water to move
		float stt = w0*w0; 
		float godo = 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;
				}
			}
		} // 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);

	
// 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 ( 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 );
	
	fflush(file); 
	fflush(fileN); 
}


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