/**************************************************************************/ // Hydrologic functions used to distribute surface water horizontally //@ Alexey Voinov -- voinov@cbl.umces.edu /**************************************************************************/ #include "SWater.h" //************************************************************************* // Simplest fluxing procedure over steep areas, assuming that fluxes are // already known (say defined by a Stella algorithm). Also moves 'Stuff' together // with water. A certain amount of water is taken from the donor cell based on the // predefined flux amount and added to the recepient cell. This is repeated over // the whole area several times. // 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 SWTransport1( CVariable& AvailWater, CVariable& HydMap, CVariable& Fluxx, CVariable& Stuff ) { 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; float 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); // SME: Sets grid ordering to default // Read control Data from file SWData if (fileIn == NULL) { CPathString pathName(Env::DataPath()); // SME: get path to Data directory 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 ); //dummy parm here fgets ( ccc, 1000, fileIn); sscanf ( ccc, "%f", &WATLEV ); //dummy parm here fgets ( ccc, 1000, fileIn); sscanf ( ccc, "%f", &NDIF ); fgets ( ccc, 1000, fileIn); sscanf ( ccc, "%d", &OPENW ); //dummy parm here 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()); // SME: get path for DriverOutput directory pathName.Add("NOut"); 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"); 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) ) // SME: loop over studyarea { const OrderedPoint& pt = grid.GetPoint(p); if( HydMap(pt) == GSTAT ) { int ir0 = pt.row(), ic0 = pt.col(); //SME: define cell coordinates 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 ", 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 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(); // make MAXDRUN iterations per function call to run water further downhill for (int iter=0; iter < MAXDRUN; iter++) { swFlux->Set(0.0); NFlux->Set(0.0); // loop over the whole area for ( p = grid.first(); p; grid.next(p) ) // SME: loop over studyarea { const OrderedPoint& pt = grid.GetPoint(p); // SME: pt same as p, but ordered float& w0 = AvailWater(pt); float nf = 0.; df = Fluxx(pt)*flux_parm; // define amount of water to be removed from cell if (df > w0) df = w0; (*swFlux)(pt) -= df; if ( (*swFlux)(pt) < -w0 ) (*swFlux)(pt) = -w0; if ( w0 > 0 ) nf = NDIF*Stuff(pt)*df/w0; // define amount of stuff to be removed if (nf > Stuff(pt)) nf = Stuff(pt); (*NFlux)(pt) -= nf; if ( -(*NFlux)(pt) > Stuff(pt) ) (*NFlux)(pt) = -Stuff(pt); rp = grid.LinkedPix( p ); // SME: gets downstream link rp if ( rp ) // if next cell is in studyarea { const OrderedPoint& rpt = grid.GetPoint(rp); (*swFlux)(rpt) += df; // add water to recipient cell (*NFlux)(rpt) += nf; // add stuff to recipient cell if (HydMap(rpt) == GSTAT) // GSTAT marks the gaging stations on River map. // This is to sum up all the water and all the material // that is moved through the cells with gaging stations { (*swTot)(rpt) += df; (*NTot)(rpt) += nf; } if (HydMap(rpt) == MOUTH) // MOUTH marks the mouth on River map { TotNout += nf; // amount removed through the MOUTH point TotalOut += df; } } } //end loop over all area AvailWater.AddData(*swFlux); // SME: based on calculated fluxes change values in cells Stuff.AddData(*NFlux); } //end loop over number of iterations // 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); } /***********************************************************************/