WatMgmt.c File Reference

Dynamic equations: Horizontal solutions of spatial raster/vector water management in canal network. More...

#include "watmgmt.h"

Include dependency graph for WatMgmt.c:

Go to the source code of this file.

Functions
void	Canal_Network_Init (float baseDatum, float *elev)
	Initialize the water managment network topology.
void	CanalReInit ()
	Re-initialize the depths and concentrations in canals under the Positional Analysis mode.
void	Run_Canal_Network (float *SWH, float *Elevation, float *MC, float *GW, float *poros, float *GWcond, double *NA, double *PA, double *SA, double *GNA, double *GPA, double *GSA, float *Unsat, float *sp_yield)
	Runs the water management network.
struct Chan *	ReadChanStruct (char *filename)
	Reads the information about canal reaches from data file.
void	ReadStructures (char *filename, float BASE_DATUM)
	Read attribute data on water control structures.
void	ReadStructFlow_PtSer (char *filename)
	Read time series data of flow at water control structures.
void	ReadStructTP_PtSer (char *filename)
	Read time series data of TP concentration at water control structures.
void	ReadStructTS_PtSer (char *filename)
	Read time series data of TS (salt/tracer) concentration at water control structures.
struct Schedule *	Read_schedule (char *sched_name, char *filename, float BASE_DATUM)
	Read target stage schedule as a recurring time-series graph.
void	Channel_configure (float *ElevMap, struct Chan *channel_first)
	Configure attributes of grid cells interacting with canal vectors.
void	getCanalElev (int chan_n)
	Assign elevation to canal using slope & distance from start point .
struct Cells *	MarkCell (int x, int y, int c_ind, float length, struct Cells *cell, int direction, int levee, int xm, int ym, int c_num, int *marked, float distTot)
	Mark cell function.
void	FluxChannel (int chan_n, float *SWH, float *ElevMap, float *MC, float *GWH, float *poros, float *GWcond, double *NA, double *PA, double *SA, double *GNA, double *GPA, double *GSA, float *Unsat, float *sp_yield)
	Runs the iterative algorithm over the network of canals.
float	f_Manning (float delta, float Water, float SW_coef)
	Surface water exchange between cell and canal.
float	f_Ground (float dh, float height, float GW_coef, float I_Length)
	Subsurface ground water exchange between cell and canal.
void	Flows_in_Structures (float *SWH, float *GW, float *poros, float *Elev, double *NA, double *PA, double *SA)
	Calling function to functions that determine flows through water control structures.
void	flowData_CanCan (struct Structure *structs, struct Structure *struct_hist_start, float *SWH, double *NA, double *PA, double *SA)
	Canal-to-Canal water control structure flows: Data-driven (historical/sfwmm).
void	flowCalc_CanCan (struct Structure *structs, float *SWH, float *GW, float *poros, float *Elev, double *NA, double *PA, double *SA)
	Canal-to-Canal water control structure flows: Calculated.
void	flowData_CelCel (struct Structure *structs, float *SWH, double *NA, double *PA, double *SA)
	Cell-to-Cell water control structure flows: Data-driven (historical/sfwmm).
void	flowCalc_CelCel (struct Structure *structs, float *SWH, float *Elev, double *NA, double *PA, double *SA)
	Cell-to-Cell water control structure flows: Calculated (rule-based).
void	flowData_CanCel (struct Structure *structs, struct Structure *struct_hist_start, float *SWH, double *NA, double *PA, double *SA)
	Canal-to-Cell water control structure flows: Data-driven (historical/sfwmm).
void	flowCalc_CanCel (struct Structure *structs, float *SWH, float *GW, float *poros, float *Elev, double *NA, double *PA, double *SA)
	Canal-to-Cell water control structure flows: Calculated (rule-based).
void	flowData_CelCan (struct Structure *structs, float *SWH, double *NA, double *PA, double *SA)
	Cell-to-Canal water control structure flows: Data-driven (historical/sfwmm).
void	flowCalc_CelCan (struct Structure *structs, float *SWH, float *GW, float *poros, float *Elev, double *NA, double *PA, double *SA)
	Cell-to-Canal water water control structure flows: Calculated (rule-based).
float	GetGraph (struct Schedule *graph, float x)
	Graph to read the time dependent schedules.
float	UTM2kmy (double UTM)
	Converter from meters in UTM to grid cell column location (zero origin).
float	UTM2kmx (double UTM)
	Converter from meters in UTM to grid cell row location (zero origin).
int	Wrdcmp (char *s, char *t)
	Compare two words.

---

Detailed Description

Dynamic equations: Horizontal solutions of spatial raster/vector water management in canal network.

This source file contains the raster-vector solutions associated w/ management infrastructure, fluxing water and constituents in the canal, levee, and water control structure network.
The major functions in the sequence in which they are found in this file:
Canal_Network_Init (main calling function to initialize canal network)
CanalReInit (position analysis mode, re-initializing canals depth/concentrations)
Run_Canal_Network (main calling function for the canal-levee-structure dynamics)
ReadChanStruct (read in data to define canal/levee network)
ReadStructures (read in data to define water control structure attributes and structure flows)
ReadStructFlow_PtSer (read in time series data of water flows at water control structures)
ReadStructTP_PtSer (read in time series data of TP concentrations at water control structures)
ReadStructTS_PtSer (read in time series data of TS (salt) concentrations at water control structures)
Read_schedule (read time series data for regulation schedules)
Channel_configure (configure the arrays of cells associated with canals/levees)
MarkCell (store canal-interacting cell coordinates and their flow attributes relative to levees)
FluxChannel (dynamic fluxes of water and solutes between canal vectors and cells)
f_Manning (overland flow equation as modified for canal-cell flows etc)
f_Ground (groundwater flow equation as modified for canal-cell flows etc)
Flows_in_Structures (calling function of the dynamic fluxes of water and solutes through structures)
flowData_CanCan (Canal-to-Canal water control structure flows: Data-driven (historical/sfwmm) )
flowData_CelCel (Cell-to-Cell water control structure flows: Data-driven (historical/sfwmm) )
flowData_CanCel (Canal-to-Cell water control structure flows: Data-driven (historical/sfwmm) )
flowData_CelCan (Cell-to-Canal water control structure flows: Data-driven (historical/sfwmm) )
flowCalc_CanCan (Canal-to-Canal water control structure flows: Calculated (rule-based) )
flowCalc_CelCel (Cell-to-Cell water control structure flows: Calculated (rule-based) )
flowCalc_CanCel (Canal-to-Cell water control structure flows: Calculated (rule-based) )
flowCalc_CelCan (Cell-to-Canal water water control structure flows: Calculated (rule-based) )
misc small utility functions

Note: documented with Doxygen, which expects specific syntax within special comments.

The Everglades Landscape Model (ELM).
last updated: July 2011

Definition in file WatMgmt.c.

---

Function Documentation

void Canal_Network_Init	(	float	baseDatum,
		float *	elev
	)

Initialize the water managment network topology.

Initialize/configure the network of canals/levees and and water control structures. The canal information stored in the global array Chan_list

Parameters:

	baseDatum	GP_DATUM_DISTANCE
	elev	SED_ELEV

Definition at line 96 of file WatMgmt.c.

References Chan::area, CanalOutFile, CanalOutFile_P, CanalOutFile_S, canDebugFile, Chan_list, ChanInFile, Channel_configure(), debug, H_OPSYS, Chan::ic_depth, Chan::ic_N_con, Chan::ic_P_con, Chan::ic_S_con, init_pvar(), MCopen, modelFileName, modelName, modelVers, Chan::N, nalloc(), Structure::next_in_list, num_chan, ON_MAP, OutputPath, Chan::P, ProjName, ReadChanStruct(), ReadStructures(), Chan::S, s0, s1, Structure::S_nam, SimAlt, SimModif, struct_first, structs, T, UNIX, usrErr(), WstructOutFile, WstructOutFile_P, and WstructOutFile_S.

Referenced by init_canals().

{ struct Structure *structs;
  int i,j;
  char filename[30];

  sprintf(filename,"CanalData");

/* Input the canals geometry file, return pointer to first vector               */
/* Configure the canals, marking the cells that interact with them              */ 
/* and storing the pointers to the cells for each of the canals in Chan_list array */
  Channel_configure(elev, ReadChanStruct (filename));  

  for( i = 0; i < num_chan; i++ )
  {
      Chan_list[i]->P = Chan_list[i]->ic_P_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
      Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
      Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
  }

  usrErr ("\tCanal geometry configured OK");
  
/* Read data on the water control structures, return pointer to the first structure */
  ReadStructures  (filename, baseDatum);

  ON_MAP[T(2,2)] = 0;
  
  MCopen = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MCopen"); /* open-water manning's coef, used in structure cell-cell bridge flow */
  init_pvar(MCopen,NULL,'f',0.05); /* this is a hard-wired parameter (value set to last argument, 0.05) - but do not do anything much with it now for the bridge flows */
  
 
 if ( debug > 0 )       /* Output the canal/levee modifications to ON_MAP */
 {
        sprintf( modelFileName, "%s/%s/Output/Debug/ON_MAP_CANAL.txt", OutputPath, ProjName );
   

  if ( ( ChanInFile = fopen( modelFileName, "w" ) ) ==  NULL )
  {
     printf( "Can't open the %s canal/levee debug file!\n", modelFileName );
     exit(-1) ;
  }
  fprintf ( ChanInFile, "ROWS=%d\nCOLUMNS=%d\nFORMAT=text\nINFO=\"Debug data: Overland flow restrictions due to levees. \"\nMISSING=0\n", s0, s1 );

  for ( i = 1; i <= s0 ; i++ ) 
    { for ( j = 1 ; j <= s1 ; j++ )
         fprintf( ChanInFile, "%4d\t", ON_MAP[T(i,j)]) ;
      fprintf ( ChanInFile, "\n" );
    }   
  fclose( ChanInFile ) ;
 }


/* Open file for writing structure flows */
    { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Canal/structsOut", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:structsOut", OutputPath, ProjName );
    }

        /* Open file to print structure flows to */
    if ( ( WstructOutFile = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }
/* Print The Header Line For This File */
    fprintf ( WstructOutFile,
              "%s %s %s %s scenario: Sum of water flows at structures (ft^3/sec)\n     Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
    structs = struct_first;
    while ( structs != NULL )
    {   fprintf( WstructOutFile, "%7s\t", structs->S_nam);       
            structs = structs->next_in_list;
    }

/* Open file for writing structure P concs */
    { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Canal/structsOut_P", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:structsOut_P", OutputPath, ProjName );
    }

        /*  Open file to print structure flows to */
    if ( ( WstructOutFile_P = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }
/* Print The Header Line For This File */
    fprintf ( WstructOutFile_P,
              "%s %s %s %s scenario: Flow weighted mean TP conc at structures ( mg/L) \n     Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
    structs = struct_first;
    while ( structs != NULL ) 
    {   fprintf( WstructOutFile_P, "%7s\t", structs->S_nam);      
            structs = structs->next_in_list; 
    }

/* Open file for writing structure Salin concs */
    { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Canal/structsOut_S", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:structsOut_S", OutputPath, ProjName );
    }

        /*  Open file to print structure flows to */
    if ( ( WstructOutFile_S = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }
/* Print The Header Line For This File */
    fprintf ( WstructOutFile_S,
              "%s %s %s %s scenario: Flow weighted mean tracer concentration at structures ( g/L) \n     Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
    structs = struct_first;
    while ( structs != NULL ) 
    {   fprintf( WstructOutFile_S, "%7s\t", structs->S_nam);      
            structs = structs->next_in_list; 
    }

/****/
    
        /* Open file to print canal stage info to */
    { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:CanalOut", OutputPath, ProjName );
    }

    if ( ( CanalOutFile = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }

    fprintf ( CanalOutFile, "%s %s %s %s scenario: Instantaneous water depth (meters, from bottom of canal) in canal Reaches\n    Date\t", &modelName, &modelVers, &SimAlt, &SimModif ); /* print the header line for this file */
    /*  channels with negative widths are reserved for new canals */
    /* v2.2 put a "R_" prefix in front of the canal ID# */
    for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile, "R_%d\t",Chan_list[i]->number);}
    

        /* Open file to print canal phosph info to */
    { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut_P", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:CanalOut_P", OutputPath, ProjName );
    }

    if ( ( CanalOutFile_P = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }

    fprintf ( CanalOutFile_P, "%s %s %s %s scenario: Instantaneous TP conc (mg/L) in canal Reaches\n    Date\t", &modelName, &modelVers, &SimAlt, &SimModif ); /* print the header line for this file */
    /* v2.2 put a "R_" prefix in front of the canal ID# */
    for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile_P, "R_%d\t",Chan_list[i]->number);}
    
        /* Open file to print canal salinity/tracer info to */
    { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut_S", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:CanalOut_S", OutputPath, ProjName );
    }

    if ( ( CanalOutFile_S = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }

    fprintf ( CanalOutFile_S, "%s %s %s %s scenario: Instantaneous tracer conc (g/L) in canal Reaches\n    Date\t", &modelName, &modelVers, &SimAlt, &SimModif ); /* print the header line for this file */
    /* v2.2 put a "R_" prefix in front of the canal ID# */
    for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile_S, "R_%d\t",Chan_list[i]->number);}
    
 

/* Open file for canal-cell flux debugging purposes */
  if( debug > 0 ) 
  { 
     if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Output/Debug/Can_debug", OutputPath, ProjName );
      else sprintf( modelFileName, "%s%s:Output:Can_debug", OutputPath, ProjName );
    

        /* Open file for debugging and error/warnings on canal-cell flows, structure flows */
    if ( ( canDebugFile = fopen( modelFileName, "w" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName );
       exit(-1) ;
    }
  
     fprintf ( canDebugFile, "Depending on level of the debug requested, printed here are data on canal-cell flows and excessive demands on data-driven structure flows. \n" ); /* v2.5 made use of this file pointer */
     fflush(canDebugFile); 
    
  }
  
  return;

}

Here is the call graph for this function:

void CanalReInit

(

)

Re-initialize the depths and concentrations in canals under the Positional Analysis mode.

Definition at line 292 of file WatMgmt.c.

References Chan::area, Chan_list, Chan::ic_depth, Chan::ic_N_con, Chan::ic_P_con, Chan::ic_S_con, Chan::N, num_chan, Chan::P, Chan::S, and Chan::wat_depth.

Referenced by reinitCanals().

{
    int i;
    for( i = 0; i < num_chan; i++ )

    {
        Chan_list[i]->wat_depth = Chan_list[i]->ic_depth;
        Chan_list[i]->P = Chan_list[i]->ic_P_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
        Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
        Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
         
    }
}

void Run_Canal_Network	(	float *	SWH,
		float *	Elevation,
		float *	MC,
		float *	GW,
		float *	poros,
		float *	GWcond,
		double *	NA,
		double *	PA,
		double *	SA,
		double *	GNA,
		double *	GPA,
		double *	GSA,
		float *	Unsat,
		float *	sp_yield
	)

Runs the water management network.

Invoke the calls to water control structure flows and the canal-cell interactions.

Parameters:

	SWH	SURFACE_WAT
	Elevation	SED_ELEV
	MC	HYD_MANNINGS_N
	GW	SAT_WATER
	poros	HP_HYD_POROSITY
	GWcond	HYD_RCCONDUCT
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT
	GNA	DINdummy
	GPA	TP_SED_WT
	GSA	SALT_SED_WT
	Unsat	UNSAT_WATER
	sp_yield	HP_HYD_SPEC_YIELD

Definition at line 330 of file WatMgmt.c.

References Chan::area, Chan::basin, CanalOutFile, CanalOutFile_P, CanalOutFile_S, canDebugFile, canPrint, canstep, Chan_list, simTime::da, debug, Flows_in_Structures(), FluxChannel(), FMOD(), hyd_iter, simTime::IsBudgEnd, simTime::IsSimulationEnd, simTime::mo, N_c_iter, num_chan, P, Chan::P_con, printchan, Chan::S_con, SensiOn, SimTime, Chan::sumHistIn, Chan::sumHistOut, Chan::sumRuleIn, Chan::sumRuleOut, TOT_P_CAN, TOT_S_CAN, TOT_VOL_CAN, Chan::wat_depth, WstructOutFile, WstructOutFile_P, WstructOutFile_S, and simTime::yr.

Referenced by horizFlow().

