Everglades Landscape Model (ELM) Home Page
Main Page | Data Structures | Directories | File List | Data Fields | Globals

WatMgmt.c

Go to the documentation of this file.
00001 
00030 /* General NOTES on revisions to this source file:
00031         april 00 VERSION 2.1 - output available for application for water quality performance
00032         july 02 VERSIOIN 2.1a - first widely-available public release of code -
00033               added some more code documentation, etc.
00034         July 03 VERSION 2.2.0 - moved declarations into header file, other similar changes
00035               to help reviewers better understand organization of codes.
00036               -Some minor code changes that are found tagged with the comment "v2.2" 
00037               -verified that "re-organized" v2.2.0 code is functionally the
00038                   same as ELMv2.1a code
00039         June 04 VERSION 2.2.1 - some additions to code documentation
00040               -Tagged some new comments w/ "v2.2.1note" 
00041         Aug 2004 VERSION 2.3.0 - added dynamic boundary conditions 
00042                 - used new spatial time series data input for dynamic boundary conditions
00043         Oct 2004 VERSION 2.3.1 - internal release 
00044                 - checked that results reasonable, not fully checked
00045         Nov 2004 VERSION 2.3.2 - documentation upgrade
00046                 - added Doxygen tags in comments
00047                 - fmod, fabs
00048         Apr 2005 v2.4.4: added (actually redid) flexible canal output
00049         Aug 2005 v2.5.0: added getCanalElev function 
00050                 -  calculate slope of elev from the two end points of a reach (instead
00051                 of local elev following raster by raster elev)
00052         May 2006 v2.5.1: added edgeMann parameter to struct Chan
00053                 - allows Manning's n associated w/ edge of canal, to accomodate topographic
00054                 lip/berm and/or denser veg along canal length
00055 
00056 */
00057    
00058 
00059 #include "watmgmt.h"
00060 
00061 
00062 /****************************************************************************************/
00072 void Canal_Network_Init(float baseDatum, float *elev )
00073 
00074 { struct Structure *structs;
00075   int i,j;
00076   char filename[30];
00077 
00078   sprintf(filename,"CanalData");
00079 
00080 /* Input the canals geometry file, return pointer to first vector               */
00081 /* Configure the canals, marking the cells that interact with them              */ 
00082 /* and storing the pointers to the cells for each of the canals in Chan_list array */
00083   Channel_configure(elev, ReadChanStruct (filename));  
00084 
00085   for( i = 0; i < num_chan; i++ )
00086   {
00087       Chan_list[i]->P = Chan_list[i]->ic_P_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
00088       Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
00089       Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
00090   }
00091 
00092   usrErr ("\tCanal geometry configured OK");
00093   
00094 /* Read data on the water control structures, return pointer to the first structure */
00095   ReadStructures  (filename, baseDatum);
00096 
00097   ON_MAP[T(2,2)] = 0;
00098   
00099   MCopen = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MCopen"); /* open-water manning's coef, used in structure cell-cell bridge flow */
00100   init_pvar(MCopen,NULL,'f',0.05); /* this is a hard-wired parameter - but do not do anything much with it now for the bridge flows */
00101   
00102  
00103  if ( debug > 0 )       /* Output the canal/levee modifications to ON_MAP */
00104  {
00105         sprintf( modelFileName, "%s/%s/Output/Debug/ON_MAP_CANAL.txt", OutputPath, ProjName );
00106    
00107 
00108   if ( ( ChanInFile = fopen( modelFileName, "w" ) ) ==  NULL )
00109   {
00110      printf( "Can't open the %s canal/levee debug file!\n", modelFileName );
00111      exit(-1) ;
00112   }
00113   fprintf ( ChanInFile, "ROWS=%d\nCOLUMNS=%d\nFORMAT=text\nINFO=\"Debug data: Overland flow restrictions due to levees. \"\nMISSING=0\n", s0, s1 );
00114 
00115   for ( i = 1; i <= s0 ; i++ ) 
00116     { for ( j = 1 ; j <= s1 ; j++ )
00117          fprintf( ChanInFile, "%4d\t", ON_MAP[T(i,j)]) ;
00118       fprintf ( ChanInFile, "\n" );
00119     }   
00120   fclose( ChanInFile ) ;
00121  }
00122 
00123 
00124 /* Open file for writing structure flows */
00125     { if(H_OPSYS==UNIX) 
00126                 sprintf( modelFileName, "%s/%s/Output/Canal/structsOut", OutputPath, ProjName );
00127       else sprintf( modelFileName, "%s%s:Output:structsOut", OutputPath, ProjName );
00128     }
00129 
00130         /* Open file to print structure flows to */
00131     if ( ( WstructOutFile = fopen( modelFileName, "w" ) ) ==  NULL )
00132         {
00133        printf( "Can't open %s file!\n", modelFileName );
00134        exit(-1) ;
00135     }
00136 /* Print The Header Line For This File */
00137     fprintf ( WstructOutFile,
00138               "%s %s %s %s scenario: Sum of water flows at structures (ft^3/sec)\n     Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00139     structs = struct_first;
00140     while ( structs != NULL )
00141     {   fprintf( WstructOutFile, "%7s\t", structs->S_nam);       
00142             structs = structs->next_in_list;
00143     }
00144 
00145 /* Open file for writing structure P concs */
00146     { if(H_OPSYS==UNIX) 
00147                 sprintf( modelFileName, "%s/%s/Output/Canal/structsOut_P", OutputPath, ProjName );
00148       else sprintf( modelFileName, "%s%s:Output:structsOut_P", OutputPath, ProjName );
00149     }
00150 
00151         /*  Open file to print structure flows to */
00152     if ( ( WstructOutFile_P = fopen( modelFileName, "w" ) ) ==  NULL )
00153         {
00154        printf( "Can't open %s file!\n", modelFileName );
00155        exit(-1) ;
00156     }
00157 /* Print The Header Line For This File */
00158     fprintf ( WstructOutFile_P,
00159               "%s %s %s %s scenario: Flow weighted mean TP conc at structures ( mg/L) \n     Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00160     structs = struct_first;
00161     while ( structs != NULL ) 
00162     {   fprintf( WstructOutFile_P, "%7s\t", structs->S_nam);      
00163             structs = structs->next_in_list; 
00164     }
00165 
00166 /* Open file for writing structure Salin concs */
00167     { if(H_OPSYS==UNIX) 
00168                 sprintf( modelFileName, "%s/%s/Output/Canal/structsOut_S", OutputPath, ProjName );
00169       else sprintf( modelFileName, "%s%s:Output:structsOut_S", OutputPath, ProjName );
00170     }
00171 
00172         /*  Open file to print structure flows to */
00173     if ( ( WstructOutFile_S = fopen( modelFileName, "w" ) ) ==  NULL )
00174         {
00175        printf( "Can't open %s file!\n", modelFileName );
00176        exit(-1) ;
00177     }
00178 /* Print The Header Line For This File */
00179     fprintf ( WstructOutFile_S,
00180               "%s %s %s %s scenario: Flow weighted mean tracer concentration at structures ( g/L) \n     Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00181     structs = struct_first;
00182     while ( structs != NULL ) 
00183     {   fprintf( WstructOutFile_S, "%7s\t", structs->S_nam);      
00184             structs = structs->next_in_list; 
00185     }
00186 
00187 /****/
00188     
00189         /* Open file to print canal stage info to */
00190     { if(H_OPSYS==UNIX) 
00191                 sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut", OutputPath, ProjName );
00192       else sprintf( modelFileName, "%s%s:Output:CanalOut", OutputPath, ProjName );
00193     }
00194 
00195     if ( ( CanalOutFile = fopen( modelFileName, "w" ) ) ==  NULL )
00196         {
00197        printf( "Can't open %s file!\n", modelFileName );
00198        exit(-1) ;
00199     }
00200 
00201     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 */
00202     /*  channels with negative widths are reserved for new canals */
00203     /* v2.2 put a "R_" prefix in front of the canal ID# */
00204     for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile, "R_%d\t",Chan_list[i]->number);}
00205     
00206 
00207         /* Open file to print canal phosph info to */
00208     { if(H_OPSYS==UNIX) 
00209                 sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut_P", OutputPath, ProjName );
00210       else sprintf( modelFileName, "%s%s:Output:CanalOut_P", OutputPath, ProjName );
00211     }
00212 
00213     if ( ( CanalOutFile_P = fopen( modelFileName, "w" ) ) ==  NULL )
00214         {
00215        printf( "Can't open %s file!\n", modelFileName );
00216        exit(-1) ;
00217     }
00218 
00219     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 */
00220     /* v2.2 put a "R_" prefix in front of the canal ID# */
00221     for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile_P, "R_%d\t",Chan_list[i]->number);}
00222     
00223         /* Open file to print canal salinity/tracer info to */
00224     { if(H_OPSYS==UNIX) 
00225                 sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut_S", OutputPath, ProjName );
00226       else sprintf( modelFileName, "%s%s:Output:CanalOut_S", OutputPath, ProjName );
00227     }
00228 
00229     if ( ( CanalOutFile_S = fopen( modelFileName, "w" ) ) ==  NULL )
00230         {
00231        printf( "Can't open %s file!\n", modelFileName );
00232        exit(-1) ;
00233     }
00234 
00235     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 */
00236     /* v2.2 put a "R_" prefix in front of the canal ID# */
00237     for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile_S, "R_%d\t",Chan_list[i]->number);}
00238     
00239  
00240 
00241 /* Open file for canal-cell flux debugging purposes */
00242   if( debug > 0 ) 
00243   { 
00244      if(H_OPSYS==UNIX) 
00245                 sprintf( modelFileName, "%s/%s/Output/Debug/Can_debug", OutputPath, ProjName );
00246       else sprintf( modelFileName, "%s%s:Output:Can_debug", OutputPath, ProjName );
00247     
00248 
00249         /* Open file for debug */
00250     if ( ( canDebugFile = fopen( modelFileName, "w" ) ) ==  NULL )
00251         {
00252        printf( "Can't open %s file!\n", modelFileName );
00253        exit(-1) ;
00254     }
00255   
00256      fprintf ( canDebugFile, "Depending on level of the debug parameter, data on canal-cell flows printed here. \n" ); /* v2.5 made use of this file pointer */
00257      fflush(canDebugFile); 
00258     
00259   }
00260   
00261   return;
00262 
00263 }
00264 
00265 /********************************************************************************************/
00268 void CanalReInit()
00269 {
00270     int i;
00271     for( i = 0; i < num_chan; i++ )
00272 
00273     {
00274         Chan_list[i]->wat_depth = Chan_list[i]->ic_depth;
00275         Chan_list[i]->P = Chan_list[i]->ic_P_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
00276         Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
00277         Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth; /* kg/m3 * m2 * m */
00278          
00279     }
00280 }
00281 
00282 
00283 
00284 
00285 /********************************************************************************************/
00306 void Run_Canal_Network ( float *SWH, float *Elevation, float *MC, float *GW, float *poros, 
00307                          float *GWcond, double *NA, double *PA, double *SA, double *GNA, double *GPA, double *GSA,
00308                          float *Unsat, float *sp_yield) 
00309 
00310 { int i;
00311  float CH_vol;
00312 
00313  if ( canPrint && (( N_c_iter++ % hyd_iter ) == 0.0) ) {
00314         if (SensiOn) { 
00315             printchan = (SimTime.IsSimulationEnd) ? (1) : (0); /* flag to indicate to print canal/struct data */
00316          }
00317          else printchan = 1;
00318  }
00319  else  printchan = 0;
00320 
00321  if (printchan == 1) {
00322      
00323      fprintf( WstructOutFile, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00324      fprintf( WstructOutFile_P, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00325 /*      fprintf( WstructOutFile_N, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] ); */
00326      fprintf( WstructOutFile_S, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00327 
00328      fprintf (CanalOutFile, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00329      fprintf (CanalOutFile_P, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00330 /*      fprintf (CanalOutFile_N, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] ); */
00331      fprintf (CanalOutFile_S, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00332  }
00333     
00334   
00335  for ( i = 0; i < num_chan; i++ )
00336  {   /* set the sum of historical/other-data flows to/from a canal per iteration to 0 */        
00337      Chan_list[i]->sumHistIn = Chan_list[i]->sumHistOut = 0.0;
00338      Chan_list[i]->sumRuleIn = Chan_list[i]->sumRuleOut = 0.0;
00339  }
00340 
00341 /* calculate flows through water control structures */
00342 Flows_in_Structures ( SWH, Elevation, MC, NA, PA, SA ); 
00343   
00344 /* output  header for canal-cell flux output canal debug file, at relatively high debug levels (writes a LOT of stuff) */
00345  if( debug > 3 ) 
00346  {
00347      fprintf( canDebugFile, "N%3d  WatDepth     Pconc     Sconc        flux_S        flux_G        flux_L           Qin          Qout #_iter      Area    error \n", N_c_iter );
00348       
00349  }
00350 
00351 /* Flux the water along a canal using the relaxation procedure */
00352   
00353  for( i = 0; i < num_chan; i++ )
00354      if (Chan_list[i]->width > 0) /* a neg width canal is skipped */
00355 
00356      { 
00357          FluxChannel ( i, SWH, Elevation, MC, GW, poros, GWcond, NA, PA, SA, GNA, GPA, GSA, Unsat, sp_yield );
00358              /* calculate total canal volume after canal fluxes */
00359          CH_vol = Chan_list[i]->area*Chan_list[i]->wat_depth; 
00360 
00361              /* sum the storages on iteration that performs budget checks */
00362          if ( !( FMOD(N_c_iter,(1/canstep) )) && (SimTime.IsBudgEnd ) ) { 
00363              TOT_VOL_CAN[Chan_list[i]->basin] += CH_vol;
00364              TOT_VOL_CAN[0] += CH_vol;
00365              TOT_P_CAN[Chan_list[i]->basin] += Chan_list[i]->P;
00366              TOT_P_CAN[0] += Chan_list[i]->P;
00367              TOT_S_CAN[Chan_list[i]->basin] += Chan_list[i]->S;
00368              TOT_S_CAN[0] += Chan_list[i]->S;
00369          }
00370                                                                                                
00371      
00372          if ( printchan )
00373          {
00374                  /* report N&P concentration in mg/L (kg/m3==g/L * 1000 = mg/L) */
00375              Chan_list[i]->P_con = (CH_vol > 0) ? (Chan_list[i]->P/CH_vol*1000.0) : 0.0;
00376 /*              Chan_list[i]->N_con = (CH_vol > 0) ? (Chan_list[i]->N/CH_vol*1000.0) : 0.0; */
00377              Chan_list[i]->S_con = (CH_vol > 0) ? (Chan_list[i]->S/CH_vol       ) : 0.0;
00378  
00379              /* v2.2 increased decimal places in output (from 6.3f to 7.4f) */
00380              fprintf( CanalOutFile, "%6.2f\t", Chan_list[i]->wat_depth );
00381              fprintf( CanalOutFile_P, "%7.4f\t", Chan_list[i]->P_con );
00382 /*          fprintf( CanalOutFile_N, "%7.4f\t", Chan_list[i]->N_con ); */
00383              fprintf( CanalOutFile_S, "%7.4f\t", Chan_list[i]->S_con ); 
00384          }
00385      } 
00386  if ( printchan )
00387  {
00388      fflush(CanalOutFile);
00389      fflush(CanalOutFile_P);
00390 /*       fflush(CanalOutFile_N); */
00391      fflush(CanalOutFile_S);
00392  }
00393   
00394  return;
00395 }
00396 
00397                   
00398 
00399 /********************************************************************************************/
00406 struct Chan *ReadChanStruct ( char* filename )
00407 { 
00408  struct Chan *channels, *chan_first, *chan_last;
00409  struct Chan_reach *C_reach, *last_reach, *next_reach;
00410  int i, index, firstReach, firstCanal = 1;
00411  char ss[422], scenario[20], modnam[20], bas_txt[8];
00412  float C_x, C_y, C_x1, C_y1;
00413 
00414  
00415  int ibas, foundBasn;
00416  
00417 
00418  { if(H_OPSYS==UNIX) 
00419      sprintf( modelFileName, "%s/%s/Data/%s.chan", ModelPath, ProjName, filename );
00420  else sprintf( modelFileName, "%s%s:Data:%s.chan", ModelPath, ProjName, filename );
00421  }
00422 
00423 /* Open file with canals data */
00424  if ( ( ChanInFile = fopen( modelFileName, "r" ) ) ==  NULL )
00425  {
00426      sprintf( msgStr,"Can't open %s file! ",modelFileName ) ; usrErr(msgStr);
00427      exit(-1) ;
00428  }
00429   
00430 /* Allocate memory for canal reach structure. A canal reach is the straight canal segment */
00431 /* defined in the canal data file by the coordinates of its firstReaching and ending points.   */
00432 /* Each canal can be made of several canal reaches              */
00433  if ( (C_reach = (struct Chan_reach *) malloc( (size_t) sizeof( struct Chan_reach ))) == NULL )
00434  {
00435      usrErr( "No memory for first canal reach\n " ) ;
00436      exit( -2 ) ;
00437  };
00438 
00439  num_chan = 0;
00440    
00441 /* Read the general canal network information */ 
00442  fgets( ss, 420, ChanInFile ); sscanf( ss, "%d", &UTM_EOrigin );
00443  fgets( ss, 420, ChanInFile ); sscanf( ss, "%d", &UTM_NOrigin );
00444  fgets( ss, 420, ChanInFile );  /* skip comment line */
00445  fgets( ss, 420, ChanInFile ); sscanf( ss,"%s %s %d %d", &modnam, &scenario, &CHo, &CHANNEL_MAX_ITER );
00446  if (strcmp(scenario,SimAlt) != 0) {
00447      sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
00448              scenario, modelFileName, SimAlt); usrErr(msgStr);
00449              exit(-1);
00450  }
00451  if (strcmp(modnam,modelName) != 0) {
00452      sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
00453              modnam, modelFileName, modelName); usrErr(msgStr);
00454              exit(-1);
00455  }
00456   
00457  fgets( ss, 420, ChanInFile );   /* skip comment line */
00458  fgets( ss, 420, ChanInFile ); sscanf ( ss, "%f %f %f ", &F_ERROR, &MINDEPTH, &C_F );
00459      /* the C_F is used only for sensitivity experiments */
00460  fgets( ss, 420, ChanInFile );   /* skip comment line */
00461    
00462 /* Read canal descriptors until none left */   
00463  while ( fgets( ss, 420, ChanInFile ) != NULL && !feof( ChanInFile ) )
00464  { 
00465      if ( *ss == '#' )  /* "#" identifies each new canal in the data file */
00466      { /* Allocate memory for canal structure */
00467          if ( (channels = (struct Chan *) malloc( (size_t) sizeof( struct Chan ))) == NULL )
00468          {
00469              printf( "No memory for first canal\n " ) ;
00470              exit( -2 ) ;
00471          };
00472                 
00473              /* Read canal information:  - canal number, 
00474                 - levee marker (1 - left levee, -1 - right, 0 - none, 2 - both sides),
00475                 - ranges of cell interaction to the left and to the right - functionality removed, should ALWAYS=1,
00476                 - canal depth and width,
00477                 - seepage coefficient through the levee,
00478                 - initial S and P concentrations (doubles) in the canal water,
00479                 - initial water depth in canal (can be >, = or < than canal depth 
00480                 - Manning's n associated w/ edge of canal, to accomodate topographic lip/berm and/or denser veg along canal length
00481                 - hydrologic basin that canal exchanges surface water with */
00482          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",
00483                      &channels->number, 
00484                      &channels->levee,
00485                      &channels->roil,   &channels->roir, 
00486                      &channels->depth,  &channels->width,  
00487                      &channels->cond,   
00488                      &channels->ic_S_con,  &channels->ic_P_con,  
00489                      &channels->wat_depth, &channels->edgeMann, &bas_txt);
00490 
00491          foundBasn=0;
00492          for (ibas=numBasn;ibas>0;ibas--) {
00493              basins = basn_list[ibas];
00494              if(strcmp(bas_txt,basins->basnTxt)==0) {
00495                  channels->basin = ibas;
00496                  foundBasn=1;
00497                  channels->family = basn_list[ibas]->family;
00498                  channels->parent = basn_list[ibas]->parent;
00499                      /* canals may ONLY belong to parent hydro basins, not to indicator regions */
00500                  if (!channels->parent) {
00501                      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",
00502                             modelName,channels->number); exit (-1);}     
00503              }
00504          }
00505          if (foundBasn==0) { printf ("Sorry - %s's canal %d's basin does not exist! Please define the basin in the basinIR file.\n",
00506                             modelName,channels->number); exit (-1);}                
00507 
00508          channels->ic_P_con *= 0.001; /* initial input value in mg/L, convert to g/L==kg/m3 */
00509          channels->ic_N_con = 0.00; /* v2.2 "removed" N */ /* initial input value in mg/L, convert to g/L==kg/m3 */
00510          channels->ic_S_con *= 1.0; /* initial input value in g/L, no conversion */
00511              /* store the initial depth for re-initializing under Positional Analysis mode */
00512          channels->ic_depth = channels->wat_depth;
00513          
00514 
00515          channels->reaches = C_reach; 
00516                         
00517          if (i < 10)    /* check for the number of input fields read */
00518          {printf ( " Error in CanalData.chan: reach input %d\n", (num_chan+1) ); exit (-1);}
00519                         
00520                             
00521          num_chan++;            /*count the number of canals (not canal reaches) */
00522          firstReach = 1;
00523                         
00524          if (firstCanal) 
00525          {
00526              firstCanal = 0;  
00527              chan_first = channels;
00528              chan_last = channels;
00529              continue; 
00530          }
00531 
00532          last_reach->next_reach = NULL;
00533          chan_last->next_in_list = channels;
00534          chan_last = channels;     
00535      }
00536      else  
00537      { /* read the pairs of coordinates for the consecutive canal reaches */
00538          i = sscanf ( ss, "%f %f", &C_x1, &C_y1 );
00539          if (i < 2)     /* check for the number of input fields read */
00540          {printf ( " Error in CanalData.chan coords: %d'th reach input \n", (num_chan+1) ); exit (-1);}
00541                 
00542              /* Convert from UTM to km units with 0,0 origin */
00543          C_x = UTM2kmx (C_x1);  C_y = UTM2kmy (C_y1);           
00544          C_reach->x1 = C_x;  C_reach->y1 = C_y;
00545                         
00546          if ( firstReach )                   
00547          {  C_reach->x0 = C_x; C_reach->y0 = C_y; 
00548          firstReach = 0;
00549          }
00550          else 
00551          {  
00552                                 /* Allocate memory for canal reach structure */
00553              if ( (next_reach = (struct Chan_reach *) malloc( (size_t) sizeof( struct Chan_reach ))) == NULL )
00554              {  printf( "No memory for first canal reach\n " ) ;
00555              exit( -2 ) ;
00556              };
00557              C_reach->next_reach = next_reach;
00558              last_reach = C_reach;
00559              C_reach = next_reach;
00560              C_reach->x0 = C_x; C_reach->y0 = C_y;
00561          }                       
00562      }  
00563  }/* end loop for reading */
00564   
00565  chan_last->next_in_list = NULL;
00566  last_reach->next_reach = NULL;
00567  free ((char *) C_reach);
00568  fclose (ChanInFile);
00569    
00570  usrErr( "\tCanal locs/attributes read OK" );
00571    
00572  return (chan_first);
00573 }
00574 
00575 /****************************************************************************************/
00582 void ReadStructures  ( char* filename, float BASE_DATUM )
00583 {
00593 struct Structure *structs, *struct_next, *struct_last;
00594 
00595  char ss[322], *s;
00596  char sched_name[20];
00597  int i, j, k;
00598  int disAggNum=0; 
00599  int IDhist[MAX_H_STRUCT]; 
00600  int lincnt; 
00601  double Jdate_read;
00602  char lline[2001], *line;
00603  int yyyy, mm, dd;
00604  char structName[20]; 
00605  char scenario[20], s_modname[20];
00606  char histTP[5]; 
00607  char histTS[5]; 
00608  
00609  num_struct_hist = 0;
00610  numTPser = numTSser = 0;
00611  
00612  { if(H_OPSYS==UNIX) 
00613      sprintf( modelFileName, "%s/%s/Data/%s.struct", ModelPath, ProjName, filename );
00614  else sprintf( modelFileName, "%s%s:Data:%s.struct", ModelPath, ProjName, filename );
00615  }
00616 
00617 /* Open file with structures data */
00618  if ( ( ChanInFile = fopen( modelFileName, "r" ) ) ==  NULL )
00619  {
00620      printf( "Can't open %s file!\n", modelFileName ) ;
00621      exit(-1) ;
00622  }
00623     
00624 /* Allocate memory for structures */
00625  if ( (structs = (struct Structure *) malloc( (size_t) sizeof( struct Structure ))) == NULL )
00626  {
00627      printf( "No memory for first structure\n " ) ;
00628      exit( -2 ) ;
00629  }
00630  struct_first = structs;
00631 
00632  fgets( ss, 320, ChanInFile );
00633  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 */
00634  if (strcmp(scenario,SimAlt) != 0) {
00635          sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
00636                  scenario, modelFileName, SimAlt); usrErr(msgStr);
00637      exit(-1);
00638   }
00639  if (strcmp(s_modname,modelName) != 0) {
00640          sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) in the canal data file!",
00641                  s_modname, modelFileName, modelName); usrErr(msgStr);
00642      exit(-1);
00643   }
00644  
00645 /* Skip 2'nd header line */
00646  fgets( ss, 320, ChanInFile );
00647 
00648 
00649 /* Read structure descriptors until none left */   
00650  while ( fgets( ss, 320, ChanInFile ) != NULL && !feof( ChanInFile ) )
00651  { 
00652      s = ss; 
00653          
00654          
00655      sscanf( s, "%d", &structs->flag ); 
00656 /* structure flag: <0 - skip structure info, structure not operating; */
00657 /*                 >0 - historical/other-data data array, flag=1 normally in this case, but is >1 for aggregated SFWMM structures; */
00658 /*                  0 - rule-based structure   */
00659      if ( structs->flag < 0 ) continue;
00660      if ( structs->flag > 0 && structs->flag <10 ) {
00661          structs->histID = ++num_struct_hist;
00662          Hist_data[structs->histID-1].flag = structs->flag;
00663      }
00664      
00665          
00666 /*    assign integer to the historical data ID, and increment the counter for structures with historical/other-data data to be read from datafile. */
00667 /*    flags >= 20 represent structures (such as the S-11s) that were aggregated in SFWMM output */
00668 /*    and need to be assigned disaggregated flows for separate ELM structures */
00669 /*    Ex: SFWMM structure name S11 may have flag = 2 in CanalData.struct.  That S11 flow holds the sum of S-11A-C flows. */
00670 /*    The flags on S-11A, S-11B, and S-11C would all have been assigned 20, 20, and 20 in the CanalData.struct file */
00671 
00672 /*           NO LONGER NECESSARY: as a temporary measure until the SFWMM output .dss file is fixed, the structures that have calculated "other" partitions */
00673 /*             (S140,S7,S8) are assigned historical flags of 9 (they are in the dss file read by a SAS routine, the flag */
00674 /*             of 9 allows them to remain in the file, but the structure is turned off because it has flag >1)  */
00675          
00676      
00677      s = Scip( s, TAB );
00678          
00679      if ( *s != TAB ) {
00680              sscanf ( s, "%s", &structs->S_nam );
00681              if ( structs->flag >0 && structs->flag <10) 
00682                  sscanf ( s, "%s", &Hist_data[structs->histID-1].S_nam );
00683              
00684                  /* for solute concentration of flows thru structures, */
00685                  /* data use flag XXser: -1== simulated, 0==static historical, 1= time-varying historical */
00686              s = Scip( s, TAB );
00687              if ( *s == TAB ) structs->TPser = -1; 
00688              else {
00689                  sscanf ( s, "%s", &histTP ); /* read the field that either has TP conc (ug/L=ppb) that is constant across time
00690                                                    or has (any) string to indicate the use of TimeSeries data from a datafile (read later)*/
00691                  if (isInteger(histTP) ) {
00692                       structs->TPser = 0; 
00693                       structs->TP = atof(histTP)*1.0e-6;  /* conc (ug/L=>g/L) that is constant across time */
00694                  }
00695                  else { structs->TPser = 1; numTPser++;} /* otherwise, we will be reading time series data later */
00696              }
00697              
00698              s = Scip( s, TAB );
00699              if ( *s == TAB ) structs->TNser = -1; /* value to indicate use of simulated, not historical, data */
00700              else {
00701                  sscanf ( s, "%f", &structs->TN );
00702                  structs->TNser = 0;
00703                  structs->TN *= 1.0e-6;/* read the TN conc (ug/L=>g/L) that is constant across time */
00704              }
00705              
00706 
00707              s = Scip( s, TAB );
00708              if ( *s == TAB ) structs->TSser = -1; /* value to indicate use of simulated, not historical, data */
00709              /* v2.2 added code to allow use of TimeSeries input data in addition to old method of a constant concentration */
00710              else {
00711                   sscanf ( s, "%s", &histTS ); /* read the field that either has TS (salt) conc (g/L=ppt) that is constant across time
00712                                                    or has (any) string to indicate the use of TimeSeries data from a datafile (read later)*/
00713                  if (isFloat(histTS) ) {
00714                       structs->TSser = 0; 
00715                       structs->TS = atof(histTS)*1.0;  /* conc (no conversion, g/L) that is constant across time */
00716                  }
00717                  else { structs->TSser = 1; numTSser++;} /* otherwise, we will be reading time series data later */
00718 
00719              }
00720              
00721              s = Scip( s, TAB );
00722              if ( *s == TAB ) structs->str_cell_i = 0;  /* non 0 if elevation at the structure location is desired */
00723              else {
00724                  sscanf( s, "%d", &structs->str_cell_i );
00725                  structs->str_cell_i++; /* this converts the data map coords to model coords (map+1) */
00726              }
00727              
00728              s = Scip( s, TAB );
00729 
00730              if ( *s == TAB ) structs->str_cell_j = 0;  /* non 0 if elevation at the structure location is desired */
00731              else {
00732                  sscanf( s, "%d", &structs->str_cell_j );
00733                  structs->str_cell_j++;
00734              }
00735 
00736              if ( structs->str_cell_i != 0 && !ON_MAP[T(structs->str_cell_i, structs->str_cell_j)] )
00737              {  printf ( "Location of water control structure %s is off-map.\n", structs->S_nam );
00738              exit (-2);                 /*check that cells obtained are within on-map boundaries */
00739              }
00740              
00741              
00742      }
00743          
00744      s = Scip( s, TAB );
00745 
00746      if ( *s == TAB ) structs->canal_fr  = 0;   /* non 0 if the structure involves a donor canal */
00747      else sscanf( s, "%d", &structs->canal_fr  );
00748      s = Scip( s, TAB );
00749      if ( *s == TAB) structs->canal_to = 0;             /* non 0 if the structure involves a receiving canal */
00750      else sscanf( s, "%d", &structs->canal_to );
00751      s = Scip( s, TAB );
00752 
00753      if ( structs->canal_fr  - CHo >= num_chan || structs->canal_to - CHo >= num_chan )
00754      {  printf ( "Canal from/to doesn't exist for water control structure %s\n", structs->S_nam );
00755      exit (-2);                 /* check that canals read are consistent with canals data */
00756      }           
00757 /*  check to ensure canals associated with structures are valid */
00758      if ( structs->canal_fr != 0 && Chan_list[structs->canal_fr- CHo]->width <= 0.1 )
00759      { printf ( "The donor canal %d has been turned off for water control structure %s\n", structs->canal_fr, structs->S_nam );
00760      exit (-2);                 
00761      }           
00762      if ( structs->canal_to != 0 && Chan_list[structs->canal_to- CHo]->width <= 0.1 )
00763      {  printf ( "The recipient canal %d has been turned off for water control structure %s\n", structs->canal_to, structs->S_nam );
00764      exit (-2);                 
00765      }           
00766          
00767      if (*s == TAB) structs->cell_i_fr = 0;  /* non 0 if the structure involves a donor cell */
00768      else 
00769      { sscanf( s, "%d", &structs->cell_i_fr );
00770      structs->cell_i_fr++;
00771      }   
00772      s = Scip( s, TAB );
00773      if (*s == TAB) structs->cell_j_fr = 0;
00774      else 
00775      { sscanf( s, "%d", &structs->cell_j_fr );
00776      structs->cell_j_fr++;
00777      }
00778      s = Scip( s, TAB );
00779          
00780 
00781      if ( structs->cell_i_fr > 2 && !ON_MAP[T(structs->cell_i_fr, structs->cell_j_fr)] )
00782      {  printf ( "Donor cell for water control structure %s is off-map.\n", structs->S_nam );
00783      exit (-2);                 /*check that cells obtained are within on-map boundaries */
00784      }
00785 
00786          /* if valid donor cell, mark it on ON_MAP */
00787      if ( structs->cell_j_fr ) ON_MAP[T(structs->cell_i_fr, structs->cell_j_fr)] += 100;
00788          
00789      if (*s == TAB) structs->cell_i_to = 0;  /* non 0 if the structure involves a recipient cell */
00790      else 
00791      { sscanf( s, "%d", &structs->cell_i_to );
00792      structs->cell_i_to++;
00793      }   
00794          
00795      s = Scip( s, TAB );
00796      if (*s == TAB) structs->cell_j_to = 0;
00797      else
00798      { sscanf( s, "%d", &structs->cell_j_to );
00799      structs->cell_j_to++;
00800      }   
00801 
00802 
00803      if ( structs->cell_i_to > 2 && !ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] )
00804      {  printf ( "Recipient cell for water control structure %s is off-map.\n", structs->S_nam );
00805      exit (-2);                 /*check that cells obtained are within on-map boundaries */
00806      }
00807 
00808          /* if valid recieving cell, mark it on ON_MAP */
00809          /* this used to check for cell_j_to > 2, not putting recipient LEC cells ON_MAP */
00810      if ( structs->cell_j_to ) ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] += 100;
00811          
00812      s = Scip( s, TAB );                /* read the targeted water levels */
00813      if (*s == TAB) structs->HW_stage = 0;
00814      else if ( isalpha (*s) ) {  /* if there is a graph for the schedule */
00815          sscanf( s, "%s", sched_name );
00816          structs->HW_graph = Read_schedule ( sched_name, filename, BASE_DATUM );
00817          structs->HW_stage = -99;  
00818      } 
00819      else   { sscanf( s, "%f", &structs->HW_stage );
00820      structs->HW_stage += BASE_DATUM;
00821      }
00822                         
00823      s = Scip( s, TAB );
00824      if (*s == TAB) structs->TW_stage = 0;
00825      else if ( isalpha (*s) ) {  /* if there is a graph for the schedule */
00826          sscanf( s, "%s", sched_name );
00827          structs->TW_graph = Read_schedule ( sched_name, filename, BASE_DATUM );
00828          structs->TW_stage = -99;
00829      } 
00830      else       { sscanf( s, "%f", &structs->TW_stage );
00831      structs->TW_stage -= BASE_DATUM; 
00832      }
00833      s = Scip( s, TAB );
00834          
00835      if (*s == TAB) structs->cell_i_HW = 0;
00836      else               /* read coordinates of the reference cell for head water */
00837      { sscanf( s, "%d", &structs->cell_i_HW );
00838      structs->cell_i_HW++;
00839      }   
00840      s = Scip( s, TAB );
00841      if (*s == TAB) structs->cell_j_HW = 0;
00842      else 
00843      { sscanf( s, "%d", &structs->cell_j_HW );
00844      structs->cell_j_HW++;
00845      }   
00846      s = Scip( s, TAB );
00847 
00848      if ( structs->cell_i_HW > s0 || structs->cell_j_HW > s1 )
00849      {  printf ( "Error in definition of water control structure %s", structs->S_nam );
00850      exit (-2); /*check that cell obtained is within array boundaries */
00851      }
00852 
00853      if (*s == TAB) structs->cell_i_TW = 0;
00854      else               /* read coordinates of the reference cell for tail water */
00855      { sscanf( s, "%d", &structs->cell_i_TW );
00856      structs->cell_i_TW++;
00857      }   
00858      s = Scip( s, TAB );
00859      if (*s == TAB) structs->cell_j_TW = 0;
00860      else 
00861      { sscanf( s, "%d", &structs->cell_j_TW );
00862      structs->cell_j_TW++;
00863      }   
00864      s = Scip( s, TAB );
00865 
00866      if ( structs->cell_i_TW > s0 || structs->cell_j_TW > s1 )
00867      {  printf ( "Error in definition of water control structure %s", structs->S_nam );
00868      exit (-2); /*check that cell obtained is within array boundaries */
00869      }
00870 
00871      if (*s == TAB || *s == '?') structs->w_coef = 0;
00872      else sscanf( s, "%f", &structs->w_coef );
00873          /* the next field after the coef is the date the file was exported, ignored here */
00874         
00875      struct_last = structs;
00876           
00877          /* Allocate memory for next structure */
00878      if ( (struct_next = (struct Structure *) malloc( (size_t) sizeof( struct Structure ))) == NULL )
00879      {
00880          printf( "Failed to allocate memory for next structure (%s)\n ", structs->S_nam ) ;
00881          exit( -2 ) ;
00882      };
00883      
00884      /* initialize the sum of structure  flows and constit solute in flows */
00885      structs->SumFlow = 0.0; 
00886      structs->Sum_P = 0.0; 
00887      structs->Sum_N = 0.0; 
00888      structs->Sum_S = 0.0; 
00889 
00890      structs->next_in_list = struct_next;
00891      structs = struct_next;
00892  }
00893    
00894  struct_last->next_in_list = NULL;
00895  free ((char *)structs);
00896  fclose (ChanInFile);
00897    
00898 
00899  usrErr( "\tWater control structure locs/attributes read OK" );
00900  usrErr0 ( "\tControl structures' water flow time series...");
00901  
00902    
00903 /**********************/
00904 /* read structure flow input file(s) */
00905      { if(H_OPSYS==UNIX) 
00906          sprintf( modelFileName, "%s/%s/Data/%s.struct_wat", ModelPath, ProjName, filename );
00907      else sprintf( modelFileName, "%s%s:Data:%s.struct_wat", ModelPath, ProjName, filename );
00908      }
00909 
00910 /* Open file with structure flow data */
00911      if ( ( F_struct_wat = fopen( modelFileName, "r" ) ) ==  NULL )
00912      {
00913          printf( "Can't open %s file!\n", modelFileName ) ;
00914          exit(-1) ;
00915      }
00916 
00917      fgets(lline, 2000, F_struct_wat); /* first header line with the alternative's name */
00918      line=lline;
00919      sscanf(line, "%s",&scenario);
00920      if (strcmp(scenario,SimAlt) != 0) {
00921          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);
00922          exit(-1);
00923      }
00924 
00925      fgets(lline, 2000, F_struct_wat); /* read 2'nd header line with column names */
00926      line=lline;
00927      line = Scip( line, TAB );
00928      line = Scip( line, TAB );
00929      line = Scip( line, TAB ); /*read past the three date columns, ready to read the variable names */
00930 
00931 
00932      /* loop through all the control structure columns */
00933      for ( i = 0; i < num_struct_hist; i++ ) { 
00934 
00935          Hist_data[i].arrayWat = (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
00936          Hist_data[i].arrayP =   (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
00937           Hist_data[i].arrayN =  (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 
00938           Hist_data[i].arrayS =  (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 
00939 
00940 /*  allocate memory for dissaggregated structures and provide a couple of attributes */
00941          if ( Hist_data[i].flag > 1 ) {
00942              structs = struct_first;  /* go through all structures to pull out the disaggregated ones */
00943              while ( structs != NULL ) 
00944              {
00945                  /* Ex: the SFWMM aggregated S10's flag is 2, S10AB&C's flags are 20 */
00946                  if ( structs->flag == 10 * Hist_data[i].flag  )  
00947                  { 
00948                         /* disAggNum counts the total number of disaggregated structures for simulation */
00949                      disAggNum++;
00950                      Hist_data[num_struct_hist+disAggNum-1].arrayWat = 
00951                              (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
00952                      Hist_data[num_struct_hist+disAggNum-1].arrayP    = 
00953                              (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
00954                      Hist_data[num_struct_hist+disAggNum-1].arrayN    = 
00955                              (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 
00956                      Hist_data[num_struct_hist+disAggNum-1].arrayS    = 
00957                              (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" ); 
00958 
00959                      Hist_data[num_struct_hist+disAggNum-1].aggID = i; /* the aggregated structure ID is attached to the dissagregated one */
00960                      strcpy( Hist_data[num_struct_hist+disAggNum-1].S_nam, structs->S_nam ); /* the disaggregated structure's name is put into the HistData structure */
00961                      Hist_data[i].aggCnt++; /* increment a count of the number of disagg structs associated w/ each agg struct */
00962                      structs->flag = 1; /* set the disagg struct's processing flag to 1, the only flag used for simulation structure flows (new input style) */
00963                      structs->histID = num_struct_hist+disAggNum; /* store the historical ID value for the disagg structures */
00964                      
00965                  } 
00966    
00967                  structs = structs->next_in_list;
00968              }             
00969          }
00970          
00971          sscanf ( line,"%s\t",&structName); /* read the structure names in the header */
00972          if (strcmp(Hist_data[i].S_nam , structName) != 0) {
00973              printf( "variable name  %s in input file does not match order of CanalData.structs.\n", structName);
00974              exit(-1); 
00975              }
00976          line = Scip( line, TAB );
00977      }
00978 
00979 
00980      lincnt = 0; /*  historical/other-data data input file record #  starting at the point where the model reads values */
00981      while ( fgets(lline, 2000, F_struct_wat) != NULL && !feof( F_struct_wat ) )
00982      {
00983          line = lline;
00984          sscanf( line, "%d %d %d",&yyyy,&mm,&dd); /* read the date of the current record */
00985          
00986         /* julian day, returns a double (args = mon, day, year, hours, minutes, seconds) */
00987          Jdate_read = julday( mm,dd,yyyy,0,0,0.0 ); 
00988          if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
00989              printf (" \n***Error\nStructure flow data (%s) init date > simulation start date\n", modelFileName); exit (-1);
00990          }
00991          if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
00992          {
00993              line = Scip( line, TAB );
00994              line = Scip( line, TAB );
00995 
00996              for ( i = 0; i < num_struct_hist; i++ ) { /* loop through all the control structure columns in the input file */
00997                  line = Scip( line, TAB );
00998                  if (*line == TAB) Hist_data[i].arrayWat[((int)(lincnt))] = 0;
00999                  else sscanf( line, "%f", &Hist_data[i].arrayWat[((int)(lincnt))] );
01000                  Hist_data[i].arrayWat[((int)(lincnt))] *= 2446.848; /* cfs => m^3/d (0.02832*60*60*24)  */
01001              }
01002          
01003          for (i = num_struct_hist; i < num_struct_hist+disAggNum; i++) { /* dissaggregate all agg structs */
01004              k = Hist_data[i].aggID; /* use the aggregated historical data WCstructure */
01005              Hist_data[i].arrayWat[((int)(lincnt))] =
01006                  Hist_data[k].arrayWat[((int)(lincnt))] / (Hist_data[k].aggCnt) ;
01007 
01008          }
01009          
01010               
01011          lincnt++; /* increment the number of daily data records */
01012      }
01013 
01014 } /* end of reading structure flow records */
01015      
01016 if (lincnt == 0) {
01017     printf (" \n***Error\nStructure flow data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
01018 }
01019 else if (lincnt != PORnumday ) {
01020     printf (" \n***Error\nStructure flow data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
01021 }
01022      
01023 sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
01024 fclose (F_struct_wat);
01025 
01026 /**********************/
01027 /* read time series structure P concentration input file(s) */
01028 /* there can be situations when we have no time series TP data (numTPser==0) */
01029  if (numTPser>0) {
01030      usrErr0 ( "\tControl structures' TP concen. time series...");
01031 
01032      { if(H_OPSYS==UNIX) 
01033          sprintf( modelFileName, "%s/%s/Data/%s.struct_TP", ModelPath, ProjName, filename );
01034      else sprintf( modelFileName, "%s%s:Data:%s.struct_TP", ModelPath, ProjName, filename );
01035      }
01036 
01037 /* Open file with structure concentration data */
01038      if ( ( F_struct_TP = fopen( modelFileName, "r" ) ) ==  NULL )
01039      {
01040          printf( "Can't open %s file!\n", modelFileName ) ;
01041          exit(-1) ;
01042      }
01043 
01044      fgets(lline, 2000, F_struct_TP); /* first header line with the alternative's name */
01045      line=lline;
01046      sscanf(line, "%s",&scenario);
01047      if (strcmp(scenario,SimAlt) != 0) {
01048          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);
01049          exit(-1);
01050      }
01051      fgets(lline, 2000, F_struct_TP); /* read 2'nd header line with column names */
01052      line=lline;
01053      line = Scip( line, TAB );
01054      line = Scip( line, TAB );
01055      line = Scip( line, TAB );
01056 
01057 /* now we're past the three date columns, ready to read the variable names */
01058      for ( j = 0; j < numTPser; j++ ) { /* loop through all the data columns */
01059          IDhist[j] = -1;
01060          strcpy(structName,"NULL"); /* a string structName read previously is retained - this NULL assignment will take care of case
01061                                    where nothing was found/read in the time series file, and thus cause an error */
01062          sscanf ( line,"%s\t",&structName); /* read the structure names in the header */
01063          for ( i = 0; i < num_struct_hist+disAggNum; i++ ) { 
01064              if (strcmp(Hist_data[i].S_nam , structName) == 0) {
01065                  IDhist[j] = i;
01066                  line = Scip( line, TAB );
01067              }
01068          }
01069          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);
01070          exit(-1); 
01071          }
01072      }
01073 
01074      lincnt = 0; /*  input file record #  starting at the point where the model reads values */
01075 
01076      while ( fgets(lline, 2000, F_struct_TP) != NULL && !feof( F_struct_TP ) )
01077      {
01078          line = lline;
01079          sscanf( line, "%d %d %d",&yyyy,&mm,&dd); /* read the date of the current record */
01080          
01081         /* julian day, returns a double (args = mon, day, year, hours, minutes, seconds) */
01082          Jdate_read = julday( mm,dd,yyyy,0,0,0.0 ); 
01083          if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
01084              printf (" \n***Error\nStructure TP data (%s) init date > simulation start date\n", modelFileName); exit (-1);
01085          }
01086          if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
01087          {
01088              line = Scip( line, TAB );
01089              line = Scip( line, TAB );
01090              for ( j = 0; j < numTPser; j++ ) { /* loop through all the control structure columns in the input file */
01091                  line = Scip( line, TAB );
01092                  if (*line == TAB) Hist_data[IDhist[j]].arrayP[((int)(lincnt))] = 0;
01093                  else sscanf( line, "%f", &Hist_data[IDhist[j]].arrayP[((int)(lincnt))] );
01094                  Hist_data[IDhist[j]].arrayP[((int)(lincnt))] *= 1.0E-6; /* data read in ug/L, convert to g/L=kg/m3 */
01095                  
01096              }
01097               
01098              lincnt++; /* increment the number of daily data records */
01099          }
01100 
01101      } /* end of reading TP conc data records */
01102      
01103      if (lincnt == 0 ) {
01104              printf (" \n***Error\nTP Time series data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
01105      }
01106      else if (lincnt != PORnumday ) {
01107              printf (" \n***Error\nTP Time series data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
01108      }
01109      sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
01110      
01111      fclose (F_struct_TP);
01112  } /* end of reading TP time series data (if needed) */
01113  
01114  
01115  
01116 /**********************/
01117 /* read time series structure tracer (salt) concentration input file(s) */
01118 /* there can be situations when we have no time series TS (salt) data (numTSser==0) */
01119  if (numTSser>0) {
01120      usrErr0 ( "\tControl structures' Tracer concen. time series...");
01121 
01122      { if(H_OPSYS==UNIX) 
01123          sprintf( modelFileName, "%s/%s/Data/%s.struct_TS", ModelPath, ProjName, filename );
01124      else sprintf( modelFileName, "%s%s:Data:%s.struct_TS", ModelPath, ProjName, filename );
01125      }
01126 
01127 /* Open file with structure concentration data */
01128      if ( ( F_struct_TS = fopen( modelFileName, "r" ) ) ==  NULL )
01129      {
01130          printf( "Can't open %s file!\n", modelFileName ) ;
01131          exit(-1) ;
01132      }
01133 
01134      fgets(lline, 2000, F_struct_TS); /* first header line with the alternative's name */
01135      line=lline;
01136      sscanf(line, "%s",&scenario);
01137      if (strcmp(scenario,SimAlt) != 0) {
01138          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);
01139          exit(-1);
01140      }
01141      fgets(lline, 2000, F_struct_TS); /* read 2'nd header line with column names */
01142      line=lline;
01143      line = Scip( line, TAB );
01144      line = Scip( line, TAB );
01145      line = Scip( line, TAB );
01146 
01147 /* now we're past the three date columns, ready to read the variable names */
01148      for ( j = 0; j < numTSser; j++ ) { /* loop through all the data columns */
01149          IDhist[j] = -1;
01150          strcpy(structName,"NULL"); /* a string structName read previously is retained - this NULL assignment will take care of case
01151                                    where nothing was found/read in the time series file, and thus cause an error */
01152          sscanf ( line,"%s\t",&structName); /* read the structure names in the header */
01153          for ( i = 0; i < num_struct_hist+disAggNum; i++ ) { 
01154              if (strcmp(Hist_data[i].S_nam , structName) == 0) {
01155                  IDhist[j] = i;
01156                  line = Scip( line, TAB );
01157              }
01158          }
01159          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);
01160          exit(-1); 
01161          }
01162      }
01163 
01164      lincnt = 0; /*  input file record #  starting at the point where the model reads values */
01165 
01166      while ( fgets(lline, 2000, F_struct_TS) != NULL && !feof( F_struct_TS ) )
01167      {
01168          line = lline;
01169          sscanf( line, "%d %d %d",&yyyy,&mm,&dd); /* read the date of the current record */
01170          
01171         /* julian day, returns a double (args = mon, day, year, hours, minutes, seconds) */
01172          Jdate_read = julday( mm,dd,yyyy,0,0,0.0 ); 
01173          if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
01174              printf (" \n***Error\nStructure Tracer data (%s) init date > simulation start date\n", modelFileName); exit (-1);
01175          }
01176          if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
01177          {
01178              line = Scip( line, TAB );
01179              line = Scip( line, TAB );
01180              for ( j = 0; j < numTSser; j++ ) { /* loop through all the control structure columns in the input file */
01181                  line = Scip( line, TAB );
01182                  if (*line == TAB) Hist_data[IDhist[j]].arrayS[((int)(lincnt))] = 0;
01183                  else sscanf( line, "%f", &Hist_data[IDhist[j]].arrayS[((int)(lincnt))] );
01184                  Hist_data[IDhist[j]].arrayS[((int)(lincnt))] *= 1.0; /* data read in g/L=kg/m3, no conversion*/
01185                  
01186              }
01187               
01188              lincnt++; /* increment the number of daily data records */
01189          }
01190 
01191      } /* end of reading tracer (salt) conc data records */
01192      
01193      if (lincnt == 0 ) {
01194              printf (" \n***Error\nTracer Time series data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
01195      }
01196      else if (lincnt != PORnumday ) {
01197              printf (" \n***Error\nTracer Time series data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
01198      }
01199      sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
01200      
01201      fclose (F_struct_TS);
01202  } /* end of reading trace (salt) time series data (if needed) */
01203  
01204 
01205  
01206  sprintf(msgStr, "\t**** Evaluate scenario: %s %s ****", &SimAlt, &SimModif ); usrErr(msgStr); WriteMsg(msgStr,1);
01207 
01208  return;
01209 }
01210 
01211 
01212 /****************************************************************************************/
01222 struct Schedule *Read_schedule ( char *sched_name, char *filename, float BASE_DATUM )
01223 { 
01224 char ss[301], *s, scenario[20], modnam[20];
01225 int i, count = 1, sched_found=0;
01226 struct Schedule *Sce_Graph;
01227 
01228    { if(H_OPSYS==UNIX) 
01229                 sprintf( modelFileName, "%s/%s/Data/%s.graph", ModelPath, ProjName, filename );
01230      else sprintf( modelFileName, "%s%s:Data:%s.graph", ModelPath, ProjName, filename );
01231    }
01232 
01233 /* Open file with schedule graphs */
01234   if ( ( schedFile = fopen( modelFileName, "r" ) ) ==  NULL )
01235         {
01236        printf( "Can't open %s file!\n", modelFileName ) ;
01237        exit(-1) ;
01238     }
01239 
01240  fgets( ss, 300, schedFile );  /* skip comment line */
01241  fgets( ss, 300, schedFile );  /* skip comment line */
01242  fgets( ss, 300, schedFile );  /* skip comment line */
01243  fgets( ss, 300, schedFile ); sscanf( ss,"%s %s", &modnam, &scenario );
01244  if (strcmp(scenario,SimAlt) != 0) {
01245      sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
01246              scenario, modelFileName, SimAlt); usrErr(msgStr);
01247              exit(-1);
01248  }
01249  if (strcmp(modnam,modelName) != 0) {
01250      sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
01251              modnam, modelFileName, modelName); usrErr(msgStr);
01252              exit(-1);
01253  }
01254 
01255 /* Read schedule descriptors until none left */   
01256   while ( fgets( ss, 300, schedFile ) != NULL && !feof( schedFile ) )
01257         { if ( !Wrdcmp ( sched_name, ss ) ) continue;
01258           sched_found = True;
01259           
01260           fgets( ss, 300, schedFile );  
01261           s = ss;
01262           /* count the number of pairs for the graph */
01263           while ( *s++ != EOL ) if ( *s == '(' ) count++;    
01264 
01265           /* Allocate memory for next schedule */
01266       if ( (Sce_Graph = (struct Schedule *) malloc( (size_t) sizeof( struct Schedule ))) == NULL )
01267        {
01268          printf( "Failed to allocate memory for Managed flow Schedule graph in (%s)\n ", sched_name ) ;
01269          exit( -2 ) ;
01270        };
01271      
01272           /* Allocate memory for next schedule graph containing 'count' points*/
01273       if ( (Sce_Graph->graph_points = 
01274                 (struct Points *) malloc( (size_t) sizeof( struct Points ) * count)) == NULL )
01275        {
01276          printf( "Failed to allocate memory for Managed flow Schedule graph in (%s)\n ", sched_name ) ;
01277          exit( -2 ) ;
01278        };
01279           
01280           Sce_Graph->num_point = count;
01281           
01282           /* Read the graph points */
01283           s = ss;
01284           for ( i = 0; i < count && *s!=EOS; i++ )
01285           { 
01286                 s = Scip( s,'(');       
01287                 sscanf( s, "%f", &Sce_Graph->graph_points[i].time ); 
01288                 s = Scip ( s, ',');
01289                 sscanf( s, "%f", &Sce_Graph->graph_points[i].value );
01290                 Sce_Graph->graph_points[i].value += BASE_DATUM;
01291           }
01292         } /* end of while */
01293         /* v2.5 included this error-trap */
01294         if (!sched_found) {
01295              sprintf(msgStr, "The schedule name (%s) was not found in the %s schedule file! (Was named in the CanalData.struct file.)",
01296              sched_name, modelFileName); usrErr(msgStr);
01297              exit(-1);
01298     }
01299 
01300         
01301   fclose (schedFile);
01302   return Sce_Graph;
01303 }
01304 
01305 /********************************************************************************************/
01317 void Channel_configure (  float *ElevMap, struct Chan *channel_first ) 
01318 { int i, j;                                                       /*coord. of current cell crossed by the canal */
01319  int  cellLoc, cellLoc_i, cellLoc_j;
01320  int HV, k, L_range, R_range;
01321  float T_length;                                                /* length of a straight segment of a canal reach (m) */
01322  float C_length;                                          /*  total length of all of the segments of a canal reach (m) */
01323  float distTot;                                       /*  cumulative (along grid cells) distance along a canal reach, from the starting point (m) */
01324  float A1, B1, A2, B2;                                                   /* parameters of the canal equation */
01325  float L, p_middle, m, i45, x, y, x1, y1, length, RoIL, RoIR;
01326  float cell_step = 1.0;
01327    
01328  float a0, b0, a1, b1;
01329  int c_num, z;
01330  struct Cells  *cell; /* addresses for adjacent cell structures */
01331  struct Chan *channel;
01332  struct Chan_reach *Reach;
01333  
01334  int numChek; /* placeholder for debugging cell elev along canals */
01335    
01336 
01337  int *marked; /* a "temporary" fix used to determine if a cell belonging to a canal has already been marked */
01338  marked = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"marked");
01339  init_pvar(marked,NULL,'d',0);
01340 
01341 /* Open file to print canal-cell interaction info to */
01342  sprintf( modelFileName, "%s/%s/Output/Debug/CanalCells_interaction.txt", OutputPath, ProjName );
01343  if ( ( CanalCellInteractDebug = fopen( modelFileName, "w" ) ) ==  NULL )
01344  {
01345      printf( "Can't open %s file!\n", modelFileName ) ;
01346      exit(-1) ;
01347  }
01348 
01349      /* allocate memory for array of pointers to Chan structs */
01350  if ( (Chan_list = 
01351        (struct Chan **) malloc( (size_t) sizeof( struct Chan *) * num_chan)) == NULL )
01352  {
01353      printf( "Failed to allocate memory for next canal (%d)\n ", channel->number ) ;
01354      exit( -2 ) ;
01355  };
01356 
01357      /* allocate memory for first cell structure */
01358  if ( (cell = (struct Cells *) malloc( (size_t) sizeof( struct Cells ))) == NULL )
01359  { 
01360      printf( "Failed to allocate memory for cell structure\n " ) ;
01361      exit( -2 ) ;
01362  }
01363 
01364      /* store the pointer to first cell in the array */
01365  channel = channel_first;
01366 
01367      /* loop over all the canal reaches (== channel, defined as an individual canal reach with an ID number) */
01368  while (channel != NULL)
01369  {
01370      Reach = channel->reaches;
01371      c_num = channel->number;
01372      z  = channel->levee;
01373      RoIL = channel->roil;
01374      RoIR = channel->roir;
01375         
01376      Chan_list[c_num - CHo] = channel;
01377         
01378      channel->cells = cell;
01379      distTot = 0;
01380      C_length = 0;
01381      C_Mark = 0;
01382     
01383          /* a neg width canal is skipped, to SKIPCHAN near end of this cycle along canals */
01384      if (channel->width < 0)   goto SKIPCHAN; 
01385     
01386      fprintf ( CanalCellInteractDebug, "N = %d  levee = %d\n", c_num, z );
01387      
01388 /* get the elevation at the starting reach segment (may be up or down stream) of the canal */
01389      cellLoc_i=(int)(Reach->x0+1);
01390      cellLoc_j=(int)(Reach->y0+1);
01391      cellLoc = T(cellLoc_i,cellLoc_j );
01392      Chan_list[c_num - CHo]->elev_start = ElevMap[cellLoc]; 
01393      if ( !ON_MAP[cellLoc] ) {
01394         printf( "Starting elevation of Reach %d is out of the active domain.  Fix the reach location.\n",c_num);
01395         exit(-2);
01396      }
01397      fprintf ( CanalCellInteractDebug, " %d, %d, StartElev = %f\n", cellLoc_i,cellLoc_j,  Chan_list[c_num - CHo]->elev_start );
01398 
01399                 
01400      while (Reach != NULL)  /* loop through canal reach segments */
01401      {
01402          a0 = Reach->x0;
01403          b0 = Reach->y0;
01404          a1 = Reach->x1;
01405          b1 = Reach->y1;
01406 
01407          x1 = 0;  y1 = 0;
01408 
01409          T_length = sqrt( (b1 - b0)*(b1 - b0) + (a1 - a0)*(a1 - a0) );  /* total reach segment length */
01410          C_length += T_length;                                                                              /* total canal reach length */
01411 
01412          if ( a1 != a0 && b1 != b0 ) 
01413          {
01414              A1 = (b1 - b0)/(a1 - a0);                      /* define the equation: y = A1x + B1 */
01415              B1 = b0 - A1*a0;
01416              A2 = 1/A1;  B2 = - B1/A1;                      /* define the equation: x = A2y + B2 */
01417          } 
01418     
01419          x = a0; y = b0;
01420          i = Round (x);    j = Round (y);   
01421               
01422              /* loop along cells of the canal reach segment */ 
01423          while ( x1 != a1 || y1 != b1 )
01424          {
01425              if ( a1 == a0 )                                     /* if goes straight horizontally */
01426              {  x1 = a0; 
01427              y1 = (j == Round (b1) ) ? b1 : j + STEP(b1 - b0)*cell_step; 
01428              length = y1 - y;
01429              }
01430        
01431              else if ( b1 == b0 )                                   /* if goes straight vertically */
01432              { x1 = (i == Round (a1) ) ? a1 : i + STEP(a1 - a0)*cell_step; 
01433              y1 = b0; 
01434              length = x1 - x;
01435              }
01436 
01437              else if ( i == Round (a1) && j == Round (b1) )
01438              { x1 = a1;  y1 = b1;                                               /* if last cell */ 
01439              length = T_length * (y1 - y)/(b1 - b0);
01440              C_Mark = 0;
01441              }
01442           
01443              else
01444              { y1 = A1*(i + STEP(a1 - a0)*cell_step) + B1;  length = T_length * (y1 - y)/(b1 - b0);
01445              x1 = A2*(j + STEP(b1 - b0)*cell_step) + B2;  L = T_length * (x1 - x)/(a1 - a0);
01446     
01447              if ( length < L )                       /* define the end coordinates and length */
01448                  x1 = i + STEP(a1 - a0)*cell_step;  /* (the next coordinate is the closest)*/
01449              else
01450              {  y1 = j + STEP(b1 - b0)*cell_step;  length = L;
01451              }
01452              }
01453 
01454              if ( Abs (a1 - a0) < Abs (b1 - b0) )
01455                      /* define the direction of canal reach and  */
01456              { p_middle = (x1 + x)/2;      /* calculate the middle point on the canal reach in cell */
01457              m = i + cell_step/2;                                                               /* m - the center of the cell */
01458              HV = 1;  /* to indicate that the canal reach goes closer to the horizontal direction */
01459              }
01460              else
01461              { p_middle = (y1 + y)/2;
01462              m = j + cell_step/2;
01463              HV = 0;
01464              }
01465              
01466                  distTot = distTot + length*celWid;  /*  cumulative (along grid cells) distance along a canal reach, from the starting point (m) */
01467          
01468                  /* mark cells depending on where the center of canal reach is rel. to the cell center*/
01469 
01470              if ( Round (a1) == Round (a0) ) /* if the canal reach goes within the range of one */
01471                      /*cell horizontally, mark the two vertical cells */
01472              { if (b0 < b1)     /* if going from left to right */
01473              {          /* if the middle point of the canal reach is lower than the */
01474                      /* cell center, mark the cell as LEFT else as RIGHT */
01475                  if ( p_middle > m ) 
01476                  { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01477                  L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01478                  }  
01479                  else 
01480                  { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i - 1, j, c_num, marked, distTot );
01481                  L_range = (int) RoIL;  R_range = (int)(RoIR - 1);
01482                  }  
01483                  for ( k = 0; k < L_range; k++ )
01484                      cell = MarkCell ( i - k - 1, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01485                  for ( k = 0; k < R_range; k++ )
01486                      cell = MarkCell ( i + k + 1, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01487              }
01488              else  /* vice versa if we go from right to left */
01489              {          /* if the middle point of the canal reach is lower than the */
01490                      /* cell center, mark the cell as RIGHT else as LEFT */
01491                  if ( p_middle > m ) 
01492                  { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01493                  L_range = (int)RoIL; R_range = (int)(RoIR - 1);
01494                  } 
01495                  else 
01496                  { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i - 1, j, c_num, marked, distTot );
01497                  L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01498                  }
01499                  for ( k = 0; k < L_range; k++ )
01500                      cell = MarkCell ( i + k + 1, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01501                  for ( k = 0; k < R_range; k++ )
01502                      cell = MarkCell ( i - k - 1, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01503              }
01504              }
01505              else if ( Round (b1) == Round (b0) )       /* if the canal reach goes vertically */
01506              { if (a0 < a1) /* if going from top to bottom */
01507              {      /* if the middle point of the canal reach is further to the right than the */
01508                      /* cell center, mark the cell as RIGHT else as LEFT */
01509                  if ( p_middle > m ) 
01510                  { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01511                  L_range = (int)RoIL; R_range = (int)(RoIR -1);
01512                  }
01513                  else 
01514                  { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j - 1, c_num, marked, distTot );
01515                  L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01516                  }
01517                  for ( k = 0; k < L_range; k++ )
01518                      cell = MarkCell ( i, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01519                  for ( k = 0; k < R_range; k++ )
01520                      cell = MarkCell ( i, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01521              }
01522              else  /* vice versa if we go from bottom to top */
01523              {      /* if the middle point of the canal reach is further to the left than the */
01524                      /* cell center, mark the cell as LEFT else as RIGHT */
01525                  if ( p_middle > m ) 
01526                  { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01527                  L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01528                  }
01529                  else 
01530                  { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j - 1, c_num, marked, distTot );
01531                  L_range = (int)RoIL; R_range = (int)(RoIR - 1);
01532                  }
01533                  for ( k = 0; k < L_range; k++ )
01534                      cell = MarkCell ( i, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01535                  for ( k = 0; k < R_range; k++ )
01536                      cell = MarkCell ( i, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01537              }
01538              }
01539              else
01540                      /* if the canal reach crosses more than one cell (in both dimensions), mark diagonal cells */
01541              { if ( A1>0 )                 /* this is an extremely confusing check, but in this */
01542                  i45 = ( HV ) ? p_middle - m : m - p_middle;  /*form it seems to be working OK */
01543              else 
01544                  i45 = p_middle - m;
01545                  /*  b0>b1, a0>a1 \/ b0<b1, a0>a1 */
01546                  /*  b0>b1, a0<a1 /\ b0<b1, a0<a1 */
01547              if ( b0 < b1 )
01548              { if ( a0 < a1 )     /* when a0<a1 the canal reach goes down in the \ direction */
01549              { if ( i45 > 0 ) 
01550              { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01551              cell = MarkCell ( i + 1, j - 1, RIGHT, length, cell, SOUTH, z, i + 1, j - 1, c_num, marked, distTot );
01552              }
01553              else 
01554              { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01555              cell = MarkCell ( i - 1, j + 1, LEFT, length, cell, EAST, z, i - 1, j + 1, c_num, marked, distTot );
01556              }
01557              for ( k = 0; k < RoIR-1; k++ )
01558                  cell = MarkCell ( i + k + 1, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01559              for ( k = 0; k < RoIL-1; k++ )
01560                  cell = MarkCell ( i - k - 1, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01561              }
01562              else                                 /* in this case the canal reach goes up in the / direction */
01563              { if ( i45 > 0 ) 
01564              { cell = MarkCell ( i, j, LEFT, length, cell, NONE, z, i, j, c_num, marked, distTot );
01565              cell = MarkCell ( i + 1, j + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01566              }
01567              else 
01568              { cell = MarkCell ( i, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01569              cell = MarkCell ( i - 1, j - 1, LEFT, length, cell, NONE, z, i - 1, j - 1, c_num, marked, distTot );
01570              }
01571              for ( k = 0; k < RoIR-1; k++ )
01572                  cell = MarkCell ( i + k + 1, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01573              for ( k = 0; k < RoIL-1; k++ )
01574                  cell = MarkCell ( i - k - 1, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01575              }
01576              }
01577                  /* if going in the opposite direction change left and right */            
01578              else
01579              { if ( a0 > a1 )       /* when a0<a1 the canal reach goes up in the \ direction */
01580              { if ( i45 < 0 ) 
01581              { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01582              cell = MarkCell ( i - 1, j + 1, RIGHT, length, cell, EAST, z, i - 1, j + 1, c_num, marked, distTot );
01583              }
01584              else 
01585              { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01586              cell = MarkCell ( i + 1, j - 1, LEFT, length, cell, SOUTH, z, i + 1, j - 1, c_num, marked, distTot );
01587              }
01588 
01589              for ( k = 0; k < RoIL-1; k++ )
01590                  cell = MarkCell ( i + k + 1, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01591              for ( k = 0; k < RoIR-1; k++ )
01592                  cell = MarkCell ( i - k - 1, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01593              }
01594              else                               /* in this case the canal reach goes down in the / direction */
01595              { if ( i45 < 0 ) 
01596              { cell = MarkCell ( i, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01597              cell = MarkCell ( i - 1, j - 1, RIGHT, length, cell, NONE, z, i - 1, j - 1, c_num, marked, distTot );
01598              }
01599              else 
01600              { cell = MarkCell ( i, j, RIGHT, length, cell, NONE, z, i, j, c_num, marked, distTot );
01601              cell = MarkCell ( i + 1, j + 1, LEFT, length, cell, ALL, z, i + 1, j + 1, c_num, marked, distTot );
01602              }
01603 
01604              for ( k = 0; k < (int)(RoIL-1); k++ )
01605                  cell = MarkCell ( i + k + 1, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01606              for ( k = 0; k < (int)(RoIR-1); k++ )
01607                  cell = MarkCell ( i - k - 1, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01608              }
01609              }
01610              }
01611          
01612              num_cell++;        /*at this point we're counting the cells only once thinking that all 
01613                                   cells have the same length of interaction irrelevant of their distance from the canal 
01614                                   and that this length will be equal to the canal length / number of cells it crosses */
01615 
01616                  /* prepare for the next step in cycle */
01617              x = x1;  y = y1;
01618              i = Round (x);     if ( a1 < a0 && i == x ) i--;         /* if we're going upwards, we */
01619                  /* need to change the rule of Rounding so that we move ahead along cells */
01620              j = Round (y);     if ( b1 < b0 && j == y ) j--;
01621              C_Mark = 1;
01622     
01623 
01624          } /* end loop along cells of the canal reach segment  */
01625          
01626          Reach = Reach->next_reach;
01627      } /* end loop over all reach segments of a canal */
01628 
01629 /* get the elevation at the ending reach segment (may be up or down stream) of the canal */
01630      cellLoc_i=(int)(a1+1);
01631      cellLoc_j=(int)(b1+1);
01632      cellLoc = T(cellLoc_i,cellLoc_j );
01633      Chan_list[c_num - CHo]->elev_end = ElevMap[cellLoc];  
01634 
01635      Chan_list[c_num - CHo]->num_of_cells = num_cell;
01636      num_cell = 0;
01637          /* store the canal attributes of:
01638             length of entire canal, 
01639             area of entire canal, 
01640             minimum volume allowed in canal,
01641             avg length of cell along reach segment,
01642             overland flow coefficient (not incl. dynamic manning's n, C_F is for sensitivity experiments only (normally=1) ),
01643             seepage flow coefficient (including GP_calibGWat groundwater flow calibration coef),
01644             land surface elevation difference between starting and ending points of canal 
01645             slope of the elevation gradient from start to end points */
01646      Chan_list[c_num - CHo]->length = C_length*celWid;
01647      Chan_list[c_num - CHo]->area = Chan_list[c_num - CHo]->width*Chan_list[c_num - CHo]->length;
01648      Chan_list[c_num - CHo]->minVol = Chan_list[c_num - CHo]->area*MINDEPTH;
01649      Chan_list[c_num - CHo]->seg_len = Chan_list[c_num - CHo]->length/Chan_list[c_num - CHo]->num_of_cells;
01650      Chan_list[c_num - CHo]->SW_flow_coef = sqrt(Chan_list[c_num - CHo]->seg_len) * sec_per_day * C_F;
01651      Chan_list[c_num - CHo]->SPG_flow_coef = Chan_list[c_num - CHo]->cond * GP_calibGWat;
01652      Chan_list[c_num - CHo]->elev_drop = Chan_list[c_num - CHo]->elev_start - Chan_list[c_num - CHo]->elev_end;
01653      Chan_list[c_num - CHo]->elev_slope = Chan_list[c_num - CHo]->elev_drop / Chan_list[c_num - CHo]->length ; /* v2.5 calculated slope */
01654      
01655      if ( !ON_MAP[cellLoc] ) {
01656         printf( "Ending elevation of Reach %d is out of the active domain.  Fix the reach location.\n",c_num);
01657         exit(-2);
01658      }
01659      fprintf ( CanalCellInteractDebug, " %d, %d, EndingElev = %f, ElevSlope= %f, ReachLength= %f\n", cellLoc_i,cellLoc_j,  
01660                                         Chan_list[c_num - CHo]->elev_end, Chan_list[c_num - CHo]->elev_slope, Chan_list[c_num - CHo]->length );
01661      
01662      getCanalElev(c_num - CHo); /* v2.5 */
01663      
01664    SKIPCHAN: /* a single goto comes here, skipping a canal with negative width */
01665      channel = channel->next_in_list;
01666 
01667      cell_last->next_cell = NULL;                       /* close the cell list */
01668 
01669  } /* end loop loop over all the canal reaches */
01670  
01671  fclose ( CanalCellInteractDebug );
01672  
01673  return; 
01674 }
01675 
01676 /********************************************************************************************/
01687 void getCanalElev (int chan_n)
01688 {
01689      int i;
01690          /* cycle along canal cells */
01691      cell = Chan_list[chan_n]->cells;
01692      i = T(cell->x, cell->y);
01693      
01694      while ( cell != NULL )
01695      { 
01696         cell->reachElev = Chan_list[chan_n]->elev_start - cell->reachDistDown * Chan_list[chan_n]->elev_slope;
01697         fprintf ( CanalCellInteractDebug, " Reach= %d, Cell (%d, %d), ReachElev= %f\n", chan_n+CHo, cell->x, cell->y, cell->reachElev );
01698         cell = cell->next_cell;
01699      }
01700 }
01701 
01702 
01703 /********************************************************************************************/
01728 struct Cells *MarkCell (int x, int y, int c_ind, float length, struct Cells *cell, int direction, int levee, 
01729                                         int xm, int ym, int c_num, int *marked, float distTot )
01730 { struct Cells *cell_next;
01731   int N, cellLoc;
01732   
01733    if ( x >= s0 || y >= s1 )
01734      {  printf ( "Error in definition of canal in cell x=%d  y=%d ", x, y );
01735         exit (-2);                      /*check that cells obtained are within array boundaries */
01736      }
01737 
01738    N = T(xm+1, ym+1);
01739    if ( !ON_MAP[N] ) return cell;
01740    
01741    if ( levee && direction != ALL )
01742    ON_MAP[N] = ((ON_MAP[N]-1) & (direction+100)) + 1;   /* modifying ON_MAP to take into account the 
01743    allowed direction of fluxing: 00 - none, 01=1 - to the East, 10=2 - to the South, 11=3 - all 
01744    ON_MAP they naturally become 1 - none, 2 -to the East, 3 - to the South, 4 - all directions 
01745    0 is reserved for cells out of the ON_MAP area*/
01746    
01747    if ( !C_Mark ) return cell; /* C_Mark is to make gaps between canals, so that there will be no 
01748                 canals being connected through cells if it turns out that they are linked to the same cells */
01749              
01750        /* a temporary fix used to determine if a cell belonging to a canal has already been marked */
01751        /***********/
01752    cellLoc = T(x+1,y+1);
01753    if (marked[cellLoc] == c_num ) return cell;
01754    marked[cellLoc] = c_num;
01755        /***********/
01756              
01757 
01758    
01759    cell->x = x+1;                                               /* shift coord to match the other maps */
01760    cell->y = y+1;
01761    cell->ind = c_ind;   
01762    cell->length = length*celWid;                /* length of canal associated with this cell * cell width m */
01763    cell->reachDistDown = distTot; /*  cumulative (by grid cells) distance along a canal reach, from the starting point (m) (v2.5) */
01764    
01765    fprintf ( CanalCellInteractDebug, " %d, %d, l/r = %d, distance= %7.1f\n", cell->x, cell->y, c_ind, cell->reachDistDown );
01766    
01767    /* allocate memory for next cell structure */
01768    if ( (cell_next = (struct Cells *) malloc((size_t) sizeof( struct Cells ))) == NULL )
01769        {
01770          printf( "Failed to allocate memory for cell structure\n " ) ;
01771          exit( -2 ) ;
01772        }
01773    cell->next_cell = cell_next;
01774    cell_last = cell;
01775 
01776    return cell_next;
01777 }
01778 
01779 /********************************************************************************************/
01805 void  FluxChannel ( int chan_n, float *SWH, float *ElevMap, float *MC, 
01806                     float *GWH, float *poros, float *GWcond,
01807                     double *NA, double *PA, double *SA, double *GNA, double *GPA, double *GSA,
01808                     float *Unsat, float *sp_yield)
01809 { 
01810                 
01842 int converged = 0;                                              
01843  float error = 1.0;                                     
01844  float prev_error = 1.0;                                        
01845  int iter = 0;                                                  
01846  int max_iter = CHANNEL_MAX_ITER;               
01847  float init_factor = 0.5;                               
01848  float factor = init_factor;                    
01849 
01850  int i;
01851  int attempt=0;                                                 
01852  double CanWatDep = Chan_list[chan_n]->wat_depth;
01853  double T_flux_S, T_flux_G, T_flux_L;           
01854  double fluxO, fluxCk;
01855  double fluxS, fluxG, fluxL; 
01856  double SW_coef, SPG_coef, GW_coef;     
01857  double N_fl, P_fl, S_fl;                               
01858  /* double SNfl, GNfl, SPfl, GPfl, SSfl, GSfl; UNUSED */
01859  double I_Length;                                               
01860  double seg_area;                                               
01861  double MIN_VOL;                                                
01862  double GW_head;                                                
01863  double h_GWflow, h_SPflow;                             
01864  double dh;                                                     
01865  double Qin, Qout, Qout1;                               
01866  double SWater;                                                 
01867  double tot_head;                                               
01868  double CH_vol;                                         
01869  double CH_vol_R;                                       
01870  double CH_bottElev;                                    
01871  double H_rad_ch, H_rad_cell;                   
01872  double can_cell_M, can_cell_A;                 
01873  float errorConv;                                               
01874 
01875 /* UNUSED below, unimplemented gwater/sfwater integration vars */ 
01876     int equilib, revers=1;
01877     float H_pot_don, H_pot_rec;
01878     
01879     float fluxTOunsat_don, fluxHoriz; /* vertical and horiz components of the calculated total groundwater Flux */
01880     float UnsatZ_don, UnsatCap_don, UnsatPot_don, Sat_vs_unsat; /* unsaturated zone capacities, potentials, in the donor cell */
01881     float UnsatZ_rec, UnsatCap_rec, UnsatPot_rec; /* unsaturated zone capacities, potentials, in the recipient cell */
01882     float sfTOsat_don, unsatTOsat_don, unsatTOsat_rec, satTOsf_rec, horizTOsat_rec; /* fluxes for vertical integration */
01883     
01884     float sedwat_don, sedwat_rec;
01885     float sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec; /* vertical flux of solutes */
01886     float sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don; /* vertical flux of solutes */
01887     float sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz; /* horiz flux of solutes */
01888 /* UNUSED above, unimplemented gwater/sfwater integration vars */ 
01889  
01890  
01891 /* look at all the structures connected to the canal and add/subtract flow from/to them */
01892  Qin = 0; Qout = 0;
01893  structs = struct_first;
01894  while ( structs != NULL )
01895  {
01896      
01897      if ( structs->canal_fr  == chan_n+CHo )  
01898      {
01899          if ( structs->flow > 0 ) Qout += structs->flow;
01900          else Qin -= structs->flow;
01901      } 
01902      if ( structs->canal_to == chan_n+CHo ) 
01903      {
01904          if ( structs->flow > 0 ) Qin += structs->flow;
01905          else Qout -= structs->flow;
01906      } 
01907     
01908      structs = structs->next_in_list;
01909  }
01910   
01911  CH_vol_R = Chan_list[chan_n]->area * CanWatDep;  /* total canal volume prior to relaxation procedure */
01912  MIN_VOL =  Chan_list[chan_n]->minVol;                          /* minimum canal volume allowed */
01913  SPG_coef = Chan_list[chan_n]->SPG_flow_coef;   /* GP_calibGWat already in this */
01914 
01915  I_Length =   Chan_list[chan_n]->seg_len;          /* avg length of cell along reach segment */
01916  seg_area = I_Length * Chan_list[chan_n]->width;  /* avg area of reach segment */
01917  can_cell_M = Chan_list[chan_n]->area*CELL_SIZE;  /* used in reducing excess flows between canal<->cell */
01918  can_cell_A = Chan_list[chan_n]->area+CELL_SIZE;  /* used in reducing excess flows between canal<->cell */
01919  
01920 
01921  
01922 /* start relaxation procedure */  
01923  do
01924  {
01925      if   ( Abs( error) < F_ERROR ) 
01926      {
01927              converged = 1; /* when we've converged, we'll go through the cell list one more time w/o modifying the depth estimate */
01928                         /* occasionally, we get a modified total flux value, and diff error, a second time with the same depth estimate! */
01929              errorConv=error;
01930              factor = 0;
01931      }
01932      else if  ( iter >= (max_iter-1) )
01933      {
01934          if (attempt == 4) {
01935              converged = 1; /* haven't actually converged, but stop now */
01936              sprintf(msgStr,"Day %6.1f: ERROR - couldn't converge on Canal %d, HALTING the iterative relaxation after 4 tries, error = %f.", 
01937                         SimTime.TIME, Chan_list[chan_n]->number, error); 
01938              WriteMsg( msgStr,True ); dynERRORnum++;
01939          }
01940          else {
01941              if (debug > 3) { sprintf(msgStr,"Day %6.1f: Warning - couldn't converge on Canal %d, redoing the iterative relaxation, error = %f.", 
01942                         SimTime.TIME, Chan_list[chan_n]->number, error); 
01943              WriteMsg( msgStr,True ); }
01944 
01945              attempt++; /* the attempt counter increment */
01946              iter = 0; /* reset iterations to 0, and double factor & max # iterations on first try */ 
01947              if (attempt == 1) { factor = init_factor*2.0;  max_iter = CHANNEL_MAX_ITER * 2; }
01948              else if (attempt == 2) { factor = init_factor*0.5; max_iter = CHANNEL_MAX_ITER * 2; }
01949              else if (attempt == 3) { factor = init_factor*4.0; max_iter = CHANNEL_MAX_ITER * 4; }
01950              else if (attempt == 4) { factor = init_factor*0.05; max_iter = CHANNEL_MAX_ITER * 4; }
01951          }
01952      }
01953 
01954      else
01955      {
01956              /* if error changed sign, we have overshot the target, therefore reverse direction */
01957          if (error * prev_error < 0 && iter != 1) 
01958              factor = - sgn(error) *  Abs( factor ) * 0.5;      /* and reduce the step */
01959          
01960              /* if error contiues to grow after first iter.(#0), reverse direction */   
01961          if (iter == 1 && error > 0) factor = - factor;
01962          
01963      }
01964      
01965      if (iter != 0 && !converged ) 
01966      {
01967          CanWatDep += factor; 
01968          prev_error = error;
01969          if ( CanWatDep < MINDEPTH ) 
01970          { 
01971              if (debug > 3) { sprintf(msgStr,"Day %6.1f: Warning - estimate is below minimum depth, revising estimate. (Canal %d, depth %f)", 
01972                         SimTime.TIME, Chan_list[chan_n]->number,CanWatDep); 
01973              WriteMsg( msgStr,True ); }            
01974 
01975              CanWatDep = MINDEPTH; 
01976             /* converged = 1; */
01977             factor = -factor*0.1; /* reverse the direction and decrease for next estimate*/
01978          }      
01979      }     
01980     
01981      T_flux_G = T_flux_S = T_flux_L = 0.0;           /* at each convergence iter, set the sum flows to 0 */
01982      CH_vol = Chan_list[chan_n]->area * CanWatDep;  /* volume of canal at the estimated water depth */
01983     
01984          /* cycle along canal cells */
01985      cell = Chan_list[chan_n]->cells;
01986      
01987      while ( cell != NULL )
01988      { 
01989          i = T(cell->x, cell->y);
01990          if ( !ON_MAP[i] )       /* check that cell is ON_MAP, else we do nothing */
01991                     /* return to beginning of do loop                */
01992          {  cell = cell->next_cell;  continue;
01993          }
01994          if (basn[i] == 0) {
01995              sprintf(msgStr,"ERROR - Cell (%d, %d) is out-of-system.", cell->x, cell->y); 
01996              WriteMsg( msgStr,True );  dynERRORnum++; }
01997 
01998          fluxS = 0.;  fluxG = 0.; fluxL = 0.; /* overland, subsurface, and seepage cell-canal fluxes set to 0 */
01999            /* v2.5 In the case of some canal reaches, there is a lip or berm, and often associated vegetation, that is 
02000            along the edge of a canal, decreasing flow rates between canal & marsh.  If user provides a non-zero roughness
02001            parameter in the CanalData.chan input file, the the mean roughness of this berm and the cell is used for 
02002            surface exchanges */
02003          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] );
02004          GW_coef = GWcond[i] * poros[HAB[i]];   
02005          SWater = SWH[i];
02006 
02007          if (debug>0 && (SWater < -GP_MinCheck) ) {
02008          sprintf(msgStr,"Day %6.1f: capacityERR - neg sfwater along canal %d, candepth %f, cell (%d, %d), celldepth %f", 
02009                  SimTime.TIME,Chan_list[chan_n]->number,CanWatDep, cell->x, cell->y, SWater); 
02010          WriteMsg( msgStr,True );  dynERRORnum++; }
02011          
02012          GW_head = GWH[i]/poros[HAB[i]];
02013          tot_head = SWater + ElevMap[i]; 
02014          CH_bottElev = cell->reachElev - Chan_list[chan_n]->depth;  /* elev of bottom of canal at cell location (v2.5 is from fixed vector slope instead of actual grid cell elev) */
02015              
02016          dh = ( CH_bottElev + CanWatDep ) - tot_head;
02017 
02018              /* hydraulic radius of canal for overland flow out of canal */
02019          H_rad_ch =   Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
02020                         SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);   /* v2.5 contrain to >= 0 */
02021             /* hydraulic radius of cell for overland flow into canal (same eqn right now) */
02022          H_rad_cell = Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
02023                         SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);     /* v2.5 contrain to >= 0 */
02024 
02025                 /* determine the heights of the water cross sections associated with the subsurf and seep flows */
02026          if (dh > 0.0) { /* from canal */
02027                 h_GWflow = Min(Chan_list[chan_n]->depth, CanWatDep);
02028                 h_SPflow = Max(CH_bottElev + CanWatDep - ElevMap[i], 0.0);
02029          }
02030          else if (dh < 0.0) { /* from cell */
02031                 h_GWflow = Max(GW_head-CH_bottElev, 0.0);
02032                 h_SPflow = Max(tot_head-ElevMap[i], 0.0);
02033          }
02034 
02035                  
02036          switch( Chan_list[chan_n]->levee )
02037          {      
02038                 
02039              case 2:    /* if levees are on both sides of the canal */
02040                          /*seepage and subsurface flow on both sides */
02041                  fluxL =  (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02042                  fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
02043                  break;
02044              case 0: /* if there are no levees */
02045                      /*overland flow either side */
02046                  if ( dh > 0 ) {
02047                          /* flux from canal to cell  */
02048                      fluxS = f_Manning ( dh, H_rad_ch, SW_coef);
02049                  }
02050                  else if (SWater>GP_DetentZ) {
02051                      fluxS =  f_Manning ( dh, H_rad_cell, SW_coef) ; 
02052                     /* fluxS =  ( (-fluxS) > (SWater-GP_DetentZ) *CELL_SIZE )  ? ( -(SWater-GP_DetentZ)*CELL_SIZE ) : (fluxS);*/
02053                     if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02054                         fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02055                  }
02056           
02057                                 /* subsurface groundwater flow both sides */
02058                  fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
02059                  break;
02060              case 1:  /* if levee is on the left */ 
02061                  if ( cell->ind < 0  )
02062                  { /*overland flow */
02063                      if ( dh > 0 ) {
02064                              /* flux from canal to cell  */
02065                          fluxS = f_Manning ( dh, H_rad_ch, SW_coef);
02066                      }
02067                      else if (SWater>GP_DetentZ) {
02068                          fluxS = f_Manning ( dh, H_rad_cell, SW_coef) ; 
02069                         /* fluxS =  ( (-fluxS) > (SWater-GP_DetentZ) *CELL_SIZE )  ? ( -(SWater-GP_DetentZ)*CELL_SIZE ) : (fluxS);*/
02070                     if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02071                         fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02072                      }
02073           
02074                  }      
02075                  else 
02076                          /*seepage on the levee side */
02077                      fluxL =  (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02078                      
02079                                 /* subsurface groundwater flow both sides */
02080                  fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
02081                  break;
02082              case -1: /* if levee is on the right */
02083                  if ( cell->ind > 0 )
02084                  { /*overland flow */
02085                      if ( dh > 0 ) {
02086                              /* flux from canal to cell */
02087                          fluxS = f_Manning( dh, H_rad_ch, SW_coef);
02088                      }
02089                      else if (SWater>GP_DetentZ) {
02090                          fluxS = f_Manning ( dh, H_rad_cell, SW_coef); 
02091                          /*fluxS =  ( (-fluxS) > (SWater-GP_DetentZ) *CELL_SIZE )  ? ( -(SWater-GP_DetentZ)*CELL_SIZE ) : (fluxS);*/
02092                     if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02093                         fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02094                      }
02095                  }      
02096                  else
02097                          /*seepage on the levee side */
02098                      fluxL =  (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02099 
02100                                 /* subsurface groundwater flow both sides */
02101                  fluxG =  (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0); 
02102                  break;  
02103              default:
02104                  printf ("Error in levee location: Channel N %d", chan_n+CHo );
02105                  exit (-1);
02106                  break;
02107          } /*end switch */
02108           
02109 
02110 
02111              /* can_cell_M and can_cell_A are the areas of the canal and the cell multiplied (M) and added (A) together, respectively */
02112              /* the flux eqn is the same as the virtual weir eqn, equilibrating the heads across two unequal area regions */
02113 
02114              /* Surface flows: making sure that the cell head is not higher than the canal head */
02115          if (fluxS>0 && ( (fluxCk = tot_head + fluxS/CELL_SIZE - (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) ) > 0.0 ) ) {
02116              fluxS = (can_cell_M*(CH_bottElev + CanWatDep)/can_cell_A - can_cell_M*tot_head/can_cell_A );
02117              fluxCk = tot_head + fluxS/CELL_SIZE - (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area);
02118              if (fluxCk > F_ERROR && debug>3) {
02119                  sprintf(msgStr,"Day %6.1f: Warning - excessive head in cell after flux from Canal %d (%f m diff)", 
02120                          SimTime.TIME, Chan_list[chan_n]->number,fluxCk); 
02121                  WriteMsg( msgStr,True );  
02122              }
02123          }
02124          
02125              /* Surface flows: making sure that the canal water level is not flooded much higher than the water level in the cell */
02126          if ( fluxS<0 && ( (fluxCk = (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) - (tot_head + fluxS/CELL_SIZE) ) > 0.0 ) ) {
02127              fluxS = (can_cell_M*(CH_bottElev + CanWatDep)/can_cell_A - can_cell_M*tot_head/can_cell_A );
02128              fluxCk =  (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) - (tot_head + fluxS/CELL_SIZE);
02129              if (fluxCk > F_ERROR && debug>3) {
02130                  sprintf(msgStr,"Day %6.1f: Warning - excessive head in Canal %d after flux from cell (%f m diff)", 
02131                          SimTime.TIME, Chan_list[chan_n]->number,fluxCk); 
02132                  WriteMsg( msgStr,True );  
02133              }
02134          }
02135 
02136              /* now ensure the canal is not drained below its minimum level */ 
02137          if ( fluxS > 0 || fluxG > 0 || fluxL > 0) /* check gw, seep, and sfwat */
02138          {
02139              fluxO  = (fluxS > 0) ? fluxS : 0;
02140              fluxO += (fluxG > 0) ? fluxG : 0;
02141              fluxO += (fluxL > 0) ? fluxL : 0;
02142     
02143                  /* using vol associated with new estimate (before subtracting losses) */
02144              if ( (fluxCk = CH_vol - MIN_VOL - fluxS - fluxG - fluxL) < 0.0 )
02145              { fluxCk =  (fluxO>0) ? ( ( fluxO + fluxCk )/fluxO) : (0.0);       /* v2.5 removed a divide by zero, but shouldn't be needed */
02146              if (fluxS>0) fluxS *= /* (double) */fluxCk; /* reduce all pos outflows by this proportion */
02147              if (fluxG>0) fluxG *= fluxCk;
02148              if (fluxL>0) fluxL *= fluxCk;
02149              if ( (fluxCk = (CH_vol - MIN_VOL - fluxS - fluxG - fluxL)/Chan_list[chan_n]->area) < -GP_MinCheck ) {
02150                  sprintf(msgStr,"Day %6.1f: capacityERR - excessive draining from Canal %d (%f m below min)", SimTime.TIME,Chan_list[chan_n]->number,fluxCk); 
02151                  WriteMsg( msgStr,True ); dynERRORnum++; }
02152              }
02153          }
02154 
02155 /******** UNUSED *********************/
02156     if (debug>99 && fluxG != 0.0) do {
02157         /* This is a part of the routine used to integrate surface/groundwater
02158          fluxes in the cell-cell Fluxes.c.  It only modifies the allocations if
02159          a cell is a recipient of groundwater flux (does not deal with cells as
02160          donor, which is handled in Fluxes.c).
02161          NOTE:****This routine is NOT fully integrated into the canal iterative solution,
02162          and thus is NOT used for now - get to it some day if it is determined to be needed.  */
02163        fluxHoriz =  fluxG;
02164 
02165 
02166 /**** recipient cell's **pre-flux** capacities */
02167         UnsatZ_rec  =  (fluxG>0) ? ( ElevMap[i] - GWH[i] / poros[HAB[i]] ) : (0.0) ; /* unsat zone depth */
02168         if (debug>2 && (UnsatZ_rec < -0.001 ) ) {
02169             sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in recipient cell (%d,%d) in canal-cell gw fluxes!", 
02170                     SimTime.TIME, UnsatZ_rec, cell->x, cell->y ); WriteMsg( msgStr,True ); dynERRORnum++; }
02171             
02172         UnsatCap_rec =  (fluxG>0) ? (UnsatZ_rec * poros[HAB[i]]) : (0.0); /* maximum pore space capacity (UnsatCap)
02173                                                               in the depth (UnsatZ) of new unsaturated zone */
02174         UnsatPot_rec  = (fluxG>0) ? (UnsatCap_rec - Unsat[i] ) : (0.0); /* (height of) the volume of pore space (soil "removed")
02175                                                            that is unoccupied in the unsat zone */
02176      /* sf-unsat-sat fluxes */
02177         horizTOsat_rec = (fluxG>0) ? (fluxHoriz) : (0.0); /* lateral inflow to soil into saturated storage */
02178         satTOsf_rec = (fluxG>0) ? (Max(fluxHoriz - UnsatPot_rec, 0.0) ) : (0.0); /* upwelling beyond soil capacity */
02179         if (fluxG>0) unsatTOsat_rec = (UnsatZ_rec > 0.0) ? /* incorporation of unsat moisture into sat storage with
02180                                                  rising water table due to horiz inflow */
02181             ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[i]] ) / UnsatZ_rec * Unsat[i] ) :
02182             (0.0);
02183         else
02184             unsatTOsat_rec = 0.0;
02185         
02186 /**** potential new head in recipient cell will be ... */
02187         if (fluxG>0)  H_pot_rec = (GWH[i] + horizTOsat_rec + unsatTOsat_rec - satTOsf_rec)
02188             / poros[HAB[i]] + (SWater + satTOsf_rec);
02189 
02190 
02191         equilib = 1;
02192         
02193     } while  (!equilib); /* end of ***incomplete** routine for integrating groundwater-surface water and preventing head reversals */
02194 /******** UNUSED *********************/
02195 
02196    /*  if (fluxG>0) fluxG=fluxHoriz; */ /* unused, part of unimplemented sf/gwater routine */
02197     
02198          CH_vol -= ( fluxS + fluxG + fluxL); /* subtract flows (pos == loss) from volume associated with new estimate */
02199          
02200          if ( converged )
02201          {      /* the nutrient mass fluxes (pre-flux volumes) (mass=kg, vol=m3)*/
02202                  /* remember pos flux is from canal */
02203                  /* (CH_vol_R is vol prior to in/out fluxes) */
02204         /* surface water canal-cell flux */               
02205              if ( fluxS != 0 ) {
02206              if ( fluxS > 0 ) /* surface water flux from canal to cell */
02207              {  
02208                  /* 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) */
02209                  N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxS) : (0.0);
02210                  P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxS) : (0.0);
02211                  S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxS) : (0.0);
02212 
02213                  if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02214                      if ( (Chan_list[chan_n]->family !=  basn_list[basn[i]]->family) && (debug > 0) ) {        /* if the flow is not within the family... */
02215                          sprintf(msgStr,"Day %6.1f: ERROR - no basin-basin canal-cell overland flow! (Reach %d, basn %d, with cell %d,%d, basin %d)", 
02216                                  SimTime.TIME,Chan_list[chan_n]->number, Chan_list[chan_n]->basin, cell->x, cell->y, basn[i]); 
02217                          WriteMsg( msgStr,True );  dynERRORnum++;
02218                      }
02219                      else {    /* if it's flow within a family, just keep
02220                                   track of what the children do among themselves */            
02221                          if  ( !Chan_list[chan_n]->parent  ) { 
02222                          /* then find out about the flow for the family's sake */
02223                              VOL_OUT_OVL[Chan_list[chan_n]->basin] += fluxS;
02224                              P_OUT_OVL[Chan_list[chan_n]->basin] += P_fl;
02225                              S_OUT_OVL[Chan_list[chan_n]->basin] += S_fl;
02226                          }
02227                          if ( !basn_list[basn[i]]->parent ) {                 
02228                                  /* then find out about the flow for the family's sake */
02229                              VOL_IN_OVL[basn[i]] += fluxS; 
02230                              P_IN_OVL[basn[i]] += P_fl; 
02231                              S_IN_OVL[basn[i]] += S_fl; 
02232                          }
02233                      }
02234                  }                
02235              }
02236              else /* surface water flux from cell to canal */
02237              {
02238                  /* this is where Disp_Calc comes in */
02239                  N_fl = ( SWH[i] > 0 ) ? NA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
02240                  P_fl = ( SWH[i] > 0 ) ? PA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
02241                  S_fl = ( SWH[i] > 0 ) ? SA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
02242                  if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02243 
02244                      if ( (Chan_list[chan_n]->family !=  basn_list[basn[i]]->family) && (debug > 0) ) {        /* if the flow is not within the family... */
02245                          sprintf(msgStr,"Day %6.1f: ERROR - no basin-basin canal-cell overland flow! (Reach %d, basn %d, with cell %d,%d, basin %d)", 
02246                                         SimTime.TIME,Chan_list[chan_n]->number, Chan_list[chan_n]->basin, cell->x, cell->y, basn[i]); 
02247                          WriteMsg( msgStr,True );  dynERRORnum++;
02248                  }
02249                  else {    /* if it's flow within a family, just keep
02250                           track of what the children do among themselves */            
02251                      if  ( !basn_list[basn[i]]->parent  ) { 
02252                          /* then find out about the flow for the family's sake */
02253                          VOL_OUT_OVL[basn[i]] -= fluxS;
02254                          P_OUT_OVL[basn[i]] -= P_fl;
02255                          S_OUT_OVL[basn[i]] -= S_fl;
02256                      }
02257                      if ( !Chan_list[chan_n]->parent ) {                 
02258                          /* then find out about the flow for the family's sake */
02259                          VOL_IN_OVL[Chan_list[chan_n]->basin] -= fluxS; 
02260                          P_IN_OVL[Chan_list[chan_n]->basin] -= P_fl; 
02261                          S_IN_OVL[Chan_list[chan_n]->basin] -= S_fl; 
02262                      }
02263                  }
02264                  }
02265 
02266 
02267              } 
02268                  /* now pass that mass flux to/from the cell/canal */
02269              NA[i] += N_fl;
02270              PA[i] += P_fl;
02271              SA[i] += S_fl;
02272              Chan_list[chan_n]->N -= N_fl;
02273              Chan_list[chan_n]->P -= P_fl;
02274              Chan_list[chan_n]->S -= S_fl;
02275              }
02276 
02277         /* groundwater canal-cell flux */                 
02278              if ( fluxG != 0 ) {
02279              if ( fluxG > 0 ) /* ground water flux from canal to cell */
02280              {  
02281                  N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxG) : (0.0);
02282                  P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxG) : (0.0);
02283                  S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxG) : (0.0);
02284 
02285                  if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02286                      if (Chan_list[chan_n]->family !=  
02287                          basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
02288                          if  ( !Chan_list[chan_n]->parent  ) { /* and if the donor or recipient is a child... */
02289                          /* then find out about the flow for the family's sake */
02290                              VOL_OUT_GW[Chan_list[chan_n]->family] += fluxG;
02291                              P_OUT_GW[Chan_list[chan_n]->family] += P_fl;
02292                              S_OUT_GW[Chan_list[chan_n]->family] += S_fl;
02293                          }
02294                          if ( !basn_list[basn[i]]->parent ) {                 
02295                            /* then find out about the flow for the family's sake */
02296                              VOL_IN_GW[basn_list[basn[i]]->family] += fluxG; 
02297                              P_IN_GW[basn_list[basn[i]]->family] += P_fl; 
02298                              S_IN_GW[basn_list[basn[i]]->family] += S_fl; 
02299                          }
02300                      VOL_OUT_GW[Chan_list[chan_n]->basin] += fluxG;
02301                      VOL_IN_GW[basn[i]] += fluxG; 
02302                      P_OUT_GW[Chan_list[chan_n]->basin] += P_fl;
02303                      P_IN_GW[basn[i]] += P_fl; 
02304                      S_OUT_GW[Chan_list[chan_n]->basin] += S_fl;
02305                      S_IN_GW[basn[i]] += S_fl; 
02306                  }
02307                  else {    /* if it's flow within a family, just keep
02308                           track of what the children do among themselves */            
02309                      if  ( !Chan_list[chan_n]->parent  ) { 
02310                          /* then find out about the flow for the family's sake */
02311                          VOL_OUT_GW[Chan_list[chan_n]->basin] += fluxG;
02312                          P_OUT_GW[Chan_list[chan_n]->basin] += P_fl;
02313                          S_OUT_GW[Chan_list[chan_n]->basin] += S_fl;
02314                      }
02315                      if ( !basn_list[basn[i]]->parent ) {                 
02316                          /* then find out about the flow for the family's sake */
02317                          VOL_IN_GW[basn[i]] += fluxG; 
02318                          P_IN_GW[basn[i]] += P_fl; 
02319                          S_IN_GW[basn[i]] += S_fl; 
02320                      }
02321                  }
02322                  }
02323              } /* end of positive flow evaluations */
02324              
02325              
02326              else  /*  ground water flux from cell to canal */
02327              {
02328                  N_fl = ( GWH[i] > 0 ) ? GNA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
02329                  P_fl = ( GWH[i] > 0 ) ? GPA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
02330                  S_fl = ( GWH[i] > 0 ) ? GSA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
02331                  if (Chan_list[chan_n]->basin != basn[i] && debug > -999) {
02332                      if (Chan_list[chan_n]->family !=  
02333                          basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
02334                          if  ( !basn_list[basn[i]]->parent  ) { /* and if the donor or recipient is a child... */
02335                          /* then find out about the flow for the family's sake */
02336                              VOL_OUT_GW[basn_list[basn[i]]->family] -= fluxG;
02337                              P_OUT_GW[basn_list[basn[i]]->family] -= P_fl;
02338                              S_OUT_GW[basn_list[basn[i]]->family] -= S_fl;
02339                          }
02340                          if ( !Chan_list[chan_n]->parent ) {                 
02341                            /* then find out about the flow for the family's sake */
02342                              VOL_IN_GW[Chan_list[chan_n]->family] -= fluxG; 
02343                              P_IN_GW[Chan_list[chan_n]->family] -= P_fl; 
02344                              S_IN_GW[Chan_list[chan_n]->family] -= S_fl; 
02345                          }
02346                      VOL_IN_GW[Chan_list[chan_n]->basin] -= fluxG;
02347                      VOL_OUT_GW[basn[i]] -= fluxG;
02348                      P_IN_GW[Chan_list[chan_n]->basin] -= P_fl;
02349                      P_OUT_GW[basn[i]] -= P_fl;
02350                      S_IN_GW[Chan_list[chan_n]->basin] -= S_fl;
02351                      S_OUT_GW[basn[i]] -= S_fl;
02352                  }
02353                  else {    /* if it's flow within a family, just keep
02354                           track of what the children do among themselves */            
02355                      if  ( !basn_list[basn[i]]->parent  ) { 
02356                          /* then find out about the flow for the family's sake */
02357                          VOL_OUT_GW[basn[i]] -= fluxG;
02358                          P_OUT_GW[basn[i]] -= P_fl;
02359                          S_OUT_GW[basn[i]] -= S_fl;
02360                      }
02361                      if ( !Chan_list[chan_n]->parent ) {                 
02362                          /* then find out about the flow for the family's sake */
02363                          VOL_IN_GW[Chan_list[chan_n]->basin] -= fluxG; 
02364                          P_IN_GW[Chan_list[chan_n]->basin] -= P_fl; 
02365                          S_IN_GW[Chan_list[chan_n]->basin] -= S_fl; 
02366                      }
02367                  }
02368                  }
02369              } 
02370                  /* now pass that mass flux to/from the cell/canal
02371                   this includes allocating nuts to surface water thru upflow (unimplemented) */
02372              /* unused, part of unimplemented sf/gwater routine */
02373 /*              SNfl = Chan_list[chan_n]->N/CH_vol_R*satTOsf_rec ; */
02374 /*              GNfl = N_fl - SNfl; */
02375 /*              SPfl = Chan_list[chan_n]->P/CH_vol_R*satTOsf_rec ; */
02376 /*              GPfl = P_fl - SPfl; */
02377 /*              SSfl = Chan_list[chan_n]->S/CH_vol_R*satTOsf_rec ; */
02378 /*              GSfl = S_fl - SSfl; */
02379 /*              GNA[i] += GNfl; */
02380 /*              GPA[i] += GPfl; */
02381 /*              GSA[i] += GSfl; */
02382              GNA[i] += N_fl;
02383              GPA[i] += P_fl;
02384              GSA[i] += S_fl;
02385              Chan_list[chan_n]->N -= N_fl;
02386              Chan_list[chan_n]->P -= P_fl;
02387              Chan_list[chan_n]->S -= S_fl;
02388                  /* upflow into surface water nuts (this is actually canal conc., but what the heck) */
02389              /* unused, part of unimplemented sf/gwater routine */
02390 /*              NA[i] += SNfl; */
02391 /*              PA[i] += SPfl; */
02392 /*              SA[i] += SSfl; */
02393              }
02394 
02395         /* seepage canal-cell flux */             
02396              if ( fluxL != 0 ) {
02397              if ( fluxL > 0 ) /* seepage flux from canal to cell */
02398              {  
02399                  N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxL) : (0.0);
02400                  P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxL) : (0.0);
02401                  S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxL) : (0.0);
02402                  if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02403                      if (Chan_list[chan_n]->family !=  
02404                          basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
02405                          if  ( !Chan_list[chan_n]->parent  ) { /* and if the donor or recipient is a child... */
02406                          /* then find out about the flow for the family's sake */
02407                              VOL_OUT_SPG[Chan_list[chan_n]->family] += fluxL;
02408                              P_OUT_SPG[Chan_list[chan_n]->family] += P_fl;
02409                              S_OUT_SPG[Chan_list[chan_n]->family] += S_fl;
02410                          }
02411                          if ( !basn_list[basn[i]]->parent ) {                 
02412                            /* then find out about the flow for the family's sake */
02413                              VOL_IN_SPG[basn_list[basn[i]]->family] += fluxL; 
02414                              P_IN_SPG[basn_list[basn[i]]->family] += P_fl; 
02415                              S_IN_SPG[basn_list[basn[i]]->family] += S_fl; 
02416                          }
02417                      VOL_OUT_SPG[Chan_list[chan_n]->basin] += fluxL;
02418                      VOL_IN_SPG[basn[i]] += fluxL; 
02419                      P_OUT_SPG[Chan_list[chan_n]->basin] += P_fl;
02420                      P_IN_SPG[basn[i]] += P_fl; 
02421                      S_OUT_SPG[Chan_list[chan_n]->basin] += S_fl;
02422                      S_IN_SPG[basn[i]] += S_fl; 
02423                  }
02424                  else {    /* if it's flow within a family, just keep
02425                           track of what the children do among themselves */            
02426                      if  ( !Chan_list[chan_n]->parent  ) { 
02427                          /* then find out about the flow for the family's sake */
02428                          VOL_OUT_SPG[Chan_list[chan_n]->basin] += fluxL;
02429                          P_OUT_SPG[Chan_list[chan_n]->basin] += P_fl;
02430                          S_OUT_SPG[Chan_list[chan_n]->basin] += S_fl;
02431                      }
02432                      if ( !basn_list[basn[i]]->parent ) {                 
02433                          /* then find out about the flow for the family's sake */
02434                          VOL_IN_SPG[basn[i]] += fluxL; 
02435                          P_IN_SPG[basn[i]] += P_fl; 
02436                          S_IN_SPG[basn[i]] += S_fl; 
02437                      }
02438                  }
02439                  }
02440                      if ( tot_head > GW_head ) { /* to cell surface water */
02441                                 NA[i] += N_fl;
02442                                 PA[i] += P_fl;
02443                                 SA[i] += S_fl;
02444                                 }
02445                          else { /* to cell groundwater */
02446                                 GNA[i] += N_fl;
02447                                 GPA[i] += P_fl;
02448                                 GSA[i] += S_fl;
02449                         }
02450              }
02451              else if (fluxL < 0) /*  seepage flux from cell to canal */
02452              {
02453                  if ( tot_head > GW_head ) { /* cell's tot head > its groundwater head */
02454                         N_fl = ( SWH[i] > 0 ) ? NA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
02455                         P_fl = ( SWH[i] > 0 ) ? PA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
02456                         S_fl = ( SWH[i] > 0 ) ? SA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
02457                         NA[i] += N_fl;
02458                         PA[i] += P_fl;
02459                         SA[i] += S_fl;
02460                  }
02461                  else {
02462                         N_fl = ( GWH[i] > 0 ) ? GNA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
02463                         P_fl = ( GWH[i] > 0 ) ? GPA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
02464                         S_fl = ( GWH[i] > 0 ) ? GSA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
02465                         GNA[i] += N_fl;
02466                         GPA[i] += P_fl;
02467                         GSA[i] += S_fl;
02468                 }
02469                  if (Chan_list[chan_n]->basin != basn[i] && debug > -999) {
02470                      if (Chan_list[chan_n]->family !=  
02471                          basn_list[basn[i]]->family ) {        /* if the flow is not within the family... */
02472                          if  ( !basn_list[basn[i]]->parent  ) { /* and if the donor or recipient is a child... */
02473                          /* then find out about the flow for the family's sake */
02474                              VOL_OUT_SPG[basn_list[basn[i]]->family] -= fluxL;
02475                              P_OUT_SPG[basn_list[basn[i]]->family] -= P_fl;
02476                              S_OUT_SPG[basn_list[basn[i]]->family] -= S_fl;
02477                          }
02478                          if ( !Chan_list[chan_n]->parent ) {                 
02479                            /* then find out about the flow for the family's sake */
02480                              VOL_IN_SPG[Chan_list[chan_n]->family] -= fluxL; 
02481                              P_IN_SPG[Chan_list[chan_n]->family] -= P_fl; 
02482                              S_IN_SPG[Chan_list[chan_n]->family] -= S_fl; 
02483                          }
02484                      VOL_IN_SPG[Chan_list[chan_n]->basin] -= fluxL;
02485                      VOL_OUT_SPG[basn[i]] -= fluxL;
02486                      P_IN_SPG[Chan_list[chan_n]->basin] -= P_fl;
02487                      P_OUT_SPG[basn[i]] -= P_fl;
02488                      S_IN_SPG[Chan_list[chan_n]->basin] -= S_fl;
02489                      S_OUT_SPG[basn[i]] -= S_fl;
02490                  }
02491                  else {    /* if it's flow within a family, just keep
02492                           track of what the children do among themselves */            
02493                      if  ( !basn_list[basn[i]]->parent  ) { 
02494                          /* then find out about the flow for the family's sake */
02495                          VOL_OUT_SPG[basn[i]] -= fluxL;
02496                          P_OUT_SPG[basn[i]] -= P_fl;
02497                          S_OUT_SPG[basn[i]] -= S_fl;
02498                      }
02499                      if ( !Chan_list[chan_n]->parent ) {                 
02500                          /* then find out about the flow for the family's sake */
02501                          VOL_IN_SPG[Chan_list[chan_n]->basin] -= fluxL; 
02502                          P_IN_SPG[Chan_list[chan_n]->basin] -= P_fl; 
02503                          S_IN_SPG[Chan_list[chan_n]->basin] -= S_fl; 
02504                      }
02505                  }
02506                  }
02507              } 
02508                  /* now pass that mass flux to/from the canal */
02509              Chan_list[chan_n]->N -= N_fl;
02510              Chan_list[chan_n]->P -= P_fl;
02511              Chan_list[chan_n]->S -= S_fl;
02512              }
02513                                 /* now we recalculate the depths of SFWater and Gwater in cell */
02514 /*              SWH[i] += (fluxS+satTOsf_rec)/CELL_SIZE;  */
02515 /*              GWH[i] += (fluxG-satTOsf_rec)/CELL_SIZE;  */
02516              SWH[i] += (fluxS)/CELL_SIZE; 
02517              GWH[i] += (fluxG)/CELL_SIZE; 
02518              if ( tot_head > GW_head ) SWH[i] += fluxL/CELL_SIZE; /* whether to add seepage to surface or */
02519              else GWH[i] += fluxL/CELL_SIZE;                      /* groundwater is not critical, as vertical updates will follow */
02520                                 
02521 
02522          } /* that's all we do when finally converged */
02523                     
02524          T_flux_S += fluxS;
02525          T_flux_G += fluxG;
02526          T_flux_L += fluxL;
02527                   
02528          cell = cell->next_cell;
02529      } /*end cycle along canal cells */
02530 
02531          
02532      error = (CanWatDep - Chan_list[chan_n]->wat_depth) - (Qin - Qout - T_flux_S - T_flux_G - T_flux_L)/ Chan_list[chan_n]->area;
02533      if (converged && error != errorConv ) {
02534          /* error the second time around should be close to that of the first */
02535          if (fabs(error-errorConv)>F_ERROR ) { 
02536              sprintf(msgStr,"Day %6.1f: ERROR - converg error != to error after cell-canal fluxes (Chan %d, est= %5.3f, last depth= %5.3f m)", 
02537                         SimTime.TIME, Chan_list[chan_n]->number, CanWatDep, Chan_list[chan_n]->wat_depth); 
02538              WriteMsg( msgStr,True );  dynERRORnum++;
02539          }
02540          
02541          if (CanWatDep>Chan_list[chan_n]->depth*3.0 || CanWatDep < 0.0) {
02542                 CanWatDep = MINDEPTH; /* lost the mass balance, but let it go for now */
02543                 sprintf(msgStr,"           Reset Chan %d depth to minimum allowed (%5.3f m)", 
02544                                 Chan_list[chan_n]->number, MINDEPTH); 
02545                 WriteMsg( msgStr,True );  }
02546                  
02547                  /*errorConv is the error upon first reaching convergence criteria; */
02548              /* we then recalc the fluxes in the same loop thru cells (w/ same depth est.), */
02549              /* and should have the same fluxes and error as before */
02550          if (debug>4) { /* this here to halt program with high debug */
02551              printf("\nDay %6.1f: ERROR - converg error != to error after cell-canal fluxes (Chan %d, est= %5.3f m)\n", 
02552                         SimTime.TIME, Chan_list[chan_n]->number, CanWatDep);
02553              exit (-1);
02554 
02555          }
02556          
02557      }
02558      
02559      ++iter;
02560 
02561  } while (!converged );  /* end relaxation procedure */
02562 
02563  if ( factor != 0 && Qout > 0) /* if relaxation was abnormally stopped because of hitting the bottom */
02564  { 
02565      sprintf(msgStr,"Day %6.1f: ERROR - hit bottom of Canal %d, stopped the iterative relaxation", SimTime.TIME, Chan_list[chan_n]->number); 
02566      WriteMsg( msgStr,True );     
02567      dynERRORnum++;         
02568             /*#######*******########  this old code should not be reached
02569              the code is questionable - an ERROR message printed
02570              to file will indicate signif problem due to entrance to this */
02571      Qout1 = (Chan_list[chan_n]->wat_depth - CanWatDep)*Chan_list[chan_n]->area + Qin - T_flux_S - T_flux_G - T_flux_L;
02572      structs = struct_first;
02573      while ( structs != NULL )   /* recalculate the outflow in the data for possible use in the next canal inflow */
02574      {  
02575          if ( structs->canal_fr  == chan_n+CHo )   structs->flow *= Qout1/Qout;
02576     
02577          structs = structs->next_in_list;
02578      }
02579      Qout = Qout1;
02580  } 
02581 
02582  Chan_list[chan_n]->wat_depth = CanWatDep;      /* new canal water depth */
02583  CH_vol =  CanWatDep * Chan_list[chan_n]->area; /* new canal volume */
02584 
02585  if( debug > 3 ) 
02586  {
02587          /* concentration in mg/L (kg/m3==g/L, convert P (not S) to mg/L) */
02588      Chan_list[chan_n]->P_con = (CH_vol > 0) ? (Chan_list[chan_n]->P/CH_vol*1000.0) : 0.0;
02589      Chan_list[chan_n]->S_con = (CH_vol > 0) ? (Chan_list[chan_n]->S/CH_vol) : 0.0;
02590      fprintf( canDebugFile, 
02591               "%4d  %8.3f  %8.4f  %8.4f  %12.1f  %12.1f  %12.1f  %12.1f  %12.1f  %5d  %9.1f  %7.5f \n",
02592               Chan_list[chan_n]->number, Chan_list[chan_n]->wat_depth, Chan_list[chan_n]->P_con,  
02593               Chan_list[chan_n]->S_con, T_flux_S, T_flux_G, T_flux_L, Qin, Qout, iter, Chan_list[chan_n]->area, error );
02594  
02595  }
02596 
02597  return;
02598 }
02599 
02600 /***************************************************************************************/
02609 float f_Manning( float delta, float Water, float SW_coef )
02610 {       float a_delta;
02611     
02612         a_delta = Abs (delta);
02613         return ( sgn( delta ) * SW_coef * pow(Water, GP_mannDepthPow) * sqrt(a_delta) * canstep);
02614 
02615 }
02616 /****************************************************************************************/
02626 float f_Ground ( float dh, float height, float GW_coef, float I_Length )
02627 {       
02628         float  w;
02629 
02630 /* porosity already part of GW_coef */
02631         w = dh * I_Length * GW_coef / (0.5*celWid) * height * canstep; 
02632 
02633         return ( w);
02634 
02635 }
02636 
02637 /****************************************************************************/
02675 void Flows_in_Structures (  float *SWH, float *Elev, float *MC, double *NA, double *PA, double *SA )
02676 {
02691         struct Structure *structs, *struct_last, *struct_hist_start;
02692     int cell_to, cell_fr, can_to, can_fr;
02693     int cell_struct;
02694     int TWf, HWf;
02695     double HeadH, HeadT, HeadH2, HeadT2;
02696     double CH_vol, cell_vol;
02697     double Nconc, Pconc, Sconc;
02698     double N_fl, P_fl, S_fl; 
02699     double chan_over;
02700     double flow_rev; 
02701     double ChanHistOut;
02702     double HeadH_drop, HeadT_drop;
02703     float FFlux; 
02704 
02705     /* from unitmod.h */
02706     extern float* boundcond_depth; /* depth in external cells */
02707 
02708 
02709 
02710     TWf = 0; HWf = 0; /* 0 when there is no schedule graphs for tail or headwater */
02711     
02712     structs = struct_first; 
02713 
02714     while ( structs != NULL )
02715     {           
02716         if ( structs->flag < 0 || (structs->flag >1) ) {
02717             structs->flow = 0.0; /* used in summing up flows for Qin and Qout later */
02718             structs->conc_P = 0.0;
02719             structs->conc_N = 0.0;
02720             structs->conc_S = 0.0;
02721             goto OUT; 
02722         } /* nothing to do when structure inactive (neg or >1 flag) */
02723         
02725         if ( structs->TW_stage == -99 ) /* this means that there is a graph for TW schedules */
02726         { 
02727         structs->TW_stage = GetGraph ( structs->TW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 0.0) ); 
02728         /* add sea level rise */
02729         structs->TW_stage = structs->TW_stage + GP_SLRise * (SimTime.TIME/365);
02730         TWf = 1; 
02731         }
02732           
02733         if ( structs->HW_stage == -99 ) /* this means that there is a graph for HW schedules */
02734         { 
02735         structs->HW_stage = GetGraph ( structs->HW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 365.0) ); 
02736         /* add sea level rise */
02737         structs->HW_stage = structs->HW_stage + GP_SLRise * (SimTime.TIME/365);
02738         HWf = 1; }
02739         
02740           
02741 /* CASE 1 */
02742 /*check that this is a canal-to-canal structure */
02743 /*  these flows are always internal to model domain */
02744         if ( structs->canal_fr  != 0 && structs->canal_to != 0 ) {
02745             can_fr = structs->canal_fr - CHo;
02746             can_to = structs->canal_to - CHo;
02747             cell_struct = T(structs->str_cell_i,structs->str_cell_j); 
02748                 
02749                 /* first check for historical/sfwmm flow */
02750             if ( structs->flag == 1 ) {
02751                 if (structs->multiOut == 1) goto OUT; /* this structure has already been processed as a multiple outflow from a canal */
02752                     /* FLOW */
02753                     /* don't do much with zero or negative flows */           
02754                 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
02755                 if (structs->flow == 0.0) { 
02756                     structs->conc_P = 0.0;
02757                     structs->conc_N = 0.0;
02758                     structs->conc_S = 0.0;
02759                     goto OUT; 
02760                 }   /* zero the reported struct flow conc (for output info only) */
02761                 else if (structs->flow < 0.0) {
02762                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
02763                              SimTime.TIME, Chan_list[can_to]->number );
02764                     WriteMsg( msgStr,True );
02765                     dynERRORnum++;
02766                     goto OUT;
02767                 }
02768 
02769                     /* for multiple historical/other-data outflows from a canal, ensure mass balance (sum of outflows !> avail volume) */
02770                     /* all flows from the canal associated with each set must be non-negative, which is currently valid */
02771                     /* the below cycle thru the structures will also process a canal->cell flow structure */
02772 
02773                 struct_hist_start = structs; /* this is the first structure w/ historical info for a particular canal */
02774                 ChanHistOut = 0.0; /* sum of all historical/other-data outflows from a canal */
02775 
02776                 while ( structs != NULL )   
02777                 {  
02778                     if  ( structs->flag == 1 && struct_hist_start->canal_fr  == structs->canal_fr   ) { /* look for outflows from same canal */
02779                         structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
02780                         ChanHistOut += structs->flow; /* sum the multiple historical/other-data outflows from canal */
02781                     }
02782                         
02783                     structs = structs->next_in_list; /* increment to next structure */
02784                 } /* end cycle thru list of remaining structures in total list */
02785                     
02786                 structs = struct_hist_start; /* now we are pointing back to the first (or only) water control structure draining the canal */
02787                          
02788                     /* concentration in donor canal for structure flow, g/L = kg/m^3 */
02789                     /* this is based on previous depth/volume calc'd in FluxChannel */
02790                     /* or use historical conc */
02791 
02792                     /* CONC */
02793                                 /* canal TOTAL volume before new structure flows */
02794                 CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth; 
02795                     /* concentrations applied to all structures draining this canal */
02796                 Pconc = (structs->TPser == -1) ? 
02797                     ( (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0) ) :
02798                     ( (structs->TPser == 0) ? (structs->TP) : ( Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
02799                 Nconc = (structs->TNser == -1) ? 
02800                     ( (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0) ) :
02801                     ( (structs->TNser == 0) ? (structs->TN) : ( Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
02802                 Sconc = (structs->TSser == -1) ? 
02803                     ( (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0) ) :
02804                     ( (structs->TSser == 0) ? (structs->TS) : ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
02805 
02806 
02807                     /* vol avail for fluxing */
02808                 CH_vol =  Chan_list[can_fr]->area*(ramp(Chan_list[can_fr]->wat_depth-MINDEPTH) );
02809                 chan_over = (ChanHistOut > 0.0 ) ? ( CH_vol/ChanHistOut ) : 1.0;
02810 
02811                 do /* loop again through the remaining sequence of structures to correct any outflows for the canal */
02812                 {   
02813                     if (structs->flag == 1 &&  struct_hist_start->canal_fr  == structs->canal_fr   ) { 
02814 
02815                         structs->multiOut = 1; /* set the flag to indicate this structure has been processed */
02816                             
02817                             /* revised FLOW */
02818                         if ( chan_over < 1.0 ) { /* if hist/other-data demand > canal vol */
02819                             if( debug > 0 ) { 
02820                                 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",
02821                                          SimTime.TIME, structs->S_nam, ChanHistOut, Chan_list[can_fr]->number,
02822                                          CH_vol, (ChanHistOut - CH_vol), Chan_list[can_fr]->wat_depth );
02823                                 
02824                             }  
02825                                 /* decrement the individual structure outflows by it's proportional contrib. */
02826                             structs->flow = (structs->flow) * chan_over ;
02827                         } 
02828                             /* sum of all realized (corrected for excess demand) historical/other-data outflows from a canal during an iteration */
02829                             /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
02830                         Chan_list[can_fr]->sumHistOut += structs->flow;         
02831 
02832                             /* the recipient changes as we cycle thru the structures (may be another canal or a cell) */
02833                         can_to = structs->canal_to - CHo;
02834                         cell_to = T(structs->cell_i_to,structs->cell_j_to); 
02835 
02836                         structs->conc_P = Pconc;/* concentrations applied to all structures draining this canal */
02837                         structs->conc_N = Nconc;
02838                         structs->conc_S = (structs->TSser == -1) ?
02839                             (Sconc) :
02840                             ( (structs->TSser == 0) ? (structs->TS) :
02841                               ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ); /* see note below on tracer conc. application */
02842                         
02843                         P_fl = structs->flow * structs->conc_P; 
02844                         N_fl = structs->flow * structs->conc_N;
02845                         S_fl = structs->flow * structs->conc_S;
02846                                 /* change in nut mass in fr canal (kg) */
02847                         Chan_list[can_fr]->N -=  N_fl ;
02848                         Chan_list[can_fr]->P -=  P_fl ;
02849                         Chan_list[can_fr]->S -=  (structs->TSser == -1) ? (S_fl) : (0);
02850                             /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
02851                                Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
02852                                that was not calculated from the available water and salt (i.e., the conc. was from the
02853                                CanalData.struct data file, either fixed or a time series) */
02854                             
02855                         if (can_to > 0) {    /* IF this is a canal recip, change in nut mass in to canal (kg) */
02856                                 /* sum of all realized (corrected for excess demand) historical/other-data outflows into a canal during an iteration */
02857                                 /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
02858                             Chan_list[can_to]->sumHistIn += structs->flow;              
02859 
02860                             Chan_list[can_to]->N +=  N_fl ;
02861                             Chan_list[can_to]->P +=  P_fl ;
02862                             Chan_list[can_to]->S +=  S_fl ;
02863                                 /* mass balance in fr and to canals */
02864                                 /* canals themselves always belong to a parent hydro basin,
02865                                    thus canal-canal struct flows do not keep track of child-parent, child-child flows */
02866                             if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > -999 ) {
02867                                 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
02868                                 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; 
02869                                 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
02870                                 P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
02871                                 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
02872                                 S_IN_STR[Chan_list [can_to]->basin] += S_fl; 
02873                             }
02874 
02875 
02876                                              
02877                         }
02878                         else if (cell_to > 0) {
02879                                 /* OR change in nut mass in downstream internal cell (kg) */
02880  
02881                             if (structs->cell_i_to > 2 ) { /* to an on-map cell */
02882                                 SWH[cell_to] += structs->flow/CELL_SIZE;
02883                                 PA[cell_to] += P_fl; 
02884                                 NA[cell_to] += N_fl; 
02885                                 SA[cell_to] += S_fl; 
02886                                     /* mass balance in fr canal and to cell */
02887                                 if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
02888                                     VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
02889                                     VOL_IN_STR[basn[cell_to]] += structs->flow; 
02890                                     P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
02891                                     P_IN_STR[basn[cell_to]] += P_fl; 
02892                                     S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
02893                                     S_IN_STR[basn[cell_to]] += S_fl; 
02894 
02895                                         /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
02896                                     if ( (Chan_list[can_fr]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
02897                                         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",
02898                                                  SimTime.TIME, Chan_list[can_fr]->number,structs->cell_i_to,structs->cell_j_to );
02899                                         WriteMsg( msgStr,True );
02900                                         dynERRORnum++;
02901                                     }
02902                                 }
02903                                 
02904 
02905 
02906                             }
02907                             else if (structs->cell_i_to == 2 && debug > -999 ) { /* to off-map cell, mass balance check */
02908                                 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow; 
02909                                 VOL_OUT_STR[0] += structs->flow; 
02910                                 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl; 
02911                                 P_OUT_STR[0] += P_fl; 
02912                                 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
02913                                 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0); /* see note above on tracer conc. application */
02914                             }
02915                         }
02916                             /* don't do much with zero or negative flows */           
02917                         if (structs->flow == 0.0) {
02918                             structs->conc_P = 0.0;
02919                             structs->conc_N = 0.0;
02920                             structs->conc_S = 0.0;
02921                         } /* for output purposes only, zero the reported conc */
02922                         else if (structs->flow < 0.0) {
02923                             sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from another cell/chan into chan %d",
02924                                      SimTime.TIME,Chan_list[can_fr]->number );
02925                             WriteMsg( msgStr,True );
02926                             dynERRORnum++;
02927                         }
02928                     }
02929 
02930                     structs = structs->next_in_list;
02931                 } while (structs != NULL); /*  cycling to the end of the list again !!!! */
02932 
02933                 structs = struct_hist_start; /* now pointing back to the first historical outflow structure that we dealt with */
02934 
02935                 goto OUT;
02936             } /* finished check for multiple (positive) outflows exceeding canal volume */
02937 
02938             else  { /* Rule- based flow calcs */
02939 
02940                 HeadH_drop = 0.5 * Abs(Chan_list[can_fr]->elev_drop); /* elevation difference between begin&end of canal */
02941                 HeadT_drop = 0.5 * Abs(Chan_list[can_to]->elev_drop);
02942                 
02943                 HeadH = HeadH_drop + /* add portion of the elevation drop along the canal length */
02944                     Elev[cell_struct] - Chan_list[can_fr]->depth + Chan_list[can_fr]->wat_depth 
02945                     + (Chan_list[can_fr]->sumHistIn - Chan_list[can_fr]->sumHistOut)/Chan_list[can_fr]->area; 
02946                                 /* in both head and tail, add net historical/other-data flows before determining hydraulic potential */
02947                     /* in tailwater only, check to see if other virtual struct has added water already, add to head */
02948                 HeadT = -HeadT_drop + /* subtract portion of the elevation drop along the canal length */
02949                     Elev [cell_struct] - Chan_list[can_to]->depth + Chan_list[can_to]->wat_depth
02950                     + ( Chan_list[can_to]->sumRuleIn + Chan_list[can_to]->sumHistIn - Chan_list[can_to]->sumHistOut)/Chan_list[can_to]->area;
02951 
02952                     /* FLOW */
02953                 if ( structs->w_coef == 0.0)  {/* check for a virtual structure, flow only in pos direction */
02954                     if ( HeadH>HeadT ) {         
02955                         structs->flow =
02956                             Chan_list[can_fr]->area*Chan_list[can_to]->area/
02957                             (Chan_list[can_fr]->area+Chan_list[can_to]->area) * (HeadH-HeadT) ;
02958                         if (debug >3) {
02959                             HeadH2 = HeadH - structs->flow/Chan_list[can_fr]->area; 
02960                             HeadT2 = HeadT + structs->flow/Chan_list[can_to]->area;
02961                             sprintf( msgStr, "Day %6.1f:  virtual flow from chan %d (%f m) to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
02962                                      SimTime.TIME,Chan_list[can_fr]->number,HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth,
02963                                      Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
02964                             WriteMsg( msgStr,True ); 
02965                             if (Abs(HeadH2 - HeadT2) > GP_MinCheck) /* virtual flows should exactly equilibrate the two canals */
02966                             { sprintf( msgStr, "Day %6.1f: ERROR virtual flow from chan %d to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
02967                                        SimTime.TIME,Chan_list[can_fr]->number,Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
02968                             WriteMsg( msgStr,True ); dynERRORnum++; }
02969 
02970                         }
02971                         
02972                         CH_vol =  Chan_list[can_fr]->area* ramp((HeadH-HeadH_drop)-Elev[cell_struct]+Chan_list[can_fr]->depth-MINDEPTH); /* avail flow volume */
02973                         if (structs->flow > CH_vol) {
02974                             structs->flow =  CH_vol;
02975                             HeadH2 = HeadH - structs->flow/Chan_list[can_fr]->area; 
02976                             HeadT2 = HeadT + structs->flow/Chan_list[can_to]->area;
02977                             if (debug>3 && Abs(HeadH2 - HeadT2) > GP_MinCheck) {
02978                                 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",
02979                                          SimTime.TIME,Chan_list[can_fr]->number,CH_vol/Chan_list[can_fr]->area,Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
02980                                 WriteMsg( msgStr,True ); }
02981                         }
02982                         
02983                             /* revise volume to be TOTAL volume for conc calc */
02984                         CH_vol = CH_vol + MINDEPTH*Chan_list[can_fr]->area;
02985                             /* conc in donor (fr) canal (kg/m3) */
02986                         structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
02987                         structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
02988                         structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
02989                     }
02990                     else  structs->flow = 0.0;
02991                     
02992                 }
02993                 
02994                 else
02995                 {
02996                     printf ("\nVirtual structs are the only rule-based canal-canal flows allowed in this version!\n"); exit(-1);
02997                         /* the code that was here is in older version,
02998                            not fully implemented (ELM driven by historical/other-data flows, only canal-canal virtual structs allowed) */
02999                     
03000                 } /* end of rule based flow calcs */
03001  
03002                     /* don't do much with zero or negative rule-based flows  */           
03003                 if (structs->flow == 0.0) {
03004                     structs->conc_P = 0.0;
03005                     structs->conc_N = 0.0;
03006                     structs->conc_S = 0.0;
03007                     goto OUT;
03008                 } /* for output purposes only, zero the reported conc */
03009                 else if (structs->flow < 0.0 && structs->w_coef != 0) {                 
03010                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG RULE-BASED FLOW from chan %d",
03011                              SimTime.TIME,Chan_list[can_to]->number);
03012                     WriteMsg( msgStr,True );
03013                     goto OUT;
03014                 } 
03015                                                                 
03016 
03017                     /* increment the sum of rule-based inflows/outflows to canals
03018                        for potential multiple virtual inflows into some canals */
03019                 Chan_list[can_fr]->sumRuleOut += structs->flow;         
03020                 Chan_list[can_to]->sumRuleIn += structs->flow;          
03021 
03022                     /* nut flows */
03023                 P_fl = structs->flow * structs->conc_P; 
03024                 N_fl = structs->flow * structs->conc_N;
03025                 S_fl = structs->flow * structs->conc_S;
03026                     /* change in nut mass in fr canal (kg) */
03027                 Chan_list[can_fr]->N -=  N_fl ;
03028                 Chan_list[can_fr]->P -=  P_fl ;
03029                 Chan_list[can_fr]->S -=  S_fl ;
03030                     /* change in nut mass in to canal (kg) */
03031                 Chan_list[can_to]->N +=  N_fl ;
03032                 Chan_list[can_to]->P +=  P_fl ;
03033                 Chan_list[can_to]->S +=  S_fl ;
03034                     /* mass balance in fr and to canals */
03035                     /* virtual structure flows can not (by definition) establish flow among basins */
03036                 if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > 0 ) {
03037                     sprintf( msgStr, "Day %6.1f: ERROR - no hydrologic basin-basin flows in virtual structure, canals %d - %d",
03038                              SimTime.TIME, Chan_list[can_to]->number, Chan_list[can_fr]->number);
03039                     WriteMsg( msgStr,True ); dynERRORnum++;
03040                 }
03041          
03042                 goto OUT;
03043             }
03044         }
03045 
03046 /* CASE 2 */
03047 /*check that this is a structure between a cell and a cell*/
03048         if ( structs->cell_i_to != 0 && structs->cell_j_to != 0  
03049              && structs->cell_i_fr != 0 && structs->cell_j_fr != 0 ) { 
03050 
03051             cell_to = T(structs->cell_i_to,structs->cell_j_to); 
03052             cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);
03053 
03054                 /* first check for historical (or other data) flow */
03055             if ( structs->flag == 1 ) {
03056                     /* FLOW */
03057                     /* don't do much with zero or negative flows */
03058                 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03059                 if (structs->flow == 0.0) {
03060                     structs->conc_P = 0.0;
03061                     structs->conc_N = 0.0;
03062                     structs->conc_S = 0.0;
03063                     goto OUT;    /* for output purposes only, zero the reported conc */
03064                 }
03065                 else if (structs->flow < 0.0) { /* all flows should be positive, but check */
03066                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from cell (%d,%d)",
03067                              SimTime.TIME, structs->cell_i_to, structs->cell_j_to);
03068                     WriteMsg( msgStr,True );
03069                     goto OUT;
03070                 }
03071                     /* for internal cells, check for excessive historical/model demand */
03072                 if (structs->cell_i_fr != 2 ) 
03073                 {
03074 
03075                     cell_vol = SWH[cell_fr] * CELL_SIZE;
03076                     if (structs->flow > cell_vol) { /* if hist/other-data demand > cell vol */
03077                         if( debug > 0 ) { 
03078                             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",
03079                                      SimTime.TIME, structs->S_nam, structs->flow, structs->cell_i_fr, structs->cell_j_fr, cell_vol, (structs->flow - cell_vol) );
03080                            
03081                         }
03082                             
03083                         structs->flow = cell_vol;
03084                     }}
03085                     /* CONCENTRATION */
03086                     /* now calculate the nutrient conc in the donor cell */
03087                     /* nut conc in flow from outside system using historical data or 0.0 */
03088                 if (structs->cell_i_fr == 2) { 
03089                     structs->conc_P = (structs->TPser == 1)  ?
03090                         Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
03091                         ((structs->TPser == 0) ? (structs->TP) : (0.0) );
03092                     structs->conc_N = (structs->TNser == 1)  ?
03093                         Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
03094                         ((structs->TNser == 0) ? (structs->TN) : (0.0) );
03095                     structs->conc_S = (structs->TSser == 1)  ?
03096                         Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
03097                         ((structs->TSser == 0) ? (structs->TS) : (0.0) );
03098                 }
03099                 else
03100                 { /* nut conc in flow within system - check to see if there are historical concs with the historical flow */
03101                     structs->conc_P = (structs->TPser == -1)  ?
03102                         ( (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03103                         ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
03104                     structs->conc_N = ( structs->TNser == -1)  ?
03105                         ( (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03106                         ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
03107                     structs->conc_S = ( structs->TSser == -1)  ?
03108                         ( (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03109                         ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
03110                 }
03111                     /* the nut and water update done after the rule-based flow calcs */
03112 /* determine mass transfer */
03113                     
03114                 P_fl = structs->flow * structs->conc_P;
03115                 N_fl = structs->flow * structs->conc_N;
03116                 S_fl = structs->flow * structs->conc_S;
03117                 
03118                 
03119                 if (structs->cell_i_to == 2) { /* "to" cell is outside system */
03120                     if (debug > -999) {/* m3 and kg nuts outflow out from the system */
03121                         VOL_OUT_STR[basn[cell_fr]] += structs->flow; /* mass balance */
03122                         VOL_OUT_STR[0] += structs->flow; 
03123                         P_OUT_STR[basn[cell_fr]] += P_fl; 
03124                         P_OUT_STR[0] += P_fl; 
03125                         S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0); /* see below note on tracer application */
03126                         S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0); /* see below note on tracer application */
03127                     } 
03128                     PA[cell_fr] -= P_fl;
03129                     NA[cell_fr] -= N_fl; 
03130                     SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03131                         /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
03132                            Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
03133                            that was not calculated from the available water and salt (i.e., the conc. was from the
03134                            CanalData.struct data file, either fixed or a time series) */
03135                         /* recalculate the surface water depth in the "from" cell */
03136                     SWH[cell_fr] -= structs->flow/CELL_SIZE;
03137                 }
03138                 else if (structs->cell_i_fr == 2) {  /* "from" cell is outside system */
03139                     if (debug > -999) {/* m3 and kg nuts inflow to the system */
03140                         VOL_IN_STR[basn[cell_to]] += structs->flow; 
03141                         VOL_IN_STR[0] += structs->flow; 
03142                         P_IN_STR[basn[cell_to]] += P_fl; 
03143                         P_IN_STR[0] += P_fl; 
03144                         S_IN_STR[basn[cell_to]] += S_fl; 
03145                         S_IN_STR[0] += S_fl; 
03146                     } 
03147                     PA[cell_to] += P_fl;
03148                     NA[cell_to] += N_fl;
03149                     SA[cell_to] += S_fl;
03150                         /* recalculate the surface water depth in the "to" cell */
03151                     SWH[cell_to] += structs->flow/CELL_SIZE;
03152                 }
03153                 else { /* internal cell-cell m3 and kg nuts flow */
03154                     PA[cell_fr] -= P_fl;
03155                     NA[cell_fr] -= N_fl;        
03156                     SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03157                         /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
03158                            Thus, in that case, we don't subtract salt from the cell because this is "introduced" salt
03159                            that was not calculated from the available water and salt (i.e., the conc. was from the
03160                            CanalData.struct data file, either fixed or a time series) */
03161                     
03162                     PA[cell_to] += P_fl;
03163                     NA[cell_to] += N_fl;
03164                     SA[cell_to] += S_fl;
03165                         /* recalculate the surface water depth in the interacting cells */
03166                     SWH[cell_fr] -= structs->flow/CELL_SIZE;
03167                     SWH[cell_to] += structs->flow/CELL_SIZE;
03168                     if (basn[cell_to] != basn[cell_fr] && debug > -999 ) {
03169                         VOL_OUT_STR[basn[cell_fr]] += structs->flow; 
03170                         VOL_IN_STR[basn[cell_to]] += structs->flow; 
03171                         P_OUT_STR[basn[cell_fr]] += P_fl; 
03172                         P_IN_STR[basn[cell_to]] += P_fl; 
03173                         S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
03174                         S_IN_STR[basn[cell_to]] += S_fl; 
03175 
03176                             /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
03177                         if ( (basn_list[basn[cell_fr]]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
03178                             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",
03179                                      SimTime.TIME, structs->cell_i_fr,structs->cell_j_fr,structs->cell_i_to,structs->cell_j_to );
03180                             WriteMsg( msgStr,True ); dynERRORnum++;
03181                         }
03182                     }
03183                     
03184                 }
03185                 goto OUT;
03186             }
03187             
03188                 /**** Rule- based flows */
03189             else {  
03190                 if ( structs->w_coef != 0.0)  {/* check for a virtual structure  */
03191                     printf ("\nVirtual structs are the only rule-based cell-cell flows in this version!\n"); exit(-1);
03192                 }
03193                 else {
03194                         /* calc'd flow, cell-cell, using the functions defined in Fluxes.c for overland flow, bridge-flow only app right now */
03195                         /* MCopen is a (whole array) of Manning's coefs that are initialized at an open-water value */
03196                     FFlux = Flux_SWcells(structs->cell_i_fr,structs->cell_j_fr,
03197                                       structs->cell_i_to,structs->cell_j_to,SWH,Elev,MCopen);
03198                         /* flow may be negative */
03199                                 /* FFlux units = m */
03200                     structs->flow = FFlux*CELL_SIZE;
03201                     structs->conc_P = (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ;
03202 /*                     structs->conc_N = (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ; */
03203                     structs->conc_S = (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ;
03204 
03205                     if (FFlux != 0.0) Flux_SWstuff ( structs->cell_i_fr,structs->cell_j_fr,
03206                                                   structs->cell_i_to,structs->cell_j_to,FFlux,SWH,SA,NA,PA);
03207                     
03208                     
03209                 }
03210                 goto OUT;
03211             } /* end of rule-based flow calc */
03212             
03213         }
03214 
03215 /* CASE 3 */
03216 /*check that this is a structure between a canal and a cell*/
03217         if ( structs->canal_fr  != 0 && structs->cell_i_to != 0 && structs->cell_j_to != 0 ) { 
03218             can_fr = structs->canal_fr  - CHo;
03219             cell_to = T(structs->cell_i_to,structs->cell_j_to);           
03220             cell_struct = T(structs->str_cell_i,structs->str_cell_j);
03221 
03222                 /* first check for historical (or other data) flow */
03223             if ( structs->flag == 1 ) {
03224 
03225                 if (structs->multiOut == 1) goto OUT; /* this structure has already been processed as a multiple outflow from a canal */
03226                 
03227                 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03228                     /* don't do much with zero or negative flows */           
03229                 if (structs->flow == 0.0) {/* zero the reported struct flow conc (for output info only) */
03230                     structs->conc_P = 0.0;
03231                     structs->conc_N = 0.0;
03232                     structs->conc_S = 0.0;
03233                     goto OUT; 
03234                 }
03235                 else if (structs->flow < 0.0) { 
03236                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from cell (%d,%d)",
03237                              SimTime.TIME, structs->cell_i_to, structs->cell_j_to );
03238                     WriteMsg( msgStr,True ); dynERRORnum++;
03239                     goto OUT;
03240                 }
03241                 
03242                     /* for multiple historical/other-data outflows from a canal, ensure mass balance (sum of outflows !> avail volume) */
03243                     /* all flows from the canal associated with each set must be non-negative, which is currently valid */
03244 
03245                 struct_hist_start = structs; /* this is the first structure w/ historical info for a particular canal */
03246                 ChanHistOut = 0.0; /* sum of all historical outflows from a canal */
03247 
03248                 while ( structs != NULL )   
03249                 {  
03250                     if  ( structs->flag == 1 && struct_hist_start->canal_fr  == structs->canal_fr   ) { /* look for outflows from same canal */
03251                         structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03252                         ChanHistOut += structs->flow; /* sum the multiple historical outflows from canal */
03253                     }
03254                         
03255                     structs = structs->next_in_list; /* increment to next structure */
03256                 } /* end cycle thru list of remaining structures in total list */
03257                     
03258                 structs = struct_hist_start; /* now we are pointing back to the first (or only) water control structure draining the canal */                 
03259 
03260                     /* concentration in donor canal for structure flow, g/L = kg/m^3 */
03261                     /* this is based on previous depth/volume calc'd in FluxChannel */
03262                     /* or use historical conc */
03263                         /* canal TOTAL volume before new structure flows */
03264                 CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth; 
03265                     /* concentrations applied to all structures draining this canal */
03266                 Pconc = (structs->TPser == -1) ? 
03267                     ( (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0) ) :
03268                     ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) ) ;
03269                 Nconc = (structs->TNser == -1) ? 
03270                     ( (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0) ) :
03271                     ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) ) ;
03272                 Sconc = (structs->TSser == -1) ? 
03273                     ( (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0) ) :
03274                     ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ) ;
03275 
03276                    /* vol avail for fluxing */
03277                 CH_vol =  Chan_list[can_fr]->area*(ramp(Chan_list[can_fr]->wat_depth-MINDEPTH) );
03278                 chan_over = (ChanHistOut > 0.0 ) ? ( CH_vol/ChanHistOut ) : 1.0;
03279 
03280                 do /* loop again through the remaining sequence of structures to correct any outflows for the canal */
03281                 {   
03282                     if ( structs->flag == 1 && struct_hist_start->canal_fr  == structs->canal_fr   ) { 
03283 
03284                         structs->multiOut = 1; /* flag to indicate this structure has been processed (below) */
03285 
03286                             /* revised FLOW */
03287                         if ( chan_over < 1.0 ) { /* if hist/other-data demand > canal vol */
03288                             if( debug > 0 ) { 
03289                                 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",
03290                                          SimTime.TIME, structs->S_nam, ChanHistOut, Chan_list[can_fr]->number,
03291                                          CH_vol, (ChanHistOut - CH_vol), Chan_list[can_fr]->wat_depth );
03292                                 
03293                             }  
03294                                 /* decrement the individual structure outflows by it's proportional contrib. */
03295                             structs->flow = (structs->flow) * chan_over ;
03296                         } 
03297                             /* sum of all realized (corrected for excess demand) historical/other-data outflows from a canal during an iteration */
03298                             /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
03299                         Chan_list[can_fr]->sumHistOut += structs->flow;         
03300 
03301                             /* the recipient changes as we cycle thru the structures (may be another cell or a canal) */
03302                         can_to = structs->canal_to - CHo;
03303                         cell_to = T(structs->cell_i_to,structs->cell_j_to); 
03304 
03305                         structs->conc_P = Pconc;/* concentrations applied to all structures draining this canal */
03306                         structs->conc_N = Nconc;
03307                         structs->conc_S = (structs->TSser == -1) ?
03308                             (Sconc) :
03309                             ( (structs->TSser == 0) ? (structs->TS) :
03310                             ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ); /* see note below on tracer conc. application */
03311                          
03312                         P_fl = structs->flow * structs->conc_P; 
03313                         N_fl = structs->flow * structs->conc_N;
03314                         S_fl = structs->flow * structs->conc_S;
03315                                 /* change in nut mass in fr canal (kg) */
03316                         Chan_list[can_fr]->N -=  N_fl ;
03317                         Chan_list[can_fr]->P -=  P_fl ;
03318                         Chan_list[can_fr]->S -=  (structs->TSser == -1) ? (S_fl) : (0);
03319                             /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
03320                                Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
03321                                that was not calculated from the available water and salt (i.e., the conc. was from the
03322                                CanalData.struct data file, either fixed or a time series) */
03323   
03324                         if (can_to > 0) {    /* IF this is a canal recip, change in nut mass in to canal (kg) */
03325                             /* sum of all realized (corrected for excess demand) historical/other-data outflows into a canal during an iteration */
03326                             /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
03327                                 Chan_list[can_to]->sumHistIn += structs->flow;          
03328 
03329                             Chan_list[can_to]->N +=  N_fl ;
03330                             Chan_list[can_to]->P +=  P_fl ;
03331                             Chan_list[can_to]->S +=  S_fl ;
03332                                 /* mass balance in fr and to canals */
03333                             if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > -999 ) {
03334                                 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03335                                 VOL_IN_STR[Chan_list [can_to]->basin] += structs->flow; 
03336                                 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03337                                 P_IN_STR[Chan_list [can_to]->basin] += P_fl; 
03338                                 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
03339                                 S_IN_STR[Chan_list [can_to]->basin] += S_fl; 
03340                             }
03341                         }
03342                         else if (cell_to > 0) {
03343                                 /* OR change in nut mass in downstream cell (kg) */
03344  
03345                             if (structs->cell_i_to > 2 ) { /* to an on-map cell */
03346                                 SWH[cell_to] += structs->flow/CELL_SIZE;
03347                                 PA[cell_to] += P_fl; 
03348                                 NA[cell_to] += N_fl; 
03349                                 SA[cell_to] += S_fl; 
03350                                     /* mass balance in fr canal and to cell */
03351                                 if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
03352                                     VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03353                                     VOL_IN_STR[basn[cell_to]] += structs->flow; 
03354                                     P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03355                                     P_IN_STR[basn[cell_to]] += P_fl; 
03356                                     S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
03357                                     S_IN_STR[basn[cell_to]] += S_fl; 
03358 
03359                                         /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
03360                                     if ( (Chan_list[can_fr]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
03361                                         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",
03362                                                  SimTime.TIME, Chan_list[can_fr]->number,structs->cell_i_to,structs->cell_j_to );
03363                                         WriteMsg( msgStr,True ); dynERRORnum++;
03364                                     }
03365                                 }
03366                             }
03367                             else if (structs->cell_i_to == 2 && debug > -999 ) { /* to off-map cell, mass balance check */
03368                                 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow; 
03369                                 VOL_OUT_STR[0] += structs->flow; 
03370                                 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl; 
03371                                 P_OUT_STR[0] += P_fl; 
03372                                 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
03373                                 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
03374                             }
03375                         }
03376                             /* don't do much with zero or negative flows */           
03377                         if (structs->flow == 0.0) {
03378                             structs->conc_P = 0.0;
03379                             structs->conc_N = 0.0;
03380                             structs->conc_S = 0.0;
03381                         } /* for output purposes only, zero the reported conc */
03382                         else if (structs->flow < 0.0) {
03383                             sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from another chan/cell into chan %d",
03384                                      SimTime.TIME, Chan_list[can_fr]->number );
03385                             WriteMsg( msgStr,True ); dynERRORnum++;
03386                         }
03387                     }
03388 
03389                     structs = structs->next_in_list;
03390                 } while (structs != NULL); /*  cycling to the end of the list again !!!! */
03391 
03392                 structs = struct_hist_start; /* now pointing back to the first historical outflow structure that we dealt with */
03393     
03394 
03395                 goto OUT;
03396                     /* finished check for multiple (positive) outflows exceeding canal volume */
03397             }
03398                 /* FLOW */
03399                 /**** Rule- based flows */
03400             else { 
03401                 HeadH = Elev [cell_struct] - Chan_list[can_fr]->depth + Chan_list[can_fr]->wat_depth
03402                                 + (Chan_list[can_fr]->sumHistIn - Chan_list[can_fr]->sumHistOut)/Chan_list[can_fr]->area; 
03403           /* in v2.3, the rule-based canal->cell flows use either tide data or SFWMM/other model data for boundary conditions */
03404                 
03405                /* v2.2 and previous eqn: HeadT = ( ( structs->cell_i_to == 2 ) ? (Elev[cell_struct] - 0.1) : (SWH[cell_to]+Elev[cell_to]) ); */
03406                 /* Head at Tailwater location: check that destination is out of model domain (always in v2.3 and prior versions)
03407                         The TWf (TailWater flag) indicates use of annually-recurring tidal data, with Sea Level Rise added.
03408                         If the TWf is false, use relative depth data from the SFWMM or other model output (boundcond_depth data is NOT stage, but depth relative to elev NGVD'29) */
03409                         /* v2.5 use cell_i_TW, cell_j_TW in cell_to location */
03410                 HeadT = ( ( structs->cell_i_to == 2 ) ? 
03411                         ( (TWf) ? (structs->TW_stage) : ( boundcond_depth[T(structs->cell_i_TW,structs->cell_j_TW)]  + Elev[cell_struct]) ): 
03412                         (SWH[cell_to]+Elev[cell_to]) );
03413                                 /* volume available for fluxing - right now, the only structures here drain border seepage canals, leave 30 cm water */
03414                 CH_vol =  Chan_list[can_fr]->area* ramp(HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth-0.3); /* avail flow volume */
03415             
03416                 if (HeadH>HeadT) 
03417                 {
03418                     /* this ****DOES NOT**** work when the HW and TW checks are being made at a remote cell */
03419                     /* ok for now, just fix when driving ELM with rules and checking remote cells */
03420                     structs->flow =  Min(structs->w_coef * ( HeadH - HeadT ) * canstep, CH_vol );                   
03421                         /* all rule based instances of canal->cell flow  go
03422                            to external cells, so this reversal check not invoked */
03423                     flow_rev = ( structs->cell_i_to != 2 ) ? 
03424                         (HeadH - structs->flow/Chan_list[can_fr]->area) - (HeadT + structs->flow/CELL_SIZE) :
03425                         (0.0);
03426                         /* the flux eqn is the same as the (new Jul98) virtual weir eqn, equilibrating the heads across two unequal area regions */
03427                     /* flow is only in positive directioin */
03428                     if (flow_rev < 0.0) {  structs->flow = 
03429                         ramp( Chan_list[can_fr]->area*CELL_SIZE*HeadH/(Chan_list[can_fr]->area+CELL_SIZE) -
03430                         Chan_list[can_fr]->area*CELL_SIZE*HeadT/(Chan_list[can_fr]->area+CELL_SIZE) );
03431                           }
03432                }
03433                 else structs->flow = 0.0; /* undefined flow if heads are equal, thus this stmt added */
03434                 
03435             
03436             
03437                     /* don't do much with zero or negative flows */           
03438                 if (structs->flow == 0.0) {
03439                     structs->conc_P = 0.0;
03440                     structs->conc_N = 0.0;
03441                     structs->conc_S = 0.0;
03442                     goto OUT;
03443                 } /* for output purposes only, zero the reported conc */
03444                 else if (structs->flow < 0) {  
03445                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG FLOW from cell (%d,%d)",
03446                              SimTime.TIME, structs->cell_i_to, structs->cell_j_to );
03447                     WriteMsg( msgStr,True ); dynERRORnum++;
03448                     goto OUT;
03449                 }
03450           
03451                     /* rule-based flows:  calculate the nutrient concentrations */
03452                     /* conc in canal -> cell flow  */
03453                 CH_vol =  Chan_list[can_fr]->area* ramp(HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth); /* TOTAL volume */
03454                     /* conc in donor (fr) canal (kg/m3) */
03455                 structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
03456                 structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
03457                 structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
03458 
03459                 P_fl = structs->flow * structs->conc_P; 
03460                 N_fl = structs->flow * structs->conc_N;
03461                 S_fl = structs->flow * structs->conc_S;
03462                   
03463                     /* increment the sum of rule-based flows to from canals
03464                        for potential multiple virtual flows associated with some canals */
03465                 Chan_list[can_fr]->sumRuleOut += structs->flow;         
03466                     /* change in nut mass in downstream cell (kg) */
03467                 if (structs->cell_i_to > 2 ) { /* to an on-map cell */
03468                     SWH[cell_to] += structs->flow/CELL_SIZE;
03469                     PA[cell_to] += P_fl; 
03470                     NA[cell_to] += N_fl; 
03471                     SA[cell_to] += S_fl; 
03472                         /* mass balance in fr canal and to cell */
03473                     if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
03474                         VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03475                         VOL_IN_STR[basn[cell_to]] += structs->flow; 
03476                         P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03477                         P_IN_STR[basn[cell_to]] += P_fl; 
03478                         S_OUT_STR[Chan_list[can_fr]->basin] += S_fl;
03479                         S_IN_STR[basn[cell_to]] += S_fl; 
03480 
03481                             /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
03482                         if ( (Chan_list[can_fr]->family ==  basn_list[basn[cell_to]]->family) && (debug > 0) ) {
03483                             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",
03484                                      SimTime.TIME, Chan_list[can_fr]->number,structs->cell_i_to,structs->cell_j_to );
03485                             WriteMsg( msgStr,True ); dynERRORnum++;
03486                         }
03487                     }
03488                 }
03489                 else if (structs->cell_i_to == 2 && debug > -999 ) { /* to off-map cell, mass balance check */
03490                     VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow; 
03491                     VOL_OUT_STR[0] += structs->flow; 
03492                     P_OUT_STR[Chan_list[can_fr]->basin] += P_fl; 
03493                     P_OUT_STR[0] += P_fl; 
03494                     S_OUT_STR[Chan_list[can_fr]->basin] += S_fl; 
03495                     S_OUT_STR[0] += S_fl; 
03496                 }
03497 
03498             
03499                     /* calc the change in mass in  canal (kg) */
03500                 Chan_list[can_fr]->N -=  N_fl ;
03501                 Chan_list[can_fr]->P -=  P_fl ;
03502                 Chan_list[can_fr]->S -=  S_fl ;
03503 
03504                 goto OUT;
03505             }
03506         }
03507         
03508     
03509 
03510 /* CASE 4 */
03511 /*check that this is a structure between a cell and a canal*/
03512             /* if from 1,1 to canal, always will require historical data */
03513         if ( structs->canal_to != 0 && structs->cell_i_fr != 0 && structs->cell_j_fr != 0 ) { 
03514             can_to = structs->canal_to - CHo;
03515             cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);           
03516             cell_struct = T(structs->str_cell_i,structs->str_cell_j);
03517 
03518                 /* first check for historical (or other data) flow */
03519             if ( structs->flag == 1 ) {
03520   
03521                 if (structs->multiOut == 1) {
03522                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
03523                              SimTime.TIME, Chan_list[can_to]->number );
03524                     WriteMsg( msgStr,True ); dynERRORnum++;
03525                     goto OUT; /* this structure has already been processed as a multiple outflow from a canal */
03526                 }
03527                 
03528                 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03529                     /* don't do much with zero or negative flows */           
03530                 if (structs->flow == 0.0) {
03531                     structs->conc_P = 0.0;
03532                     structs->conc_N = 0.0;
03533                     structs->conc_S = 0.0;
03534                     goto OUT; 
03535                 }
03536                 else if (structs->flow <0) { 
03537                     sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
03538                              SimTime.TIME, Chan_list[can_to]->number );
03539                     WriteMsg( msgStr,True ); dynERRORnum++;
03540                     goto OUT;
03541                 } 
03542 
03543                 if (structs->cell_i_fr != 2 ) /* cell (2,2) == (1,1) in maps, not in model domain */
03544                 {
03545 
03546                     cell_vol = SWH[cell_fr] * CELL_SIZE;
03547                     if (structs->flow > cell_vol) {
03548                         if( debug > 0 ) { 
03549                             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",
03550                                      SimTime.TIME, structs->S_nam, structs->flow, structs->cell_i_fr, structs->cell_j_fr, cell_vol, (structs->flow - cell_vol) );
03551                             
03552                         }
03553                         structs->flow = cell_vol;
03554                     }}
03555                     /* calculate the nutrient conc in the donor cell */
03556                     /* nut conc in flow from outside system using historical data or 0.0 */
03557                 if (structs->cell_i_fr == 2) { 
03558                     structs->conc_P = (structs->TPser == 1)  ?
03559                         Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
03560                         ((structs->TPser == 0) ? (structs->TP) : (0.0) );
03561                     structs->conc_N = (structs->TNser == 1)  ?
03562                         Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
03563                         ((structs->TNser == 0) ? (structs->TN) : (0.0) );
03564                     structs->conc_S = (structs->TSser == 1)  ?
03565                         Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
03566                         ((structs->TSser == 0) ? (structs->TS) : (0.0) );
03567                 }
03568                 else
03569                 { /* nut conc in flow within system - check to see if there are historical concs with the historical flow */
03570                     structs->conc_P = (structs->TPser == -1)  ?
03571                         ( (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03572                         ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
03573                     structs->conc_N = ( structs->TNser == -1)  ?
03574                         ( (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03575                         ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
03576                     structs->conc_S = ( structs->TSser == -1)  ?
03577                         ( (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03578                         ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
03579                 }
03580                     /* now determine the nutrient mass (kg) flow from cell to canal */
03581                     /* 1000 is to convert g/m3 to kg */ 
03582                 P_fl = structs->flow * structs->conc_P;
03583                 N_fl = structs->flow * structs->conc_N;
03584                 S_fl = structs->flow * structs->conc_S;
03585                 if (structs->cell_i_fr == 2 ) {/* "from" cell is outside system */
03586                     if (debug > -999) {  
03587                         VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; /* mass balance */
03588                         VOL_IN_STR[0] += structs->flow; 
03589                         P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
03590                         P_IN_STR[0] += P_fl; 
03591                         S_IN_STR[Chan_list[can_to]->basin] += S_fl; 
03592                         S_IN_STR[0] += S_fl; 
03593                     } }
03594                 else { /* internal cell-> canal flow */
03595                     SWH[cell_fr] -= structs->flow/CELL_SIZE;
03596                     PA[cell_fr] -= P_fl; 
03597                     NA[cell_fr] -= N_fl; 
03598                     SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03599                             /* for the salt (S), we sometimes want to apply a tracer to internal canals/cells.
03600                                Thus, in that case, we don't subtract salt from the canal because this is "introduced" salt
03601                                that was not calculated from the available water and salt (i.e., the conc. was from the
03602                                CanalData.struct data file, either fixed or a time series) */
03603                         /* mass balance in fr cell and to canal */
03604                     if (basn[cell_fr] != Chan_list[can_to]->basin && debug > -999 ) {
03605                         VOL_OUT_STR[basn[cell_fr]] += structs->flow;
03606                         VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; 
03607                         P_OUT_STR[basn[cell_fr]] += P_fl;
03608                         P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
03609                         S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0); /* see above note on tracer application */
03610                         S_IN_STR[Chan_list[can_to]->basin] += S_fl; 
03611 
03612                             /* all structure flows, EXCEPT canal-canal, must be between parent (family) hydro basins */
03613                         if ( (Chan_list[can_to]->family ==  basn_list[basn[cell_fr]]->family) && (debug > 0) ) {
03614                             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",
03615                                      SimTime.TIME, structs->cell_i_fr, structs->cell_j_fr, Chan_list[can_to]->number );
03616                             WriteMsg( msgStr,True ); dynERRORnum++;
03617                         }
03618                     }
03619                 }
03620                 
03621                   /* sum of all realized (corrected for excess demand) historical/other-data outflows into a canal during an iteration */
03622                   /* used to determine any allowable rule-based canal flows after the historical/other-data demands */
03623                 Chan_list[can_to]->sumHistIn += structs->flow;          
03624 
03625                    /* change in nut mass in "to" canal (kg) */
03626                 Chan_list[can_to]->N +=  N_fl ;
03627                 Chan_list[can_to]->P +=  P_fl ;
03628                 Chan_list[can_to]->S +=  S_fl ;
03629 
03630                 goto OUT;
03631             }
03632 
03633             else { /* rule based flow calc */
03634                 if (structs->cell_i_fr != 2  ) {
03635                     sprintf( msgStr, "Day %6.1f: ERROR in flow configuration for structure from cell (%d,%d) to chan %d: this rule-based flow is restricted to flow from outside model domain.",
03636                              SimTime.TIME, structs->cell_i_fr, structs->cell_j_fr, Chan_list[can_to]->number );
03637                     WriteMsg( msgStr,True );   dynERRORnum++;    
03638                 }
03639                 HeadH = (structs->HW_stage) ; 
03640           /* in v2.3, the rule-based cell->canal flows use only tide data for boundary conditions */
03641                 HeadT = Elev [cell_struct] - Chan_list[can_to]->depth + Chan_list[can_to]->wat_depth
03642                     + (Chan_list[can_to]->sumHistIn - Chan_list[can_to]->sumHistOut)/Chan_list[can_to]->area;
03643             
03644                 if (HeadH>HeadT) 
03645                 {
03646                     structs->flow =  Max(structs->w_coef * ( HeadH - HeadT ) * canstep, 0 );                   
03647                 }
03648                 else structs->flow = 0.0; /* undefined flow if heads are equal, thus this stmt added */
03649                 
03650                 /* don't do much with zero or negative flows */           
03651             if (structs->flow == 0.0) {
03652                 structs->conc_P = 0.0;
03653                 structs->conc_N = 0.0;
03654                 structs->conc_S = 0.0;
03655                 goto OUT;
03656             } /* for output purposes only, zero the reported conc */
03657 
03658                 /*  rule-based: calculate the nutrients - this is restricted to calculated tidal boundary condition flows  from outside of system */
03659                 /* concentration in structure flows, g/L = kg/m^3; mass of stock in kg */
03660 
03661                 if (structs->cell_i_fr == 2) { 
03662                     structs->conc_P = (structs->TPser == 1)  ?
03663                         Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
03664                         ((structs->TPser == 0) ? (structs->TP) : (0.0) );
03665                     structs->conc_N = (structs->TNser == 1)  ?
03666                         Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
03667                         ((structs->TNser == 0) ? (structs->TN) : (0.0) );
03668                     structs->conc_S = (structs->TSser == 1)  ?
03669                         Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
03670                         ((structs->TSser == 0) ? (structs->TS) : (0.0) );
03671                 }
03672             P_fl = structs->flow * structs->conc_P; 
03673             N_fl = structs->flow * structs->conc_N;
03674             S_fl = structs->flow * structs->conc_S;
03675              
03676                 /* this calculated flow is restricted to an outside-system cell, which does not need updating */
03677                 /* these are nutrient stocks, not concentrations */  
03678                 /* mass balance in fr cell and to canal */
03679             if (structs->cell_i_fr == 2  ) {
03680                 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow; 
03681                 VOL_IN_STR[0] += structs->flow; 
03682                 P_IN_STR[Chan_list[can_to]->basin] += P_fl; 
03683                 P_IN_STR[0] += P_fl; 
03684                 S_IN_STR[Chan_list[can_to]->basin] += S_fl; 
03685                 S_IN_STR[0] += S_fl; 
03686             }
03687             else        {  
03688                 sprintf( msgStr, "Day %6.1f: ERROR in flow to chan %d: cannot have calculated cell->canal flows except from outside of model domain. ",
03689                          SimTime.TIME, Chan_list[can_to]->number);
03690                 WriteMsg( msgStr,True ); dynERRORnum++;
03691                 goto OUT; 
03692             }
03693             
03694 
03695                 /* change in nut mass in canal (kg) */
03696             Chan_list[can_to]->N +=  N_fl ;
03697             Chan_list[can_to]->P +=  P_fl ;
03698             Chan_list[can_to]->S +=  S_fl ;
03699             
03700             goto OUT;
03701             }
03702         }
03703         
03704  
03705 /* CASE 5 */
03706 /* if none of the above structure flow conditions */
03707         printf ("\nError in data for structure N %s: impossible condition! \n", structs->S_nam );
03708         exit (-1);
03709         
03710       OUT:
03711         if (TWf) { structs->TW_stage = -99; TWf = 0; }
03712         if (HWf) { structs->HW_stage = -99; HWf = 0; }
03713 
03714             /* sum of flows and solute masses thru structs over the can_Intvl canal summary interval */
03715         structs->SumFlow += structs->flow; 
03716         structs->Sum_P += structs->conc_P * structs->flow; 
03717 /*        structs->Sum_N += structs->conc_N * structs->flow; */
03718         structs->Sum_S += structs->conc_S * structs->flow; 
03719 
03720         structs->multiOut = 0; /* reset the multi-outflow flag to zero after done */
03721         
03722             /* print sum of flows and the flow weighted mean conc thru structures */
03723         if ( printchan ) { 
03724             fprintf( WstructOutFile, "%7.0f\t", (structs->SumFlow)/(2446.848) ); /* convert to cfs (0.02832*60*60*24) */
03725             fprintf( WstructOutFile_P, "%7.4f\t", 
03726                      (structs->SumFlow > 0.0) ? (structs->Sum_P/structs->SumFlow*1000.0) : (0.0) ); /* convert g/L to mg/L */
03727 /*            fprintf( WstructOutFile_N, "%7.4f\t", 
03728               (structs->SumFlow > 0.0) ? (structs->Sum_N/structs->SumFlow*1000.0) : (0.0) );*/ /* convert g/L to mg/L */
03729             fprintf( WstructOutFile_S, "%7.4f\t", 
03730                      (structs->SumFlow > 0.0) ? (structs->Sum_S/structs->SumFlow) : (0.0) ); /* g/L, no conversion */
03731 
03732             structs->SumFlow = 0.0; 
03733             structs->Sum_P = 0.0;
03734 /*            structs->Sum_N = 0.0; */
03735             structs->Sum_S = 0.0;
03736         }
03737         
03738             /* increment to the next structure in the list */
03739         structs = structs->next_in_list;
03740     }
03741 
03742     if (printchan) {
03743         fflush( WstructOutFile );
03744         fflush( WstructOutFile_P );
03745 /*         fflush( WstructOutFile_N ); */
03746         fflush( WstructOutFile_S );
03747     }
03748  
03749     return;
03750     
03751 }
03752 
03753 /****************************************************************************************/
03762 float GetGraph (struct Schedule *graph, float x)
03763 {
03764     float s;
03765     int ig = 0, Np;
03766         
03767     Np = graph->num_point;
03768 
03769     while(1) 
03770     {
03771         if (x <= graph->graph_points[ig].time) 
03772         { if ( ig > 0 ) ig--; 
03773         else return ( graph->graph_points[0].value );
03774         }
03775         else if (x > graph->graph_points[ig].time && x <= graph->graph_points[ig+1].time) 
03776         {
03777             s =         ( graph->graph_points[ig+1].value - graph->graph_points[ig].value )/
03778                 ( graph->graph_points[ig+1].time - graph->graph_points[ig].time );
03779             return ( s * ( x - graph->graph_points[ig].time ) + graph->graph_points[ig].value ); 
03780         }
03781         else if (x > graph->graph_points[ig+1].time) 
03782         { if ( ig < Np ) ig++; 
03783         else return ( graph->graph_points[Np].value );
03784         }
03785     }
03786 }
03787 
03788 
03789 /********************************************************************************************/
03795 float UTM2kmy ( double UTM )
03796 { return ( Abs(UTM - UTM_EOrigin)/celWid );  
03797 }
03798 
03804 float UTM2kmx ( double UTM )
03805 { return ( Abs(UTM - UTM_NOrigin)/celWid );              
03806 }
03807 
03808 /****************************************************************************************/
03815 int Wrdcmp ( char *s, char *t )
03816 { 
03817   for ( ; *s == *t; s++, t++ )
03818          if ( isspace (*(s+1)) || *(s+1) == EOS ) return 1;
03819         
03820   return 0;
03821   
03822 }
03823 
03824 /******/

Generated on Thu Jul 6 11:17:57 2006 for ELM source code by  doxygen 1.3.9.1