00001
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
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
00081
00082
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;
00088 Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00089 Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00090 }
00091
00092 usrErr ("\tCanal geometry configured OK");
00093
00094
00095 ReadStructures (filename, baseDatum);
00096
00097 ON_MAP[T(2,2)] = 0;
00098
00099 MCopen = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MCopen");
00100 init_pvar(MCopen,NULL,'f',0.05);
00101
00102
00103 if ( debug > 0 )
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
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
00131 if ( ( WstructOutFile = fopen( modelFileName, "w" ) ) == NULL )
00132 {
00133 printf( "Can't open %s file!\n", modelFileName );
00134 exit(-1) ;
00135 }
00136
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
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
00152 if ( ( WstructOutFile_P = fopen( modelFileName, "w" ) ) == NULL )
00153 {
00154 printf( "Can't open %s file!\n", modelFileName );
00155 exit(-1) ;
00156 }
00157
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
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
00173 if ( ( WstructOutFile_S = fopen( modelFileName, "w" ) ) == NULL )
00174 {
00175 printf( "Can't open %s file!\n", modelFileName );
00176 exit(-1) ;
00177 }
00178
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
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 );
00202
00203
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
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 );
00220
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
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 );
00236
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
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
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" );
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;
00276 Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00277 Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
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);
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
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
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 {
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
00342 Flows_in_Structures ( SWH, Elevation, MC, NA, PA, SA );
00343
00344
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
00352
00353 for( i = 0; i < num_chan; i++ )
00354 if (Chan_list[i]->width > 0)
00355
00356 {
00357 FluxChannel ( i, SWH, Elevation, MC, GW, poros, GWcond, NA, PA, SA, GNA, GPA, GSA, Unsat, sp_yield );
00358
00359 CH_vol = Chan_list[i]->area*Chan_list[i]->wat_depth;
00360
00361
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
00375 Chan_list[i]->P_con = (CH_vol > 0) ? (Chan_list[i]->P/CH_vol*1000.0) : 0.0;
00376
00377 Chan_list[i]->S_con = (CH_vol > 0) ? (Chan_list[i]->S/CH_vol ) : 0.0;
00378
00379
00380 fprintf( CanalOutFile, "%6.2f\t", Chan_list[i]->wat_depth );
00381 fprintf( CanalOutFile_P, "%7.4f\t", Chan_list[i]->P_con );
00382
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
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
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
00431
00432
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
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 );
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 );
00458 fgets( ss, 420, ChanInFile ); sscanf ( ss, "%f %f %f ", &F_ERROR, &MINDEPTH, &C_F );
00459
00460 fgets( ss, 420, ChanInFile );
00461
00462
00463 while ( fgets( ss, 420, ChanInFile ) != NULL && !feof( ChanInFile ) )
00464 {
00465 if ( *ss == '#' )
00466 {
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
00474
00475
00476
00477
00478
00479
00480
00481
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
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;
00509 channels->ic_N_con = 0.00;
00510 channels->ic_S_con *= 1.0;
00511
00512 channels->ic_depth = channels->wat_depth;
00513
00514
00515 channels->reaches = C_reach;
00516
00517 if (i < 10)
00518 {printf ( " Error in CanalData.chan: reach input %d\n", (num_chan+1) ); exit (-1);}
00519
00520
00521 num_chan++;
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 {
00538 i = sscanf ( ss, "%f %f", &C_x1, &C_y1 );
00539 if (i < 2)
00540 {printf ( " Error in CanalData.chan coords: %d'th reach input \n", (num_chan+1) ); exit (-1);}
00541
00542
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
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 }
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
00618 if ( ( ChanInFile = fopen( modelFileName, "r" ) ) == NULL )
00619 {
00620 printf( "Can't open %s file!\n", modelFileName ) ;
00621 exit(-1) ;
00622 }
00623
00624
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);
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
00646 fgets( ss, 320, ChanInFile );
00647
00648
00649
00650 while ( fgets( ss, 320, ChanInFile ) != NULL && !feof( ChanInFile ) )
00651 {
00652 s = ss;
00653
00654
00655 sscanf( s, "%d", &structs->flag );
00656
00657
00658
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
00667
00668
00669
00670
00671
00672
00673
00674
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
00685
00686 s = Scip( s, TAB );
00687 if ( *s == TAB ) structs->TPser = -1;
00688 else {
00689 sscanf ( s, "%s", &histTP );
00690
00691 if (isInteger(histTP) ) {
00692 structs->TPser = 0;
00693 structs->TP = atof(histTP)*1.0e-6;
00694 }
00695 else { structs->TPser = 1; numTPser++;}
00696 }
00697
00698 s = Scip( s, TAB );
00699 if ( *s == TAB ) structs->TNser = -1;
00700 else {
00701 sscanf ( s, "%f", &structs->TN );
00702 structs->TNser = 0;
00703 structs->TN *= 1.0e-6;
00704 }
00705
00706
00707 s = Scip( s, TAB );
00708 if ( *s == TAB ) structs->TSser = -1;
00709
00710 else {
00711 sscanf ( s, "%s", &histTS );
00712
00713 if (isFloat(histTS) ) {
00714 structs->TSser = 0;
00715 structs->TS = atof(histTS)*1.0;
00716 }
00717 else { structs->TSser = 1; numTSser++;}
00718
00719 }
00720
00721 s = Scip( s, TAB );
00722 if ( *s == TAB ) structs->str_cell_i = 0;
00723 else {
00724 sscanf( s, "%d", &structs->str_cell_i );
00725 structs->str_cell_i++;
00726 }
00727
00728 s = Scip( s, TAB );
00729
00730 if ( *s == TAB ) structs->str_cell_j = 0;
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);
00739 }
00740
00741
00742 }
00743
00744 s = Scip( s, TAB );
00745
00746 if ( *s == TAB ) structs->canal_fr = 0;
00747 else sscanf( s, "%d", &structs->canal_fr );
00748 s = Scip( s, TAB );
00749 if ( *s == TAB) structs->canal_to = 0;
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);
00756 }
00757
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;
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);
00784 }
00785
00786
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;
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);
00806 }
00807
00808
00809
00810 if ( structs->cell_j_to ) ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] += 100;
00811
00812 s = Scip( s, TAB );
00813 if (*s == TAB) structs->HW_stage = 0;
00814 else if ( isalpha (*s) ) {
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) ) {
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
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);
00851 }
00852
00853 if (*s == TAB) structs->cell_i_TW = 0;
00854 else
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);
00869 }
00870
00871 if (*s == TAB || *s == '?') structs->w_coef = 0;
00872 else sscanf( s, "%f", &structs->w_coef );
00873
00874
00875 struct_last = structs;
00876
00877
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
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
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
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);
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);
00926 line=lline;
00927 line = Scip( line, TAB );
00928 line = Scip( line, TAB );
00929 line = Scip( line, TAB );
00930
00931
00932
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
00941 if ( Hist_data[i].flag > 1 ) {
00942 structs = struct_first;
00943 while ( structs != NULL )
00944 {
00945
00946 if ( structs->flag == 10 * Hist_data[i].flag )
00947 {
00948
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;
00960 strcpy( Hist_data[num_struct_hist+disAggNum-1].S_nam, structs->S_nam );
00961 Hist_data[i].aggCnt++;
00962 structs->flag = 1;
00963 structs->histID = num_struct_hist+disAggNum;
00964
00965 }
00966
00967 structs = structs->next_in_list;
00968 }
00969 }
00970
00971 sscanf ( line,"%s\t",&structName);
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;
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);
00985
00986
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++ ) {
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;
01001 }
01002
01003 for (i = num_struct_hist; i < num_struct_hist+disAggNum; i++) {
01004 k = Hist_data[i].aggID;
01005 Hist_data[i].arrayWat[((int)(lincnt))] =
01006 Hist_data[k].arrayWat[((int)(lincnt))] / (Hist_data[k].aggCnt) ;
01007
01008 }
01009
01010
01011 lincnt++;
01012 }
01013
01014 }
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
01028
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
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);
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);
01052 line=lline;
01053 line = Scip( line, TAB );
01054 line = Scip( line, TAB );
01055 line = Scip( line, TAB );
01056
01057
01058 for ( j = 0; j < numTPser; j++ ) {
01059 IDhist[j] = -1;
01060 strcpy(structName,"NULL");
01061
01062 sscanf ( line,"%s\t",&structName);
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;
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);
01080
01081
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++ ) {
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;
01095
01096 }
01097
01098 lincnt++;
01099 }
01100
01101 }
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 }
01113
01114
01115
01116
01117
01118
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
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);
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);
01142 line=lline;
01143 line = Scip( line, TAB );
01144 line = Scip( line, TAB );
01145 line = Scip( line, TAB );
01146
01147
01148 for ( j = 0; j < numTSser; j++ ) {
01149 IDhist[j] = -1;
01150 strcpy(structName,"NULL");
01151
01152 sscanf ( line,"%s\t",&structName);
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;
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);
01170
01171
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++ ) {
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;
01185
01186 }
01187
01188 lincnt++;
01189 }
01190
01191 }
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 }
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
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 );
01241 fgets( ss, 300, schedFile );
01242 fgets( ss, 300, schedFile );
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
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
01263 while ( *s++ != EOL ) if ( *s == '(' ) count++;
01264
01265
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
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
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 }
01293
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;
01319 int cellLoc, cellLoc_i, cellLoc_j;
01320 int HV, k, L_range, R_range;
01321 float T_length;
01322 float C_length;
01323 float distTot;
01324 float A1, B1, A2, B2;
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;
01331 struct Chan *channel;
01332 struct Chan_reach *Reach;
01333
01334 int numChek;
01335
01336
01337 int *marked;
01338 marked = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"marked");
01339 init_pvar(marked,NULL,'d',0);
01340
01341
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
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
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
01365 channel = channel_first;
01366
01367
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
01384 if (channel->width < 0) goto SKIPCHAN;
01385
01386 fprintf ( CanalCellInteractDebug, "N = %d levee = %d\n", c_num, z );
01387
01388
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)
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) );
01410 C_length += T_length;
01411
01412 if ( a1 != a0 && b1 != b0 )
01413 {
01414 A1 = (b1 - b0)/(a1 - a0);
01415 B1 = b0 - A1*a0;
01416 A2 = 1/A1; B2 = - B1/A1;
01417 }
01418
01419 x = a0; y = b0;
01420 i = Round (x); j = Round (y);
01421
01422
01423 while ( x1 != a1 || y1 != b1 )
01424 {
01425 if ( a1 == a0 )
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 )
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;
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 )
01448 x1 = i + STEP(a1 - a0)*cell_step;
01449 else
01450 { y1 = j + STEP(b1 - b0)*cell_step; length = L;
01451 }
01452 }
01453
01454 if ( Abs (a1 - a0) < Abs (b1 - b0) )
01455
01456 { p_middle = (x1 + x)/2;
01457 m = i + cell_step/2;
01458 HV = 1;
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;
01467
01468
01469
01470 if ( Round (a1) == Round (a0) )
01471
01472 { if (b0 < b1)
01473 {
01474
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
01489 {
01490
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) )
01506 { if (a0 < a1)
01507 {
01508
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
01523 {
01524
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
01541 { if ( A1>0 )
01542 i45 = ( HV ) ? p_middle - m : m - p_middle;
01543 else
01544 i45 = p_middle - m;
01545
01546
01547 if ( b0 < b1 )
01548 { if ( a0 < a1 )
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
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
01578 else
01579 { if ( a0 > a1 )
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
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++;
01613
01614
01615
01616
01617 x = x1; y = y1;
01618 i = Round (x); if ( a1 < a0 && i == x ) i--;
01619
01620 j = Round (y); if ( b1 < b0 && j == y ) j--;
01621 C_Mark = 1;
01622
01623
01624 }
01625
01626 Reach = Reach->next_reach;
01627 }
01628
01629
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
01638
01639
01640
01641
01642
01643
01644
01645
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 ;
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);
01663
01664 SKIPCHAN:
01665 channel = channel->next_in_list;
01666
01667 cell_last->next_cell = NULL;
01668
01669 }
01670
01671 fclose ( CanalCellInteractDebug );
01672
01673 return;
01674 }
01675
01676
01687 void getCanalElev (int chan_n)
01688 {
01689 int i;
01690
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);
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;
01743
01744
01745
01746
01747 if ( !C_Mark ) return cell;
01748
01749
01750
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;
01760 cell->y = y+1;
01761 cell->ind = c_ind;
01762 cell->length = length*celWid;
01763 cell->reachDistDown = distTot;
01764
01765 fprintf ( CanalCellInteractDebug, " %d, %d, l/r = %d, distance= %7.1f\n", cell->x, cell->y, c_ind, cell->reachDistDown );
01766
01767
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
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
01876 int equilib, revers=1;
01877 float H_pot_don, H_pot_rec;
01878
01879 float fluxTOunsat_don, fluxHoriz;
01880 float UnsatZ_don, UnsatCap_don, UnsatPot_don, Sat_vs_unsat;
01881 float UnsatZ_rec, UnsatCap_rec, UnsatPot_rec;
01882 float sfTOsat_don, unsatTOsat_don, unsatTOsat_rec, satTOsf_rec, horizTOsat_rec;
01883
01884 float sedwat_don, sedwat_rec;
01885 float sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec;
01886 float sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don;
01887 float sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz;
01888
01889
01890
01891
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;
01912 MIN_VOL = Chan_list[chan_n]->minVol;
01913 SPG_coef = Chan_list[chan_n]->SPG_flow_coef;
01914
01915 I_Length = Chan_list[chan_n]->seg_len;
01916 seg_area = I_Length * Chan_list[chan_n]->width;
01917 can_cell_M = Chan_list[chan_n]->area*CELL_SIZE;
01918 can_cell_A = Chan_list[chan_n]->area+CELL_SIZE;
01919
01920
01921
01922
01923 do
01924 {
01925 if ( Abs( error) < F_ERROR )
01926 {
01927 converged = 1;
01928
01929 errorConv=error;
01930 factor = 0;
01931 }
01932 else if ( iter >= (max_iter-1) )
01933 {
01934 if (attempt == 4) {
01935 converged = 1;
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++;
01946 iter = 0;
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
01957 if (error * prev_error < 0 && iter != 1)
01958 factor = - sgn(error) * Abs( factor ) * 0.5;
01959
01960
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
01977 factor = -factor*0.1;
01978 }
01979 }
01980
01981 T_flux_G = T_flux_S = T_flux_L = 0.0;
01982 CH_vol = Chan_list[chan_n]->area * CanWatDep;
01983
01984
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] )
01991
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.;
01999
02000
02001
02002
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;
02015
02016 dh = ( CH_bottElev + CanWatDep ) - tot_head;
02017
02018
02019 H_rad_ch = Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
02020 SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);
02021
02022 H_rad_cell = Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
02023 SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);
02024
02025
02026 if (dh > 0.0) {
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) {
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:
02040
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:
02045
02046 if ( dh > 0 ) {
02047
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
02053 if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02054 fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02055 }
02056
02057
02058 fluxG = (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0);
02059 break;
02060 case 1:
02061 if ( cell->ind < 0 )
02062 {
02063 if ( dh > 0 ) {
02064
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
02070 if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02071 fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02072 }
02073
02074 }
02075 else
02076
02077 fluxL = (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02078
02079
02080 fluxG = (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0);
02081 break;
02082 case -1:
02083 if ( cell->ind > 0 )
02084 {
02085 if ( dh > 0 ) {
02086
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
02092 if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02093 fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02094 }
02095 }
02096 else
02097
02098 fluxL = (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02099
02100
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 }
02108
02109
02110
02111
02112
02113
02114
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
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
02137 if ( fluxS > 0 || fluxG > 0 || fluxL > 0)
02138 {
02139 fluxO = (fluxS > 0) ? fluxS : 0;
02140 fluxO += (fluxG > 0) ? fluxG : 0;
02141 fluxO += (fluxL > 0) ? fluxL : 0;
02142
02143
02144 if ( (fluxCk = CH_vol - MIN_VOL - fluxS - fluxG - fluxL) < 0.0 )
02145 { fluxCk = (fluxO>0) ? ( ( fluxO + fluxCk )/fluxO) : (0.0);
02146 if (fluxS>0) fluxS *= fluxCk;
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
02156 if (debug>99 && fluxG != 0.0) do {
02157
02158
02159
02160
02161
02162
02163 fluxHoriz = fluxG;
02164
02165
02166
02167 UnsatZ_rec = (fluxG>0) ? ( ElevMap[i] - GWH[i] / poros[HAB[i]] ) : (0.0) ;
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);
02173
02174 UnsatPot_rec = (fluxG>0) ? (UnsatCap_rec - Unsat[i] ) : (0.0);
02175
02176
02177 horizTOsat_rec = (fluxG>0) ? (fluxHoriz) : (0.0);
02178 satTOsf_rec = (fluxG>0) ? (Max(fluxHoriz - UnsatPot_rec, 0.0) ) : (0.0);
02179 if (fluxG>0) unsatTOsat_rec = (UnsatZ_rec > 0.0) ?
02180
02181 ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[i]] ) / UnsatZ_rec * Unsat[i] ) :
02182 (0.0);
02183 else
02184 unsatTOsat_rec = 0.0;
02185
02186
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);
02194
02195
02196
02197
02198 CH_vol -= ( fluxS + fluxG + fluxL);
02199
02200 if ( converged )
02201 {
02202
02203
02204
02205 if ( fluxS != 0 ) {
02206 if ( fluxS > 0 )
02207 {
02208
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) ) {
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 {
02220
02221 if ( !Chan_list[chan_n]->parent ) {
02222
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
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
02237 {
02238
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) ) {
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 {
02250
02251 if ( !basn_list[basn[i]]->parent ) {
02252
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
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
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
02278 if ( fluxG != 0 ) {
02279 if ( fluxG > 0 )
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 ) {
02288 if ( !Chan_list[chan_n]->parent ) {
02289
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
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 {
02308
02309 if ( !Chan_list[chan_n]->parent ) {
02310
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
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 }
02324
02325
02326 else
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 ) {
02334 if ( !basn_list[basn[i]]->parent ) {
02335
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
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 {
02354
02355 if ( !basn_list[basn[i]]->parent ) {
02356
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
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
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
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
02389
02390
02391
02392
02393 }
02394
02395
02396 if ( fluxL != 0 ) {
02397 if ( fluxL > 0 )
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 ) {
02405 if ( !Chan_list[chan_n]->parent ) {
02406
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
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 {
02425
02426 if ( !Chan_list[chan_n]->parent ) {
02427
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
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 ) {
02441 NA[i] += N_fl;
02442 PA[i] += P_fl;
02443 SA[i] += S_fl;
02444 }
02445 else {
02446 GNA[i] += N_fl;
02447 GPA[i] += P_fl;
02448 GSA[i] += S_fl;
02449 }
02450 }
02451 else if (fluxL < 0)
02452 {
02453 if ( tot_head > GW_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 ) {
02472 if ( !basn_list[basn[i]]->parent ) {
02473
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
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 {
02492
02493 if ( !basn_list[basn[i]]->parent ) {
02494
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
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
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
02514
02515
02516 SWH[i] += (fluxS)/CELL_SIZE;
02517 GWH[i] += (fluxG)/CELL_SIZE;
02518 if ( tot_head > GW_head ) SWH[i] += fluxL/CELL_SIZE;
02519 else GWH[i] += fluxL/CELL_SIZE;
02520
02521
02522 }
02523
02524 T_flux_S += fluxS;
02525 T_flux_G += fluxG;
02526 T_flux_L += fluxL;
02527
02528 cell = cell->next_cell;
02529 }
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
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;
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
02548
02549
02550 if (debug>4) {
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 );
02562
02563 if ( factor != 0 && Qout > 0)
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
02569
02570
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 )
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;
02583 CH_vol = CanWatDep * Chan_list[chan_n]->area;
02584
02585 if( debug > 3 )
02586 {
02587
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
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
02706 extern float* boundcond_depth;
02707
02708
02709
02710 TWf = 0; HWf = 0;
02711
02712 structs = struct_first;
02713
02714 while ( structs != NULL )
02715 {
02716 if ( structs->flag < 0 || (structs->flag >1) ) {
02717 structs->flow = 0.0;
02718 structs->conc_P = 0.0;
02719 structs->conc_N = 0.0;
02720 structs->conc_S = 0.0;
02721 goto OUT;
02722 }
02723
02725 if ( structs->TW_stage == -99 )
02726 {
02727 structs->TW_stage = GetGraph ( structs->TW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 0.0) );
02728
02729 structs->TW_stage = structs->TW_stage + GP_SLRise * (SimTime.TIME/365);
02730 TWf = 1;
02731 }
02732
02733 if ( structs->HW_stage == -99 )
02734 {
02735 structs->HW_stage = GetGraph ( structs->HW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 365.0) );
02736
02737 structs->HW_stage = structs->HW_stage + GP_SLRise * (SimTime.TIME/365);
02738 HWf = 1; }
02739
02740
02741
02742
02743
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
02750 if ( structs->flag == 1 ) {
02751 if (structs->multiOut == 1) goto OUT;
02752
02753
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 }
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
02770
02771
02772
02773 struct_hist_start = structs;
02774 ChanHistOut = 0.0;
02775
02776 while ( structs != NULL )
02777 {
02778 if ( structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
02779 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
02780 ChanHistOut += structs->flow;
02781 }
02782
02783 structs = structs->next_in_list;
02784 }
02785
02786 structs = struct_hist_start;
02787
02788
02789
02790
02791
02792
02793
02794 CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth;
02795
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
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
02812 {
02813 if (structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
02814
02815 structs->multiOut = 1;
02816
02817
02818 if ( chan_over < 1.0 ) {
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
02826 structs->flow = (structs->flow) * chan_over ;
02827 }
02828
02829
02830 Chan_list[can_fr]->sumHistOut += structs->flow;
02831
02832
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;
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))] ) );
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
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
02851
02852
02853
02854
02855 if (can_to > 0) {
02856
02857
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
02864
02865
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);
02872 S_IN_STR[Chan_list [can_to]->basin] += S_fl;
02873 }
02874
02875
02876
02877 }
02878 else if (cell_to > 0) {
02879
02880
02881 if (structs->cell_i_to > 2 ) {
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
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);
02893 S_IN_STR[basn[cell_to]] += S_fl;
02894
02895
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 ) {
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);
02913 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0);
02914 }
02915 }
02916
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 }
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);
02932
02933 structs = struct_hist_start;
02934
02935 goto OUT;
02936 }
02937
02938 else {
02939
02940 HeadH_drop = 0.5 * Abs(Chan_list[can_fr]->elev_drop);
02941 HeadT_drop = 0.5 * Abs(Chan_list[can_to]->elev_drop);
02942
02943 HeadH = HeadH_drop +
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
02947
02948 HeadT = -HeadT_drop +
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
02953 if ( structs->w_coef == 0.0) {
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)
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);
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
02984 CH_vol = CH_vol + MINDEPTH*Chan_list[can_fr]->area;
02985
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
02998
02999
03000 }
03001
03002
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 }
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
03018
03019 Chan_list[can_fr]->sumRuleOut += structs->flow;
03020 Chan_list[can_to]->sumRuleIn += structs->flow;
03021
03022
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
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
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
03035
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
03047
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
03055 if ( structs->flag == 1 ) {
03056
03057
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;
03064 }
03065 else if (structs->flow < 0.0) {
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
03072 if (structs->cell_i_fr != 2 )
03073 {
03074
03075 cell_vol = SWH[cell_fr] * CELL_SIZE;
03076 if (structs->flow > 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
03086
03087
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 {
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
03112
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) {
03120 if (debug > -999) {
03121 VOL_OUT_STR[basn[cell_fr]] += structs->flow;
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);
03126 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0);
03127 }
03128 PA[cell_fr] -= P_fl;
03129 NA[cell_fr] -= N_fl;
03130 SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03131
03132
03133
03134
03135
03136 SWH[cell_fr] -= structs->flow/CELL_SIZE;
03137 }
03138 else if (structs->cell_i_fr == 2) {
03139 if (debug > -999) {
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
03151 SWH[cell_to] += structs->flow/CELL_SIZE;
03152 }
03153 else {
03154 PA[cell_fr] -= P_fl;
03155 NA[cell_fr] -= N_fl;
03156 SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03157
03158
03159
03160
03161
03162 PA[cell_to] += P_fl;
03163 NA[cell_to] += N_fl;
03164 SA[cell_to] += S_fl;
03165
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);
03174 S_IN_STR[basn[cell_to]] += S_fl;
03175
03176
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
03189 else {
03190 if ( structs->w_coef != 0.0) {
03191 printf ("\nVirtual structs are the only rule-based cell-cell flows in this version!\n"); exit(-1);
03192 }
03193 else {
03194
03195
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
03199
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
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 }
03212
03213 }
03214
03215
03216
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
03223 if ( structs->flag == 1 ) {
03224
03225 if (structs->multiOut == 1) goto OUT;
03226
03227 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03228
03229 if (structs->flow == 0.0) {
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
03243
03244
03245 struct_hist_start = structs;
03246 ChanHistOut = 0.0;
03247
03248 while ( structs != NULL )
03249 {
03250 if ( structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
03251 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03252 ChanHistOut += structs->flow;
03253 }
03254
03255 structs = structs->next_in_list;
03256 }
03257
03258 structs = struct_hist_start;
03259
03260
03261
03262
03263
03264 CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth;
03265
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
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
03281 {
03282 if ( structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
03283
03284 structs->multiOut = 1;
03285
03286
03287 if ( chan_over < 1.0 ) {
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
03295 structs->flow = (structs->flow) * chan_over ;
03296 }
03297
03298
03299 Chan_list[can_fr]->sumHistOut += structs->flow;
03300
03301
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;
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))] ) );
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
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
03320
03321
03322
03323
03324 if (can_to > 0) {
03325
03326
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
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);
03339 S_IN_STR[Chan_list [can_to]->basin] += S_fl;
03340 }
03341 }
03342 else if (cell_to > 0) {
03343
03344
03345 if (structs->cell_i_to > 2 ) {
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
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);
03357 S_IN_STR[basn[cell_to]] += S_fl;
03358
03359
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 ) {
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);
03373 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0);
03374 }
03375 }
03376
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 }
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);
03391
03392 structs = struct_hist_start;
03393
03394
03395 goto OUT;
03396
03397 }
03398
03399
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
03404
03405
03406
03407
03408
03409
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
03414 CH_vol = Chan_list[can_fr]->area* ramp(HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth-0.3);
03415
03416 if (HeadH>HeadT)
03417 {
03418
03419
03420 structs->flow = Min(structs->w_coef * ( HeadH - HeadT ) * canstep, CH_vol );
03421
03422
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
03427
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;
03434
03435
03436
03437
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 }
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
03452
03453 CH_vol = Chan_list[can_fr]->area* ramp(HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth);
03454
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
03464
03465 Chan_list[can_fr]->sumRuleOut += structs->flow;
03466
03467 if (structs->cell_i_to > 2 ) {
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
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
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 ) {
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
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
03511
03512
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
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;
03526 }
03527
03528 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03529
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 )
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
03556
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 {
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
03581
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 ) {
03586 if (debug > -999) {
03587 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow;
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 {
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
03600
03601
03602
03603
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);
03610 S_IN_STR[Chan_list[can_to]->basin] += S_fl;
03611
03612
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
03622
03623 Chan_list[can_to]->sumHistIn += structs->flow;
03624
03625
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 {
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
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;
03649
03650
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 }
03657
03658
03659
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
03677
03678
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
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
03706
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
03715 structs->SumFlow += structs->flow;
03716 structs->Sum_P += structs->conc_P * structs->flow;
03717
03718 structs->Sum_S += structs->conc_S * structs->flow;
03719
03720 structs->multiOut = 0;
03721
03722
03723 if ( printchan ) {
03724 fprintf( WstructOutFile, "%7.0f\t", (structs->SumFlow)/(2446.848) );
03725 fprintf( WstructOutFile_P, "%7.4f\t",
03726 (structs->SumFlow > 0.0) ? (structs->Sum_P/structs->SumFlow*1000.0) : (0.0) );
03727
03728
03729 fprintf( WstructOutFile_S, "%7.4f\t",
03730 (structs->SumFlow > 0.0) ? (structs->Sum_S/structs->SumFlow) : (0.0) );
03731
03732 structs->SumFlow = 0.0;
03733 structs->Sum_P = 0.0;
03734
03735 structs->Sum_S = 0.0;
03736 }
03737
03738
03739 structs = structs->next_in_list;
03740 }
03741
03742 if (printchan) {
03743 fflush( WstructOutFile );
03744 fflush( WstructOutFile_P );
03745
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