{ int i;
 float CH_vol;

 if ( canPrint && (( N_c_iter++ % hyd_iter ) == 0.0) ) {
        if (SensiOn) { 
            printchan = (SimTime.IsSimulationEnd) ? (1) : (0); /* flag to indicate to print canal/struct data */
         }
         else printchan = 1;
 }
 else  printchan = 0;

 if (printchan == 1) {
     
     fprintf( WstructOutFile, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
     fprintf( WstructOutFile_P, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
/*      fprintf( WstructOutFile_N, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] ); */
     fprintf( WstructOutFile_S, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );

     fprintf (CanalOutFile, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
     fprintf (CanalOutFile_P, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
/*      fprintf (CanalOutFile_N, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] ); */
     fprintf (CanalOutFile_S, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
 }
    
  
 for ( i = 0; i < num_chan; i++ )
 {   /* set the sum of historical (or other-model) and Rule-based flows to/from a canal per iteration to 0 */   
     Chan_list[i]->sumHistIn = Chan_list[i]->sumHistOut = 0.0;
     Chan_list[i]->sumRuleIn = Chan_list[i]->sumRuleOut = 0.0;
 }

/* calculate flows through water control structures */
 Flows_in_Structures ( SWH, GW, poros, Elevation, NA, PA, SA ); 
  
/* output  header for canal-cell flux output canal debug file, at relatively high debug levels (writes a LOT of stuff) */
 if( debug > 3 ) 
 {
     fprintf( canDebugFile, "N%3d  WatDepth     Pconc     Sconc        flux_S        flux_G        flux_L           Qin          Qout #_iter      Area    error \n", N_c_iter );
      
 }

/* Flux the water along a canal using the relaxation procedure */
  
 for( i = 0; i < num_chan; i++ )
     if (Chan_list[i]->width > 0) /* a neg width canal is skipped */

     { 
         FluxChannel ( i, SWH, Elevation, MC, GW, poros, GWcond, NA, PA, SA, GNA, GPA, GSA, Unsat, sp_yield );
             /* calculate total canal volume after canal fluxes */
         CH_vol = Chan_list[i]->area*Chan_list[i]->wat_depth; 

             /* sum the storages on iteration that performs budget checks */
         if ( !( FMOD(N_c_iter,(1/canstep) )) && (SimTime.IsBudgEnd ) ) { 
             TOT_VOL_CAN[Chan_list[i]->basin] += CH_vol;
             TOT_VOL_CAN[0] += CH_vol;
             TOT_P_CAN[Chan_list[i]->basin] += Chan_list[i]->P;
             TOT_P_CAN[0] += Chan_list[i]->P;
             TOT_S_CAN[Chan_list[i]->basin] += Chan_list[i]->S;
             TOT_S_CAN[0] += Chan_list[i]->S;
         }
                                                                                               
     
         if ( printchan )
         {
                 /* report N&P concentration in mg/L (kg/m3==g/L * 1000 = mg/L) */
             Chan_list[i]->P_con = (CH_vol > 0) ? (Chan_list[i]->P/CH_vol*1000.0) : 0.0;
/*              Chan_list[i]->N_con = (CH_vol > 0) ? (Chan_list[i]->N/CH_vol*1000.0) : 0.0; */
             Chan_list[i]->S_con = (CH_vol > 0) ? (Chan_list[i]->S/CH_vol       ) : 0.0;
 
             /* v2.2 increased decimal places in output (from 6.3f to 7.4f) */
             fprintf( CanalOutFile, "%6.2f\t", Chan_list[i]->wat_depth );
             fprintf( CanalOutFile_P, "%7.4f\t", Chan_list[i]->P_con );
/*          fprintf( CanalOutFile_N, "%7.4f\t", Chan_list[i]->N_con ); */
             fprintf( CanalOutFile_S, "%7.4f\t", Chan_list[i]->S_con ); 
         }
     } 
 if ( printchan )
 {
     fflush(CanalOutFile);
     fflush(CanalOutFile_P);
/*       fflush(CanalOutFile_N); */
     fflush(CanalOutFile_S);
 }
  
 return;
}

Here is the call graph for this function:

struct Chan* ReadChanStruct

(

char *

filename

)

[read]

Reads the information about canal reaches from data file.

A canal in the model is termed a canal reach. A canal reach is defined by one or more straight segments. A segment is defined by two points with geographic coordinates. Thus, a curved canal reach is comprised of multiple, small segments.

Parameters:

filename

Canal/levee attributes data file.

Returns:

chan_first Pointer to first canal data struct

Definition at line 435 of file WatMgmt.c.

References Chan::basin, basins, basn_list, basndef::basnTxt, C_F, ChanInFile, CHANNEL_MAX_ITER, CHo, Chan::cond, Chan::depth, Chan::edgeMann, F_ERROR, basndef::family, Chan::family, H_OPSYS, Chan::ic_depth, Chan::ic_N_con, Chan::ic_P_con, Chan::ic_S_con, Chan::levee, MINDEPTH, modelFileName, modelName, ModelPath, msgStr, Chan::next_in_list, Chan_segm::next_segm, num_chan, numBasn, Chan::number, basndef::parent, Chan::parent, ProjName, Chan::roil, Chan::roir, Chan::segments, SimAlt, UNIX, usrErr(), UTM2kmx(), UTM2kmy(), UTM_EOrigin, UTM_NOrigin, Chan::wat_depth, Chan::width, Chan_segm::x0, Chan_segm::x1, Chan_segm::y0, and Chan_segm::y1.

Referenced by Canal_Network_Init().

{ 
 struct Chan *channels, *chan_first, *chan_last;
 struct Chan_segm *C_segm, *last_segm, *next_segm;
 int i, index, firstSegm, firstCanal = 1;
 char ss[622], scenario[20], modnam[20], bas_txt[8];
 float C_x, C_y, C_x1, C_y1;

 
 int ibas, foundBasn;
 

 { if(H_OPSYS==UNIX) 
     sprintf( modelFileName, "%s/%s/Data/%s.chan", ModelPath, ProjName, filename );
 else sprintf( modelFileName, "%s%s:Data:%s.chan", ModelPath, ProjName, filename );
 }

/* Open file with canals data */
 if ( ( ChanInFile = fopen( modelFileName, "r" ) ) ==  NULL )
 {
     sprintf( msgStr,"Can't open %s file! ",modelFileName ) ; usrErr(msgStr);
     exit(-1) ;
 }
  
/* Allocate memory for canal reach segment structure. A canal segment is the straight canal segment */
/* defined in the canal data file by the coordinates of its beginning and ending points.   */
/* Each canal reach is comprised of one or more canal segments          */
 if ( (C_segm = (struct Chan_segm *) malloc( (size_t) sizeof( struct Chan_segm ))) == NULL )
 {
     usrErr( "No memory for first canal reach segment\n " ) ;
     exit( -2 ) ;
 };

 num_chan = 0;
   
/* Read the general canal network information */ 
 fgets( ss, 620, ChanInFile ); sscanf( ss, "%d", &UTM_EOrigin );
 fgets( ss, 620, ChanInFile ); sscanf( ss, "%d", &UTM_NOrigin );
 fgets( ss, 620, ChanInFile );  /* skip comment line */
 fgets( ss, 620, ChanInFile ); sscanf( ss,"%s %s %d %d", &modnam, &scenario, &CHo, &CHANNEL_MAX_ITER );
 if (strcmp(scenario,SimAlt) != 0) {
     sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
             scenario, modelFileName, SimAlt); usrErr(msgStr);
             exit(-1);
 }
 if (strcmp(modnam,modelName) != 0) {
     sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
             modnam, modelFileName, modelName); usrErr(msgStr);
             exit(-1);
 }
  
 fgets( ss, 620, ChanInFile );   /* skip comment line */
 fgets( ss, 620, ChanInFile ); sscanf ( ss, "%f %f %f ", &F_ERROR, &MINDEPTH, &C_F );
     /* the C_F is used only for sensitivity experiments */
 fgets( ss, 620, ChanInFile );   /* skip comment line */
   
/* Read canal descriptors until none left */   
 while ( fgets( ss, 620, ChanInFile ) != NULL && !feof( ChanInFile ) )
 { 
     if ( *ss == '#' )  /* "#" identifies each new canal in the data file */
     { /* Allocate memory for canal reach structure */
         if ( (channels = (struct Chan *) malloc( (size_t) sizeof( struct Chan ))) == NULL )
         {
             printf( "No memory for first canal\n " ) ;
             exit( -2 ) ;
         };
                
             /* Read canal information:  - canal number, 
                - levee marker (1 - left levee, -1 - right, 0 - none, 2 - both sides),
                - ranges of cell interaction to the left and to the right - functionality removed, should ALWAYS=1,
                - canal depth and width,
                - seepage coefficient through the levee,
                - initial S and P concentrations (doubles) in the canal water,
                - initial water depth in canal (can be >, = or < than canal depth 
                - Manning's n associated w/ edge of canal, to accomodate topographic lip/berm and/or denser veg along canal length
                - hydrologic basin that canal exchanges surface water with */
         i = sscanf( ss+1, "%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%s",
                     &channels->number, 
                     &channels->levee,
                     &channels->roil,   &channels->roir, 
                     &channels->depth,  &channels->width,  
                     &channels->cond,   
                     &channels->ic_S_con,  &channels->ic_P_con,  
                     &channels->wat_depth, &channels->edgeMann, &bas_txt);

         foundBasn=0;
         for (ibas=numBasn;ibas>0;ibas--) {
             basins = basn_list[ibas];
             if(strcmp(bas_txt,basins->basnTxt)==0) {
                 channels->basin = ibas;
                 foundBasn=1;
                 channels->family = basn_list[ibas]->family;
                 channels->parent = basn_list[ibas]->parent;
                     /* canals may ONLY belong to parent hydro basins, not to indicator regions */
                 if (!channels->parent) {
                     printf ("Sorry - %s's canal %d's basin is not a parent hydrologic basin! Please assign to a hydro basin in the CanalData.chan file.\n",
                            modelName,channels->number); exit (-1);}     
             }
         }
         if (foundBasn==0) { printf ("Sorry - %s's canal %d's basin does not exist! Please define the basin in the basinIR file.\n",
                            modelName,channels->number); exit (-1);}                

         channels->ic_P_con *= 0.001; /* initial input value in mg/L, convert to g/L==kg/m3 */
         channels->ic_N_con = 0.00; /* v2.2 "removed" N */ /* initial input value in mg/L, convert to g/L==kg/m3 */
         channels->ic_S_con *= 1.0; /* initial input value in g/L, no conversion */
             /* store the initial depth for re-initializing under Positional Analysis mode */
         channels->ic_depth = channels->wat_depth;
         

         channels->segments = C_segm; 
                        
         if (i < 10)    /* check for the number of input fields read */
         {printf ( " Error in CanalData.chan: reach input %d\n", (num_chan+1) ); exit (-1);}
                        
                            
         num_chan++;            /*count the number of canal reaches (not canal segments) */
         firstSegm = 1;
                        
         if (firstCanal) 
         {
             firstCanal = 0;  
             chan_first = channels;
             chan_last = channels;
             continue; 
         }

         last_segm->next_segm = NULL;
         chan_last->next_in_list = channels;
         chan_last = channels;     
     }
     else  
     { /* read the pairs of coordinates for the consecutive canal segments */
         i = sscanf ( ss, "%f %f", &C_x1, &C_y1 );
         if (i < 2)     /* check for the number of input fields read */
         {printf ( " Error in CanalData.chan coords: Reach %d  \n", (num_chan+1) ); exit (-1);}
                
             /* Convert from UTM to km units with 0,0 origin */
         C_x = UTM2kmx (C_x1);  C_y = UTM2kmy (C_y1);           
         C_segm->x1 = C_x;  C_segm->y1 = C_y;
                        
         if ( firstSegm )                    
         {  C_segm->x0 = C_x; C_segm->y0 = C_y; 
         firstSegm = 0;
         }
         else 
         {  
                                /* Allocate memory for canal reach segment structure */
             if ( (next_segm = (struct Chan_segm *) malloc( (size_t) sizeof( struct Chan_segm ))) == NULL )
             {  printf( "No memory for canal reach segment\n " ) ;
             exit( -2 ) ;
             };
             C_segm->next_segm = next_segm;
             last_segm = C_segm;
             C_segm = next_segm;
             C_segm->x0 = C_x; C_segm->y0 = C_y;
         }                       
     }  
 }/* end loop for reading */
  
 chan_last->next_in_list = NULL;
 last_segm->next_segm = NULL;
 free ((char *) C_segm);
 fclose (ChanInFile);
   
 usrErr( "\tCanal locs/attributes read OK" );
   
 return (chan_first);
}

Here is the call graph for this function:

void ReadStructures	(	char *	filename,
		float	BASE_DATUM
	)

Read attribute data on water control structures.

Parameters:

	filename	Water control structures attributes file name
	BASE_DATUM	GP_DATUM_DISTANCE

Variables local to function

disAggNum Total # of disaggregated water control structures
IDhist ID of historical data set
lincnt Flow data input file record # starting at the point where the model reads values
structName Water control structure name in time series file
histTP Input data field that either has TP conc (ug/L=ppb) that is constant across time or has (any) string to indicate the use of TimeSeries data from a datafile
histTS Input data field that either has TS (salt) conc (g/L=ppt) that is constant across time or has (any) string to indicate the use of TimeSeries data from a datafile

Definition at line 611 of file WatMgmt.c.

References Structure::canal_fr, Structure::canal_to, Structure::cell_HW, Structure::cell_i_fr, Structure::cell_i_HW, Structure::cell_i_to, Structure::cell_i_TW, Structure::cell_j_fr, Structure::cell_j_HW, Structure::cell_j_to, Structure::cell_j_TW, Structure::cell_TW, Chan_list, ChanInFile, CHo, disAggNum, HistData::flag, Structure::flag, Structure::flagHWsched, Structure::flagHWtide, Structure::flagTWsched, Structure::flagTWtide, H_OPSYS, Hist_data, Structure::histID, Structure::HW_graph, Structure::HW_stage, IDhist, isFloat(), isInteger(), MAX_H_STRUCT, modelFileName, modelName, ModelPath, msgStr, Structure::next_in_list, num_chan, num_struct_hist, NumScheds, numTPser, numTSser, ON_MAP, ProjName, Read_schedule(), ReadStructFlow_PtSer(), ReadStructTP_PtSer(), ReadStructTS_PtSer(), s0, s1, HistData::S_nam, Structure::S_nam, Scip(), SimAlt, SimModif, Structure::SrcDest, Structure::str_cell_i, Structure::str_cell_j, struct_first, structs, Structure::Sum_N, Structure::Sum_P, Structure::Sum_S, Structure::SumFlow, T, TAB, Structure::TN, Structure::TNser, Structure::TP, Structure::TPser, Structure::TS, Structure::TSser, Structure::TW_graph, Structure::TW_stage, UNIX, usrErr(), Structure::w_coef, Chan::width, and WriteMsg().

Referenced by Canal_Network_Init().

{
struct Structure *structs, *struct_next, *struct_last;

 char ss[622], *s;
 char sched_name[20];
 int i, j, k;
 int IDhist[MAX_H_STRUCT]; 
 int lincnt; 
 double Jdate_read;
 char lline[2001], *line;
 int yyyy, mm, dd;
 char structName[20]; 
 char scenario[20], s_modname[20];
 char histTP[5]; 
 char histTS[5]; 
 char srcType[5], destType[5]; /* v2.8 */
 char prefixTide[]="tid_"; /* v2.8 - the Headwater and Tailwater schedules have the 'tid_' prefix (in CanalData.graph) if they are tidally based, this will be used to check for such a schedule */
 
 num_struct_hist = 0;
 numTPser = numTSser = disAggNum = 0; /* v2.6 */
 
 { if(H_OPSYS==UNIX) 
     sprintf( modelFileName, "%s/%s/Data/%s.struct", ModelPath, ProjName, filename );
 else sprintf( modelFileName, "%s%s:Data:%s.struct", ModelPath, ProjName, filename );
 }

/* Open file with structures data */
 if ( ( ChanInFile = fopen( modelFileName, "r" ) ) ==  NULL )
 {
     printf( "Can't open %s file!\n", modelFileName ) ;
     exit(-1) ;
 }
    
/* Allocate memory for structures */
 if ( (structs = (struct Structure *) malloc( (size_t) sizeof( struct Structure ))) == NULL )
 {
     printf( "No memory for first structure\n " ) ;
     exit( -2 ) ;
 }
 struct_first = structs;

 fgets( ss, 620, ChanInFile );
 sscanf(ss, "%s\t%s\t%s\t%s",&scenario,&scenario,&s_modname,&scenario); /*  2 junk strings followed by the model name and actual scenario ID */
 if (strcmp(scenario,SimAlt) != 0) {
         sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
                 scenario, modelFileName, SimAlt); usrErr(msgStr);
     exit(-1);
  }
 if (strcmp(s_modname,modelName) != 0) {
         sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) in the canal data file!",
                 s_modname, modelFileName, modelName); usrErr(msgStr);
     exit(-1);
  }
 
/* Skip 2'nd header line */
 fgets( ss, 620, ChanInFile );


/* Read structure descriptors until none left */   
 while ( fgets( ss, 620, ChanInFile ) != NULL && !feof( ChanInFile ) )
 { 
     s = ss; 
         
         
     sscanf( s, "%d", &structs->flag ); 
/* structure flag: <0 - skip structure info, structure not operating; */
/*                 >0 - historical/other-data data array, flag=1 normally in this case, but is >1 for aggregated SFWMM structures; */
/*                  0 - rule-based structure   */
     if ( structs->flag < 0 ) continue;
     if ( structs->flag > 0 && structs->flag <10 ) {
         structs->histID = ++num_struct_hist;
         Hist_data[structs->histID-1].flag = structs->flag;
     }
     
         
/*    assign integer to the historical data ID, and increment the counter for structures with historical/other-data data to be read from datafile. */
/*    flags >= 20 represent structures (such as the S-11s) that were aggregated in SFWMM output */
/*    and need to be assigned disaggregated flows for separate ELM structures */
/*    Ex: SFWMM structure name S11 may have flag = 2 in CanalData.struct.  That S11 flow holds the sum of S-11A-C flows. */
/*    The flags on S-11A, S-11B, and S-11C would all have been assigned 20, 20, and 20 in the CanalData.struct file */

/*           NO LONGER NECESSARY (old note, kept here in v2.2+): as a temporary measure until the SFWMM output .dss file is fixed, the structures that have calculated "other" partitions */
/*             (S140,S7,S8) are assigned historical flags of 9 (they are in the dss file read by a SAS routine, the flag */
/*             of 9 allows them to remain in the file, but the structure is turned off because it has flag >1)  */
         
     
     s = Scip( s, TAB );
         
     if ( *s != TAB ) {
             sscanf ( s, "%s", &structs->S_nam );
             if ( structs->flag >0 && structs->flag <10) 
                 sscanf ( s, "%s", &Hist_data[structs->histID-1].S_nam );
             
                 /* for solute concentration of flows thru structures, */
                 /* data use flag XXser: -1== simulated, 0==static historical, 1= time-varying historical */
             s = Scip( s, TAB );
             if ( *s == TAB ) structs->TPser = -1; 
             else {
                 sscanf ( s, "%s", &histTP ); /* read the field that either has TP conc (ug/L=ppb) that is constant across time
                                                   or has (any) string to indicate the use of TimeSeries data from a datafile (read later)*/
                 if (isInteger(histTP) ) {
                      structs->TPser = 0; 
                      structs->TP = atof(histTP)*1.0e-6;  /* conc (ug/L=>g/L) that is constant across time */
                 }
                 else { structs->TPser = 1; numTPser++;} /* otherwise, we will be reading time series data later */
             }
             
             s = Scip( s, TAB );
             if ( *s == TAB ) structs->TNser = -1; /* value to indicate use of simulated, not historical, data */
             else {
                 sscanf ( s, "%f", &structs->TN );
                 structs->TNser = 0;
                 structs->TN *= 1.0e-6;/* read the TN conc (ug/L=>g/L) that is constant across time */
             }
             

             s = Scip( s, TAB );
             if ( *s == TAB ) structs->TSser = -1; /* value to indicate use of simulated, not historical, data */
             /* v2.2 added code to allow use of TimeSeries input data in addition to old method of a constant concentration */
             else {
                  sscanf ( s, "%s", &histTS ); /* read the field that either has TS (salt) conc (g/L=ppt) that is constant across time
                                                   or has (any) string to indicate the use of TimeSeries data from a datafile (read later)*/
                 if (isFloat(histTS) ) {
                      structs->TSser = 0; 
                      structs->TS = atof(histTS)*1.0;  /* conc (no conversion, g/L) that is constant across time */
                 }
                 else { structs->TSser = 1; numTSser++;} /* otherwise, we will be reading time series data later */

             }
             
             s = Scip( s, TAB );
             if ( *s == TAB ) structs->str_cell_i = 0;  /* non 0 if elevation at the structure location is desired */
             else {
                 sscanf( s, "%d", &structs->str_cell_i );
                 structs->str_cell_i++; /* this converts the data map coords to model coords (map+1) */
             }
             
             s = Scip( s, TAB );

             if ( *s == TAB ) structs->str_cell_j = 0;  /* non 0 if elevation at the structure location is desired */
             else {
                 sscanf( s, "%d", &structs->str_cell_j );
                 structs->str_cell_j++;
             }

             if ( structs->str_cell_i != 0 && !ON_MAP[T(structs->str_cell_i, structs->str_cell_j)] )
             {  printf ( "Location of water control structure %s is off-map.\n", structs->S_nam );
             exit (-2);                 /*check that cells obtained are within on-map boundaries */
             }
             
             
     }
         
     s = Scip( s, TAB );

     if ( *s == TAB ) structs->canal_fr  = 0;   /* non 0 if the structure involves a donor canal */
     else sscanf( s, "%d", &structs->canal_fr  );
     s = Scip( s, TAB );
     if ( *s == TAB) structs->canal_to = 0;             /* non 0 if the structure involves a receiving canal */
     else sscanf( s, "%d", &structs->canal_to );
     s = Scip( s, TAB );

     if ( structs->canal_fr  - CHo >= num_chan || structs->canal_to - CHo >= num_chan )
     {  printf ( "Canal from/to doesn't exist for water control structure %s\n", structs->S_nam );
     exit (-2);                 /* check that canals read are consistent with canals data */
     }           
/*  check to ensure canals associated with structures are valid */
     if ( structs->canal_fr != 0 && Chan_list[structs->canal_fr- CHo]->width <= 0.1 )
     { printf ( "The donor canal %d has been turned off for water control structure %s\n", structs->canal_fr, structs->S_nam );
     exit (-2);                 
     }           
     if ( structs->canal_to != 0 && Chan_list[structs->canal_to- CHo]->width <= 0.1 )
     {  printf ( "The recipient canal %d has been turned off for water control structure %s\n", structs->canal_to, structs->S_nam );
     exit (-2);                 
     }           
         
     if (*s == TAB) structs->cell_i_fr = 0;  /* non 0 if the structure involves a donor cell */
     else 
     { sscanf( s, "%d", &structs->cell_i_fr );
     structs->cell_i_fr++;
     }   
     s = Scip( s, TAB );
     if (*s == TAB) structs->cell_j_fr = 0;
     else 
     { sscanf( s, "%d", &structs->cell_j_fr );
     structs->cell_j_fr++;
     }
     s = Scip( s, TAB );
         

     if ( structs->cell_i_fr > 2 && !ON_MAP[T(structs->cell_i_fr, structs->cell_j_fr)] )
     {  printf ( "Donor cell for water control structure %s is off-map.\n", structs->S_nam );
     exit (-2);                 /*check that cells obtained are within on-map boundaries */
     }

         /* if valid donor cell, mark it on ON_MAP */
     if ( structs->cell_j_fr ) ON_MAP[T(structs->cell_i_fr, structs->cell_j_fr)] += 100;
         
     if (*s == TAB) structs->cell_i_to = 0;  /* non 0 if the structure involves a recipient cell */
     else 
     { sscanf( s, "%d", &structs->cell_i_to );
     structs->cell_i_to++;
     }   
         
     s = Scip( s, TAB );
     if (*s == TAB) structs->cell_j_to = 0;
     else
     { sscanf( s, "%d", &structs->cell_j_to );
     structs->cell_j_to++;
     }   


     if ( structs->cell_i_to > 2 && !ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] )
     {  printf ( "Recipient cell for water control structure %s is off-map.\n", structs->S_nam );
     exit (-2);                 /*check that cells obtained are within on-map boundaries */
     }

         /* if valid recieving cell, mark it on ON_MAP */
         /* this used to check for cell_j_to > 2, not putting recipient LEC cells ON_MAP */
     if ( structs->cell_j_to ) ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] += 100;

/* the following seven field of data are for rule-based water management */

     /* read the targeted HeadWater stage info */
     s = Scip( s, TAB );                
     if (*s == TAB) {
        structs->flagHWsched = -1; /* v2.8 */
     }
     else if ( isalpha (*s) ) {  /* if there is a graph for the schedule */
         sscanf( s, "%s", sched_name );
         structs->HW_graph = Read_schedule ( sched_name, filename, BASE_DATUM );
         structs->flagHWsched = 1; /* v2.8 this is a time-varying target */
         structs->flagHWtide = ( strncmp(sched_name, prefixTide ,4) == 0) ? (1) : (0); /* true if this is a tidally-based target  */
     } 
     else   { sscanf( s, "%f", &structs->HW_stage );
        structs->HW_stage += BASE_DATUM;
        structs->flagHWsched = 0; /* v2.8 this is a constant target TODO - this should not be available as an option */
     }
                
                
         
     /* read the targeted TailWater stage info */
     s = Scip( s, TAB );
     if (*s == TAB) {
        structs->flagTWsched = -1; /* v2.8 */
     }
     else if ( isalpha (*s) ) {  /* if there is a graph for the schedule */
         sscanf( s, "%s", sched_name );
         structs->TW_graph = Read_schedule ( sched_name, filename, BASE_DATUM );
         structs->flagTWsched = 1; /* v2.8 this is a time-varying target */
         structs->flagTWtide = ( strncmp(sched_name, prefixTide ,4) == 0) ? (1) : (0); /* true if this is a tidally-based target  */
     } 
     else       { sscanf( s, "%f", &structs->TW_stage );
        structs->TW_stage += BASE_DATUM; /* prior to v2.8, the input data was designed to be negative values in the dbase for tailwater - this is NOT the case any longer */
        structs->flagTWsched = 0; /* v2.8 this is a constant target TODO - this should not be available as an option */
     }

     /* read the northing (row) coordinate of the location where the HeadWater stage will be checked */
     s = Scip( s, TAB );
     if (*s == TAB) structs->cell_i_HW = 0;
     else               
     { sscanf( s, "%d", &structs->cell_i_HW );
     structs->cell_i_HW++;
     }   

     /* read the easting (col) coordinate of the location where the HeadWater stage will be checked */
     s = Scip( s, TAB );
     if (*s == TAB) structs->cell_j_HW = 0;
     else 
     { sscanf( s, "%d", &structs->cell_j_HW );
     structs->cell_j_HW++;
     }   

     if ( structs->cell_i_HW > s0 || structs->cell_j_HW > s1 )
     {  printf ( "Error in definition of HW location of water control structure %s", structs->S_nam );
     exit (-2); /*check that cell obtained is within array boundaries */
     }
     else
     { (structs->cell_i_HW != 0) ? (structs->cell_HW = T(structs->cell_i_HW, structs->cell_j_HW) ) : (0); /* v2.8 determine this location now, instead of during dynamic calcs */
     }
     
//     if ( ( structs->cell_HW  != 0) && (structs->flagHWsched < 0) ) /* v2.8 TODO will need to clean up the dbase to remove HW and TW cell locations, for structures that have such fields filled in but are data-driven */
//     {  printf ( "Error: lacking target stage or stage sched for HW location of water control structure %s", structs->S_nam );
//     exit (-2);       
//     }
     

     /* read the northing (row) coordinate of the location where the TailWater stage will be checked */
     s = Scip( s, TAB );
     if (*s == TAB) structs->cell_i_TW = 0;
     else               
     { sscanf( s, "%d", &structs->cell_i_TW );
     structs->cell_i_TW++;
     }   

     /* read the easting (col) coordinate of the location where the TailWater stage will be checked */
     s = Scip( s, TAB );
     if (*s == TAB) structs->cell_j_TW = 0;
     else 
     { sscanf( s, "%d", &structs->cell_j_TW );
     structs->cell_j_TW++;
     }   

     if ( structs->cell_i_TW > s0 || structs->cell_j_TW > s1 )
     {  printf ( "Error in definition of TW location of water control structure %s", structs->S_nam );
     exit (-2); /*check that cell obtained is within array boundaries */
     }
     else
     { (structs->cell_i_TW != 0) ? (structs->cell_TW = T(structs->cell_i_TW, structs->cell_j_TW) ) : (0); /* v2.8 determine this location now, instead of during dynamic calcs */
     }

//     if ( ( structs->cell_TW  != 0) && (structs->flagTWsched < 0) ) /* v2.8 TODO will need to clean up the dbase to remove HW and TW cell locations, for structures that have such fields filled in but are data-driven */
//     {  printf ( "Error: lacking target stage or stage sched for TW location of water control structure %s", structs->S_nam );
//     exit (-2);       
//     }
     
     /* read the flow coefficient, that was derived in dbase to allow flow calcs via a simple weir-flow eqn */
     s = Scip( s, TAB );
     if (*s == TAB || *s == '?') structs->w_coef = 0;
     else sscanf( s, "%f", &structs->w_coef );
        
     struct_last = structs;
          
         /* Allocate memory for next structure */
     if ( (struct_next = (struct Structure *) malloc( (size_t) sizeof( struct Structure ))) == NULL )
     {
         printf( "Failed to allocate memory for next structure (%s)\n ", structs->S_nam ) ;
         exit( -2 ) ;
     };
     
/* v2.8 redesign of structure flow choices */
        /* assemble the attribute of the Source and Destination raster-vector type ('Can'al or 'Cel'l) */
    if ( structs->canal_fr  != 0 ) sprintf( srcType, "Can"); 
    if ( structs->canal_to  != 0 ) sprintf( destType, "Can"); 
    if ( structs->cell_i_fr != 0 && structs->cell_j_fr != 0 ) sprintf( srcType, "Cel"); 
    if ( structs->cell_i_to != 0 && structs->cell_j_to != 0 ) sprintf( destType, "Cel"); 
    sprintf( structs->SrcDest, "%s%s", srcType, destType ); 


     /* initialize the sum of structure  flows and constit solute in flows */
     structs->SumFlow = 0.0; 
     structs->Sum_P = 0.0; 
     structs->Sum_N = 0.0; 
     structs->Sum_S = 0.0; 

     structs->next_in_list = struct_next;
     structs = struct_next;
 }
   
 struct_last->next_in_list = NULL;
 free ((char *)structs);
 fclose (ChanInFile);
   

 sprintf(msgStr,"\tUsing %d target stage schedules for rule-based flows", NumScheds);  usrErr(msgStr);
 usrErr( "\tWater control structure locs/attributes read OK" );
 
 /* v2.6 - encoded the reading of these data into separate functions - no change in data or results */
 if (num_struct_hist>0) ReadStructFlow_PtSer(filename);
 if (numTPser>0) ReadStructTP_PtSer(filename);
 if (numTSser>0) ReadStructTS_PtSer(filename);
 
 sprintf(msgStr, "\t**** Evaluate scenario: %s %s ****", &SimAlt, &SimModif ); usrErr(msgStr); WriteMsg(msgStr,1);

 return;
}

Here is the call graph for this function:

void ReadStructFlow_PtSer

(

char *

filename

)

Read time series data of flow at water control structures.

v2.6 - made this a separate function instead of the old nested if - then statements

Definition at line 994 of file WatMgmt.c.

References HistData::aggCnt, HistData::aggID, HistData::arrayN, HistData::arrayP, HistData::arrayS, HistData::arrayWat, disAggNum, F_struct_wat, HistData::flag, Structure::flag, H_OPSYS, Hist_data, Structure::histID, Jdate_end, Jdate_init, julday(), modelFileName, ModelPath, msgStr, nalloc(), Structure::next_in_list, num_struct_hist, PORnumday, ProjName, Structure::S_nam, Scip(), SimAlt, struct_first, structs, TAB, UNIX, usrErr(), and usrErr0().

Referenced by ReadStructures().

{
     char lline[2001], *line;
     char scenario[20];
     int i, k; 
     int lincnt; 
     char structName[20]; 
     int yyyy, mm, dd;
     double Jdate_read;
     usrErr0 ( "\tControl structures' water flow time series...");
     if(H_OPSYS==UNIX) 
         sprintf( modelFileName, "%s/%s/Data/%s.struct_wat", ModelPath, ProjName, filename );
     else sprintf( modelFileName, "%s%s:Data:%s.struct_wat", ModelPath, ProjName, filename );
     

/* Open file with structure flow data */
     if ( ( F_struct_wat = fopen( modelFileName, "r" ) ) ==  NULL )
     {
         printf( "Can't open %s file!\n", modelFileName ) ;
         exit(-1) ;
     }

     fgets(lline, 2000, F_struct_wat); /* first header line with the alternative's name */
     line=lline;
     sscanf(line, "%s",&scenario);
     if (strcmp(scenario,SimAlt) != 0) {
         sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!", scenario, modelFileName, SimAlt); usrErr(msgStr);
         exit(-1);
     }

     fgets(lline, 2000, F_struct_wat); /* read 2'nd header line with column names */
     line=lline;
     line = Scip( line, TAB );
     line = Scip( line, TAB );
     line = Scip( line, TAB ); /*read past the three date columns, ready to read the variable names */


     /* loop through all the control structure columns */
     for ( i = 0; i < num_struct_hist; i++ ) { 

         Hist_data[i].arrayWat = (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
         Hist_data[i].arrayP =   (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
          Hist_data[i].arrayN =  (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 
          Hist_data[i].arrayS =  (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 

/*  allocate memory for dissaggregated structures and provide a couple of attributes */
         if ( Hist_data[i].flag > 1 ) {
             structs = struct_first;  /* go through all structures to pull out the disaggregated ones */
             while ( structs != NULL ) 
             {
                 /* Ex: the SFWMM aggregated S10's flag is 2, S10AB&C's flags are 20 */
                 if ( structs->flag == 10 * Hist_data[i].flag  )  
                 { 
                        /* disAggNum counts the total number of disaggregated structures for simulation */
                     disAggNum++;
                     Hist_data[num_struct_hist+disAggNum-1].arrayWat = 
                             (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
                     Hist_data[num_struct_hist+disAggNum-1].arrayP    = 
                             (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
                     Hist_data[num_struct_hist+disAggNum-1].arrayN    = 
                             (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 
                     Hist_data[num_struct_hist+disAggNum-1].arrayS    = 
                             (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 

                     Hist_data[num_struct_hist+disAggNum-1].aggID = i; /* the aggregated structure ID is attached to the dissagregated one */
                     strcpy( Hist_data[num_struct_hist+disAggNum-1].S_nam, structs->S_nam ); /* the disaggregated structure's name is put into the HistData structure */
                     Hist_data[i].aggCnt++; /* increment a count of the number of disagg structs associated w/ each agg struct */
                     structs->flag = 1; /* set the disagg struct's processing flag to 1, the only flag used for simulation structure flows (new input style) */
                     structs->histID = num_struct_hist+disAggNum; /* store the historical ID value for the disagg structures */
                     
                 } 
   
                 structs = structs->next_in_list;
             }             
         }
         
         sscanf ( line,"%s\t",&structName); /* read the structure names in the header */
         if (strcmp(Hist_data[i].S_nam , structName) != 0) {
             printf( "variable name  %s in input file does not match order of CanalData.structs.\n", structName);
             exit(-1); 
             }
         line = Scip( line, TAB );
     }


     lincnt = 0; /*  historical/other-data data input file record #  starting at the point where the model reads values */
     while ( fgets(lline, 2000, F_struct_wat) != NULL && !feof( F_struct_wat ) )
     {
         line = lline;
         sscanf( line, "%d %d %d",&yyyy,&mm,&dd); /* read the date of the current record */
         
        /* julian day, returns a double (args = mon, day, year, hours, minutes, seconds) */
         Jdate_read = julday( mm,dd,yyyy,0,0,0.0 ); 
         if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
             printf (" \n***Error\nStructure flow data (%s) init date > simulation start date\n", modelFileName); exit (-1);
         }
         if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
         {
             line = Scip( line, TAB );
             line = Scip( line, TAB );

             for ( i = 0; i < num_struct_hist; i++ ) { /* loop through all the control structure columns in the input file */
                 line = Scip( line, TAB );
                 if (*line == TAB) Hist_data[i].arrayWat[((int)(lincnt))] = 0;
                 else sscanf( line, "%f", &Hist_data[i].arrayWat[((int)(lincnt))] );
                 Hist_data[i].arrayWat[((int)(lincnt))] *= 2446.848; /* cfs => m^3/d (0.02832*60*60*24)  */
             }
         
         for (i = num_struct_hist; i < num_struct_hist+disAggNum; i++) { /* dissaggregate all agg structs */
             k = Hist_data[i].aggID; /* use the aggregated historical data WCstructure */
             Hist_data[i].arrayWat[((int)(lincnt))] =
                 Hist_data[k].arrayWat[((int)(lincnt))] / (Hist_data[k].aggCnt) ;

         }
         
              
         lincnt++; /* increment the number of daily data records */
     }

} /* end of reading structure flow records */
     
if (lincnt == 0) {
    printf (" \n***Error\nStructure flow data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
}
else if (lincnt != PORnumday ) {
    printf (" \n***Error\nStructure flow data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
}
     
sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
fclose (F_struct_wat);
} /* end of function for flow record input */

Here is the call graph for this function:

void ReadStructTP_PtSer

(

char *

filename

)

Read time series data of TP concentration at water control structures.

v2.6 - made this a separate function instead of the old nested if - then statements

Definition at line 1132 of file WatMgmt.c.

References HistData::arrayP, disAggNum, F_struct_TP, H_OPSYS, Hist_data, IDhist, Jdate_end, Jdate_init, julday(), modelFileName, ModelPath, msgStr, num_struct_hist, numTPser, PORnumday, ProjName, Structure::S_nam, Scip(), SimAlt, TAB, UNIX, usrErr(), and usrErr0().

Referenced by ReadStructures().

{
     char lline[2001], *line;
     char scenario[20];
     int i, j, k; 
     int lincnt; 
     /* int IDhist[MAX_H_STRUCT]; */ 
     char structName[20]; 
     int yyyy, mm, dd;
     double Jdate_read;

     usrErr0 ( "\tControl structures' TP concen. time series...");

     { if(H_OPSYS==UNIX) 
         sprintf( modelFileName, "%s/%s/Data/%s.struct_TP", ModelPath, ProjName, filename );
     else sprintf( modelFileName, "%s%s:Data:%s.struct_TP", ModelPath, ProjName, filename );
     }

/* Open file with structure concentration data */
     if ( ( F_struct_TP = fopen( modelFileName, "r" ) ) ==  NULL )
     {
         printf( "Can't open %s file!\n", modelFileName ) ;
         exit(-1) ;
     }

     fgets(lline, 2000, F_struct_TP); /* first header line with the alternative's name */
     line=lline;
     sscanf(line, "%s",&scenario);
     if (strcmp(scenario,SimAlt) != 0) {
         sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!", scenario, modelFileName, SimAlt); usrErr(msgStr);
         exit(-1);
     }
     fgets(lline, 2000, F_struct_TP); /* read 2'nd header line with column names */
     line=lline;
     line = Scip( line, TAB );
     line = Scip( line, TAB );
     line = Scip( line, TAB );

/* now we're past the three date columns, ready to read the variable names */
     for ( j = 0; j < numTPser; j++ ) { /* loop through all the data columns */
         IDhist[j] = -1;
         strcpy(structName,"NULL"); /* a string structName read previously is retained - this NULL assignment will take care of case
                                   where nothing was found/read in the time series file, and thus cause an error */
         sscanf ( line,"%s\t",&structName); /* read the structure names in the header */
         for ( i = 0; i < num_struct_hist+disAggNum; i++ ) { 
             if (strcmp(Hist_data[i].S_nam , structName) == 0) {
                 IDhist[j] = i;
                 line = Scip( line, TAB );
             }
         }
         if (IDhist[j] == -1 ) { printf( "variable name %s in %s not found in CanalData.structs. (if 'NULL', a variable in CanalData.structs defined as having a time series column was not found in %s)\n", structName, modelFileName, modelFileName);
         exit(-1); 
         }
     }

     lincnt = 0; /*  input file record #  starting at the point where the model reads values */

     while ( fgets(lline, 2000, F_struct_TP) != NULL && !feof( F_struct_TP ) )
     {
         line = lline;
         sscanf( line, "%d %d %d",&yyyy,&mm,&dd); /* read the date of the current record */
         
        /* julian day, returns a double (args = mon, day, year, hours, minutes, seconds) */
         Jdate_read = julday( mm,dd,yyyy,0,0,0.0 ); 
         if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
             printf (" \n***Error\nStructure TP data (%s) init date > simulation start date\n", modelFileName); exit (-1);
         }
         if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
         {
             line = Scip( line, TAB );
             line = Scip( line, TAB );
             for ( j = 0; j < numTPser; j++ ) { /* loop through all the control structure columns in the input file */
                 line = Scip( line, TAB );
                 if (*line == TAB) Hist_data[IDhist[j]].arrayP[((int)(lincnt))] = 0;
                 else sscanf( line, "%f", &Hist_data[IDhist[j]].arrayP[((int)(lincnt))] );
                 Hist_data[IDhist[j]].arrayP[((int)(lincnt))] *= 1.0E-6; /* data read in ug/L, convert to g/L=kg/m3 */
                 
             }
              
             lincnt++; /* increment the number of daily data records */
         }

     } /* end of reading TP conc data records */
     
     if (lincnt == 0 ) {
             printf (" \n***Error\nTP Time series data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
     }
     else if (lincnt != PORnumday ) {
             printf (" \n***Error\nTP Time series data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
     }
     sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
     
     fclose (F_struct_TP);
 } /* end of reading TP time series data (if needed) */

Here is the call graph for this function:

void ReadStructTS_PtSer

(

char *

filename

)

Read time series data of TS (salt/tracer) concentration at water control structures.

v2.6 - made this a separate function instead of the old nested if - then statements

Definition at line 1234 of file WatMgmt.c.

References HistData::arrayS, disAggNum, F_struct_TS, H_OPSYS, Hist_data, IDhist, Jdate_end, Jdate_init, julday(), modelFileName, ModelPath, msgStr, num_struct_hist, numTSser, PORnumday, ProjName, Structure::S_nam, Scip(), SimAlt, TAB, UNIX, usrErr(), and usrErr0().

Referenced by ReadStructures().

{
 char lline[2001], *line;
 char scenario[20];
 int i, j, k; 
 int lincnt; 
 /* int IDhist[MAX_H_STRUCT];  */
 char structName[20]; 
 int yyyy, mm, dd;
 double Jdate_read;

     usrErr0 ( "\tControl structures' Tracer concen. time series...");

     { if(H_OPSYS==UNIX) 
         sprintf( modelFileName, "%s/%s/Data/%s.struct_TS", ModelPath, ProjName, filename );
     else sprintf( modelFileName, "%s%s:Data:%s.struct_TS", ModelPath, ProjName, filename );
     }

/* Open file with structure concentration data */
     if ( ( F_struct_TS = fopen( modelFileName, "r" ) ) ==  NULL )
     {
         printf( "Can't open %s file!\n", modelFileName ) ;
         exit(-1) ;
     }

     fgets(lline, 2000, F_struct_TS); /* first header line with the alternative's name */
     line=lline;
     sscanf(line, "%s",&scenario);
     if (strcmp(scenario,SimAlt) != 0) {
         sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!", scenario, modelFileName, SimAlt); usrErr(msgStr);
         exit(-1);
     }
     fgets(lline, 2000, F_struct_TS); /* read 2'nd header line with column names */
     line=lline;
     line = Scip( line, TAB );
     line = Scip( line, TAB );
     line = Scip( line, TAB );

/* now we're past the three date columns, ready to read the variable names */
     for ( j = 0; j < numTSser; j++ ) { /* loop through all the data columns */
         IDhist[j] = -1;
         strcpy(structName,"NULL"); /* a string structName read previously is retained - this NULL assignment will take care of case
                                   where nothing was found/read in the time series file, and thus cause an error */
         sscanf ( line,"%s\t",&structName); /* read the structure names in the header */
         for ( i = 0; i < num_struct_hist+disAggNum; i++ ) { 
             if (strcmp(Hist_data[i].S_nam , structName) == 0) {
                 IDhist[j] = i;
                 line = Scip( line, TAB );
             }
         }
         if (IDhist[j] == -1 ) { printf( "variable name %s in %s not found in CanalData.structs. (if 'NULL', a variable in CanalData.structs defined as having a time series column was not found in %s)\n", structName, modelFileName, modelFileName);
         exit(-1); 
         }
     }

     lincnt = 0; /*  input file record #  starting at the point where the model reads values */

     while ( fgets(lline, 2000, F_struct_TS) != NULL && !feof( F_struct_TS ) )
     {
         line = lline;
         sscanf( line, "%d %d %d",&yyyy,&mm,&dd); /* read the date of the current record */
         
        /* julian day, returns a double (args = mon, day, year, hours, minutes, seconds) */
         Jdate_read = julday( mm,dd,yyyy,0,0,0.0 ); 
         if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
             printf (" \n***Error\nStructure Tracer data (%s) init date > simulation start date\n", modelFileName); exit (-1);
         }
         if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
         {
             line = Scip( line, TAB );
             line = Scip( line, TAB );
             for ( j = 0; j < numTSser; j++ ) { /* loop through all the control structure columns in the input file */
                 line = Scip( line, TAB );
                 if (*line == TAB) Hist_data[IDhist[j]].arrayS[((int)(lincnt))] = 0;
                 else sscanf( line, "%f", &Hist_data[IDhist[j]].arrayS[((int)(lincnt))] );
                 Hist_data[IDhist[j]].arrayS[((int)(lincnt))] *= 1.0; /* data read in g/L=kg/m3, no conversion*/
                 
             }
              
             lincnt++; /* increment the number of daily data records */
         }

     } /* end of reading tracer (salt) conc data records */
     
     if (lincnt == 0 ) {
             printf (" \n***Error\nTracer Time series data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
     }
     else if (lincnt != PORnumday ) {
             printf (" \n***Error\nTracer Time series data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
     }
     sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
     
     fclose (F_struct_TS);
 } /* end of function for reading tracer (salt) time series data (if needed) */

Here is the call graph for this function:

struct Schedule* Read_schedule	(	char *	sched_name,
		char *	filename,
		float	BASE_DATUM
	)			[read]

Read target stage schedule as a recurring time-series graph.

Parameters:

	sched_name	The name of the target stage location
	filename	Input file name for graph
	BASE_DATUM	GP_DATUM_DISTANCE

Returns:

Sce_Graph struct of the time series schedule graph

Definition at line 1342 of file WatMgmt.c.

References EOL, EOS, Schedule::graph_points, H_OPSYS, modelFileName, modelName, ModelPath, msgStr, Schedule::num_point, NumScheds, ProjName, schedFile, Schedule::schedName, Scip(), SimAlt, Points::time, True, UNIX, usrErr(), Points::value, Wrdcmp(), and WriteMsg().

Referenced by ReadStructures().

{ 
char ss[601], *s, scenario[30], modnam[30];
int i, count = 1, sched_found=0;
struct Schedule *Sce_Graph;

   { if(H_OPSYS==UNIX) 
                sprintf( modelFileName, "%s/%s/Data/%s.graph", ModelPath, ProjName, filename );
     else sprintf( modelFileName, "%s%s:Data:%s.graph", ModelPath, ProjName, filename );
   }

/* Open file with schedule graphs */
  if ( ( schedFile = fopen( modelFileName, "r" ) ) ==  NULL )
        {
       printf( "Can't open %s file!\n", modelFileName ) ;
       exit(-1) ;
    }

 fgets( ss, 600, schedFile );  /* skip comment line */
 fgets( ss, 600, schedFile );  /* skip comment line */
 fgets( ss, 600, schedFile );  /* skip comment line */
 fgets( ss, 600, schedFile ); sscanf( ss,"%s %s", &modnam, &scenario );
 if (strcmp(scenario,SimAlt) != 0) {
     sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
             scenario, modelFileName, SimAlt); usrErr(msgStr);
             exit(-1);
 }
 if (strcmp(modnam,modelName) != 0) {
     sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
             modnam, modelFileName, modelName); usrErr(msgStr);
             exit(-1);
 }

/* Read schedule descriptors until none left */   
  while ( fgets( ss, 600, schedFile ) != NULL && !feof( schedFile ) )
        { if ( !Wrdcmp ( ss, sched_name ) ) continue; /* v2.8 - bugfix - reversed arg order in this version: Wrdcmp looks to see if the next char is a space in the first arg; the ss, sched_name is correct order for this to work (old code allowed sched_names that started with the same string to be treated as the same) */
          sched_found = True;
          
          fgets( ss, 600, schedFile );  
          s = ss;
          /* count the number of pairs for the graph */
          while ( *s++ != EOL ) if ( *s == '(' ) count++;    

          /* Allocate memory for next schedule */
      if ( (Sce_Graph = (struct Schedule *) malloc( (size_t) sizeof( struct Schedule ))) == NULL )
       {
         printf( "Failed to allocate memory for Managed flow Schedule graph in (%s)\n ", sched_name ) ;
         exit( -2 ) ;
       };
     
          /* Allocate memory for next schedule graph containing 'count' points*/
      if ( (Sce_Graph->graph_points = 
                (struct Points *) malloc( (size_t) sizeof( struct Points ) * count)) == NULL )
       {
         printf( "Failed to allocate memory for Managed flow Schedule graph in (%s)\n ", sched_name ) ;
         exit( -2 ) ;
       };
          
          Sce_Graph->num_point = count;
         
          s = ss;
          for ( i = 0; i < count && *s!=EOS; i++ )
          { 
                s = Scip( s,'(');       
                sscanf( s, "%f", &Sce_Graph->graph_points[i].time ); 
                s = Scip ( s, ',');
                sscanf( s, "%f", &Sce_Graph->graph_points[i].value );
                Sce_Graph->graph_points[i].value += BASE_DATUM;
          }
        } /* end of while */
        /* v2.5 included this error-trap */
        if (!sched_found) {
             sprintf(msgStr, "The schedule name (%s) was not found in the %s schedule file! (Was named in the CanalData.struct file.)",
             sched_name, modelFileName); usrErr(msgStr);
             exit(-1);
    }

         /* v2.8 writing this to Driver1.out debug file */
         sprintf(msgStr, "Read target schedule graph: %s, which contained %d time points.",strcpy(Sce_Graph->schedName, sched_name), Sce_Graph->num_point);
     WriteMsg(msgStr,1);  
     NumScheds = NumScheds+1; /* used just for writing info to console on number of schedules you've read */
     
        
  fclose (schedFile);
  return Sce_Graph;
}

Here is the call graph for this function:

void Channel_configure	(	float *	ElevMap,
		struct Chan *	channel_first
	)

Configure attributes of grid cells interacting with canal vectors.

Identifies the cells that are crossed by the canal, calculates the length of canal within each cell and marks the cells directly interacting with the canal. Store the pointers to the adjacent cell arrays.

Parameters:

	ElevMap	SED_ELEV
	channel_first	Pointer to first canal data structure

Definition at line 1441 of file WatMgmt.c.

References Abs, ALL, Chan::area, C_F, C_Mark, CanalCellInteractDebug, cell, cell_last, Chan::cells, celWid, Chan_list, CHo, Chan::cond, EAST, Chan::elev_drop, Chan::elev_end, Chan::elev_slope, Chan::elev_start, getCanalElev(), GP_calibGWat, init_pvar(), LEFT, Chan::length, Chan::levee, MarkCell(), MINDEPTH, Chan::minVol, modelFileName, nalloc(), Cells::next_cell, Chan::next_in_list, Chan_segm::next_segm, NONE, num_cell, num_chan, Chan::num_of_cells, Chan::number, ON_MAP, OutputPath, ProjName, RIGHT, Chan::roil, Chan::roir, Round(), s0, s1, sec_per_day, Chan::seg_len, Chan::segments, SOUTH, Chan::SPG_flow_coef, STEP, Chan::SW_flow_coef, T, Chan::width, Chan_segm::x0, Chan_segm::x1, Chan_segm::y0, and Chan_segm::y1.

Referenced by Canal_Network_Init().

{ int i, j;                                                       /*coord. of current cell crossed by the canal */
 int  cellLoc, cellLoc_i, cellLoc_j;
 int HV, k, L_range, R_range;
 float T_length;                                                /* length of a straight segment of a canal reach (m) */
 float C_length;                                          /*  total length of all of the segments of a canal reach (m) */
 float distTot;                                       /*  cumulative (along grid cells) distance along a canal reach, from the starting point (m) */
 float A1, B1, A2, B2;                                                   /* parameters of the canal equation */
 float L, p_middle, m, i45, x, y, x1, y1, length, RoIL, RoIR;
 float cell_step = 1.0;
   
 float a0, b0, a1, b1;
 int c_num, z;
 struct Cells  *cell; /* addresses for adjacent cell structures */
 struct Chan *channel;
 struct Chan_segm *Segm;
 
 int numChek; /* placeholder for debugging cell elev along canals */
   

 int *marked; /* a "temporary" fix used to determine if a cell belonging to a canal has already been marked */
 marked = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"marked");
 init_pvar(marked,NULL,'d',0);

/* Open file to print canal-cell interaction info to */
 sprintf( modelFileName, "%s/%s/Output/Debug/CanalCells_interaction.txt", OutputPath, ProjName );
 if ( ( CanalCellInteractDebug = fopen( modelFileName, "w" ) ) ==  NULL )
 {
     printf( "Can't open %s file!\n", modelFileName ) ;
     exit(-1) ;
 }

     /* allocate memory for array of pointers to Chan structs */
 if ( (Chan_list = 
       (struct Chan **) malloc( (size_t) sizeof( struct Chan *) * num_chan)) == NULL )
 {
     printf( "Failed to allocate memory for next canal (%d)\n ", channel->number ) ;
     exit( -2 ) ;
 };

     /* allocate memory for first cell structure */
 if ( (cell = (struct Cells *) malloc( (size_t) sizeof( struct Cells ))) == NULL )
 { 
     printf( "Failed to allocate memory for cell structure\n " ) ;
     exit( -2 ) ;
 }

     /* store the pointer to first cell in the array */
 channel = channel_first;

     /* loop over all the canal reaches (== channel, defined as an individual canal reach with an ID number) */
 while (channel != NULL)
 {
     Segm = channel->segments;
     c_num = channel->number;
     z  = channel->levee;
     RoIL = channel->roil;
     RoIR = channel->roir;
        
     Chan_list[c_num - CHo] = channel;
        
     channel->cells = cell;
     distTot = 0;
     C_length = 0;
     C_Mark = 0;
    
         /* a negative-width canal is skipped, going to SKIPCHAN line that is located (almost) at the end of this function */
     if (channel->width < 0)   goto SKIPCHAN; 
    
     fprintf ( CanalCellInteractDebug, "N = %d  levee = %d\n", c_num, z );
     
/* get the elevation at the starting reach segment (may be up or down stream) of the canal */
     cellLoc_i=(int)(Segm->x0+1);
     cellLoc_j=(int)(Segm->y0+1);
     cellLoc = T(cellLoc_i,cellLoc_j );
     Chan_list[c_num - CHo]->elev_start = ElevMap[cellLoc]; 
     if ( !ON_MAP[cellLoc] ) {
        printf( "Starting elevation of Reach %d is out of the active domain.  Fix the reach location.\n",c_num);
        exit(-2);
     }
     fprintf ( CanalCellInteractDebug, " %d, %d, StartElev = %f\n", cellLoc_i,cellLoc_j,  Chan_list[c_num - CHo]->elev_start );

                
     while (Segm != NULL)  /* loop through canal reach segments */
     {
         a0 = Segm->x0;
         b0 = Segm->y0;
         a1 = Segm->x1;
         b1 = Segm->y1;

         x1 = 0;  y1 = 0;

         T_length = sqrt( (b1 - b0)*(b1 - b0) + (a1 - a0)*(a1 - a0) );  /* total reach segment length */
         C_length += T_length;                                                                              /* total canal reach length */

         if ( a1 != a0 && b1 != b0 ) 
         {
             A1 = (b1 - b0)/(a1 - a0);                      /* define the equation: y = A1x + B1 */
             B1 = b0 - A1*a0;
             A2 = 1/A1;  B2 = - B1/A1;                      /* define the equation: x = A2y + B2 */
         } 
    
         x = a0; y = b0;
         i = Round (x);    j = Round (y);   
              
             /* loop along cells of the canal reach segment */ 
         while ( x1 != a1 || y1 != b1 )
         {
             if ( a1 == a0 )                                     /* if goes straight horizontally */
             {  x1 = a0; 
             y1 = (j == Round (b1) ) ? b1 : j + STEP(b1 - b0)*cell_step; 
             length = y1 - y;
             }
       
             else if ( b1 == b0 )                                   /* if goes straight vertically */
             { x1 = (i == Round (a1) ) ? a1 : i + STEP(a1 - a0)*cell_step; 
             y1 = b0; 
             length = x1 - x;
             }

             else if ( i == Round (a1) && j == Round (b1) )
             { x1 = a1;  y1 = b1;                                               /* if last cell */ 
             length = T_length * (y1 - y)/(b1 - b0);
             C_Mark = 0;
             }
          
             else
             { y1 = A1*(i + STEP(a1 - a0)*cell_step) + B1;  length = T_length * (y1 - y)/(b1 - b0);
             x1 = A2*(j + STEP(b1 - b0)*cell_step) + B2;  L = T_length * (x1 - x)/(a1 - a0);
    
             if ( length < L )                       /* define the end coordinates and length */
                 x1 = i + STEP(a1 - a0)*cell_step;  /* (the next coordinate is the closest)*/
             else
             {  y1 = j + STEP(b1 - b0)*cell_step;  length = L;
             }
             }

             if ( Abs (a1 - a0) < Abs (b1 - b0) )
                     /* define the direction of canal reach and  */
             { p_middle = (x1 + x)/2;      /* calculate the middle point on the canal reach in cell */
             m = i + cell_step/2;                                                               /* m - the center of the cell */
             HV = 1;  /* to indicate that the canal reach goes closer to the horizontal direction */
             }
             else
             { p_middle = (y1 + y)/2;
             m = j + cell_step/2;
             HV = 0;
             }
             
                 distTot = distTot + length*celWid;  /*  cumulative (along grid cells) distance along a canal reach, from the starting point (m) */
         
                 /* mark cells depending on where the center of canal reach is rel. to the cell center*/

             if ( Round (a1) == Round (a0) ) /* if the canal reach goes within the range of one */
                     /*cell horizontally, mark the two vertical cells */
             { if (b0 < b1)     /* if going from left to right */
             {          /* if the middle point of the canal reach is lower than the */
                     /* cell center, mark the cell as LEFT else as RIGHT */
                 if ( p_middle > m ) 
                 { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i, j, c_num, marked, distTot );
                 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
                 }  
                 else 
                 { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i - 1, j, c_num, marked, distTot );
                 L_range = (int) RoIL;  R_range = (int)(RoIR - 1);
                 }  
                 for ( k = 0; k < L_range; k++ )
                     cell = MarkCell ( i - k - 1, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
                 for ( k = 0; k < R_range; k++ )
                     cell = MarkCell ( i + k + 1, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             else  /* vice versa if we go from right to left */
             {          /* if the middle point of the canal reach is lower than the */
                     /* cell center, mark the cell as RIGHT else as LEFT */
                 if ( p_middle > m ) 
                 { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i, j, c_num, marked, distTot );
                 L_range = (int)RoIL; R_range = (int)(RoIR - 1);
                 } 
                 else 
                 { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i - 1, j, c_num, marked, distTot );
                 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
                 }
                 for ( k = 0; k < L_range; k++ )
                     cell = MarkCell ( i + k + 1, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
                 for ( k = 0; k < R_range; k++ )
                     cell = MarkCell ( i - k - 1, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             }
             else if ( Round (b1) == Round (b0) )       /* if the canal reach goes vertically */
             { if (a0 < a1) /* if going from top to bottom */
             {      /* if the middle point of the canal reach is further to the right than the */
                     /* cell center, mark the cell as RIGHT else as LEFT */
                 if ( p_middle > m ) 
                 { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
                 L_range = (int)RoIL; R_range = (int)(RoIR -1);
                 }
                 else 
                 { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j - 1, c_num, marked, distTot );
                 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
                 }
                 for ( k = 0; k < L_range; k++ )
                     cell = MarkCell ( i, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
                 for ( k = 0; k < R_range; k++ )
                     cell = MarkCell ( i, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             else  /* vice versa if we go from bottom to top */
             {      /* if the middle point of the canal reach is further to the left than the */
                     /* cell center, mark the cell as LEFT else as RIGHT */
                 if ( p_middle > m ) 
                 { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
                 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
                 }
                 else 
                 { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j - 1, c_num, marked, distTot );
                 L_range = (int)RoIL; R_range = (int)(RoIR - 1);
                 }
                 for ( k = 0; k < L_range; k++ )
                     cell = MarkCell ( i, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
                 for ( k = 0; k < R_range; k++ )
                     cell = MarkCell ( i, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             }
             else
                     /* if the canal reach crosses more than one cell (in both dimensions), mark diagonal cells */
             { if ( A1>0 )                 /* this is an extremely confusing check, but in this */
                 i45 = ( HV ) ? p_middle - m : m - p_middle;  /*form it seems to be working OK */
             else 
                 i45 = p_middle - m;
                 /*  b0>b1, a0>a1 \/ b0<b1, a0>a1 */
                 /*  b0>b1, a0<a1 /\ b0<b1, a0<a1 */
             if ( b0 < b1 )
             { if ( a0 < a1 )     /* when a0<a1 the canal reach goes down in the \ direction */
             { if ( i45 > 0 ) 
             { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i + 1, j - 1, RIGHT, length, cell, SOUTH, z, i + 1, j - 1, c_num, marked, distTot );
             }
             else 
             { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i - 1, j + 1, LEFT, length, cell, EAST, z, i - 1, j + 1, c_num, marked, distTot );
             }
             for ( k = 0; k < RoIR-1; k++ )
                 cell = MarkCell ( i + k + 1, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             for ( k = 0; k < RoIL-1; k++ )
                 cell = MarkCell ( i - k - 1, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             else                                 /* in this case the canal reach goes up in the / direction */
             { if ( i45 > 0 ) 
             { cell = MarkCell ( i, j, LEFT, length, cell, NONE, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i + 1, j + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             else 
             { cell = MarkCell ( i, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i - 1, j - 1, LEFT, length, cell, NONE, z, i - 1, j - 1, c_num, marked, distTot );
             }
             for ( k = 0; k < RoIR-1; k++ )
                 cell = MarkCell ( i + k + 1, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             for ( k = 0; k < RoIL-1; k++ )
                 cell = MarkCell ( i - k - 1, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             }
                 /* if going in the opposite direction change left and right */            
             else
             { if ( a0 > a1 )       /* when a0<a1 the canal reach goes up in the \ direction */
             { if ( i45 < 0 ) 
             { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i - 1, j + 1, RIGHT, length, cell, EAST, z, i - 1, j + 1, c_num, marked, distTot );
             }
             else 
             { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i + 1, j - 1, LEFT, length, cell, SOUTH, z, i + 1, j - 1, c_num, marked, distTot );
             }

             for ( k = 0; k < RoIL-1; k++ )
                 cell = MarkCell ( i + k + 1, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             for ( k = 0; k < RoIR-1; k++ )
                 cell = MarkCell ( i - k - 1, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             else                               /* in this case the canal reach goes down in the / direction */
             { if ( i45 < 0 ) 
             { cell = MarkCell ( i, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i - 1, j - 1, RIGHT, length, cell, NONE, z, i - 1, j - 1, c_num, marked, distTot );
             }
             else 
             { cell = MarkCell ( i, j, RIGHT, length, cell, NONE, z, i, j, c_num, marked, distTot );
             cell = MarkCell ( i + 1, j + 1, LEFT, length, cell, ALL, z, i + 1, j + 1, c_num, marked, distTot );
             }

             for ( k = 0; k < (int)(RoIL-1); k++ )
                 cell = MarkCell ( i + k + 1, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             for ( k = 0; k < (int)(RoIR-1); k++ )
                 cell = MarkCell ( i - k - 1, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
             }
             }
             }
         
             num_cell++;        /*at this point we're counting the cells only once thinking that all 
                                  cells have the same length of interaction irrelevant of their distance from the canal 
                                  and that this length will be equal to the canal length / number of cells it crosses */

                 /* prepare for the next step in cycle */
             x = x1;  y = y1;
             i = Round (x);     if ( a1 < a0 && i == x ) i--;         /* if we're going upwards, we */
                 /* need to change the rule of Rounding so that we move ahead along cells */
             j = Round (y);     if ( b1 < b0 && j == y ) j--;
             C_Mark = 1;
    

         } /* end loop along cells of the canal reach segment  */
         
         Segm = Segm->next_segm;
     } /* end loop over all reach segments of a canal */

/* get the elevation at the ending reach segment (may be up or down stream) of the canal */
     cellLoc_i=(int)(a1+1);
     cellLoc_j=(int)(b1+1);
     cellLoc = T(cellLoc_i,cellLoc_j );
     Chan_list[c_num - CHo]->elev_end = ElevMap[cellLoc];  

     Chan_list[c_num - CHo]->num_of_cells = num_cell;
     num_cell = 0;
         /* store the canal attributes of:
            length of entire canal, 
            area of entire canal, 
            minimum volume allowed in canal,
            avg length of cell along reach segment,
            overland flow coefficient (not incl. dynamic manning's n, C_F is for sensitivity experiments only (normally=1) ),
            seepage flow coefficient (including GP_calibGWat groundwater flow calibration coef),
            land surface elevation difference between starting and ending points of canal 
            slope of the elevation gradient from start to end points */
     Chan_list[c_num - CHo]->length = C_length*celWid;
     Chan_list[c_num - CHo]->area = Chan_list[c_num - CHo]->width*Chan_list[c_num - CHo]->length;
     Chan_list[c_num - CHo]->minVol = Chan_list[c_num - CHo]->area*MINDEPTH;
     Chan_list[c_num - CHo]->seg_len = Chan_list[c_num - CHo]->length/Chan_list[c_num - CHo]->num_of_cells;
     Chan_list[c_num - CHo]->SW_flow_coef = sqrt(Chan_list[c_num - CHo]->seg_len) * sec_per_day * C_F;
     Chan_list[c_num - CHo]->SPG_flow_coef = Chan_list[c_num - CHo]->cond * GP_calibGWat;
     Chan_list[c_num - CHo]->elev_drop = Chan_list[c_num - CHo]->elev_start - Chan_list[c_num - CHo]->elev_end;
     Chan_list[c_num - CHo]->elev_slope = Chan_list[c_num - CHo]->elev_drop / Chan_list[c_num - CHo]->length ; /* v2.5 calculated slope */
     
     if ( !ON_MAP[cellLoc] ) {
        printf( "Ending elevation of Reach %d is out of the active domain.  Fix the reach location.\n",c_num);
        exit(-2);
     }
     fprintf ( CanalCellInteractDebug, " %d, %d, EndingElev = %f, ElevSlope= %f, ReachLength= %f\n", cellLoc_i,cellLoc_j,  
                                        Chan_list[c_num - CHo]->elev_end, Chan_list[c_num - CHo]->elev_slope, Chan_list[c_num - CHo]->length );
     
     getCanalElev(c_num - CHo); /* v2.5 */
     
   SKIPCHAN: /* a single goto comes here, skipping a canal with negative width */
     channel = channel->next_in_list;

     cell_last->next_cell = NULL;                       /* close the cell list */

 } /* end loop loop over all the canal reaches */
 
 fclose ( CanalCellInteractDebug );
 
 return; 
}

Here is the call graph for this function:

void getCanalElev

(

int

chan_n

)

Assign elevation to canal using slope & distance from start point
.

This function does two things:
1. For each grid cell along a canal, we determine the elevation (relative to model base datum) of the top of the canal bank, using the slope of the canal and the distance from the beginning point of the canal. (function new in v2.5)
2. Calculate the mean elevation for a reach.

Parameters:

chan_n

Canal reach number (adjusted for the beginning ID number)

Returns:

void

Definition at line 1812 of file WatMgmt.c.

References CanalCellInteractDebug, cell, Chan::cells, Chan_list, CHo, Chan::elev_avg, Chan::elev_slope, Chan::elev_start, Exit(), Cells::next_cell, Cells::reachDistDown, Cells::reachElev, T, Cells::x, and Cells::y.

Referenced by Channel_configure().

{
     int i, celcount=0;
     Chan_list[chan_n]->elev_avg = 0.0 ; /* v2.8 calculating an avg elevation for reach */
         /* cycle along canal cells */
     cell = Chan_list[chan_n]->cells;
     i = T(cell->x, cell->y);
     
     while ( cell != NULL )
     { 
        celcount++;
        cell->reachElev = Chan_list[chan_n]->elev_start - cell->reachDistDown * Chan_list[chan_n]->elev_slope;
        Chan_list[chan_n]->elev_avg = Chan_list[chan_n]->elev_avg + cell->reachElev;
        fprintf ( CanalCellInteractDebug, " Reach= %d, Cell (%d, %d), ReachElev= %f\n", chan_n+CHo, cell->x, cell->y, cell->reachElev );
        cell = cell->next_cell;
     }
     if (celcount == 0) { printf(" Error in getting cells along canal for elevation avg\n"); Exit(-1); }
     Chan_list[chan_n]->elev_avg = Chan_list[chan_n]->elev_avg/celcount;
}

Here is the call graph for this function:

struct Cells* MarkCell	(	int	x,
		int	y,
		int	c_ind,
		float	length,
		struct Cells *	cell,
		int	direction,
		int	levee,
		int	xm,
		int	ym,
		int	c_num,
		int *	marked,
		float	distTot
	)			[read]

Mark cell function.

This function does two things:
1. It stores the cell coordinates (x, y) in the interaction array, also identifying the location of this cell with respect to the canal;
2. It marks the cell (xm,ym) on the ON_MAP array to prevent horizontal overland fluxes if there is a levee associated with this canal.

Parameters:

	x	Grid cell row
	y	Grid cell column
	c_ind	Cell index = 1 for left cells and -1 for right cells
	length	Length of canal associated with this cell
	cell	A struct of Cells
	direction	Direction of interaction: 1 - none, 2 -to the East, 3 - to the South, 4 - all directions
	levee	Levee location attribute
	xm	Grid cell row marked
	ym	Grid cell column marked
	c_num	Canal ID number
	marked	Denote a marked cell
	distTot	Cumulative (along grid cells) distance along a canal reach, from the starting point (m)

Returns:

cell Pointer to struct of Cells

Definition at line 1858 of file WatMgmt.c.

References ALL, C_Mark, CanalCellInteractDebug, cell_last, celWid, Cells::ind, Cells::length, Cells::next_cell, ON_MAP, Cells::reachDistDown, s0, s1, T, Cells::x, and Cells::y.

Referenced by Channel_configure().

{ struct Cells *cell_next;
  int N, cellLoc;
  
   if ( x > s0 || y > s1 ) /* v2.6 changed evaluation to ">" (was ">=") */
     {  printf ( "Error in definition of canal in cell row=%d  col=%d ", x, y );
        exit (-2);                      /*check that cells obtained are within array boundaries */
     }

   N = T(xm+1, ym+1);
   if ( !ON_MAP[N] ) return cell;
   
   if ( levee && direction != ALL )
   ON_MAP[N] = ((ON_MAP[N]-1) & (direction+100)) + 1;   /* modifying ON_MAP to take into account the 
   allowed direction of fluxing: 00 - none, 01=1 - to the East, 10=2 - to the South, 11=3 - all 
   ON_MAP they naturally become 1 - none, 2 -to the East, 3 - to the South, 4 - all directions 
   0 is reserved for cells out of the ON_MAP area*/
   
   if ( !C_Mark ) return cell; /* C_Mark is to make gaps between canals, so that there will be no 
                canals being connected through cells if it turns out that they are linked to the same cells */
             
       /* a temporary fix used to determine if a cell belonging to a canal has already been marked */
       /***********/
   cellLoc = T(x+1,y+1);
   if (marked[cellLoc] == c_num ) return cell;
   marked[cellLoc] = c_num;
       /***********/
             

   
   cell->x = x+1;                                               /* shift coord to match the other maps */
   cell->y = y+1;
   cell->ind = c_ind;   
   cell->length = length*celWid;                /* length of canal associated with this cell * cell width m */
   cell->reachDistDown = distTot; /*  cumulative (by grid cells) distance along a canal reach, from the starting point (m) (v2.5) */
   
   fprintf ( CanalCellInteractDebug, " %d, %d, l/r = %d, distance= %7.1f\n", cell->x, cell->y, c_ind, cell->reachDistDown );
   
   /* allocate memory for next cell structure */
   if ( (cell_next = (struct Cells *) malloc((size_t) sizeof( struct Cells ))) == NULL )
       {
         printf( "Failed to allocate memory for cell structure\n " ) ;
         exit( -2 ) ;
       }
   cell->next_cell = cell_next;
   cell_last = cell;

   return cell_next;
}

void FluxChannel	(	int	chan_n,
		float *	SWH,
		float *	ElevMap,
		float *	MC,
		float *	GWH,
		float *	poros,
		float *	GWcond,
		double *	NA,
		double *	PA,
		double *	SA,
		double *	GNA,
		double *	GPA,
		double *	GSA,
		float *	Unsat,
		float *	sp_yield
	)

Runs the iterative algorithm over the network of canals.

The iterative relaxation routine calculates the new head in the canal and the water exchange with the adjacent cells, including constituent fluxes.

Parameters:

	chan_n	Canal ID number
	SWH	SURFACE_WAT
	ElevMap	SED_ELEV
	MC	HYD_MANNINGS_N
	GWH	SAT_WATER
	poros	HP_HYD_POROSITY
	GWcond	HYD_RCCONDUCT
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT
	GNA	DINdummy
	GPA	TP_SED_WT
	GSA	SALT_SED_WT
	Unsat	UNSAT_WATER
	sp_yield	HP_HYD_SPEC_YIELD

Variables local to function

converged The convergence flag
error The error calculated as the diff between the incoming volume and the absorbed volume
prev_error The error calculated in the previous iteration
iter The iteration number in the iterative relaxation procedure
max_iter The max number of iterations allowed
init_factor The depth increment for the initial guess (now same as "factor")
factor The depth increment for the new guess (now same as "init_factor")
attempt The counter of number of attempts at convergence
CanWatDep The water level height in canal
T_flux_S, T_flux_G, T_flux_L The total net fluxes among canal and cells (Surface, Groundwater, Levee seepage)
fluxS, fluxG, fluxL The flux among canal and cells (Surface, Groundwater, Levee seepage)
fluxO, fluxCk Flow-checks to ensure the canal is not drained below its minimum level
SW_coef, SPG_coef, GW_coef Aggregated flow coefficients (Surface, Groundwater, Levee seepage)
N_fl, P_fl, S_fl Flux of Nitrogen (UNUSED), Phosphorus, and Salt/tracer among canal and cells
I_Length Average length of cell along reach segment
seg_area Average area of reach segment
MIN_VOL Minimum canal volume allowed
GW_head Groundwater head
h_GWflow, h_SPflow Water heights for groundwater and seepage flows
dh Difference between hydraulic heads
Qin, Qout, Qout1 Flow volumes
SWater SURFACE_WAT
tot_head Grid cell water head
CH_vol Total volume of canal at the estimated water depth
CH_vol_R Total canal volume prior to relaxation procedure
CH_bottElev Elevation of bottom of canal at cell location
H_rad_ch, H_rad_cell Hydraulic radii, canal (channel) and cell
can_cell_M, can_cell_A Temp vars used in reducing excess flows between canal<->cell
errorConv Numerical error after convergence

Definition at line 1935 of file WatMgmt.c.

References Abs, Chan::area, Chan::basin, basn, basn_list, Structure::canal_fr, Structure::canal_to, canDebugFile, cell, CELL_SIZE, Chan::cells, Chan_list, CHANNEL_MAX_ITER, CHo, debug, Chan::depth, dynERRORnum, Chan::edgeMann, F_ERROR, f_Ground(), f_Manning(), basndef::family, Chan::family, Structure::flow, GP_DetentZ, GP_MinCheck, HAB, Cells::ind, Max, Min, MINDEPTH, Chan::minVol, msgStr, Chan::N, Cells::next_cell, Structure::next_in_list, Chan::number, ON_MAP, Chan::P, P, Chan::P_con, P_IN_GW, P_IN_OVL, P_IN_SPG, P_OUT_GW, P_OUT_OVL, P_OUT_SPG, ramp, Cells::reachElev, Chan::S, Chan::S_con, S_IN_GW, S_IN_OVL, S_IN_SPG, S_OUT_GW, S_OUT_OVL, S_OUT_SPG, Chan::seg_len, sgn, SimTime, Chan::SPG_flow_coef, struct_first, structs, T, simTime::TIME, True, VOL_IN_GW, VOL_IN_OVL, VOL_IN_SPG, VOL_OUT_GW, VOL_OUT_OVL, VOL_OUT_SPG, Chan::wat_depth, Chan::width, WriteMsg(), Cells::x, and Cells::y.

Referenced by Run_Canal_Network().

{ 
                
int converged = 0;                                              
 float error = 1.0;                                     
 float prev_error = 1.0;                                        
 int iter = 0;                                                  
 int max_iter = CHANNEL_MAX_ITER;               
 float init_factor = 0.5;                               
 float factor = init_factor;                    

 int i;
 int attempt=0;                                                 
 double CanWatDep = Chan_list[chan_n]->wat_depth;
 double T_flux_S, T_flux_G, T_flux_L;           
 double fluxO, fluxCk;
 double fluxS, fluxG, fluxL; 
 double SW_coef, SPG_coef, GW_coef;     
 double N_fl, P_fl, S_fl;                               
 /* double SNfl, GNfl, SPfl, GPfl, SSfl, GSfl; UNUSED */
 double I_Length;                                               
 double seg_area;                                               
 double MIN_VOL;                                                
 double GW_head;                                                
 double h_GWflow, h_SPflow;                             
 double dh;                                                     
 double Qin, Qout, Qout1;                               
 double SWater;                                                 
 double tot_head;                                               
 double CH_vol;                                         
 double CH_vol_R;                                       
 double CH_bottElev;                                    
 double H_rad_ch, H_rad_cell;                   
 double can_cell_M, can_cell_A;                 
 float errorConv;                                               

/* UNUSED below, unimplemented gwater/sfwater integration vars */ 
    int equilib, revers=1;
    float H_pot_don, H_pot_rec;
    
    float fluxTOunsat_don, fluxHoriz; /* vertical and horiz components of the calculated total groundwater Flux */
    float UnsatZ_don, UnsatCap_don, UnsatPot_don, Sat_vs_unsat; /* unsaturated zone capacities, potentials, in the donor cell */
    float UnsatZ_rec, UnsatCap_rec, UnsatPot_rec; /* unsaturated zone capacities, potentials, in the recipient cell */
    float sfTOsat_don, unsatTOsat_don, unsatTOsat_rec, satTOsf_rec, horizTOsat_rec; /* fluxes for vertical integration */
    
    float sedwat_don, sedwat_rec;
    float sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec; /* vertical flux of solutes */
    float sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don; /* vertical flux of solutes */
    float sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz; /* horiz flux of solutes */
/* UNUSED above, unimplemented gwater/sfwater integration vars */ 
 
 
/* look at all the structures connected to the canal and add/subtract flow from/to them */
 Qin = 0; Qout = 0;
 structs = struct_first;
 while ( structs != NULL )
 {
     
     if ( structs->canal_fr  == chan_n+CHo )  
     {
         if ( structs->flow > 0 ) Qout += structs->flow;
         else Qin -= structs->flow;
     } 
     if ( structs->canal_to == chan_n+CHo ) 
     {
         if ( structs->flow > 0 ) Qin += structs->flow;
         else Qout -= structs->flow;
     } 
    
     structs = structs->next_in_list;
 }
  
 CH_vol_R = Chan_list[chan_n]->area * CanWatDep;  /* total canal volume prior to relaxation procedure */
 MIN_VOL =  Chan_list[chan_n]->minVol;                          /* minimum canal volume allowed */
 SPG_coef = Chan_list[chan_n]->SPG_flow_coef;   /* GP_calibGWat already in this */

 I_Length =   Chan_list[chan_n]->seg_len;          /* avg length of cell along reach segment */
 seg_area = I_Length * Chan_list[chan_n]->width;  /* avg area of reach segment */
 can_cell_M = Chan_list[chan_n]->area*CELL_SIZE;  /* used in reducing excess flows between canal<->cell */
 can_cell_A = Chan_list[chan_n]->area+CELL_SIZE;  /* used in reducing excess flows between canal<->cell */
 

 
/* start relaxation procedure */  
 do
 {
     if   ( Abs( error) < F_ERROR ) 
     {
             converged = 1; /* when we've converged, we'll go through the cell list one more time w/o modifying the depth estimate */
                        /* occasionally, we get a modified total flux value, and diff error, a second time with the same depth estimate! */
             errorConv=error;
             factor = 0;
     }
     else if  ( iter >= (max_iter-1) )
     {
         if (attempt == 4) {
             converged = 1; /* haven't actually converged, but stop now */
             sprintf(msgStr,"Day %6.1f: ERROR - couldn't converge on Canal %d, HALTING the iterative relaxation after 4 tries, error = %f.", 
                        SimTime.TIME, Chan_list[chan_n]->number, error); 
             WriteMsg( msgStr,True ); dynERRORnum++;
         }
         else {
             if (debug > 3) { sprintf(msgStr,"Day %6.1f: Warning - couldn't converge on Canal %d, redoing the iterative relaxation, error = %f.", 
                        SimTime.TIME, Chan_list[chan_n]->number, error); 
             WriteMsg( msgStr,True ); }

             attempt++; /* the attempt counter increment */
             iter = 0; /* reset iterations to 0, and double factor & max # iterations on first try */ 
             if (attempt == 1) { factor = init_factor*2.0;  max_iter = CHANNEL_MAX_ITER * 2; }
             else if (attempt == 2) { factor = init_factor*0.5; max_iter = CHANNEL_MAX_ITER * 2; }
             else if (attempt == 3) { factor = init_factor*4.0; max_iter = CHANNEL_MAX_ITER * 4; }
             else if (attempt == 4) { factor = init_factor*0.05; max_iter = CHANNEL_MAX_ITER * 4; }
         }
     }

     else
     {
             /* if error changed sign, we have overshot the target, therefore reverse direction */
         if (error * prev_error < 0 && iter != 1) 
             factor = - sgn(error) *  Abs( factor ) * 0.5;      /* and reduce the step */
         
             /* if error contiues to grow after first iter.(#0), reverse direction */   
         if (iter == 1 && error > 0) factor = - factor;
         
     }
     
     if (iter != 0 && !converged ) 
     {
         CanWatDep += factor; 
         prev_error = error;
         if ( CanWatDep < MINDEPTH ) 
         { 
             if (debug > 3) { sprintf(msgStr,"Day %6.1f: Warning - estimate is below minimum depth, revising estimate. (Canal %d, depth %f)", 
                        SimTime.TIME, Chan_list[chan_n]->number,CanWatDep); 
             WriteMsg( msgStr,True ); }            

             CanWatDep = MINDEPTH; 
            /* converged = 1; */
            factor = -factor*0.1; /* reverse the direction and decrease for next estimate*/
         }      
     }     
    
     T_flux_G = T_flux_S = T_flux_L = 0.0;           /* at each convergence iter, set the sum flows to 0 */
     CH_vol = Chan_list[chan_n]->area * CanWatDep;  /* volume of canal at the estimated water depth */
    
         /* cycle along canal cells */
     cell = Chan_list[chan_n]->cells;
     
     while ( cell != NULL )
     { 
         i = T(cell->x, cell->y);
         if ( !ON_MAP[i] )       /* check that cell is ON_MAP, else we do nothing */
                    /* return to beginning of do loop                */
         {  cell = cell->next_cell;  continue;
         }
         if (basn[i] == 0) {
             sprintf(msgStr,"ERROR - Cell (%d, %d) is out-of-system.", cell->x, cell->y); 
             WriteMsg( msgStr,True );  dynERRORnum++; }

         fluxS = 0.;  fluxG = 0.; fluxL = 0.; /* overland, subsurface, and seepage cell-canal fluxes set to 0 */
           /* v2.5 In the case of some canal reaches, there is a lip or berm, and often associated vegetation, that is 
           along the edge of a canal, decreasing flow rates between canal & marsh.  If user provides a non-zero roughness
           parameter in the CanalData.chan input file, the the mean roughness of this berm and the cell is used for 
           surface exchanges */
         SW_coef = ( MC[i] == 0.0 ) ? 0 : Chan_list[chan_n]->SW_flow_coef / ( (Chan_list[chan_n]->edgeMann > 0.0) ?  (MC[i] + Chan_list[chan_n]->edgeMann)/2 : MC[i] );
         GW_coef = GWcond[i] * poros[HAB[i]];   
         SWater = SWH[i];

         if (debug>0 && (SWater < -GP_MinCheck) ) {
         sprintf(msgStr,"Day %6.1f: capacityERR - neg sfwater along canal %d, candepth %f, cell (%d, %d), celldepth %f", 
                 SimTime.TIME,Chan_list[chan_n]->number,CanWatDep, cell->x, cell->y, SWater); 
         WriteMsg( msgStr,True );  dynERRORnum++; }
         
         GW_head = Min( GWH[i]/poros[HAB[i]], ElevMap[i]); /* v2.8, eqn modified */
         tot_head = ramp(SWater) + GW_head; /* v2.8, eqn modified */
         CH_bottElev = cell->reachElev - Chan_list[chan_n]->depth;  /* elev of bottom of canal at cell location (v2.5+: this is from fixed vector slope instead of actual grid cell elev) */
             
         dh = ( CH_bottElev + CanWatDep ) - tot_head;

             /* hydraulic radius of canal for overland flow out of canal */
         H_rad_ch =   Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
                        SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);   
            /* hydraulic radius of cell for overland flow into canal (same eqn right now) */
         H_rad_cell = Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
                        SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);     

                /* determine the heights of the water cross sections associated with the subsurf and seep flows */
         if (dh > 0.0) { /* from canal */
                h_GWflow = Min(Chan_list[chan_n]->depth, CanWatDep );
                h_SPflow = Max(CH_bottElev + CanWatDep - ElevMap[i], 0.0);
         }
         else if (dh < 0.0) { /* from cell */
                h_GWflow = Max(GW_head-CH_bottElev, 0.0);
                h_SPflow = Max(tot_head-ElevMap[i], 0.0);
         }
         else {
                h_GWflow =  0.0;
                h_SPflow =  0.0;
         }
         

                 
         switch( Chan_list[chan_n]->levee )
         {      
                
             case 2:    /* if levees are on both sides of the canal */
                         /*seepage and subsurface flow on both sides */
                 fluxL =  (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
                 fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
                 break;
             case 0: /* if there are no levees */
                     /*overland flow either side */
                 if ( dh > 0 ) {
                         /* flux from canal to cell  */
                     fluxS = f_Manning ( dh, H_rad_ch, SW_coef);
                 }
                 else if (SWater>GP_DetentZ) {
                     fluxS =  f_Manning ( dh, H_rad_cell, SW_coef) ; 
                    if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
                        fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
                 }
          
                                /* subsurface groundwater flow both sides */
                 fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
                 break;
             case 1:  /* if levee is on the left */ 
                 if ( cell->ind < 0  )
                 { /*overland flow */
                     if ( dh > 0 ) {
                             /* flux from canal to cell  */
                         fluxS = f_Manning ( dh, H_rad_ch, SW_coef);
                     }
                     else if (SWater>GP_DetentZ) {
                         fluxS = f_Manning ( dh, H_rad_cell, SW_coef) ; 
                    if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
                        fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
                     }
          
                 }      
                 else 
                         /*seepage on the levee side */
                     fluxL =  (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
                     
                                /* subsurface groundwater flow both sides */
                 fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
                 break;
             case -1: /* if levee is on the right */
                 if ( cell->ind > 0 )
                 { /*overland flow */
                     if ( dh > 0 ) {
                             /* flux from canal to cell */
                         fluxS = f_Manning( dh, H_rad_ch, SW_coef);
                     }
                     else if (SWater>GP_DetentZ) {
                         fluxS = f_Manning ( dh, H_rad_cell, SW_coef); 
                    if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
                        fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
                     }
                 }      
                 else
                         /*seepage on the levee side */
                     fluxL =  (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);

                                /* subsurface groundwater flow both sides */
                 fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
                 break;  
             default:
                 printf ("Error in levee location: Channel N %d", chan_n+CHo );
                 exit (-1);
                 break;
         } /*end switch */
          


             /* can_cell_M and can_cell_A are the areas of the canal and the cell multiplied (M) and added (A) together, respectively */
             /* the flux eqn is similar to the simple weir eqn for virtual structures, equilibrating the heads across two unequal area regions */

             /* Surface flows: making sure that the cell head is not higher than the canal head */
         if (fluxS>0 && ( (fluxCk = tot_head + fluxS/CELL_SIZE - (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) ) > 0.0 ) ) {
             fluxS = (can_cell_M*(CH_bottElev + CanWatDep)/can_cell_A - can_cell_M*tot_head/can_cell_A );
             fluxCk = tot_head + fluxS/CELL_SIZE - (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area);
             if (fluxCk > F_ERROR && debug>3) {
                 sprintf(msgStr,"Day %6.1f: Warning - excessive head in cell after flux from Canal %d (%f m diff)", 
                         SimTime.TIME, Chan_list[chan_n]->number,fluxCk); 
                 WriteMsg( msgStr,True );  
             }
         }
         
             /* Surface flows: making sure that the canal water level is not flooded much higher than the water level in the cell */
         if ( fluxS<0 && ( (fluxCk = (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) - (tot_head + fluxS/CELL_SIZE) ) > 0.0 ) ) {
             fluxS = (can_cell_M*(CH_bottElev + CanWatDep)/can_cell_A - can_cell_M*tot_head/can_cell_A );
             fluxCk =  (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) - (tot_head + fluxS/CELL_SIZE);
             if (fluxCk > F_ERROR && debug>3) {
                 sprintf(msgStr,"Day %6.1f: Warning - excessive head in Canal %d after flux from cell (%f m diff)", 
                         SimTime.TIME, Chan_list[chan_n]->number,fluxCk); 
                 WriteMsg( msgStr,True );  
             }
         }

             /* now ensure the canal is not drained below its minimum level */ 
         if ( fluxS > 0 || fluxG > 0 || fluxL > 0) /* check gw, seep, and sfwat */
         {
             fluxO  = (fluxS > 0) ? fluxS : 0;
             fluxO += (fluxG > 0) ? fluxG : 0;
             fluxO += (fluxL > 0) ? fluxL : 0;
    
                 /* using vol associated with new estimate (before subtracting losses) */
             if ( (fluxCk = CH_vol - MIN_VOL - fluxS - fluxG - fluxL) < 0.0 )
             { fluxCk =  (fluxO>0) ? ( ( fluxO + fluxCk )/fluxO) : (0.0);       /* v2.5 removed a divide by zero, but shouldn't be needed */
             if (fluxS>0) fluxS *= fluxCk; /* reduce all pos outflows by this proportion */
             if (fluxG>0) fluxG *= fluxCk;
             if (fluxL>0) fluxL *= fluxCk;
             if ( (fluxCk = (CH_vol - MIN_VOL - fluxS - fluxG - fluxL)/Chan_list[chan_n]->area) < -GP_MinCheck ) {
                 sprintf(msgStr,"Day %6.1f: capacityERR - excessive draining from Canal %d (%f m below min)", SimTime.TIME,Chan_list[chan_n]->number,fluxCk); 
                 WriteMsg( msgStr,True ); dynERRORnum++; }
             }
         }

/******** Below is UNUSED *********************/
    if (debug>99 && fluxG != 0.0) do {
        /* This is a part of the routine used to integrate surface/groundwater
         fluxes in the cell-cell Fluxes.c.  It only modifies the allocations if
         a cell is a recipient of groundwater flux (does not deal with cells as
         donor, which is handled in Fluxes.c).
         NOTE:****This routine is NOT fully integrated into the canal iterative solution,
         and thus is NOT used for now - get to it some day if it is determined to be needed.  */
       fluxHoriz =  fluxG;


/**** recipient cell's **pre-flux** capacities */
        UnsatZ_rec  =  (fluxG>0) ? ( ElevMap[i] - GWH[i] / poros[HAB[i]] ) : (0.0) ; /* unsat zone depth */
        if (debug>2 && (UnsatZ_rec < -0.001 ) ) {
            sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in recipient cell (%d,%d) in canal-cell gw fluxes!", 
                    SimTime.TIME, UnsatZ_rec, cell->x, cell->y ); WriteMsg( msgStr,True ); dynERRORnum++; }
            
        UnsatCap_rec =  (fluxG>0) ? (UnsatZ_rec * poros[HAB[i]]) : (0.0); /* maximum pore space capacity (UnsatCap)
                                                              in the depth (UnsatZ) of new unsaturated zone */
        UnsatPot_rec  = (fluxG>0) ? (UnsatCap_rec - Unsat[i] ) : (0.0); /* (height of) the volume of pore space (soil "removed")
                                                           that is unoccupied in the unsat zone */
     /* sf-unsat-sat fluxes */
        horizTOsat_rec = (fluxG>0) ? (fluxHoriz) : (0.0); /* lateral inflow to soil into saturated storage */
        satTOsf_rec = (fluxG>0) ? (Max(fluxHoriz - UnsatPot_rec, 0.0) ) : (0.0); /* upwelling beyond soil capacity */
        if (fluxG>0) unsatTOsat_rec = (UnsatZ_rec > 0.0) ? /* incorporation of unsat moisture into sat storage with
                                                 rising water table due to horiz inflow */
            ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[i]] ) / UnsatZ_rec * Unsat[i] ) :
            (0.0);
        else
            unsatTOsat_rec = 0.0;
        
/**** potential new head in recipient cell will be ... */
        if (fluxG>0)  H_pot_rec = (GWH[i] + horizTOsat_rec + unsatTOsat_rec - satTOsf_rec)
            / poros[HAB[i]] + (SWater + satTOsf_rec);


        equilib = 1;
        
    } while  (!equilib); /* end of ***incomplete** routine for integrating groundwater-surface water and preventing head reversals */
/******** Above is UNUSED *********************/

   /*  if (fluxG>0) fluxG=fluxHoriz; */ /* unused, part of unimplemented sf/gwater routine */
    
         CH_vol -= ( fluxS + fluxG + fluxL); /* subtract flows (pos == loss) from volume associated with new estimate */
         
         if ( converged )
         {      /* the nutrient mass fluxes (pre-flux volumes) (mass=kg, vol=m3)*/
                 /* remember pos flux is from canal */
                 /* (CH_vol_R is vol prior to in/out fluxes) */
        /* surface water canal-cell flux */               
             if ( fluxS != 0 ) {
             if ( fluxS > 0 ) /* surface water flux from canal to cell */
             {  
                 /* this is where Disp_Calc comes in */ /* hcf TEST to see about reducing dispersion, halfing the flux of salt (i.e., ca. 500m grid dispersion) */
                 N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxS) : (0.0);
                 P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxS) : (0.0);
                 S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxS) : (0.0);

                 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
                     if ( (Chan_list[chan_n]->family !=  basn_list[basn[i]]->family) && (debug > 0) ) {        /* if the flow is not within the family... */
                         sprintf(msgStr,"Day %6.1f: ERROR - no basin-basin canal-cell overland flow! (Reach %d, basn %d, with cell %d,%d, basin %d)", 
                                 SimTime.TIME,Chan_list[chan_n]->number, Chan_list[chan_n]->basin, cell->x, cell->y, basn[i]); 
                         WriteMsg( msgStr,True );  dynERRORnum++;
                     }
                     else {    /* if it's flow within a family, just keep
                                  track of what the children do among themselves */            
                         if  ( !Chan_list[chan_n]->parent  ) { 
                         /* then find out about the flow for the family's sake */
                             VOL_OUT_OVL[Chan_list[chan_n]->basin] += fluxS;
                             P_OUT_OVL[Chan_list[chan_n]->basin] += P_fl;
                             S_OUT_OVL[Chan_list[chan_n]->basin] += S_fl;
                         }
                         if ( !basn_list[basn[i]]->parent ) {                 
                                 /* then find out about the flow for the family's sake */
                             VOL_IN_OVL[basn[i]] += fluxS; 
                             P_IN_OVL[basn[i]] += P_fl; 
                             S_IN_OVL[basn[i]] += S_fl; 
                         }
                     }
                 }                
             }
             else /* surface water flux from cell to canal */
             {
                 /* this is where Disp_Calc comes in */
                 N_fl = ( SWH[i] > 0 ) ? NA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
                 P_fl = ( SWH[i] > 0 ) ? PA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
                 S_fl = ( SWH[i] > 0 ) ? SA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
                 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {

                     if ( (Chan_list[chan_n]->family !=  basn_list[basn[i]]->family) && (debug > 0) ) {        /* if the flow is not within the family... */
                         sprintf(msgStr,"Day %6.1f: ERROR - no basin-basin canal-cell overland flow! (Reach %d, basn %d, with cell %d,%d, basin %d)", 
                                        SimTime.TIME,Chan_list[chan_n]->number, Chan_list[chan_n]->basin, cell->x, cell->y, basn[i]); 
                         WriteMsg( msgStr,True );  dynERRORnum++;
                 }
                 else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                     if  ( !basn_list[basn[i]]->parent  ) { 
                         /* then find out about the flow for the family's sake */
                         VOL_OUT_OVL[basn[i]] -= fluxS;
                         P_OUT_OVL[basn[i]] -= P_fl;
                         S_OUT_OVL[basn[i]] -= S_fl;
                     }
                     if ( !Chan_list[chan_n]->parent ) {                 
                         /* then find out about the flow for the family's sake */
                         VOL_IN_OVL[Chan_list[chan_n]->basin] -= fluxS; 
                         P_IN_OVL[Chan_list[chan_n]->basin] -= P_fl; 
                         S_IN_OVL[Chan_list[chan_n]->basin] -= S_fl; 
                     }
                 }
                 }


             } 
                 /* now pass that mass flux to/from the cell/canal */
             NA[i] += N_fl;
             PA[i] += P_fl;
             SA[i] += S_fl;
             Chan_list[chan_n]->N -= N_fl;
             Chan_list[chan_n]->P -= P_fl;
             Chan_list[chan_n]->S -= S_fl;
             }

        /* groundwater canal-cell flux */                 
             if ( fluxG != 0 ) {
             if ( fluxG > 0 ) /* ground water flux from canal to cell */
             {  
                 N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxG) : (0.0);
                 P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxG) : (0.0);
                 S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxG) : (0.0);

                 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
                     if (Chan_list[chan_n]->family !=  
                         basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
                         if  ( !Chan_list[chan_n]->parent  ) { /* and if the donor or recipient is a child... */
                         /* then find out about the flow for the family's sake */
                             VOL_OUT_GW[Chan_list[chan_n]->family] += fluxG;
                             P_OUT_GW[Chan_list[chan_n]->family] += P_fl;
                             S_OUT_GW[Chan_list[chan_n]->family] += S_fl;
                         }
                         if ( !basn_list[basn[i]]->parent ) {                 
                           /* then find out about the flow for the family's sake */
                             VOL_IN_GW[basn_list[basn[i]]->family] += fluxG; 
                             P_IN_GW[basn_list[basn[i]]->family] += P_fl; 
                             S_IN_GW[basn_list[basn[i]]->family] += S_fl; 
                         }
                     VOL_OUT_GW[Chan_list[chan_n]->basin] += fluxG;
                     VOL_IN_GW[basn[i]] += fluxG; 
                     P_OUT_GW[Chan_list[chan_n]->basin] += P_fl;
                     P_IN_GW[basn[i]] += P_fl; 
                     S_OUT_GW[Chan_list[chan_n]->basin] += S_fl;
                     S_IN_GW[basn[i]] += S_fl; 
                 }
                 else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                     if  ( !Chan_list[chan_n]->parent  ) { 
                         /* then find out about the flow for the family's sake */
                         VOL_OUT_GW[Chan_list[chan_n]->basin] += fluxG;
                         P_OUT_GW[Chan_list[chan_n]->basin] += P_fl;
                         S_OUT_GW[Chan_list[chan_n]->basin] += S_fl;
                     }
                     if ( !basn_list[basn[i]]->parent ) {                 
                         /* then find out about the flow for the family's sake */
                         VOL_IN_GW[basn[i]] += fluxG; 
                         P_IN_GW[basn[i]] += P_fl; 
                         S_IN_GW[basn[i]] += S_fl; 
                     }
                 }
                 }
             } /* end of positive flow evaluations */
             
             
             else  /*  ground water flux from cell to canal */
             {
                 N_fl = ( GWH[i] > 0 ) ? GNA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
                 P_fl = ( GWH[i] > 0 ) ? GPA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
                 S_fl = ( GWH[i] > 0 ) ? GSA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
                 if (Chan_list[chan_n]->basin != basn[i] && debug > -999) {
                     if (Chan_list[chan_n]->family !=  
                         basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
                         if  ( !basn_list[basn[i]]->parent  ) { /* and if the donor or recipient is a child... */
                         /* then find out about the flow for the family's sake */
                             VOL_OUT_GW[basn_list[basn[i]]->family] -= fluxG;
                             P_OUT_GW[basn_list[basn[i]]->family] -= P_fl;
                             S_OUT_GW[basn_list[basn[i]]->family] -= S_fl;
                         }
                         if ( !Chan_list[chan_n]->parent ) {                 
                           /* then find out about the flow for the family's sake */
                             VOL_IN_GW[Chan_list[chan_n]->family] -= fluxG; 
                             P_IN_GW[Chan_list[chan_n]->family] -= P_fl; 
                             S_IN_GW[Chan_list[chan_n]->family] -= S_fl; 
                         }
                     VOL_IN_GW[Chan_list[chan_n]->basin] -= fluxG;
                     VOL_OUT_GW[basn[i]] -= fluxG;
                     P_IN_GW[Chan_list[chan_n]->basin] -= P_fl;
                     P_OUT_GW[basn[i]] -= P_fl;
                     S_IN_GW[Chan_list[chan_n]->basin] -= S_fl;
                     S_OUT_GW[basn[i]] -= S_fl;
                 }
                 else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                     if  ( !basn_list[basn[i]]->parent  ) { 
                         /* then find out about the flow for the family's sake */
                         VOL_OUT_GW[basn[i]] -= fluxG;
                         P_OUT_GW[basn[i]] -= P_fl;
                         S_OUT_GW[basn[i]] -= S_fl;
                     }
                     if ( !Chan_list[chan_n]->parent ) {                 
                         /* then find out about the flow for the family's sake */
                         VOL_IN_GW[Chan_list[chan_n]->basin] -= fluxG; 
                         P_IN_GW[Chan_list[chan_n]->basin] -= P_fl; 
                         S_IN_GW[Chan_list[chan_n]->basin] -= S_fl; 
                     }
                 }
                 }
             } 
                 /* now pass that mass flux to/from the cell/canal
                  this includes allocating nuts to surface water thru upflow (unimplemented) */
             /* unused, part of unimplemented sf/gwater routine */
/*              SNfl = Chan_list[chan_n]->N/CH_vol_R*satTOsf_rec ; */
/*              GNfl = N_fl - SNfl; */
/*              SPfl = Chan_list[chan_n]->P/CH_vol_R*satTOsf_rec ; */
/*              GPfl = P_fl - SPfl; */
/*              SSfl = Chan_list[chan_n]->S/CH_vol_R*satTOsf_rec ; */
/*              GSfl = S_fl - SSfl; */
/*              GNA[i] += GNfl; */
/*              GPA[i] += GPfl; */
/*              GSA[i] += GSfl; */
             GNA[i] += N_fl;
             GPA[i] += P_fl;
             GSA[i] += S_fl;
             Chan_list[chan_n]->N -= N_fl;
             Chan_list[chan_n]->P -= P_fl;
             Chan_list[chan_n]->S -= S_fl;
                 /* upflow into surface water nuts (this is actually canal conc., but what the heck) */
             /* unused, part of unimplemented sf/gwater routine */
/*              NA[i] += SNfl; */
/*              PA[i] += SPfl; */
/*              SA[i] += SSfl; */
             }

        /* seepage canal-cell flux */             
             if ( fluxL != 0 ) {
             if ( fluxL > 0 ) /* seepage flux from canal to cell */
             {  
                 N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxL) : (0.0);
                 P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxL) : (0.0);
                 S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxL) : (0.0);
                 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
                     if (Chan_list[chan_n]->family !=  
                         basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
                         if  ( !Chan_list[chan_n]->parent  ) { /* and if the donor or recipient is a child... */
                         /* then find out about the flow for the family's sake */
                             VOL_OUT_SPG[Chan_list[chan_n]->family] += fluxL;
                             P_OUT_SPG[Chan_list[chan_n]->family] += P_fl;
                             S_OUT_SPG[Chan_list[chan_n]->family] += S_fl;
                         }
                         if ( !basn_list[basn[i]]->parent ) {                 
                           /* then find out about the flow for the family's sake */
                             VOL_IN_SPG[basn_list[basn[i]]->family] += fluxL; 
                             P_IN_SPG[basn_list[basn[i]]->family] += P_fl; 
                             S_IN_SPG[basn_list[basn[i]]->family] += S_fl; 
                         }
                     VOL_OUT_SPG[Chan_list[chan_n]->basin] += fluxL;
                     VOL_IN_SPG[basn[i]] += fluxL; 
                     P_OUT_SPG[Chan_list[chan_n]->basin] += P_fl;
                     P_IN_SPG[basn[i]] += P_fl; 
                     S_OUT_SPG[Chan_list[chan_n]->basin] += S_fl;
                     S_IN_SPG[basn[i]] += S_fl; 
                 }
                 else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                     if  ( !Chan_list[chan_n]->parent  ) { 
                         /* then find out about the flow for the family's sake */
                         VOL_OUT_SPG[Chan_list[chan_n]->basin] += fluxL;
                         P_OUT_SPG[Chan_list[chan_n]->basin] += P_fl;
                         S_OUT_SPG[Chan_list[chan_n]->basin] += S_fl;
                     }
                     if ( !basn_list[basn[i]]->parent ) {                 
                         /* then find out about the flow for the family's sake */
                         VOL_IN_SPG[basn[i]] += fluxL; 
                         P_IN_SPG[basn[i]] += P_fl; 
                         S_IN_SPG[basn[i]] += S_fl; 
                     }
                 }
                 }
                     if ( tot_head > GW_head ) { /* to cell surface water */
                                NA[i] += N_fl;
                                PA[i] += P_fl;
                                SA[i] += S_fl;
                                }
                         else { /* to cell groundwater */
                                GNA[i] += N_fl;
                                GPA[i] += P_fl;
                                GSA[i] += S_fl;
                        }
             }
             else if (fluxL < 0) /*  seepage flux from cell to canal */
             {
                 if ( tot_head > GW_head ) { /* cell's tot head > its groundwater head */
                        N_fl = ( SWH[i] > 0 ) ? NA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
                        P_fl = ( SWH[i] > 0 ) ? PA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
                        S_fl = ( SWH[i] > 0 ) ? SA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
                        NA[i] += N_fl;
                        PA[i] += P_fl;
                        SA[i] += S_fl;
                 }
                 else {
                        N_fl = ( GWH[i] > 0 ) ? GNA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
                        P_fl = ( GWH[i] > 0 ) ? GPA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
                        S_fl = ( GWH[i] > 0 ) ? GSA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
                        GNA[i] += N_fl;
                        GPA[i] += P_fl;
                        GSA[i] += S_fl;
                }
                 if (Chan_list[chan_n]->basin != basn[i] && debug > -999) {
                     if (Chan_list[chan_n]->family !=  
                         basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
                         if  ( !basn_list[basn[i]]->parent  ) { /* and if the donor or recipient is a child... */
                         /* then find out about the flow for the family's sake */
                             VOL_OUT_SPG[basn_list[basn[i]]->family] -= fluxL;
                             P_OUT_SPG[basn_list[basn[i]]->family] -= P_fl;
                             S_OUT_SPG[basn_list[basn[i]]->family] -= S_fl;
                         }
                         if ( !Chan_list[chan_n]->parent ) {                 
                           /* then find out about the flow for the family's sake */
                             VOL_IN_SPG[Chan_list[chan_n]->family] -= fluxL; 
                             P_IN_SPG[Chan_list[chan_n]->family] -= P_fl; 
                             S_IN_SPG[Chan_list[chan_n]->family] -= S_fl; 
                         }
                     VOL_IN_SPG[Chan_list[chan_n]->basin] -= fluxL;
                     VOL_OUT_SPG[basn[i]] -= fluxL;
                     P_IN_SPG[Chan_list[chan_n]->basin] -= P_fl;
                     P_OUT_SPG[basn[i]] -= P_fl;
                     S_IN_SPG[Chan_list[chan_n]->basin] -= S_fl;
                     S_OUT_SPG[basn[i]] -= S_fl;
                 }
                 else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                     if  ( !basn_list[basn[i]]->parent  ) { 
                         /* then find out about the flow for the family's sake */
                         VOL_OUT_SPG[basn[i]] -= fluxL;
                         P_OUT_SPG[basn[i]] -= P_fl;
                         S_OUT_SPG[basn[i]] -= S_fl;
                     }
                     if ( !Chan_list[chan_n]->parent ) {                 
                         /* then find out about the flow for the family's sake */
                         VOL_IN_SPG[Chan_list[chan_n]->basin] -= fluxL; 
                         P_IN_SPG[Chan_list[chan_n]->basin] -= P_fl; 
                         S_IN_SPG[Chan_list[chan_n]->basin] -= S_fl; 
                     }
                 }
                 }
             } 
                 /* now pass that mass flux to/from the canal */
             Chan_list[chan_n]->N -= N_fl;
             Chan_list[chan_n]->P -= P_fl;
             Chan_list[chan_n]->S -= S_fl;
             }
                                /* now we recalculate the depths of SFWater and Gwater in cell */
/*              SWH[i] += (fluxS+satTOsf_rec)/CELL_SIZE;  */
/*              GWH[i] += (fluxG-satTOsf_rec)/CELL_SIZE;  */
             SWH[i] += (fluxS)/CELL_SIZE; 
             GWH[i] += (fluxG)/CELL_SIZE; 
             if ( tot_head > GW_head ) SWH[i] += fluxL/CELL_SIZE; /* whether to add seepage to surface or */
             else GWH[i] += fluxL/CELL_SIZE;                      /* groundwater is not critical, as vertical updates will follow */
                                

         } /* that's all we do when finally converged */
                    
         T_flux_S += fluxS;
         T_flux_G += fluxG;
         T_flux_L += fluxL;
                  
         cell = cell->next_cell;
     } /*end cycle along canal cells */

         
     error = (CanWatDep - Chan_list[chan_n]->wat_depth) - (Qin - Qout - T_flux_S - T_flux_G - T_flux_L)/ Chan_list[chan_n]->area;
     if (converged && error != errorConv ) {
         /* error the second time around should be close to that of the first */
         if (fabs(error-errorConv)>F_ERROR ) { 
             sprintf(msgStr,"Day %6.1f: ERROR - converg error != to error after cell-canal fluxes (Chan %d, est= %5.3f, last depth= %5.3f m)", 
                        SimTime.TIME, Chan_list[chan_n]->number, CanWatDep, Chan_list[chan_n]->wat_depth); 
             WriteMsg( msgStr,True );  dynERRORnum++;
         }
         
         if (CanWatDep>Chan_list[chan_n]->depth*3.0 || CanWatDep < 0.0) {
                CanWatDep = MINDEPTH; /* lost the mass balance, but let it go for now */
                sprintf(msgStr,"           Reset Chan %d depth to minimum allowed (%5.3f m)", 
                                Chan_list[chan_n]->number, MINDEPTH); 
                WriteMsg( msgStr,True );  }
                 
                 /*errorConv is the error upon first reaching convergence criteria; */
             /* we then recalc the fluxes in the same loop thru cells (w/ same depth est.), */
             /* and should have the same fluxes and error as before */
         if (debug>4) { /* this here to halt program with high debug */
             printf("\nDay %6.1f: ERROR - converg error != to error after cell-canal fluxes (Chan %d, est= %5.3f m)\n", 
                        SimTime.TIME, Chan_list[chan_n]->number, CanWatDep);
             exit (-1);

         }
         
     }
     
     ++iter;

 } while (!converged );  /* end relaxation procedure */

 if ( factor != 0 && Qout > 0) /* if relaxation was abnormally stopped because of hitting the bottom */
 { 
     sprintf(msgStr,"Day %6.1f: ERROR - hit bottom of Canal %d, stopped the iterative relaxation", SimTime.TIME, Chan_list[chan_n]->number); 
     WriteMsg( msgStr,True );     
     dynERRORnum++;         
            /*#######*******########  this old code should not be reached
             the code is questionable - an ERROR message printed
             to file will indicate signif problem due to entrance to this */
     Qout1 = (Chan_list[chan_n]->wat_depth - CanWatDep)*Chan_list[chan_n]->area + Qin - T_flux_S - T_flux_G - T_flux_L;
     structs = struct_first;
     while ( structs != NULL )   /* recalculate the outflow in the data for possible use in the next canal inflow */
     {  
         if ( structs->canal_fr  == chan_n+CHo )   structs->flow *= Qout1/Qout;
    
         structs = structs->next_in_list;
     }
     Qout = Qout1;
 } 

 Chan_list[chan_n]->wat_depth = CanWatDep;      /* new canal water depth */
 CH_vol =  CanWatDep * Chan_list[chan_n]->area; /* new canal volume */

 if( debug > 3 ) 
 {
         /* concentration in mg/L (kg/m3==g/L, convert P (not S) to mg/L) */
     Chan_list[chan_n]->P_con = (CH_vol > 0) ? (Chan_list[chan_n]->P/CH_vol*1000.0) : 0.0;
     Chan_list[chan_n]->S_con = (CH_vol > 0) ? (Chan_list[chan_n]->S/CH_vol) : 0.0;
     fprintf( canDebugFile, 
              "%4d  %8.3f  %8.4f  %8.4f  %12.1f  %12.1f  %12.1f  %12.1f  %12.1f  %5d  %9.1f  %7.5f \n",
              Chan_list[chan_n]->number, Chan_list[chan_n]->wat_depth, Chan_list[chan_n]->P_con,  
              Chan_list[chan_n]->S_con, T_flux_S, T_flux_G, T_flux_L, Qin, Qout, iter, Chan_list[chan_n]->area, error );
 
 }

 return;
}

Here is the call graph for this function:

float f_Manning	(	float	delta,
		float	Water,
		float	SW_coef
	)

Surface water exchange between cell and canal.

Parameters:

	delta	Head difference between cell and canal (m)
	Water	Hydraulic radius (m)
	SW_coef	Aggregated flow coefficient (m^0.5 * sec)/m^(1/3)

Returns:

m^3

Definition at line 2741 of file WatMgmt.c.

References Abs, canstep, GP_mannDepthPow, and sgn.

Referenced by FluxChannel().

{       float a_delta;
    
        a_delta = Abs (delta);
        return ( sgn( delta ) * SW_coef * pow(Water, GP_mannDepthPow) * sqrt(a_delta) * canstep);

}

float f_Ground	(	float	dh,
		float	height,
		float	GW_coef,
		float	I_Length
	)

Subsurface ground water exchange between cell and canal.

Parameters:

	dh	Head difference between cell and canal (m)
	height	Height of water (m)
	GW_coef	Aggregated flow coefficient (m/d)
	I_Length	Average length of cell along reach segment

Returns:

m^3

Definition at line 2759 of file WatMgmt.c.

References canstep, and celWid.

Referenced by FluxChannel().

{       
        float  w;

/* porosity already part of GW_coef */
        w = dh * I_Length * GW_coef / (0.5*celWid) * height * canstep; 

        return ( w);

}

void Flows_in_Structures	(	float *	SWH,
		float *	GW,
		float *	poros,
		float *	Elev,
		double *	NA,
		double *	PA,
		double *	SA
	)

Calling function to functions that determine flows through water control structures.

Parameters:

	SWH	SURFACE_WAT
	GW	SAT_WATER
	poros	HP_HYD_POROSITY
	Elev	SED_ELEV
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Remarks:

Database (Structs_attr) formatting notes/requirements:

• 1a) For data-driven structures, no particular sequence is necessary, but the sequence of the structures in the CanalData.struct and CanalData.struct_wat files must be the same (and the prog will catch and whip you if you try otherwise).
• 1b) Virtual structures and rule-based structures (i.e., anything other than data-driven structures) MUST all come after any data-driven structures.
• 2) Due to the cycle that checks for multiple historical/other-data outflows from a canal, be sure that the database (to CanalData.struct) does not erroneously have 2 types (canal and cell) of recipients for a given structure.
• 3) There need to be 2 header lines (as described in the dbase: run-name and col IDs).
• 4) Flow is expected to ALWAYS be positive, even though there remain checks for reverse flows in the code. For a two-way structure, designate a separate structure for the reverse flow.
• 5) As labeled in dbase, the concentration of P and N in structure flows are ug/L (ppb), and S is in g/L (ppt). (The code converts and deals with kg, m^3, and kg/m^3 (==g/L) ).
• 6) The water flows in the input data are cfs, converted and used here as m^3/d.
• 7) For structures with model-calculated flows, don't forget to provide the cell location of the structure itself so that the elevation at a canal-canal flow structure is available.
• 8) ******* You MAY NOT have multiple historical/other-data demands from a particular cell - only canals are allowed such multiple demands.
• 9) If you want to apply a concentration to an internal structure that is one of multiple outflows from a canal, you need to apply that conc to all structures (actually the first in the sequence of the database) draining that canal. (This would only be for tracer studies etc., all internal model structure concentrations are calc'd).
  

Flow to/from the on-map, within-model-domain area, can only occur by a flow to/from a "cell" mapped at (1,1). That's based on our standard (0,0) upper left cell: the code uses a (1,1) origin - so that (1,1) from a map => (2,2) in the code!

Variables local to function

flow_code dimless attribute, representing the type of raster-vector flow

Note:

re-activated used of head and tail-water stage schedules, for tidal boundary conditions and targets for simple rule-based calc's (v2.8) of managed flows

Definition at line 2811 of file WatMgmt.c.

References Structure::conc_N, Structure::conc_P, Structure::conc_S, Structure::flag, Structure::flagHWsched, Structure::flagHWtide, Structure::flagTWsched, Structure::flagTWtide, Structure::flow, flowCalc_CanCan(), flowCalc_CanCel(), flowCalc_CelCan(), flowCalc_CelCel(), flowData_CanCan(), flowData_CanCel(), flowData_CelCan(), flowData_CelCel(), FMOD(), GetGraph(), GP_SLRise, Structure::HW_graph, Structure::HW_stage, Structure::multiOut, Structure::next_in_list, printchan, Structure::S_nam, SimTime, Structure::SrcDest, struct_first, structs, Structure::Sum_P, Structure::Sum_S, Structure::SumFlow, simTime::TIME, Structure::TW_graph, Structure::TW_stage, WstructOutFile, WstructOutFile_P, and WstructOutFile_S.

Referenced by Run_Canal_Network().

{
        struct Structure *structs, *struct_hist_start;
    int flow_code; /* v2.8 */
    
    structs = struct_first; 

    while ( structs != NULL )
    {           
        if ( structs->flag < 0 || (structs->flag >1) ) { /* nothing to do when structure inactive (neg or >1 flag) */
            structs->flow = 0.0; /* used in summing up flows for Qin and Qout later */
            structs->conc_P = 0.0;
            structs->conc_N = 0.0;
            structs->conc_S = 0.0;
        } 
        else {
        
        if ( structs->flagTWsched == 1 ) /* this means that there is a time series schedule graph for TailWater target */
        { 
                structs->TW_stage = GetGraph ( structs->TW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 0.0) ); 
                
                /* add sea level rise (annual rate in GlobalParms input) if this is a tidaly-driven schedule */
                if (structs->flagTWtide == 1) structs->TW_stage  = structs->TW_stage + (SimTime.TIME/365) * GP_SLRise ;
        }
          
        if ( structs->flagHWsched == 1 ) /* this means that there is a time series schedule graph for HeadWater target */
        { 
                structs->HW_stage = GetGraph ( structs->HW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 0.0) ); 
                
                /* add sea level rise (annual rate in GlobalParms input) if this is a tidaly-driven schedule */
                if (structs->flagHWtide == 1) structs->HW_stage  = structs->HW_stage + (SimTime.TIME/365) * GP_SLRise ;
        }
        
        if (strcmp(structs->SrcDest,"CanCan") == 0) flow_code = 0;
        if (strcmp(structs->SrcDest,"CelCel") == 0) flow_code = 1;
        if (strcmp(structs->SrcDest,"CanCel") == 0) flow_code = 2;
        if (strcmp(structs->SrcDest,"CelCan") == 0) flow_code = 3;
        
        /* call the modules to move water & constituents thru the specific source-destination type of structure */
        switch(flow_code) {
                case 0:
                        if (structs->flag == 1) flowData_CanCan(structs, struct_hist_start, SWH, NA, PA, SA); 
                        else flowCalc_CanCan(structs, SWH, GW, poros, Elev, NA, PA, SA); 
                        break; 
                case 1:
                        if (structs->flag == 1) flowData_CelCel(structs, SWH, NA, PA, SA); 
                        else flowCalc_CelCel(structs, SWH, Elev, NA, PA, SA);
                        break; 
                case 2:
                        if (structs->flag == 1) flowData_CanCel(structs, struct_hist_start, SWH, NA, PA, SA); 
                        else flowCalc_CanCel(structs, SWH, GW, poros, Elev, NA, PA, SA); 
                        break;
                case 3:
                        if (structs->flag == 1) flowData_CelCan(structs, SWH, NA, PA, SA); 
                        else flowCalc_CelCan(structs, SWH, GW, poros, Elev, NA, PA, SA);
                        break; 
                default: {
                        printf ("\nError in source-destination data for structure N %s: impossible condition! \n", structs->S_nam );
                        exit (-1);
                }
        }
                }       

            /* sum of flows and solute masses thru structs over the can_Intvl canal summary interval */
        structs->SumFlow += structs->flow; 
        structs->Sum_P += structs->conc_P * structs->flow; 
/*        structs->Sum_N += structs->conc_N * structs->flow; */
        structs->Sum_S += structs->conc_S * structs->flow; 

        structs->multiOut = 0; /* reset the multi-outflow flag to zero after done */
        
            /* print sum of flows and the flow weighted mean conc thru structures */
        if ( printchan ) { 
            fprintf( WstructOutFile, "%7.0f\t", (structs->SumFlow)/(2446.848) ); /* convert to cfs (0.02832*60*60*24) */
            fprintf( WstructOutFile_P, "%7.4f\t", 
                     (structs->SumFlow > 0.0) ? (structs->Sum_P/structs->SumFlow*1000.0) : (0.0) ); /* convert g/L to mg/L */
/*            fprintf( WstructOutFile_N, "%7.4f\t", 
              (structs->SumFlow > 0.0) ? (structs->Sum_N/structs->SumFlow*1000.0) : (0.0) );*/ /* convert g/L to mg/L */
            fprintf( WstructOutFile_S, "%7.4f\t", 
                     (structs->SumFlow > 0.0) ? (structs->Sum_S/structs->SumFlow) : (0.0) ); /* g/L, no conversion */

            structs->SumFlow = 0.0; 
            structs->Sum_P = 0.0;
/*            structs->Sum_N = 0.0; */
            structs->Sum_S = 0.0;
        }
        
            /* increment to the next structure in the list */
        structs = structs->next_in_list;
    }

    if (printchan) {
        fflush( WstructOutFile );
        fflush( WstructOutFile_P );
/*         fflush( WstructOutFile_N ); */
        fflush( WstructOutFile_S );
    }
 
    return;
    
}

Here is the call graph for this function:

void flowData_CanCan	(	struct Structure *	structs,
		struct Structure *	struct_hist_start,
		float *	SWH,
		double *	NA,
		double *	PA,
		double *	SA
	)

Canal-to-Canal water control structure flows: Data-driven (historical/sfwmm).

Parameters:

	SWH	SURFACE_WAT
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Definition at line 2927 of file WatMgmt.c.

References Chan::area, HistData::arrayN, HistData::arrayP, HistData::arrayS, HistData::arrayWat, Chan::basin, basn, Structure::canal_fr, Structure::canal_to, canDebugFile, canstep, Structure::cell_i_to, Structure::cell_j_to, CELL_SIZE, Chan_list, CHo, Structure::conc_N, Structure::conc_P, Structure::conc_S, debug, dynERRORnum, Structure::flag, Structure::flow, Hist_data, Structure::histID, MINDEPTH, msgStr, Structure::multiOut, Chan::N, Structure::next_in_list, Chan::number, Chan::P, P_IN_STR, P_OUT_STR, ramp, Chan::S, S_IN_STR, Structure::S_nam, S_OUT_STR, SimTime, Structure::str_cell_i, Structure::str_cell_j, Chan::sumHistIn, Chan::sumHistOut, T, simTime::TIME, Structure::TN, Structure::TNser, Structure::TP, Structure::TPser, True, Structure::TS, Structure::TSser, VOL_IN_STR, VOL_OUT_STR, Chan::wat_depth, and WriteMsg().

Referenced by Flows_in_Structures().

{
    int cell_to, cell_fr, can_to, can_fr;
    int cell_struct;
    double CH_vol;
    double Nconc, Pconc, Sconc;
    double N_fl, P_fl, S_fl; 
    double chan_over;
    double ChanHistOut;

    can_fr = structs->canal_fr - CHo;
    can_to = structs->canal_to - CHo;
    cell_struct = T(structs->str_cell_i,structs->str_cell_j); 
                
    /* DATA-DRIVEN STRUCTURE FLOW */
                if (structs->multiOut == 1) return; /* this structure has already been processed as a multiple outflow from a canal */
                    /* FLOW */
                    /* don't do much with zero or negative flows */           
                structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
                if (structs->flow == 0.0) { 
                    structs->conc_P = 0.0;
                    structs->conc_N = 0.0;
                    structs->conc_S = 0.0;
                    return; 
                }   /* zero the reported struct flow conc (for output info only - mass is used in tracking actual flows) */
                else if (structs->flow < 0.0) {
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
                             SimTime.TIME, Chan_list[can_to]->number );
                    WriteMsg( msgStr,True );
                    dynERRORnum++;
                    return;
                }

                    /* for multiple historical/other-data outflows from a canal, ensure mass balance (sum of outflows !> avail volume) */
                    /* all flows from the canal associated with each set must be non-negative, which is currently valid */
                    /* the below cycle thru the structures will also process a canal->cell flow structure */

                struct_hist_start = structs; /* this is the first structure w/ historical info for a particular canal */
                ChanHistOut = 0.0; /* sum of all historical/other-data outflows from a canal */

                while ( structs != NULL )   
                {  
                    if  ( structs->flag == 1 && struct_hist_start->canal_fr  == structs->canal_fr   ) { /* look for outflows from same canal */
                        structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
                        ChanHistOut += structs->flow; /* sum the multiple historical/other-data outflows from canal */
                    }
                        
                    structs = structs->next_in_list; /* increment to next structure */
                } /* end cycle thru list of remaining structures in total list */
                    
                structs = struct_hist_start; /* now we are pointing back to the first (or only) water control structure draining the canal */
                         
                    /* concentration in donor canal for structure flow, g/L = kg/m^3 */
                    /* this is based on previous depth/volume calc'd in FluxChannel */
                    /* or use historical conc */

                    /* CONC */
                                /* canal TOTAL volume before new structure flows */
                CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth; 
                    /* concentrations applied to all structures draining this canal */
                Pconc = (structs->TPser == -1) ? 
                    ( (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0) ) :
                    ( (structs->TPser == 0) ? (structs->TP) : ( Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
                Nconc = (structs->TNser == -1) ? 
                    ( (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0) ) :
                    ( (structs->TNser == 0) ? (structs->TN) : ( Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
                Sconc = (structs->TSser == -1) ? 
                    ( (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0) ) :
                    ( (structs->TSser == 0) ? (structs->TS) : ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );


                    /* vol avail for fluxing */
                CH_vol =  Chan_list[can_fr]->area*(ramp(Chan_list[can_fr]->wat_depth-MINDEPTH) );
                chan_over = (ChanHistOut > 0.0 ) ? ( CH_vol/ChanHistOut ) : 1.0;

                do /* loop again through the remaining sequence of structures to correct any outflows for the canal */
                {   
                    if (structs->flag == 1 &&  struct_hist_start->canal_fr  == structs->canal_fr   ) { 

                        structs->multiOut = 1; /* set the flag to indicate this structure has been processed */
                           
                            /* revised FLOW */
                        if ( chan_over < 1.0 ) { /* if hist/other-data demand > canal vol */
                            if( debug > 0 ) { 
                                fprintf( canDebugFile, "Day \t%6.1f\t %s:\tTotDemand= \t%7.0f\t; exceeded chan \t%d\t avail vol= \t%7.0f \tby \t%7.0f \tm3 (\t%5.3f \tm deep).\n",
                                         SimTime.TIME, structs->S_nam, ChanHistOut, Chan_list[can_fr]->number,
                                         CH_vol, (ChanHistOut - CH_vol), Chan_list[can_fr]->wat_depth );
                                
                            }  
                                /* decrement the individual structure outflows by it's proportional contrib. */
                            structs->flow = (structs->flow) * chan_over ;
                        } 
                            /* sum of all realized (corrected for excess demand) historical/other-data outflows from a canal during an iteration */
                            /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
                        Chan_list[can_fr]->sumHistOut += structs->flow;         

                            /* the recipient changes as we cycle thru the structures (may be another canal or a cell) */
                        can_to = structs->canal_to - CHo;
                        cell_to = T(structs->cell_i_to,structs->cell_j_to); 

                        structs->conc_P = Pconc;/* concentrations applied to all structures draining this canal */
                        structs->conc_N = Nconc;
                        structs->conc_S = (structs->TSser == -1) ?
                            (Sconc) :
                            ( (structs->TSser == 0) ? (structs->TS) :
                              ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ); /* see note below on tracer conc. application */
                        
                        P_fl = structs->flow * structs->conc_P; 
                        N_fl = structs->flow * structs->conc_N;
                        S_fl = structs->flow * structs->conc_S;
                                /* change in nut mass in fr canal (kg) */
                        Chan_list[can_fr]->N -=  N_fl ;
                        Chan_list[can_fr]->P -=  P_fl ;
                        Chan_list[can_fr]->S -=  (structs->TSser == -1) ? (S_fl) : (0);
                            /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
                               Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
                               that was not calculated from the available water and salt (i.e., the conc. was from the
                               CanalData.struct data file, either fixed or a time series) */
                            
                        if (can_to > 0) {    /* IF this is a canal recip, change in nut mass in to canal (kg) */
                                /* sum of all realized (corrected for excess demand) historical/other-data outflows into a canal during an iteration */
                                /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
                            Chan_list[can_to]->sumHistIn += structs->flow;              

                            Chan_list[can_to]->N +=  N_fl ;
                            Chan_list[can_to]->P +=  P_fl ;
                            Chan_list[can_to]->S +=  S_fl ;
                                /* mass balance in fr and to canals */
                                /* canals themselves always belong to a parent hydro basin,
                                   thus canal-canal struct flows do not keep track of child-parent, child-child flows */
                            if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > -999 ) {
                                VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
                                VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; 
                                P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
                                P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
                                S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
                                S_IN_STR[Chan_list [can_to]->basin] += S_fl; 
                            }


                                             
                        }
                        else if (cell_to > 0) {
                                /* OR change in nut mass in downstream internal cell (kg) */
 
                            if (structs->cell_i_to > 2 ) { /* to an on-map cell */
                                SWH[cell_to] += structs->flow/CELL_SIZE;
                                PA[cell_to] += P_fl; 
                                NA[cell_to] += N_fl; 
                                SA[cell_to] += S_fl; 
                                    /* mass balance in fr canal and to cell */
                                if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
                                    VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
                                    VOL_IN_STR[basn[cell_to]] += structs->flow; 
                                    P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
                                    P_IN_STR[basn[cell_to]] += P_fl; 
                                    S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
                                    S_IN_STR[basn[cell_to]] += S_fl; 

                                        /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
//           v2.8.4                         if ( (Chan_list[can_fr]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
//                                        sprintf( msgStr, "Day %6.1f: ERROR in basin configuration for structure from chan %d to cell (%d,%d): flow must be between diff hydrologic basins",
//                                                 SimTime.TIME, Chan_list[can_fr]->number,structs->cell_i_to,structs->cell_j_to );
//                                        WriteMsg( msgStr,True );
//                                        dynERRORnum++;
//                                    }
                                }
                                


                            }
                            else if (structs->cell_i_to == 2 && debug > -999 ) { /* to off-map cell, mass balance check */
                                VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow; 
                                VOL_OUT_STR[0] += structs->flow; 
                                P_OUT_STR[Chan_list[can_fr]->basin] += P_fl; 
                                P_OUT_STR[0] += P_fl; 
                                S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
                                S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
                            }
                        }
                            /* don't do much with zero or negative flows */           
                        if (structs->flow == 0.0) {
                            structs->conc_P = 0.0;
                            structs->conc_N = 0.0;
                            structs->conc_S = 0.0;
                        } /* for output purposes only, zero the reported conc */
                        else if (structs->flow < 0.0) {
                            sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from another cell/chan into chan %d",
                                     SimTime.TIME,Chan_list[can_fr]->number );
                            WriteMsg( msgStr,True );
                            dynERRORnum++;
                        }
                    }

                    structs = structs->next_in_list;
                } while (structs != NULL); /*  cycling to the end of the list again !!!! */

                structs = struct_hist_start; /* now pointing back to the first historical outflow structure that we dealt with */

               

} /* end of data-driven (historical/sfwmm) canal-canal flows */

Here is the call graph for this function:

void flowCalc_CanCan	(	struct Structure *	structs,
		float *	SWH,
		float *	GW,
		float *	poros,
		float *	Elev,
		double *	NA,
		double *	PA,
		double *	SA
	)

Canal-to-Canal water control structure flows: Calculated.

Parameters:

	SWH	SURFACE_WAT
	GW	SAT_WATER
	poros	HP_HYD_POROSITY
	Elev	SED_ELEV
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Remarks:

The following are allowable configurations:

• 1) Virtual structure: links two reaches within the same basin (which equilibrates the two separate reaches in a time step). Multiple virtual (or data-driven or rule-based) structures may add to or subtract from water to a single canal reach (but this is generally not a great practice - consider other reach-discretization options). Virtual structures will operate in the sequence found within the CanalData.struct database.

• 2) Rule-based structure: a limited case of a rule-based managed structure, with PUMPED flows between two canal reaches within the same basin. In this case, you are required to have HeadWater and TailWater Schedules that are checked in remote cells.

• Remember that all calculated-flow structures should be placed after all data-driven structures in the input database.
• For most objectives, it is preferable to give all virtual structures precedence over any rule-based structures (put earlier in sequence of the input database).

TODO - put in place the CanalData.struct reading checks to ensure the conditions described above are met.

Definition at line 3156 of file WatMgmt.c.

References Abs, Chan::area, Structure::canal_fr, Structure::canal_to, canstep, Structure::cell_HW, Structure::cell_TW, Chan_list, CHo, Structure::conc_N, Structure::conc_P, Structure::conc_S, debug, Chan::depth, dynERRORnum, Structure::flagHWsched, Structure::flagTWsched, Structure::flow, GP_MinCheck, HAB, Structure::HW_stage, Min, MINDEPTH, msgStr, Chan::N, Chan::number, Chan::P, P, ramp, Chan::S, Structure::S_nam, SimTime, Structure::str_cell_i, Structure::str_cell_j, Chan::sumHistIn, Chan::sumHistOut, Chan::sumRuleIn, Chan::sumRuleOut, T, simTime::TIME, True, Structure::TW_stage, usrErr(), Structure::w_coef, Chan::wat_depth, and WriteMsg().

Referenced by Flows_in_Structures().

{
    int can_to, can_fr;
    int cell_struct;
    float HeadH, HeadT, HeadH2, HeadT2;
    double CH_vol;
    double N_fl, P_fl, S_fl; 
    float HeadH_drop, HeadT_drop;

    float HWcheck, TWcheck, HWdiff, TWdiff; /* v2.8 */
    int struct_open;

    can_fr = structs->canal_fr - CHo;
    can_to = structs->canal_to - CHo;
    cell_struct = T(structs->str_cell_i,structs->str_cell_j); 

    /* VIRTUAL STRUCTURE FLOW */
                if ( structs->w_coef == 0.0)  {/* check for a virtual structure, flow only in pos direction */
                        HeadH_drop = 0.5 * Abs(Chan_list[can_fr]->elev_drop); /* Head_drops are midpoints of elevation difference between begin&end of canal reaches */
                        HeadT_drop = 0.5 * Abs(Chan_list[can_to]->elev_drop);
                /* v2.8 - added the sumRuleIn and sumRuleOut */
                                /* in both head and tail, add net historical/other-data flows before determining hydraulic potential */
                        HeadH = HeadH_drop + /* add portion of the elevation drop along the canal length */
                        Elev[cell_struct] - Chan_list[can_fr]->depth + Chan_list[can_fr]->wat_depth 
                        + (Chan_list[can_fr]->sumHistIn - Chan_list[can_fr]->sumHistOut)/Chan_list[can_fr]->area 
                        + (Chan_list[can_fr]->sumRuleIn - Chan_list[can_fr]->sumRuleOut)/Chan_list[can_fr]->area; 
                    /* OLD COMMENT: in tailwater only, check to see if other virtual struct has added water already, add to head */
                        HeadT = -HeadT_drop + /* subtract portion of the elevation drop along the canal length */
                        Elev [cell_struct] - Chan_list[can_to]->depth + Chan_list[can_to]->wat_depth
                        + (Chan_list[can_to]->sumHistIn - Chan_list[can_to]->sumHistOut)/Chan_list[can_to]->area
                        + (Chan_list[can_to]->sumRuleIn - Chan_list[can_to]->sumRuleOut)/Chan_list[can_to]->area;

                    if ( HeadH>HeadT ) {         
                        structs->flow =
                            Chan_list[can_fr]->area*Chan_list[can_to]->area/
                            (Chan_list[can_fr]->area+Chan_list[can_to]->area) * (HeadH-HeadT) ;
                        if (debug >3) {
                            HeadH2 = HeadH - structs->flow/Chan_list[can_fr]->area; 
                            HeadT2 = HeadT + structs->flow/Chan_list[can_to]->area;
                            sprintf( msgStr, "Day %6.1f:  virtual flow from chan %d (%f m) to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
                                     SimTime.TIME,Chan_list[can_fr]->number,HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth,
                                     Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
                            WriteMsg( msgStr,True ); 
                            if (Abs(HeadH2 - HeadT2) > GP_MinCheck) /* virtual flows should exactly equilibrate the two canals */
                            { sprintf( msgStr, "Day %6.1f: ERROR virtual flow from chan %d to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
                                       SimTime.TIME,Chan_list[can_fr]->number,Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
                            WriteMsg( msgStr,True ); dynERRORnum++; }

                        }
                        
                        CH_vol =  Chan_list[can_fr]->area* ramp((HeadH-HeadH_drop)-Elev[cell_struct]+Chan_list[can_fr]->depth-MINDEPTH); /* avail flow volume */
                        if (structs->flow > CH_vol) {
                            structs->flow =  CH_vol;
                            HeadH2 = HeadH - structs->flow/Chan_list[can_fr]->area; 
                            HeadT2 = HeadT + structs->flow/Chan_list[can_to]->area;
                            if (debug>3 && Abs(HeadH2 - HeadT2) > GP_MinCheck) {
                                sprintf( msgStr, "Day %6.1f: Warning: stage didn't equilibrate in modified virtual struct flow from chan %d (%f m) to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
                                         SimTime.TIME,Chan_list[can_fr]->number,CH_vol/Chan_list[can_fr]->area,Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
                                WriteMsg( msgStr,True ); }
                        }
                        
                            /* revise volume to be TOTAL volume for conc calc */
                        CH_vol = CH_vol + MINDEPTH*Chan_list[can_fr]->area;
                            /* conc in donor (fr) canal (kg/m3) */
                        structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
                        structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
                        structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
                    } 
                    
                    else  structs->flow = 0.0; /* when HeadH is not greater than HeadTailwater */
                    
                }
                

    /* RULE-BASED STRUCTURE FLOW */
                else  /* a non-zero w_coef indicates a flow calc based on management rules */
                {
                    if ( (structs->flagTWsched != 1) || (structs->flagHWsched != 1) ) {
                        sprintf (msgStr,"\nYou must assign both Headwater and Tailwater schedules for the rule-based canal-to-canal structure %s!\n", structs->S_nam); 
                        usrErr( msgStr );
                        exit(-1);
                        }
                    if ( (structs->cell_HW == 0) || (structs->cell_TW == 0) ) {
                        sprintf (msgStr,"\nYou must assign cells to check for both the Headwater and Tailwater schedules for the rule-based canal-to-canal structure %s!\n", structs->S_nam); 
                        usrErr( msgStr );
                        exit(-1);
                        }

                                struct_open = 0; /* the structure is closed, prior to making any checks of remote cell stages relative to their targets */

                                /* headwater and tailwater checks are not at the source and the destination of any potential flow, but are at remote cells:
                                   thus any flow from this schedule will not be head-tail based, but is considered 
                                   to be from a pump whose rate of pumping is dependent on the degree to which stage in the remote cells are approaching targets */
                                HWcheck = GW[structs->cell_HW]/poros[HAB[structs->cell_HW]] + SWH[structs->cell_HW] ;
                                TWcheck = GW[structs->cell_TW]/poros[HAB[structs->cell_TW]] + SWH[structs->cell_TW] ;
                                
                                
                                struct_open = ( (structs->TW_stage > TWcheck) && (structs->HW_stage < HWcheck) ) ? (1) : (0);

                    if  ( struct_open == 1)
                                {
                                                HWdiff =  HWcheck - structs->HW_stage;
                                                TWdiff =  structs->TW_stage - TWcheck;
                                        CH_vol =  Chan_list[can_fr]->area * ramp(Chan_list[can_fr]->wat_depth-MINDEPTH); /* avail flow volume */ 
                                structs->flow =  Min(structs->w_coef * ( HWdiff+TWdiff ) * canstep, CH_vol ); 
                        }
                    else structs->flow = 0.0;

                    
                } /* end of management rule based flow calcs */
 
                    /* don't do much with zero or negative rule-based flows  */           
                if (structs->flow == 0.0) {
                    structs->conc_P = 0.0;
                    structs->conc_N = 0.0;
                    structs->conc_S = 0.0;
                    return;
                } /* for output purposes only, zero the reported conc */
                else if (structs->flow < 0.0 && structs->w_coef != 0) {                 
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG RULE-BASED FLOW from chan %d",
                             SimTime.TIME,Chan_list[can_to]->number);
                    WriteMsg( msgStr,True );
                    return;
                } 

                            /* revise volume to be TOTAL volume for conc calc */
                        CH_vol = CH_vol + MINDEPTH*Chan_list[can_fr]->area;
                            /* conc in donor (fr) canal (kg/m3) */
                        structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
                        structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
                        structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
                                                                

                    /* increment the sum of rule-based inflows/outflows to canals
                       for potential multiple virtual inflows into some canals */
                Chan_list[can_fr]->sumRuleOut += structs->flow;         
                Chan_list[can_to]->sumRuleIn += structs->flow;          

                    /* nut flows */
                P_fl = structs->flow * structs->conc_P; 
                N_fl = structs->flow * structs->conc_N;
                S_fl = structs->flow * structs->conc_S;
                    /* change in nut mass in fr canal (kg) */
                Chan_list[can_fr]->N -=  N_fl ;
                Chan_list[can_fr]->P -=  P_fl ;
                Chan_list[can_fr]->S -=  S_fl ;
                    /* change in nut mass in to canal (kg) */
                Chan_list[can_to]->N +=  N_fl ;
                Chan_list[can_to]->P +=  P_fl ;
                Chan_list[can_to]->S +=  S_fl ;
                    /* mass balance in fr and to canals */
                    /* virtual structure flows can not (by definition) establish flow among basins */
                if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > 0 ) {
                    sprintf( msgStr, "Day %6.1f: ERROR - no hydrologic basin-basin flows in canal-to-canal structure, canals %d - %d",
                             SimTime.TIME, Chan_list[can_to]->number, Chan_list[can_fr]->number);
                    WriteMsg( msgStr,True ); dynERRORnum++;
                }
         
} /* end of calculated canal-canal flows */

Here is the call graph for this function:

void flowData_CelCel	(	struct Structure *	structs,
		float *	SWH,
		double *	NA,
		double *	PA,
		double *	SA
	)

Cell-to-Cell water control structure flows: Data-driven (historical/sfwmm).

Parameters:

	SWH	SURFACE_WAT
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Definition at line 3325 of file WatMgmt.c.

References HistData::arrayN, HistData::arrayP, HistData::arrayS, HistData::arrayWat, basn, basn_list, canDebugFile, canstep, Structure::cell_i_fr, Structure::cell_i_to, Structure::cell_j_fr, Structure::cell_j_to, CELL_SIZE, Structure::conc_N, Structure::conc_P, Structure::conc_S, debug, dynERRORnum, Structure::flow, Hist_data, Structure::histID, msgStr, P_IN_STR, P_OUT_STR, S_IN_STR, Structure::S_nam, S_OUT_STR, SimTime, T, simTime::TIME, Structure::TN, Structure::TNser, Structure::TP, Structure::TPser, True, Structure::TS, Structure::TSser, VOL_IN_STR, VOL_OUT_STR, and WriteMsg().

Referenced by Flows_in_Structures().

{ 
    int cell_to, cell_fr;

    double cell_vol;
    double N_fl, P_fl, S_fl; 

            cell_to = T(structs->cell_i_to,structs->cell_j_to); 
            cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);

    /* DATA-DRIVEN STRUCTURE FLOW */
                    /* don't do much with zero or negative flows */
                structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
                if (structs->flow == 0.0) {
                    structs->conc_P = 0.0;
                    structs->conc_N = 0.0;
                    structs->conc_S = 0.0;
                    return;    /* for output purposes only, zero the reported conc */
                }
                else if (structs->flow < 0.0) { /* all flows should be positive, but check */
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from cell (%d,%d)",
                             SimTime.TIME, structs->cell_i_to, structs->cell_j_to);
                    WriteMsg( msgStr,True );
                    return;
                }
                    /* for internal cells, check for excessive historical/model demand */
                if (structs->cell_i_fr != 2 ) 
                {

                    cell_vol = SWH[cell_fr] * CELL_SIZE;
                    if (structs->flow > cell_vol) { /* if hist/other-data demand > cell vol */
                        if( debug > 0 ) { 
                            fprintf( canDebugFile, "Day \t%6.1f\t %s:\tTotDemand= \t%7.0f\t; exceeded cell \t(%d,%d)\t vol= \t%7.0f \tby \t%7.0f\t m3.\n",
                                     SimTime.TIME, structs->S_nam, structs->flow, structs->cell_i_fr, structs->cell_j_fr, cell_vol, (structs->flow - cell_vol) );
                           
                        }
                            
                        structs->flow = cell_vol;
                    } 
                    }
                    /* CONCENTRATION */
                    /* now calculate the nutrient conc in the donor cell */
                    /* nut conc in flow from outside system using historical data or 0.0 */
                if (structs->cell_i_fr == 2) { 
                    structs->conc_P = (structs->TPser == 1)  ?
                        Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
                        ((structs->TPser == 0) ? (structs->TP) : (0.0) );
                    structs->conc_N = (structs->TNser == 1)  ?
                        Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
                        ((structs->TNser == 0) ? (structs->TN) : (0.0) );
                    structs->conc_S = (structs->TSser == 1)  ?
                        Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
                        ((structs->TSser == 0) ? (structs->TS) : (0.0) );
                }
                else
                { /* nut conc in flow within system - check to see if there are historical concs with the historical flow */
                    structs->conc_P = (structs->TPser == -1)  ?
                        ( (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
                        ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
                    structs->conc_N = ( structs->TNser == -1)  ?
                        ( (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
                        ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
                    structs->conc_S = ( structs->TSser == -1)  ?
                        ( (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
                        ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
                }
                    /* the nut and water update done after the rule-based flow calcs */
/* determine mass transfer */
                    
                P_fl = structs->flow * structs->conc_P;
                N_fl = structs->flow * structs->conc_N;
                S_fl = structs->flow * structs->conc_S;
                
                
                if (structs->cell_i_to == 2) { /* "to" cell is outside system */
                    if (debug > -999) {/* m3 and kg nuts outflow out from the system */
                        VOL_OUT_STR[basn[cell_fr]] += structs->flow; /* mass balance */
                        VOL_OUT_STR[0] += structs->flow; 
                        P_OUT_STR[basn[cell_fr]] += P_fl; 
                        P_OUT_STR[0] += P_fl; 
                        S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0); /* see below note on tracer application */
                        S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0); /* see below note on tracer application */
                    } 
                    PA[cell_fr] -= P_fl;
                    NA[cell_fr] -= N_fl; 
                    SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
                        /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
                           Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
                           that was not calculated from the available water and salt (i.e., the conc. was from the
                           CanalData.struct data file, either fixed or a time series) */
                        /* recalculate the surface water depth in the "from" cell */
                    SWH[cell_fr] -= structs->flow/CELL_SIZE;
                }
                else if (structs->cell_i_fr == 2) {  /* "from" cell is outside system */
                    if (debug > -999) {/* m3 and kg nuts inflow to the system */
                        VOL_IN_STR[basn[cell_to]] += structs->flow; 
                        VOL_IN_STR[0] += structs->flow; 
                        P_IN_STR[basn[cell_to]] += P_fl; 
                        P_IN_STR[0] += P_fl; 
                        S_IN_STR[basn[cell_to]] += S_fl; 
                        S_IN_STR[0] += S_fl; 
                    } 
                    PA[cell_to] += P_fl;
                    NA[cell_to] += N_fl;
                    SA[cell_to] += S_fl;
                        /* recalculate the surface water depth in the "to" cell */
                    SWH[cell_to] += structs->flow/CELL_SIZE;
                }
                else { /* internal cell-cell m3 and kg nuts flow */
                    PA[cell_fr] -= P_fl;
                    NA[cell_fr] -= N_fl;        
                    SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
                        /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
                           Thus, in that case, we don't subtract salt from the cell because this is "introduced" salt
                           that was not calculated from the available water and salt (i.e., the conc. was from the
                           CanalData.struct data file, either fixed or a time series) */
                    
                    PA[cell_to] += P_fl;
                    NA[cell_to] += N_fl;
                    SA[cell_to] += S_fl;
                        /* recalculate the surface water depth in the interacting cells */
                    SWH[cell_fr] -= structs->flow/CELL_SIZE;
                    SWH[cell_to] += structs->flow/CELL_SIZE;
                    if (basn[cell_to] != basn[cell_fr] && debug > -999 ) {
                        VOL_OUT_STR[basn[cell_fr]] += structs->flow; 
                        VOL_IN_STR[basn[cell_to]] += structs->flow; 
                        P_OUT_STR[basn[cell_fr]] += P_fl; 
                        P_IN_STR[basn[cell_to]] += P_fl; 
                        S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
                        S_IN_STR[basn[cell_to]] += S_fl; 

                            /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
                        if ( (basn_list[basn[cell_fr]]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
                            sprintf( msgStr, "Day %6.1f: ERROR in basin configuration for structure from cell (%d,%d) to cell (%d,%d): flow must be between diff hydrologic basins",
                                     SimTime.TIME, structs->cell_i_fr,structs->cell_j_fr,structs->cell_i_to,structs->cell_j_to );
                            WriteMsg( msgStr,True ); dynERRORnum++;
                        }
                    }
                    
                }
                return;
} /* end of data-driven (historical/sfwmm) cell-cell flows */

Here is the call graph for this function:

void flowCalc_CelCel	(	struct Structure *	structs,
		float *	SWH,
		float *	Elev,
		double *	NA,
		double *	PA,
		double *	SA
	)

Cell-to-Cell water control structure flows: Calculated (rule-based).

Parameters:

	SWH	SURFACE_WAT
	Elev	SED_ELEV
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Remarks:

This is set up only for manning's equation flow, specifically for allowing un-managed (no rules) marsh flows through a large opening in a levee (e.g., used for the multiple bridges along Alligator Alley highway in the regional ELM; use it elsewhere if you understand what it does).

• This is calculated structure flow can be either positive or negative.
• Both cells must be in same hydrologic basin.

Definition at line 3481 of file WatMgmt.c.

References Structure::cell_i_fr, Structure::cell_i_to, Structure::cell_j_fr, Structure::cell_j_to, CELL_SIZE, Structure::conc_P, Structure::conc_S, Structure::flow, Flux_SWcells(), Flux_SWstuff(), MCopen, SF_WT_VEL_mag, T, and Structure::w_coef.

Referenced by Flows_in_Structures().

{ 
    int cell_to, cell_fr;

    double cell_vol;
    double N_fl, P_fl, S_fl; 
    float FFlux; 
    extern float *SF_WT_VEL_mag; 

            cell_to = T(structs->cell_i_to,structs->cell_j_to); 
            cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);

    /* RULE-BASED STRUCTURE FLOW */
                if ( structs->w_coef != 0.0)  {/* check for a virtual structure  */
                    printf ("\nSorry - no rule-based cell-cell flows in this model version!\n"); exit(-1);
                }
                else {
    /* UN-RULEY, calculated STRUCTURE FLOW */
                        /* special case of calc'd flow, cell-cell, using the functions defined in Fluxes.c for overland flow.
                           This is used ONLY for grid-cell wide overland flows under/through bridges in this version */
                        /* MCopen is a (whole array) of a single Manning's coef that are initialized at an open-water value */
                    FFlux = Flux_SWcells(structs->cell_i_fr,structs->cell_j_fr,
                                      structs->cell_i_to,structs->cell_j_to,SWH,Elev,MCopen);
                        /* flow may be negative */
                                /* FFlux units = m */
                    structs->flow = FFlux*CELL_SIZE;
                    structs->conc_P = (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ;
/*                     structs->conc_N = (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ; */
                    structs->conc_S = (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ;

                    if (FFlux != 0.0) Flux_SWstuff ( structs->cell_i_fr,structs->cell_j_fr,
                                                  structs->cell_i_to,structs->cell_j_to,FFlux,SWH,SA,NA,PA,SF_WT_VEL_mag);
                    
                    
                }
} /* end of calculated cell-cell flows */

Here is the call graph for this function:

void flowData_CanCel	(	struct Structure *	structs,
		struct Structure *	struct_hist_start,
		float *	SWH,
		double *	NA,
		double *	PA,
		double *	SA
	)

Canal-to-Cell water control structure flows: Data-driven (historical/sfwmm).

Parameters:

	SWH	SURFACE_WAT
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Definition at line 3528 of file WatMgmt.c.

References Chan::area, HistData::arrayN, HistData::arrayP, HistData::arrayS, HistData::arrayWat, Chan::basin, basn, Structure::canal_fr, Structure::canal_to, canDebugFile, canstep, Structure::cell_i_to, Structure::cell_j_to, CELL_SIZE, Chan_list, CHo, Structure::conc_N, Structure::conc_P, Structure::conc_S, debug, dynERRORnum, Structure::flag, Structure::flow, Hist_data, Structure::histID, MINDEPTH, msgStr, Structure::multiOut, Chan::N, Structure::next_in_list, Chan::number, Chan::P, P_IN_STR, P_OUT_STR, ramp, Chan::S, S_IN_STR, Structure::S_nam, S_OUT_STR, SimTime, Structure::str_cell_i, Structure::str_cell_j, Chan::sumHistIn, Chan::sumHistOut, T, simTime::TIME, Structure::TN, Structure::TNser, Structure::TP, Structure::TPser, True, Structure::TS, Structure::TSser, VOL_IN_STR, VOL_OUT_STR, Chan::wat_depth, and WriteMsg().

Referenced by Flows_in_Structures().

{ 
    int cell_to, can_to, can_fr;
    int cell_struct;
    double CH_vol;
    double Nconc, Pconc, Sconc;
    double N_fl, P_fl, S_fl; 
    double chan_over;
    double ChanHistOut;

            can_fr = structs->canal_fr  - CHo;
            cell_to = T(structs->cell_i_to,structs->cell_j_to);           
            cell_struct = T(structs->str_cell_i,structs->str_cell_j);

 
    /* DATA-DRIVEN STRUCTURE FLOW */
                if (structs->multiOut == 1) return; /* this structure has already been processed as a multiple outflow from a canal */
                
                structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
                    /* don't do much with zero or negative flows */           
                if (structs->flow == 0.0) {/* zero the reported struct flow conc (for output info only) */
                    structs->conc_P = 0.0;
                    structs->conc_N = 0.0;
                    structs->conc_S = 0.0;
                    return; 
                }
                else if (structs->flow < 0.0) { 
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from cell (%d,%d)",
                             SimTime.TIME, structs->cell_i_to, structs->cell_j_to );
                    WriteMsg( msgStr,True ); dynERRORnum++;
                    return;
                }
                
                    /* for multiple historical/other-data outflows from a canal, ensure mass balance (sum of outflows !> avail volume) */
                    /* all flows from the canal associated with each set must be non-negative, which is currently valid */

                struct_hist_start = structs; /* this is the first structure w/ historical info for a particular canal */
                ChanHistOut = 0.0; /* sum of all historical outflows from a canal */

                while ( structs != NULL )   
                {  
                    if  ( structs->flag == 1 && struct_hist_start->canal_fr  == structs->canal_fr   ) { /* look for outflows from same canal */
                        structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
                        ChanHistOut += structs->flow; /* sum the multiple historical outflows from canal */
                    }
                        
                    structs = structs->next_in_list; /* increment to next structure */
                } /* end cycle thru list of remaining structures in total list */
                    
                structs = struct_hist_start; /* now we are pointing back to the first (or only) water control structure draining the canal */                 

                    /* concentration in donor canal for structure flow, g/L = kg/m^3 */
                    /* this is based on previous depth/volume calc'd in FluxChannel */
                    /* or use historical conc */
                        /* canal TOTAL volume before new structure flows */
                CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth; 
                    /* concentrations applied to all structures draining this canal */
                Pconc = (structs->TPser == -1) ? 
                    ( (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0) ) :
                    ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) ) ;
                Nconc = (structs->TNser == -1) ? 
                    ( (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0) ) :
                    ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) ) ;
                Sconc = (structs->TSser == -1) ? 
                    ( (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0) ) :
                    ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ) ;

                   /* vol avail for fluxing */
                CH_vol =  Chan_list[can_fr]->area*(ramp(Chan_list[can_fr]->wat_depth-MINDEPTH) );
                chan_over = (ChanHistOut > 0.0 ) ? ( CH_vol/ChanHistOut ) : 1.0;

                do /* loop again through the remaining sequence of structures to correct any outflows for the canal */
                {   
                    if ( structs->flag == 1 && struct_hist_start->canal_fr  == structs->canal_fr   ) { 

                        structs->multiOut = 1; /* flag to indicate this structure has been processed (below) */

                            /* revised FLOW */
                        if ( chan_over < 1.0 ) { /* if hist/other-data demand > canal vol */
                            if( debug > 0 ) { 
                                fprintf( canDebugFile, "Day \t%6.1f\t %s:\tTotDemand= \t%7.0f\t; exceeded chan \t%d\t avail vol= \t%7.0f\t by \t%7.0f\t m3 (\t%5.3f\t m deep).\n",
                                         SimTime.TIME, structs->S_nam, ChanHistOut, Chan_list[can_fr]->number,
                                         CH_vol, (ChanHistOut - CH_vol), Chan_list[can_fr]->wat_depth );
                                
                            }  
                                /* decrement the individual structure outflows by it's proportional contrib. */
                            structs->flow = (structs->flow) * chan_over ;
                        } 
                            /* sum of all realized (corrected for excess demand) historical/other-data outflows from a canal during an iteration */
                            /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
                        Chan_list[can_fr]->sumHistOut += structs->flow;         

                            /* the recipient changes as we cycle thru the structures (may be another cell or a canal) */
                        can_to = structs->canal_to - CHo;
                        cell_to = T(structs->cell_i_to,structs->cell_j_to); 

                        structs->conc_P = Pconc;/* concentrations applied to all structures draining this canal */
                        structs->conc_N = Nconc;
                        structs->conc_S = (structs->TSser == -1) ?
                            (Sconc) :
                            ( (structs->TSser == 0) ? (structs->TS) :
                            ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ); /* see note below on tracer conc. application */
                         
                        P_fl = structs->flow * structs->conc_P; 
                        N_fl = structs->flow * structs->conc_N;
                        S_fl = structs->flow * structs->conc_S;
                                /* change in nut mass in fr canal (kg) */
                        Chan_list[can_fr]->N -=  N_fl ;
                        Chan_list[can_fr]->P -=  P_fl ;
                        Chan_list[can_fr]->S -=  (structs->TSser == -1) ? (S_fl) : (0);
                            /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
                               Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
                               that was not calculated from the available water and salt (i.e., the conc. was from the
                               CanalData.struct data file, either fixed or a time series) */
  
                        if (can_to > 0) {    /* IF this is a canal recip, change in nut mass in to canal (kg) */
                            /* sum of all realized (corrected for excess demand) historical/other-data outflows into a canal during an iteration */
                            /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
                                Chan_list[can_to]->sumHistIn += structs->flow;          

                            Chan_list[can_to]->N +=  N_fl ;
                            Chan_list[can_to]->P +=  P_fl ;
                            Chan_list[can_to]->S +=  S_fl ;
                                /* mass balance in fr and to canals */
                            if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > -999 ) {
                                VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
                                VOL_IN_STR[Chan_list [can_to]->basin] += structs->flow; 
                                P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
                                P_IN_STR[Chan_list [can_to]->basin] += P_fl; 
                                S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
                                S_IN_STR[Chan_list [can_to]->basin] += S_fl; 
                            }
                        }
                        else if (cell_to > 0) {
                                /* OR change in nut mass in downstream cell (kg) */
 
                            if (structs->cell_i_to > 2 ) { /* to an on-map cell */
                                SWH[cell_to] += structs->flow/CELL_SIZE;
                                PA[cell_to] += P_fl; 
                                NA[cell_to] += N_fl; 
                                SA[cell_to] += S_fl; 
                                    /* mass balance in fr canal and to cell */
                                if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
                                    VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
                                    VOL_IN_STR[basn[cell_to]] += structs->flow; 
                                    P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
                                    P_IN_STR[basn[cell_to]] += P_fl; 
                                    S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
                                    S_IN_STR[basn[cell_to]] += S_fl; 

                                        /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
//   v2.8.4                                 if ( (Chan_list[can_fr]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
//                                        sprintf( msgStr, "Day %6.1f: ERROR in basin configuration for structure from chan %d to cell (%d,%d): flow must be between diff hydrologic basins",
//                                                 SimTime.TIME, Chan_list[can_fr]->number,structs->cell_i_to,structs->cell_j_to );
//                                        WriteMsg( msgStr,True ); dynERRORnum++;
//                                    }
                                }
                            }
                            else if (structs->cell_i_to == 2 && debug > -999 ) { /* to off-map cell, mass balance check */
                                VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow; 
                                VOL_OUT_STR[0] += structs->flow; 
                                P_OUT_STR[Chan_list[can_fr]->basin] += P_fl; 
                                P_OUT_STR[0] += P_fl; 
                                S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
                                S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
                            }
                        }
                            /* don't do much with zero or negative flows */           
                        if (structs->flow == 0.0) {
                            structs->conc_P = 0.0;
                            structs->conc_N = 0.0;
                            structs->conc_S = 0.0;
                        } /* for output purposes only, zero the reported conc */
                        else if (structs->flow < 0.0) {
                            sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from another chan/cell into chan %d",
                                     SimTime.TIME, Chan_list[can_fr]->number );
                            WriteMsg( msgStr,True ); dynERRORnum++;
                        }
                    }

                    structs = structs->next_in_list;
                } while (structs != NULL); /*  cycling to the end of the list again !!!! */

                structs = struct_hist_start; /* now pointing back to the first historical outflow structure that we dealt with */
    

                return;
                    
} /* end of data-driven (historical/sfwmm) canal-cell flows */

Here is the call graph for this function:

void flowCalc_CanCel	(	struct Structure *	structs,
		float *	SWH,
		float *	GW,
		float *	poros,
		float *	Elev,
		double *	NA,
		double *	PA,
		double *	SA
	)

Canal-to-Cell water control structure flows: Calculated (rule-based).

Parameters:

	SWH	SURFACE_WAT
	GW	SAT_WATER
	poros	HP_HYD_POROSITY
	Elev	SED_ELEV
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Remarks:

The following are allowable configurations, set up for (positive) LOSSES FROM THE MODEL DOMAIN only:

• 1) Recipient cell is off-map and no (head or tail) schedule target available: stage-based calculation, using other-model boundary condition for external tailwater stage estimate (must have provided a cell location in input dbase for this off-map cell) In this case, the difference between donor canal stage and recipient (external) cell stage is used to deterimine if any (positive-only) flow occurs.

• 2) Recipient cell is off-map and a TAILWATER SCHEDULE TARGET is available: stage-based calculation, using tailwater schedule for tailwater stage estimate. In this case, the difference between donor canal stage and recipient (external) cell stage is used to deterimine if any (positive-only) flow occurs.

• 3) Recipient cell is off-map and a HEADWATER SCHEDULE TARGET is available: PUMPED outflow, using the schedule for determining magnitude of difference between target and stage in remote cell. The head water target check is made at a remote on-map grid cell.

• NOTE that you MAY NOT have both head and tail water schedules for this canal-cell rule-based structure.

• Remember that all calculated-flow structures should be placed after all data-driven structures in the input database.
• For most objectives, it is preferable to give all (canal-canal) virtual structures precedence over any rule-based structures (put earlier in sequence of the input database).

TODO - put in place the CanalData.struct reading checks to ensure the conditions described above are met.

Definition at line 3747 of file WatMgmt.c.

References Chan::area, Chan::basin, basn, boundcond_depth, Structure::canal_fr, canstep, Structure::cell_HW, Structure::cell_i_to, Structure::cell_i_TW, Structure::cell_j_to, Structure::cell_j_TW, CELL_SIZE, Chan_list, CHo, Structure::conc_N, Structure::conc_P, Structure::conc_S, debug, Chan::depth, dynERRORnum, Chan::elev_avg, Structure::flagHWsched, Structure::flagTWsched, Structure::flow, HAB, Structure::HW_stage, Min, MINDEPTH, msgStr, Chan::N, Chan::P, P, P_IN_STR, P_OUT_STR, ramp, Chan::S, S_IN_STR, Structure::S_nam, S_OUT_STR, SimTime, Structure::str_cell_i, Structure::str_cell_j, Chan::sumHistIn, Chan::sumHistOut, Chan::sumRuleOut, T, simTime::TIME, True, Structure::TW_stage, usrErr(), VOL_IN_STR, VOL_OUT_STR, Structure::w_coef, Chan::wat_depth, and WriteMsg().

Referenced by Flows_in_Structures().

{ 
    int cell_to, can_fr;
    int cell_struct;
    float HeadH, HeadT;
    double CH_vol;
    double N_fl, P_fl, S_fl; 
    float flow_rev; 

    /* from unitmod.h */
    extern float* boundcond_depth; /* depth in external cells */
    
    float HWcheck, TWcheck, Hdiff; /* v2.8 */
    int struct_open;

    /* RULE-BASED STRUCTURE FLOW */
            can_fr = structs->canal_fr  - CHo;
            cell_to = T(structs->cell_i_to,structs->cell_j_to);           
            cell_struct = T(structs->str_cell_i,structs->str_cell_j);

                if (structs->cell_i_to != 2  ) {
                    sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow is restricted to outflow from the model domain.",
                             structs->S_nam );
                    usrErr( msgStr );   exit(-1);    
                }
                if (structs->flagTWsched == 1 && structs->flagHWsched == 1 ) {
                    sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow can not have both HeadWater and TailWater Schedules.",
                             structs->S_nam );
                    usrErr( msgStr);   exit(-1);    
                }
                if ( structs->flagHWsched != 1 && structs->flagTWsched != 1 && structs->cell_i_TW == 0  ) {
                    sprintf( msgStr, "ERROR in flow configuration for %s structure : without a head or tail water schedule, this rule-based flow must have a tailwater-check-cell defined (to acquire data from the boundary-condition model).",
                             structs->S_nam );
                    usrErr( msgStr );   exit(-1);    
                }
                                  /* v2.8 - added the sumRuleIn and sumRuleOut */
                                  /* v2.8 - change HeadH from "Elev [cell_struct] -" ... to "Chan_list[can_fr]->elev_avg -" ... */

                HeadH = Chan_list[can_fr]->elev_avg - Chan_list[can_fr]->depth + Chan_list[can_fr]->wat_depth
                                + (Chan_list[can_fr]->sumHistIn - Chan_list[can_fr]->sumHistOut)/Chan_list[can_fr]->area 
                                + (Chan_list[can_fr]->sumRuleIn - Chan_list[can_fr]->sumRuleOut)/Chan_list[can_fr]->area; 
                                        
                HeadT =  (structs->flagTWsched == 1) ? (structs->TW_stage) : ( boundcond_depth[T(structs->cell_i_TW,structs->cell_j_TW)]  + Elev[cell_struct]);
                
                Hdiff = HeadH - HeadT;
                                /* volume available for fluxing -  leave MINDEPTH water; using the head estimate after previously calc'd losses is */
                CH_vol =  Chan_list[can_fr]->area* ramp(HeadH-Chan_list[can_fr]->elev_avg+Chan_list[can_fr]->depth-MINDEPTH); 
            
                        struct_open = 1; /* the structure is open, prior to making any checks of remote cell stages relative to their target schedule */

                        if (structs->flagHWsched == 1 && structs->cell_HW != 0) { 
                                /* headwater check is not at the source of any potential flow, but is at a remote cell:
                                   thus any flow from this schedule will not be head-tail based, but is considered 
                                   to be from a pump whose variable rate of pumping is dependent on the degree to which stage in the remote cell is approaching target */
                                HWcheck = GW[structs->cell_HW]/poros[HAB[structs->cell_HW]] + SWH[structs->cell_HW] ;
                                struct_open = (structs->HW_stage < HWcheck) ? (1) : (0);
                                Hdiff =  HWcheck - structs->HW_stage ;  
                        }
                        

                if  ( (Hdiff > 0.0) && struct_open == 1)
                {
                    structs->flow =  Min(structs->w_coef * ( Hdiff ) * canstep, CH_vol );   
                    
                        /* all rule based instances of canal->cell flow  go
                           to external cells, so this reversal check should not be reached */
                    flow_rev = ( structs->cell_i_to != 2 ) ? 
                        (HeadH - structs->flow/Chan_list[can_fr]->area) - (HeadT + structs->flow/CELL_SIZE) :
                        (0.0);
                        /* the flux eqn is the same as the weir eqn for virtual structures, equilibrating the heads across two unequal area regions */
                    /* flow is only in positive direction */
                    if (flow_rev < 0.0) {  structs->flow = 
                        ramp( Chan_list[can_fr]->area*CELL_SIZE*HeadH/(Chan_list[can_fr]->area+CELL_SIZE) -
                        Chan_list[can_fr]->area*CELL_SIZE*HeadT/(Chan_list[can_fr]->area+CELL_SIZE) );
                          }
               }
                else structs->flow = 0.0; 
                
            
            
                    /* don't do much with zero or negative flows */           
                if (structs->flow == 0.0) {
                    structs->conc_P = 0.0;
                    structs->conc_N = 0.0;
                    structs->conc_S = 0.0;
                    return;
                } /* for output purposes only, zero the reported conc */
                else if (structs->flow < 0) {  
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG FLOW from cell (%d,%d)",
                             SimTime.TIME, structs->cell_i_to, structs->cell_j_to );
                    WriteMsg( msgStr,True ); dynERRORnum++;
                    return;
                }
          
                    /* rule-based flows:  calculate the nutrient concentrations */
                    /* conc in canal -> cell flow  */
                CH_vol =  Chan_list[can_fr]->area* ramp(HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth); /* TOTAL volume */
                    /* conc in donor (fr) canal (kg/m3) */
                structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
                structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
                structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);

                P_fl = structs->flow * structs->conc_P; 
                N_fl = structs->flow * structs->conc_N;
                S_fl = structs->flow * structs->conc_S;
                  
                    /* increment the sum of rule-based flows from canal
                       for potential multiple virtual structure flows associated with some canals */
                Chan_list[can_fr]->sumRuleOut += structs->flow;         
                
                    /* change in nut mass in downstream cell (kg) */
                if (structs->cell_i_to > 2 ) { /* to an on-map cell */
                    SWH[cell_to] += structs->flow/CELL_SIZE;
                    PA[cell_to] += P_fl; 
                    NA[cell_to] += N_fl; 
                    SA[cell_to] += S_fl; 
                        /* mass balance in fr canal and to cell */
                    if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
                        VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
                        VOL_IN_STR[basn[cell_to]] += structs->flow; 
                        P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
                        P_IN_STR[basn[cell_to]] += P_fl; 
                        S_OUT_STR[Chan_list[can_fr]->basin] += S_fl;
                        S_IN_STR[basn[cell_to]] += S_fl; 

                            /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
//     v2.8.4                   if ( (Chan_list[can_fr]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
//                            sprintf( msgStr, "Day %6.1f: ERROR in basin configuration for structure from chan %d to cell (%d,%d): flow must be between diff hydrologic basins",
//                                     SimTime.TIME, Chan_list[can_fr]->number,structs->cell_i_to,structs->cell_j_to );
//                            WriteMsg( msgStr,True ); dynERRORnum++;
//                        }
                    }
                }
                else if (structs->cell_i_to == 2 && debug > -999 ) { /* to off-map cell, mass balance check */
                    VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow; 
                    VOL_OUT_STR[0] += structs->flow; 
                    P_OUT_STR[Chan_list[can_fr]->basin] += P_fl; 
                    P_OUT_STR[0] += P_fl; 
                    S_OUT_STR[Chan_list[can_fr]->basin] += S_fl; 
                    S_OUT_STR[0] += S_fl; 
                }

            
                    /* calc the change in mass in  canal (kg) */
                Chan_list[can_fr]->N -=  N_fl ;
                Chan_list[can_fr]->P -=  P_fl ;
                Chan_list[can_fr]->S -=  S_fl ;
} /* end of calculated canal-cell flows */

Here is the call graph for this function:

void flowData_CelCan	(	struct Structure *	structs,
		float *	SWH,
		double *	NA,
		double *	PA,
		double *	SA
	)

Cell-to-Canal water control structure flows: Data-driven (historical/sfwmm).

Parameters:

	SWH	SURFACE_WAT
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Definition at line 3907 of file WatMgmt.c.

References HistData::arrayN, HistData::arrayP, HistData::arrayS, HistData::arrayWat, Chan::basin, basn, basn_list, Structure::canal_to, canDebugFile, canstep, Structure::cell_i_fr, Structure::cell_j_fr, CELL_SIZE, Chan_list, CHo, Structure::conc_N, Structure::conc_P, Structure::conc_S, debug, dynERRORnum, Structure::flow, Hist_data, Structure::histID, msgStr, Structure::multiOut, Chan::N, Chan::number, Chan::P, P_IN_STR, P_OUT_STR, Chan::S, S_IN_STR, Structure::S_nam, S_OUT_STR, SimTime, Structure::str_cell_i, Structure::str_cell_j, Chan::sumHistIn, T, simTime::TIME, Structure::TN, Structure::TNser, Structure::TP, Structure::TPser, True, Structure::TS, Structure::TSser, VOL_IN_STR, VOL_OUT_STR, and WriteMsg().

Referenced by Flows_in_Structures().

{ 
    int cell_fr, can_to;
    int cell_struct;
    double cell_vol;
    double N_fl, P_fl, S_fl; 

            can_to = structs->canal_to - CHo;
            cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);           
            cell_struct = T(structs->str_cell_i,structs->str_cell_j);

    /* DATA-DRIVEN STRUCTURE FLOW */
  
                if (structs->multiOut == 1) {
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
                             SimTime.TIME, Chan_list[can_to]->number );
                    WriteMsg( msgStr,True ); dynERRORnum++;
                    return; /* this structure has already been processed as a multiple outflow from a canal */
                }
                
                structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
                    /* don't do much with zero or negative flows */           
                if (structs->flow == 0.0) {
                    structs->conc_P = 0.0;
                    structs->conc_N = 0.0;
                    structs->conc_S = 0.0;
                    return; 
                }
                else if (structs->flow <0) { 
                    sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
                             SimTime.TIME, Chan_list[can_to]->number );
                    WriteMsg( msgStr,True ); dynERRORnum++;
                    return;
                } 

                if (structs->cell_i_fr != 2 ) /* cell (2,2) == (1,1) in maps, not in model domain */
                {

                    cell_vol = SWH[cell_fr] * CELL_SIZE;
                    if (structs->flow > cell_vol) {
                        if( debug > 0 ) { 
                            fprintf( canDebugFile, "Day \t%6.1f\t %s:\tTotDemand= \t%7.0f\t; exceeded cell \t(%d,%d)\t vol= \t%7.0f\t by \t%7.0f\t m3.\n",
                                     SimTime.TIME, structs->S_nam, structs->flow, structs->cell_i_fr, structs->cell_j_fr, cell_vol, (structs->flow - cell_vol) );
                            
                        }
                        structs->flow = cell_vol;
                    }}
                    /* calculate the nutrient conc in the donor cell */
                    /* nut conc in flow from outside system using historical data or 0.0 */
                if (structs->cell_i_fr == 2) { 
                    structs->conc_P = (structs->TPser == 1)  ?
                        Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
                        ((structs->TPser == 0) ? (structs->TP) : (0.0) );
                    structs->conc_N = (structs->TNser == 1)  ?
                        Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
                        ((structs->TNser == 0) ? (structs->TN) : (0.0) );
                    structs->conc_S = (structs->TSser == 1)  ?
                        Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
                        ((structs->TSser == 0) ? (structs->TS) : (0.0) );
                }
                else
                { /* nut conc in flow within system - check to see if there are historical concs with the historical flow */
                    structs->conc_P = (structs->TPser == -1)  ?
                        ( (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
                        ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
                    structs->conc_N = ( structs->TNser == -1)  ?
                        ( (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
                        ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
                    structs->conc_S = ( structs->TSser == -1)  ?
                        ( (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
                        ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
                }
                    /* now determine the nutrient mass (kg) flow from cell to canal */
                    /* 1000 is to convert g/m3 to kg */ 
                P_fl = structs->flow * structs->conc_P;
                N_fl = structs->flow * structs->conc_N;
                S_fl = structs->flow * structs->conc_S;
                if (structs->cell_i_fr == 2 ) {/* "from" cell is outside system */
                    if (debug > -999) {  
                        VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; /* mass balance */
                        VOL_IN_STR[0] += structs->flow; 
                        P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
                        P_IN_STR[0] += P_fl; 
                        S_IN_STR[Chan_list[can_to]->basin] += S_fl; 
                        S_IN_STR[0] += S_fl; 
                    } }
                else { /* internal cell-> canal flow */
                    SWH[cell_fr] -= structs->flow/CELL_SIZE;
                    PA[cell_fr] -= P_fl; 
                    NA[cell_fr] -= N_fl; 
                    SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
                            /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
                               Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
                               that was not calculated from the available water and salt (i.e., the conc. was from the
                               CanalData.struct data file, either fixed or a time series) */
                        /* mass balance in fr cell and to canal */
                    if (basn[cell_fr] != Chan_list[can_to]->basin && debug > -999 ) {
                        VOL_OUT_STR[basn[cell_fr]] += structs->flow;
                        VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; 
                        P_OUT_STR[basn[cell_fr]] += P_fl;
                        P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
                        S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
                        S_IN_STR[Chan_list[can_to]->basin] += S_fl; 

                            /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
                        if ( (Chan_list[can_to]->family ==  basn_list[basn[cell_fr]]->family) && (debug > 0) ) {
                            sprintf( msgStr, "Day %6.1f: ERROR in basin configuration for structure from cell (%d,%d) to chan %d: flow must be between diff hydrologic basins",
                                     SimTime.TIME, structs->cell_i_fr, structs->cell_j_fr, Chan_list[can_to]->number );
                            WriteMsg( msgStr,True ); dynERRORnum++;
                        }
                    }
                }
                
                  /* sum of all realized (corrected for excess demand) historical/other-data outflows into a canal during an iteration */
                  /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
                Chan_list[can_to]->sumHistIn += structs->flow;          

                   /* change in nut mass in "to" canal (kg) */
                Chan_list[can_to]->N +=  N_fl ;
                Chan_list[can_to]->P +=  P_fl ;
                Chan_list[can_to]->S +=  S_fl ;
                return;

} /* end of data-driven (historical/sfwmm) cell-canal flows */

Here is the call graph for this function:

void flowCalc_CelCan	(	struct Structure *	structs,
		float *	SWH,
		float *	GW,
		float *	poros,
		float *	Elev,
		double *	NA,
		double *	PA,
		double *	SA
	)

Cell-to-Canal water water control structure flows: Calculated (rule-based).

Parameters:

	SWH	SURFACE_WAT
	GW	SAT_WATER
	poros	HP_HYD_POROSITY
	Elev	SED_ELEV
	NA	DINdummy
	PA	TP_SF_WT
	SA	SALT_SURF_WT

Remarks:

The following are allowable configurations, set up for (positive) INPUTS TO THE MODEL DOMAIN only:

• 1) Donor cell is off-map and no (head or tail) schedule target available: stage-based calculation, using other-model boundary condition for external headwater stage estimate (must have provided a cell location in input dbase for this off-map cell) In this case, the difference between (external) donor cell stage and recipient canal stage is used to deterimine if any (positive-only) flow occurs. THIS HAS NEVER BEEN USED, as of ELM v2.8.

• 2) Donor cell is off-map and a HEADWATER SCHEDULE TARGET is available: stage-based calculation, using headwater schedule for headwater stage estimate. In this case, the difference between (external) donor cell stage and recipient canal stage is used to deterimine if any (positive-only) flow occurs.

• 3) Donor cell is off-map and a TAILWATER SCHEDULE TARGET is available: PUMPED inflow, using the schedule for determining magnitude of difference between target and stage in remote cell. The tail water target check is made at a remote on-map grid cell.

• NOTE that you MAY NOT have both head and tail water schedules for this cell-canl rule-based structure.

• Remember that all calculated-flow structures should be placed after all data-driven structures in the input database.
• For most objectives, it is preferable to give all (canal-canal) virtual structures precedence over any rule-based structures (put earlier in sequence of the input database).

TODO - put in place the CanalData.struct reading checks to ensure the conditions described above are met.

Definition at line 4061 of file WatMgmt.c.

References Chan::area, HistData::arrayN, HistData::arrayP, HistData::arrayS, Chan::basin, boundcond_depth, Structure::canal_to, canstep, Structure::cell_i_fr, Structure::cell_i_HW, Structure::cell_j_HW, Structure::cell_TW, Chan_list, CHo, Structure::conc_N, Structure::conc_P, Structure::conc_S, Chan::depth, dynERRORnum, Chan::elev_avg, Structure::flagHWsched, Structure::flagTWsched, Structure::flow, HAB, Hist_data, Structure::histID, Structure::HW_stage, Max, msgStr, Chan::N, Chan::number, Chan::P, P_IN_STR, Chan::S, S_IN_STR, Structure::S_nam, SimTime, Chan::sumHistIn, Chan::sumHistOut, Chan::sumRuleIn, T, simTime::TIME, Structure::TN, Structure::TNser, Structure::TP, Structure::TPser, True, Structure::TS, Structure::TSser, Structure::TW_stage, usrErr(), VOL_IN_STR, Structure::w_coef, Chan::wat_depth, and WriteMsg().

Referenced by Flows_in_Structures().

{ 
    int cell_fr, can_to;
    int cell_struct;
    float HeadH, HeadT;
    double N_fl, P_fl, S_fl; 

    /* from unitmod.h */
    extern float* boundcond_depth; /* depth in external cells */

    float HWcheck, TWcheck, Hdiff; /* v2.8 */
    int struct_open;

    /* RULE-BASED STRUCTURE FLOW */
            can_to = structs->canal_to - CHo;

                if (structs->cell_i_fr != 2  ) {
                    sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow is restricted to inflow from outside the model domain.",
                             structs->S_nam );
                    usrErr( msgStr);   exit(-1);    
                }
                if (structs->flagTWsched == 1 && structs->flagHWsched == 1 ) {
                    sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow can not have both HeadWater and TailWater Schedules.",
                             structs->S_nam );
                    usrErr( msgStr);   exit(-1);    
                }
                if ( structs->flagTWsched != 1 && structs->flagHWsched != 1 && structs->cell_i_HW == 0  ) {
                    sprintf( msgStr, "ERROR in flow configuration for %s structure : without a head or tail water schedule, this rule-based flow must have a headwater-check-cell defined (to acquire data from the boundary-condition model).",
                             structs->S_nam );
                    usrErr( msgStr );   exit(-1);    
                }
                
                HeadH = (structs->flagHWsched == 1) ? (structs->HW_stage) : ( boundcond_depth[T(structs->cell_i_HW,structs->cell_j_HW)]  + Elev[cell_struct]); 
          /* in v2.3, the rule-based cell->canal flows use only tide data for boundary conditions */
          /* v2.8 - add the sumRuleIn and sumRuleOut */
                HeadT = Chan_list[can_to]->elev_avg - Chan_list[can_to]->depth + Chan_list[can_to]->wat_depth
                    + (Chan_list[can_to]->sumHistIn - Chan_list[can_to]->sumHistOut)/Chan_list[can_to]->area
                    + (Chan_list[can_to]->sumRuleIn - Chan_list[can_to]->sumRuleOut)/Chan_list[can_to]->area;
            
                        Hdiff = HeadH - HeadT;
                        
                        struct_open = 1; /* the structure is open, prior to making any checks of remote cell stages relative to their targets */
                        
                        if (structs->flagTWsched == 1 && structs->cell_TW != 0) { 
                                /* tailwater check is not at the destination of any potential flow, but is at a remote cell:
                                   thus any flow from this schedule will not be head-tail based, but is considered 
                                   to be from a pump whose rate of pumping is dependent on the degree to which stage in the remote cell is approaching target */
                                TWcheck = GW[structs->cell_TW]/poros[HAB[structs->cell_TW]] + SWH[structs->cell_TW] ;
                                struct_open = (structs->TW_stage > TWcheck) ? (1) : (0);
                                Hdiff =  structs->TW_stage - TWcheck;  
                        }
                        

                if  ( struct_open == 1)
                {
                    structs->flow =  Max(structs->w_coef * ( Hdiff ) * canstep, 0 );     
//                    printf ("  On %s, %d %d, TWcheck= %f, target= %f, struct is %d, Hdiff=%f, flow=%f \n",structs->S_nam,structs->cell_i_TW, structs->cell_j_TW, TWcheck,structs->TW_stage, struct_open,Hdiff,structs->flow);
                }
                else structs->flow = 0.0; /* undefined flow if heads are equal, thus this stmt added */
                
                /* don't do much with zero or negative flows */           
            if (structs->flow == 0.0) {
                structs->conc_P = 0.0;
                structs->conc_N = 0.0;
                structs->conc_S = 0.0;
                return;
            } /* for output purposes only, zero the reported conc */

                /*  rule-based: calculate the nutrients - this is restricted to calculated tidal boundary condition flows  from outside of system */
                /* concentration in structure flows, g/L = kg/m^3; mass of stock in kg */

                if (structs->cell_i_fr == 2) { 
                    structs->conc_P = (structs->TPser == 1)  ?
                        Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
                        ((structs->TPser == 0) ? (structs->TP) : (0.0) );
                    structs->conc_N = (structs->TNser == 1)  ?
                        Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
                        ((structs->TNser == 0) ? (structs->TN) : (0.0) );
                    structs->conc_S = (structs->TSser == 1)  ?
                        Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
                        ((structs->TSser == 0) ? (structs->TS) : (0.0) );
                }
            P_fl = structs->flow * structs->conc_P; 
            N_fl = structs->flow * structs->conc_N;
            S_fl = structs->flow * structs->conc_S;
             
                    /* increment the sum of rule-based flows to canal
                       for potential multiple virtual structure flows associated with some canals */
            Chan_list[can_to]->sumRuleIn += structs->flow;              

                /* this calculated flow is restricted to an outside-system cell, which does not need updating */
                /* these are nutrient stocks, not concentrations */  
                /* mass balance in fr cell and to canal */
            if (structs->cell_i_fr == 2  ) {
                VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; 
                VOL_IN_STR[0] += structs->flow; 
                P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
                P_IN_STR[0] += P_fl; 
                S_IN_STR[Chan_list[can_to]->basin] += S_fl; 
                S_IN_STR[0] += S_fl; 
            }
            else        {  
                sprintf( msgStr, "Day %6.1f: ERROR in flow to chan %d: cannot have calculated cell->canal flows except from outside of model domain. ",
                         SimTime.TIME, Chan_list[can_to]->number);
                WriteMsg( msgStr,True ); dynERRORnum++;
                return; 
            }
            

                /* change in nut mass in canal (kg) */
            Chan_list[can_to]->N +=  N_fl ;
            Chan_list[can_to]->P +=  P_fl ;
            Chan_list[can_to]->S +=  S_fl ;
            
} /* end of calculated cell-canal flows */

Here is the call graph for this function:

float GetGraph	(	struct Schedule *	graph,
		float	x
	)

Graph to read the time dependent schedules.

Parameters:

	graph	A struct of Schedule
	x	Time, julian day within recurring year

Returns:

value The target stage (m)

Definition at line 4190 of file WatMgmt.c.

References Schedule::graph_points, Schedule::num_point, Points::time, and Points::value.

Referenced by Flows_in_Structures().

{
    float s;
    int ig = 0, Np;
        
    Np = graph->num_point;

    while(1) 
    {
        if (x <= graph->graph_points[ig].time) 
        { if ( ig > 0 ) ig--; 
        else return ( graph->graph_points[0].value );
        }
        else if (x > graph->graph_points[ig].time && x <= graph->graph_points[ig+1].time) 
        {
            s =         ( graph->graph_points[ig+1].value - graph->graph_points[ig].value )/
                ( graph->graph_points[ig+1].time - graph->graph_points[ig].time );
            return ( s * ( x - graph->graph_points[ig].time ) + graph->graph_points[ig].value ); 
        }
        else if (x > graph->graph_points[ig+1].time) 
        { if ( ig < Np ) ig++; 
        else return ( graph->graph_points[Np].value );
        }
    }
}

float UTM2kmy

(

double

UTM

)

Converter from meters in UTM to grid cell column location (zero origin).

Parameters:

UTM

Geographic coordinate (m, ELM uses UTM zone17, NAD'27, but not important here)

Returns:

Grid cell column number

Definition at line 4223 of file WatMgmt.c.

References Abs, celWid, and UTM_EOrigin.

Referenced by ReadChanStruct().

{ return ( Abs(UTM - UTM_EOrigin)/celWid );  
}

float UTM2kmx

(

double

UTM

)

Converter from meters in UTM to grid cell row location (zero origin).

Parameters:

UTM

Geographic coordinate (m, ELM uses UTM zone17, NAD'27, but not important here)

Returns:

Grid cell row number

Definition at line 4232 of file WatMgmt.c.

References Abs, celWid, and UTM_NOrigin.

Referenced by ReadChanStruct().

{ return ( Abs(UTM - UTM_NOrigin)/celWid );              
}

int Wrdcmp	(	char *	s,
		char *	t
	)

Compare two words.

Parameters:

	s	A string
	t	Another string

Returns:

True/false

Definition at line 4243 of file WatMgmt.c.

References EOS.

Referenced by Read_schedule().

{ 
  for ( ; *s == *t; s++, t++ )
         if ( isspace (*(s+1)) || *(s+1) == EOS ) return 1;
        
  return 0;
  
}