00001
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083 #include "watmgmt.h"
00084
00085
00086
00096 void Canal_Network_Init(float baseDatum, float *elev )
00097
00098 { struct Structure *structs;
00099 int i,j;
00100 char filename[30];
00101
00102 sprintf(filename,"CanalData");
00103
00104
00105
00106
00107 Channel_configure(elev, ReadChanStruct (filename));
00108
00109 for( i = 0; i < num_chan; i++ )
00110 {
00111 Chan_list[i]->P = Chan_list[i]->ic_P_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00112 Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00113 Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00114 }
00115
00116 usrErr ("\tCanal geometry configured OK");
00117
00118
00119 ReadStructures (filename, baseDatum);
00120
00121 ON_MAP[T(2,2)] = 0;
00122
00123 MCopen = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MCopen");
00124 init_pvar(MCopen,NULL,'f',0.05);
00125
00126
00127 if ( debug > 0 )
00128 {
00129 sprintf( modelFileName, "%s/%s/Output/Debug/ON_MAP_CANAL.txt", OutputPath, ProjName );
00130
00131
00132 if ( ( ChanInFile = fopen( modelFileName, "w" ) ) == NULL )
00133 {
00134 printf( "Can't open the %s canal/levee debug file!\n", modelFileName );
00135 exit(-1) ;
00136 }
00137 fprintf ( ChanInFile, "ROWS=%d\nCOLUMNS=%d\nFORMAT=text\nINFO=\"Debug data: Overland flow restrictions due to levees. \"\nMISSING=0\n", s0, s1 );
00138
00139 for ( i = 1; i <= s0 ; i++ )
00140 { for ( j = 1 ; j <= s1 ; j++ )
00141 fprintf( ChanInFile, "%4d\t", ON_MAP[T(i,j)]) ;
00142 fprintf ( ChanInFile, "\n" );
00143 }
00144 fclose( ChanInFile ) ;
00145 }
00146
00147
00148
00149 { if(H_OPSYS==UNIX)
00150 sprintf( modelFileName, "%s/%s/Output/Canal/structsOut", OutputPath, ProjName );
00151 else sprintf( modelFileName, "%s%s:Output:structsOut", OutputPath, ProjName );
00152 }
00153
00154
00155 if ( ( WstructOutFile = fopen( modelFileName, "w" ) ) == NULL )
00156 {
00157 printf( "Can't open %s file!\n", modelFileName );
00158 exit(-1) ;
00159 }
00160
00161 fprintf ( WstructOutFile,
00162 "%s %s %s %s scenario: Sum of water flows at structures (ft^3/sec)\n Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00163 structs = struct_first;
00164 while ( structs != NULL )
00165 { fprintf( WstructOutFile, "%7s\t", structs->S_nam);
00166 structs = structs->next_in_list;
00167 }
00168
00169
00170 { if(H_OPSYS==UNIX)
00171 sprintf( modelFileName, "%s/%s/Output/Canal/structsOut_P", OutputPath, ProjName );
00172 else sprintf( modelFileName, "%s%s:Output:structsOut_P", OutputPath, ProjName );
00173 }
00174
00175
00176 if ( ( WstructOutFile_P = fopen( modelFileName, "w" ) ) == NULL )
00177 {
00178 printf( "Can't open %s file!\n", modelFileName );
00179 exit(-1) ;
00180 }
00181
00182 fprintf ( WstructOutFile_P,
00183 "%s %s %s %s scenario: Flow weighted mean TP conc at structures ( mg/L) \n Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00184 structs = struct_first;
00185 while ( structs != NULL )
00186 { fprintf( WstructOutFile_P, "%7s\t", structs->S_nam);
00187 structs = structs->next_in_list;
00188 }
00189
00190
00191 { if(H_OPSYS==UNIX)
00192 sprintf( modelFileName, "%s/%s/Output/Canal/structsOut_S", OutputPath, ProjName );
00193 else sprintf( modelFileName, "%s%s:Output:structsOut_S", OutputPath, ProjName );
00194 }
00195
00196
00197 if ( ( WstructOutFile_S = fopen( modelFileName, "w" ) ) == NULL )
00198 {
00199 printf( "Can't open %s file!\n", modelFileName );
00200 exit(-1) ;
00201 }
00202
00203 fprintf ( WstructOutFile_S,
00204 "%s %s %s %s scenario: Flow weighted mean tracer concentration at structures ( g/L) \n Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00205 structs = struct_first;
00206 while ( structs != NULL )
00207 { fprintf( WstructOutFile_S, "%7s\t", structs->S_nam);
00208 structs = structs->next_in_list;
00209 }
00210
00211
00212
00213
00214 { if(H_OPSYS==UNIX)
00215 sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut", OutputPath, ProjName );
00216 else sprintf( modelFileName, "%s%s:Output:CanalOut", OutputPath, ProjName );
00217 }
00218
00219 if ( ( CanalOutFile = fopen( modelFileName, "w" ) ) == NULL )
00220 {
00221 printf( "Can't open %s file!\n", modelFileName );
00222 exit(-1) ;
00223 }
00224
00225 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 );
00226
00227
00228 for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile, "R_%d\t",Chan_list[i]->number);}
00229
00230
00231
00232 { if(H_OPSYS==UNIX)
00233 sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut_P", OutputPath, ProjName );
00234 else sprintf( modelFileName, "%s%s:Output:CanalOut_P", OutputPath, ProjName );
00235 }
00236
00237 if ( ( CanalOutFile_P = fopen( modelFileName, "w" ) ) == NULL )
00238 {
00239 printf( "Can't open %s file!\n", modelFileName );
00240 exit(-1) ;
00241 }
00242
00243 fprintf ( CanalOutFile_P, "%s %s %s %s scenario: Instantaneous TP conc (mg/L) in canal Reaches\n Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00244
00245 for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile_P, "R_%d\t",Chan_list[i]->number);}
00246
00247
00248 { if(H_OPSYS==UNIX)
00249 sprintf( modelFileName, "%s/%s/Output/Canal/CanalOut_S", OutputPath, ProjName );
00250 else sprintf( modelFileName, "%s%s:Output:CanalOut_S", OutputPath, ProjName );
00251 }
00252
00253 if ( ( CanalOutFile_S = fopen( modelFileName, "w" ) ) == NULL )
00254 {
00255 printf( "Can't open %s file!\n", modelFileName );
00256 exit(-1) ;
00257 }
00258
00259 fprintf ( CanalOutFile_S, "%s %s %s %s scenario: Instantaneous tracer conc (g/L) in canal Reaches\n Date\t", &modelName, &modelVers, &SimAlt, &SimModif );
00260
00261 for( i = 0; i < num_chan; i++ ) {if (Chan_list[i]->width > 0) fprintf ( CanalOutFile_S, "R_%d\t",Chan_list[i]->number);}
00262
00263
00264
00265
00266 if( debug > 0 )
00267 {
00268 if(H_OPSYS==UNIX)
00269 sprintf( modelFileName, "%s/%s/Output/Debug/Can_debug", OutputPath, ProjName );
00270 else sprintf( modelFileName, "%s%s:Output:Can_debug", OutputPath, ProjName );
00271
00272
00273
00274 if ( ( canDebugFile = fopen( modelFileName, "w" ) ) == NULL )
00275 {
00276 printf( "Can't open %s file!\n", modelFileName );
00277 exit(-1) ;
00278 }
00279
00280 fprintf ( canDebugFile, "Depending on level of the debug requested, printed here are data on canal-cell flows and excessive demands on data-driven structure flows. \n" );
00281 fflush(canDebugFile);
00282
00283 }
00284
00285 return;
00286
00287 }
00288
00289
00292 void CanalReInit()
00293 {
00294 int i;
00295 for( i = 0; i < num_chan; i++ )
00296
00297 {
00298 Chan_list[i]->wat_depth = Chan_list[i]->ic_depth;
00299 Chan_list[i]->P = Chan_list[i]->ic_P_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00300 Chan_list[i]->N = Chan_list[i]->ic_N_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00301 Chan_list[i]->S = Chan_list[i]->ic_S_con * Chan_list[i]->area * Chan_list[i]->ic_depth;
00302
00303 }
00304 }
00305
00306
00307
00308
00309
00330 void Run_Canal_Network ( float *SWH, float *Elevation, float *MC, float *GW, float *poros,
00331 float *GWcond, double *NA, double *PA, double *SA, double *GNA, double *GPA, double *GSA,
00332 float *Unsat, float *sp_yield)
00333
00334 { int i;
00335 float CH_vol;
00336
00337 if ( canPrint && (( N_c_iter++ % hyd_iter ) == 0.0) ) {
00338 if (SensiOn) {
00339 printchan = (SimTime.IsSimulationEnd) ? (1) : (0);
00340 }
00341 else printchan = 1;
00342 }
00343 else printchan = 0;
00344
00345 if (printchan == 1) {
00346
00347 fprintf( WstructOutFile, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00348 fprintf( WstructOutFile_P, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00349
00350 fprintf( WstructOutFile_S, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00351
00352 fprintf (CanalOutFile, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00353 fprintf (CanalOutFile_P, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00354
00355 fprintf (CanalOutFile_S, "\n%d/%d/%d\t",SimTime.yr[0],SimTime.mo[0],SimTime.da[0] );
00356 }
00357
00358
00359 for ( i = 0; i < num_chan; i++ )
00360 {
00361 Chan_list[i]->sumHistIn = Chan_list[i]->sumHistOut = 0.0;
00362 Chan_list[i]->sumRuleIn = Chan_list[i]->sumRuleOut = 0.0;
00363 }
00364
00365
00366 Flows_in_Structures ( SWH, GW, poros, Elevation, NA, PA, SA );
00367
00368
00369 if( debug > 3 )
00370 {
00371 fprintf( canDebugFile, "N%3d WatDepth Pconc Sconc flux_S flux_G flux_L Qin Qout #_iter Area error \n", N_c_iter );
00372
00373 }
00374
00375
00376
00377 for( i = 0; i < num_chan; i++ )
00378 if (Chan_list[i]->width > 0)
00379
00380 {
00381 FluxChannel ( i, SWH, Elevation, MC, GW, poros, GWcond, NA, PA, SA, GNA, GPA, GSA, Unsat, sp_yield );
00382
00383 CH_vol = Chan_list[i]->area*Chan_list[i]->wat_depth;
00384
00385
00386 if ( !( FMOD(N_c_iter,(1/canstep) )) && (SimTime.IsBudgEnd ) ) {
00387 TOT_VOL_CAN[Chan_list[i]->basin] += CH_vol;
00388 TOT_VOL_CAN[0] += CH_vol;
00389 TOT_P_CAN[Chan_list[i]->basin] += Chan_list[i]->P;
00390 TOT_P_CAN[0] += Chan_list[i]->P;
00391 TOT_S_CAN[Chan_list[i]->basin] += Chan_list[i]->S;
00392 TOT_S_CAN[0] += Chan_list[i]->S;
00393 }
00394
00395
00396 if ( printchan )
00397 {
00398
00399 Chan_list[i]->P_con = (CH_vol > 0) ? (Chan_list[i]->P/CH_vol*1000.0) : 0.0;
00400
00401 Chan_list[i]->S_con = (CH_vol > 0) ? (Chan_list[i]->S/CH_vol ) : 0.0;
00402
00403
00404 fprintf( CanalOutFile, "%6.2f\t", Chan_list[i]->wat_depth );
00405 fprintf( CanalOutFile_P, "%7.4f\t", Chan_list[i]->P_con );
00406
00407 fprintf( CanalOutFile_S, "%7.4f\t", Chan_list[i]->S_con );
00408 }
00409 }
00410 if ( printchan )
00411 {
00412 fflush(CanalOutFile);
00413 fflush(CanalOutFile_P);
00414
00415 fflush(CanalOutFile_S);
00416 }
00417
00418 return;
00419 }
00420
00421
00422
00423
00435 struct Chan *ReadChanStruct ( char* filename )
00436 {
00437 struct Chan *channels, *chan_first, *chan_last;
00438 struct Chan_segm *C_segm, *last_segm, *next_segm;
00439 int i, index, firstSegm, firstCanal = 1;
00440 char ss[622], scenario[20], modnam[20], bas_txt[8];
00441 float C_x, C_y, C_x1, C_y1;
00442
00443
00444 int ibas, foundBasn;
00445
00446
00447 { if(H_OPSYS==UNIX)
00448 sprintf( modelFileName, "%s/%s/Data/%s.chan", ModelPath, ProjName, filename );
00449 else sprintf( modelFileName, "%s%s:Data:%s.chan", ModelPath, ProjName, filename );
00450 }
00451
00452
00453 if ( ( ChanInFile = fopen( modelFileName, "r" ) ) == NULL )
00454 {
00455 sprintf( msgStr,"Can't open %s file! ",modelFileName ) ; usrErr(msgStr);
00456 exit(-1) ;
00457 }
00458
00459
00460
00461
00462 if ( (C_segm = (struct Chan_segm *) malloc( (size_t) sizeof( struct Chan_segm ))) == NULL )
00463 {
00464 usrErr( "No memory for first canal reach segment\n " ) ;
00465 exit( -2 ) ;
00466 };
00467
00468 num_chan = 0;
00469
00470
00471 fgets( ss, 620, ChanInFile ); sscanf( ss, "%d", &UTM_EOrigin );
00472 fgets( ss, 620, ChanInFile ); sscanf( ss, "%d", &UTM_NOrigin );
00473 fgets( ss, 620, ChanInFile );
00474 fgets( ss, 620, ChanInFile ); sscanf( ss,"%s %s %d %d", &modnam, &scenario, &CHo, &CHANNEL_MAX_ITER );
00475 if (strcmp(scenario,SimAlt) != 0) {
00476 sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
00477 scenario, modelFileName, SimAlt); usrErr(msgStr);
00478 exit(-1);
00479 }
00480 if (strcmp(modnam,modelName) != 0) {
00481 sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
00482 modnam, modelFileName, modelName); usrErr(msgStr);
00483 exit(-1);
00484 }
00485
00486 fgets( ss, 620, ChanInFile );
00487 fgets( ss, 620, ChanInFile ); sscanf ( ss, "%f %f %f ", &F_ERROR, &MINDEPTH, &C_F );
00488
00489 fgets( ss, 620, ChanInFile );
00490
00491
00492 while ( fgets( ss, 620, ChanInFile ) != NULL && !feof( ChanInFile ) )
00493 {
00494 if ( *ss == '#' )
00495 {
00496 if ( (channels = (struct Chan *) malloc( (size_t) sizeof( struct Chan ))) == NULL )
00497 {
00498 printf( "No memory for first canal\n " ) ;
00499 exit( -2 ) ;
00500 };
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 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",
00512 &channels->number,
00513 &channels->levee,
00514 &channels->roil, &channels->roir,
00515 &channels->depth, &channels->width,
00516 &channels->cond,
00517 &channels->ic_S_con, &channels->ic_P_con,
00518 &channels->wat_depth, &channels->edgeMann, &bas_txt);
00519
00520 foundBasn=0;
00521 for (ibas=numBasn;ibas>0;ibas--) {
00522 basins = basn_list[ibas];
00523 if(strcmp(bas_txt,basins->basnTxt)==0) {
00524 channels->basin = ibas;
00525 foundBasn=1;
00526 channels->family = basn_list[ibas]->family;
00527 channels->parent = basn_list[ibas]->parent;
00528
00529 if (!channels->parent) {
00530 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",
00531 modelName,channels->number); exit (-1);}
00532 }
00533 }
00534 if (foundBasn==0) { printf ("Sorry - %s's canal %d's basin does not exist! Please define the basin in the basinIR file.\n",
00535 modelName,channels->number); exit (-1);}
00536
00537 channels->ic_P_con *= 0.001;
00538 channels->ic_N_con = 0.00;
00539 channels->ic_S_con *= 1.0;
00540
00541 channels->ic_depth = channels->wat_depth;
00542
00543
00544 channels->segments = C_segm;
00545
00546 if (i < 10)
00547 {printf ( " Error in CanalData.chan: reach input %d\n", (num_chan+1) ); exit (-1);}
00548
00549
00550 num_chan++;
00551 firstSegm = 1;
00552
00553 if (firstCanal)
00554 {
00555 firstCanal = 0;
00556 chan_first = channels;
00557 chan_last = channels;
00558 continue;
00559 }
00560
00561 last_segm->next_segm = NULL;
00562 chan_last->next_in_list = channels;
00563 chan_last = channels;
00564 }
00565 else
00566 {
00567 i = sscanf ( ss, "%f %f", &C_x1, &C_y1 );
00568 if (i < 2)
00569 {printf ( " Error in CanalData.chan coords: Reach %d \n", (num_chan+1) ); exit (-1);}
00570
00571
00572 C_x = UTM2kmx (C_x1); C_y = UTM2kmy (C_y1);
00573 C_segm->x1 = C_x; C_segm->y1 = C_y;
00574
00575 if ( firstSegm )
00576 { C_segm->x0 = C_x; C_segm->y0 = C_y;
00577 firstSegm = 0;
00578 }
00579 else
00580 {
00581
00582 if ( (next_segm = (struct Chan_segm *) malloc( (size_t) sizeof( struct Chan_segm ))) == NULL )
00583 { printf( "No memory for canal reach segment\n " ) ;
00584 exit( -2 ) ;
00585 };
00586 C_segm->next_segm = next_segm;
00587 last_segm = C_segm;
00588 C_segm = next_segm;
00589 C_segm->x0 = C_x; C_segm->y0 = C_y;
00590 }
00591 }
00592 }
00593
00594 chan_last->next_in_list = NULL;
00595 last_segm->next_segm = NULL;
00596 free ((char *) C_segm);
00597 fclose (ChanInFile);
00598
00599 usrErr( "\tCanal locs/attributes read OK" );
00600
00601 return (chan_first);
00602 }
00603
00604
00611 void ReadStructures ( char* filename, float BASE_DATUM )
00612 {
00622 struct Structure *structs, *struct_next, *struct_last;
00623
00624 char ss[622], *s;
00625 char sched_name[20];
00626 int i, j, k;
00627 int IDhist[MAX_H_STRUCT];
00628 int lincnt;
00629 double Jdate_read;
00630 char lline[2001], *line;
00631 int yyyy, mm, dd;
00632 char structName[20];
00633 char scenario[20], s_modname[20];
00634 char histTP[5];
00635 char histTS[5];
00636 char srcType[5], destType[5];
00637 char prefixTide[]="tid_";
00638
00639 num_struct_hist = 0;
00640 numTPser = numTSser = disAggNum = 0;
00641
00642 { if(H_OPSYS==UNIX)
00643 sprintf( modelFileName, "%s/%s/Data/%s.struct", ModelPath, ProjName, filename );
00644 else sprintf( modelFileName, "%s%s:Data:%s.struct", ModelPath, ProjName, filename );
00645 }
00646
00647
00648 if ( ( ChanInFile = fopen( modelFileName, "r" ) ) == NULL )
00649 {
00650 printf( "Can't open %s file!\n", modelFileName ) ;
00651 exit(-1) ;
00652 }
00653
00654
00655 if ( (structs = (struct Structure *) malloc( (size_t) sizeof( struct Structure ))) == NULL )
00656 {
00657 printf( "No memory for first structure\n " ) ;
00658 exit( -2 ) ;
00659 }
00660 struct_first = structs;
00661
00662 fgets( ss, 620, ChanInFile );
00663 sscanf(ss, "%s\t%s\t%s\t%s",&scenario,&scenario,&s_modname,&scenario);
00664 if (strcmp(scenario,SimAlt) != 0) {
00665 sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
00666 scenario, modelFileName, SimAlt); usrErr(msgStr);
00667 exit(-1);
00668 }
00669 if (strcmp(s_modname,modelName) != 0) {
00670 sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) in the canal data file!",
00671 s_modname, modelFileName, modelName); usrErr(msgStr);
00672 exit(-1);
00673 }
00674
00675
00676 fgets( ss, 620, ChanInFile );
00677
00678
00679
00680 while ( fgets( ss, 620, ChanInFile ) != NULL && !feof( ChanInFile ) )
00681 {
00682 s = ss;
00683
00684
00685 sscanf( s, "%d", &structs->flag );
00686
00687
00688
00689 if ( structs->flag < 0 ) continue;
00690 if ( structs->flag > 0 && structs->flag <10 ) {
00691 structs->histID = ++num_struct_hist;
00692 Hist_data[structs->histID-1].flag = structs->flag;
00693 }
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707 s = Scip( s, TAB );
00708
00709 if ( *s != TAB ) {
00710 sscanf ( s, "%s", &structs->S_nam );
00711 if ( structs->flag >0 && structs->flag <10)
00712 sscanf ( s, "%s", &Hist_data[structs->histID-1].S_nam );
00713
00714
00715
00716 s = Scip( s, TAB );
00717 if ( *s == TAB ) structs->TPser = -1;
00718 else {
00719 sscanf ( s, "%s", &histTP );
00720
00721 if (isInteger(histTP) ) {
00722 structs->TPser = 0;
00723 structs->TP = atof(histTP)*1.0e-6;
00724 }
00725 else { structs->TPser = 1; numTPser++;}
00726 }
00727
00728 s = Scip( s, TAB );
00729 if ( *s == TAB ) structs->TNser = -1;
00730 else {
00731 sscanf ( s, "%f", &structs->TN );
00732 structs->TNser = 0;
00733 structs->TN *= 1.0e-6;
00734 }
00735
00736
00737 s = Scip( s, TAB );
00738 if ( *s == TAB ) structs->TSser = -1;
00739
00740 else {
00741 sscanf ( s, "%s", &histTS );
00742
00743 if (isFloat(histTS) ) {
00744 structs->TSser = 0;
00745 structs->TS = atof(histTS)*1.0;
00746 }
00747 else { structs->TSser = 1; numTSser++;}
00748
00749 }
00750
00751 s = Scip( s, TAB );
00752 if ( *s == TAB ) structs->str_cell_i = 0;
00753 else {
00754 sscanf( s, "%d", &structs->str_cell_i );
00755 structs->str_cell_i++;
00756 }
00757
00758 s = Scip( s, TAB );
00759
00760 if ( *s == TAB ) structs->str_cell_j = 0;
00761 else {
00762 sscanf( s, "%d", &structs->str_cell_j );
00763 structs->str_cell_j++;
00764 }
00765
00766 if ( structs->str_cell_i != 0 && !ON_MAP[T(structs->str_cell_i, structs->str_cell_j)] )
00767 { printf ( "Location of water control structure %s is off-map.\n", structs->S_nam );
00768 exit (-2);
00769 }
00770
00771
00772 }
00773
00774 s = Scip( s, TAB );
00775
00776 if ( *s == TAB ) structs->canal_fr = 0;
00777 else sscanf( s, "%d", &structs->canal_fr );
00778 s = Scip( s, TAB );
00779 if ( *s == TAB) structs->canal_to = 0;
00780 else sscanf( s, "%d", &structs->canal_to );
00781 s = Scip( s, TAB );
00782
00783 if ( structs->canal_fr - CHo >= num_chan || structs->canal_to - CHo >= num_chan )
00784 { printf ( "Canal from/to doesn't exist for water control structure %s\n", structs->S_nam );
00785 exit (-2);
00786 }
00787
00788 if ( structs->canal_fr != 0 && Chan_list[structs->canal_fr- CHo]->width <= 0.1 )
00789 { printf ( "The donor canal %d has been turned off for water control structure %s\n", structs->canal_fr, structs->S_nam );
00790 exit (-2);
00791 }
00792 if ( structs->canal_to != 0 && Chan_list[structs->canal_to- CHo]->width <= 0.1 )
00793 { printf ( "The recipient canal %d has been turned off for water control structure %s\n", structs->canal_to, structs->S_nam );
00794 exit (-2);
00795 }
00796
00797 if (*s == TAB) structs->cell_i_fr = 0;
00798 else
00799 { sscanf( s, "%d", &structs->cell_i_fr );
00800 structs->cell_i_fr++;
00801 }
00802 s = Scip( s, TAB );
00803 if (*s == TAB) structs->cell_j_fr = 0;
00804 else
00805 { sscanf( s, "%d", &structs->cell_j_fr );
00806 structs->cell_j_fr++;
00807 }
00808 s = Scip( s, TAB );
00809
00810
00811 if ( structs->cell_i_fr > 2 && !ON_MAP[T(structs->cell_i_fr, structs->cell_j_fr)] )
00812 { printf ( "Donor cell for water control structure %s is off-map.\n", structs->S_nam );
00813 exit (-2);
00814 }
00815
00816
00817 if ( structs->cell_j_fr ) ON_MAP[T(structs->cell_i_fr, structs->cell_j_fr)] += 100;
00818
00819 if (*s == TAB) structs->cell_i_to = 0;
00820 else
00821 { sscanf( s, "%d", &structs->cell_i_to );
00822 structs->cell_i_to++;
00823 }
00824
00825 s = Scip( s, TAB );
00826 if (*s == TAB) structs->cell_j_to = 0;
00827 else
00828 { sscanf( s, "%d", &structs->cell_j_to );
00829 structs->cell_j_to++;
00830 }
00831
00832
00833 if ( structs->cell_i_to > 2 && !ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] )
00834 { printf ( "Recipient cell for water control structure %s is off-map.\n", structs->S_nam );
00835 exit (-2);
00836 }
00837
00838
00839
00840 if ( structs->cell_j_to ) ON_MAP[T(structs->cell_i_to, structs->cell_j_to)] += 100;
00841
00842
00843
00844
00845 s = Scip( s, TAB );
00846 if (*s == TAB) {
00847 structs->flagHWsched = -1;
00848 }
00849 else if ( isalpha (*s) ) {
00850 sscanf( s, "%s", sched_name );
00851 structs->HW_graph = Read_schedule ( sched_name, filename, BASE_DATUM );
00852 structs->flagHWsched = 1;
00853 structs->flagHWtide = ( strncmp(sched_name, prefixTide ,4) == 0) ? (1) : (0);
00854 }
00855 else { sscanf( s, "%f", &structs->HW_stage );
00856 structs->HW_stage += BASE_DATUM;
00857 structs->flagHWsched = 0;
00858 }
00859
00860
00861
00862
00863 s = Scip( s, TAB );
00864 if (*s == TAB) {
00865 structs->flagTWsched = -1;
00866 }
00867 else if ( isalpha (*s) ) {
00868 sscanf( s, "%s", sched_name );
00869 structs->TW_graph = Read_schedule ( sched_name, filename, BASE_DATUM );
00870 structs->flagTWsched = 1;
00871 structs->flagTWtide = ( strncmp(sched_name, prefixTide ,4) == 0) ? (1) : (0);
00872 }
00873 else { sscanf( s, "%f", &structs->TW_stage );
00874 structs->TW_stage += BASE_DATUM;
00875 structs->flagTWsched = 0;
00876 }
00877
00878
00879 s = Scip( s, TAB );
00880 if (*s == TAB) structs->cell_i_HW = 0;
00881 else
00882 { sscanf( s, "%d", &structs->cell_i_HW );
00883 structs->cell_i_HW++;
00884 }
00885
00886
00887 s = Scip( s, TAB );
00888 if (*s == TAB) structs->cell_j_HW = 0;
00889 else
00890 { sscanf( s, "%d", &structs->cell_j_HW );
00891 structs->cell_j_HW++;
00892 }
00893
00894 if ( structs->cell_i_HW > s0 || structs->cell_j_HW > s1 )
00895 { printf ( "Error in definition of HW location of water control structure %s", structs->S_nam );
00896 exit (-2);
00897 }
00898 else
00899 { (structs->cell_i_HW != 0) ? (structs->cell_HW = T(structs->cell_i_HW, structs->cell_j_HW) ) : (0);
00900 }
00901
00902
00903
00904
00905
00906
00907
00908
00909 s = Scip( s, TAB );
00910 if (*s == TAB) structs->cell_i_TW = 0;
00911 else
00912 { sscanf( s, "%d", &structs->cell_i_TW );
00913 structs->cell_i_TW++;
00914 }
00915
00916
00917 s = Scip( s, TAB );
00918 if (*s == TAB) structs->cell_j_TW = 0;
00919 else
00920 { sscanf( s, "%d", &structs->cell_j_TW );
00921 structs->cell_j_TW++;
00922 }
00923
00924 if ( structs->cell_i_TW > s0 || structs->cell_j_TW > s1 )
00925 { printf ( "Error in definition of TW location of water control structure %s", structs->S_nam );
00926 exit (-2);
00927 }
00928 else
00929 { (structs->cell_i_TW != 0) ? (structs->cell_TW = T(structs->cell_i_TW, structs->cell_j_TW) ) : (0);
00930 }
00931
00932
00933
00934
00935
00936
00937
00938 s = Scip( s, TAB );
00939 if (*s == TAB || *s == '?') structs->w_coef = 0;
00940 else sscanf( s, "%f", &structs->w_coef );
00941
00942 struct_last = structs;
00943
00944
00945 if ( (struct_next = (struct Structure *) malloc( (size_t) sizeof( struct Structure ))) == NULL )
00946 {
00947 printf( "Failed to allocate memory for next structure (%s)\n ", structs->S_nam ) ;
00948 exit( -2 ) ;
00949 };
00950
00951
00952
00953 if ( structs->canal_fr != 0 ) sprintf( srcType, "Can");
00954 if ( structs->canal_to != 0 ) sprintf( destType, "Can");
00955 if ( structs->cell_i_fr != 0 && structs->cell_j_fr != 0 ) sprintf( srcType, "Cel");
00956 if ( structs->cell_i_to != 0 && structs->cell_j_to != 0 ) sprintf( destType, "Cel");
00957 sprintf( structs->SrcDest, "%s%s", srcType, destType );
00958
00959
00960
00961 structs->SumFlow = 0.0;
00962 structs->Sum_P = 0.0;
00963 structs->Sum_N = 0.0;
00964 structs->Sum_S = 0.0;
00965
00966 structs->next_in_list = struct_next;
00967 structs = struct_next;
00968 }
00969
00970 struct_last->next_in_list = NULL;
00971 free ((char *)structs);
00972 fclose (ChanInFile);
00973
00974
00975 sprintf(msgStr,"\tUsing %d target stage schedules for rule-based flows", NumScheds); usrErr(msgStr);
00976 usrErr( "\tWater control structure locs/attributes read OK" );
00977
00978
00979 if (num_struct_hist>0) ReadStructFlow_PtSer(filename);
00980 if (numTPser>0) ReadStructTP_PtSer(filename);
00981 if (numTSser>0) ReadStructTS_PtSer(filename);
00982
00983 sprintf(msgStr, "\t**** Evaluate scenario: %s %s ****", &SimAlt, &SimModif ); usrErr(msgStr); WriteMsg(msgStr,1);
00984
00985 return;
00986 }
00987
00988
00989
00994 void ReadStructFlow_PtSer( char* filename)
00995 {
00996 char lline[2001], *line;
00997 char scenario[20];
00998 int i, k;
00999 int lincnt;
01000 char structName[20];
01001 int yyyy, mm, dd;
01002 double Jdate_read;
01003 usrErr0 ( "\tControl structures' water flow time series...");
01004 if(H_OPSYS==UNIX)
01005 sprintf( modelFileName, "%s/%s/Data/%s.struct_wat", ModelPath, ProjName, filename );
01006 else sprintf( modelFileName, "%s%s:Data:%s.struct_wat", ModelPath, ProjName, filename );
01007
01008
01009
01010 if ( ( F_struct_wat = fopen( modelFileName, "r" ) ) == NULL )
01011 {
01012 printf( "Can't open %s file!\n", modelFileName ) ;
01013 exit(-1) ;
01014 }
01015
01016 fgets(lline, 2000, F_struct_wat);
01017 line=lline;
01018 sscanf(line, "%s",&scenario);
01019 if (strcmp(scenario,SimAlt) != 0) {
01020 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);
01021 exit(-1);
01022 }
01023
01024 fgets(lline, 2000, F_struct_wat);
01025 line=lline;
01026 line = Scip( line, TAB );
01027 line = Scip( line, TAB );
01028 line = Scip( line, TAB );
01029
01030
01031
01032 for ( i = 0; i < num_struct_hist; i++ ) {
01033
01034 Hist_data[i].arrayWat = (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01035 Hist_data[i].arrayP = (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01036 Hist_data[i].arrayN = (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01037 Hist_data[i].arrayS = (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01038
01039
01040 if ( Hist_data[i].flag > 1 ) {
01041 structs = struct_first;
01042 while ( structs != NULL )
01043 {
01044
01045 if ( structs->flag == 10 * Hist_data[i].flag )
01046 {
01047
01048 disAggNum++;
01049 Hist_data[num_struct_hist+disAggNum-1].arrayWat =
01050 (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01051 Hist_data[num_struct_hist+disAggNum-1].arrayP =
01052 (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01053 Hist_data[num_struct_hist+disAggNum-1].arrayN =
01054 (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01055 Hist_data[num_struct_hist+disAggNum-1].arrayS =
01056 (float*) nalloc( ( (PORnumday+2)*sizeof(float) ), "timeseries temp" );
01057
01058 Hist_data[num_struct_hist+disAggNum-1].aggID = i;
01059 strcpy( Hist_data[num_struct_hist+disAggNum-1].S_nam, structs->S_nam );
01060 Hist_data[i].aggCnt++;
01061 structs->flag = 1;
01062 structs->histID = num_struct_hist+disAggNum;
01063
01064 }
01065
01066 structs = structs->next_in_list;
01067 }
01068 }
01069
01070 sscanf ( line,"%s\t",&structName);
01071 if (strcmp(Hist_data[i].S_nam , structName) != 0) {
01072 printf( "variable name %s in input file does not match order of CanalData.structs.\n", structName);
01073 exit(-1);
01074 }
01075 line = Scip( line, TAB );
01076 }
01077
01078
01079 lincnt = 0;
01080 while ( fgets(lline, 2000, F_struct_wat) != NULL && !feof( F_struct_wat ) )
01081 {
01082 line = lline;
01083 sscanf( line, "%d %d %d",&yyyy,&mm,&dd);
01084
01085
01086 Jdate_read = julday( mm,dd,yyyy,0,0,0.0 );
01087 if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
01088 printf (" \n***Error\nStructure flow data (%s) init date > simulation start date\n", modelFileName); exit (-1);
01089 }
01090 if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
01091 {
01092 line = Scip( line, TAB );
01093 line = Scip( line, TAB );
01094
01095 for ( i = 0; i < num_struct_hist; i++ ) {
01096 line = Scip( line, TAB );
01097 if (*line == TAB) Hist_data[i].arrayWat[((int)(lincnt))] = 0;
01098 else sscanf( line, "%f", &Hist_data[i].arrayWat[((int)(lincnt))] );
01099 Hist_data[i].arrayWat[((int)(lincnt))] *= 2446.848;
01100 }
01101
01102 for (i = num_struct_hist; i < num_struct_hist+disAggNum; i++) {
01103 k = Hist_data[i].aggID;
01104 Hist_data[i].arrayWat[((int)(lincnt))] =
01105 Hist_data[k].arrayWat[((int)(lincnt))] / (Hist_data[k].aggCnt) ;
01106
01107 }
01108
01109
01110 lincnt++;
01111 }
01112
01113 }
01114
01115 if (lincnt == 0) {
01116 printf (" \n***Error\nStructure flow data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
01117 }
01118 else if (lincnt != PORnumday ) {
01119 printf (" \n***Error\nStructure flow data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
01120 }
01121
01122 sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
01123 fclose (F_struct_wat);
01124 }
01125
01126
01127
01132 void ReadStructTP_PtSer( char* filename )
01133 {
01134 char lline[2001], *line;
01135 char scenario[20];
01136 int i, j, k;
01137 int lincnt;
01138
01139 char structName[20];
01140 int yyyy, mm, dd;
01141 double Jdate_read;
01142
01143 usrErr0 ( "\tControl structures' TP concen. time series...");
01144
01145 { if(H_OPSYS==UNIX)
01146 sprintf( modelFileName, "%s/%s/Data/%s.struct_TP", ModelPath, ProjName, filename );
01147 else sprintf( modelFileName, "%s%s:Data:%s.struct_TP", ModelPath, ProjName, filename );
01148 }
01149
01150
01151 if ( ( F_struct_TP = fopen( modelFileName, "r" ) ) == NULL )
01152 {
01153 printf( "Can't open %s file!\n", modelFileName ) ;
01154 exit(-1) ;
01155 }
01156
01157 fgets(lline, 2000, F_struct_TP);
01158 line=lline;
01159 sscanf(line, "%s",&scenario);
01160 if (strcmp(scenario,SimAlt) != 0) {
01161 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);
01162 exit(-1);
01163 }
01164 fgets(lline, 2000, F_struct_TP);
01165 line=lline;
01166 line = Scip( line, TAB );
01167 line = Scip( line, TAB );
01168 line = Scip( line, TAB );
01169
01170
01171 for ( j = 0; j < numTPser; j++ ) {
01172 IDhist[j] = -1;
01173 strcpy(structName,"NULL");
01174
01175 sscanf ( line,"%s\t",&structName);
01176 for ( i = 0; i < num_struct_hist+disAggNum; i++ ) {
01177 if (strcmp(Hist_data[i].S_nam , structName) == 0) {
01178 IDhist[j] = i;
01179 line = Scip( line, TAB );
01180 }
01181 }
01182 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);
01183 exit(-1);
01184 }
01185 }
01186
01187 lincnt = 0;
01188
01189 while ( fgets(lline, 2000, F_struct_TP) != NULL && !feof( F_struct_TP ) )
01190 {
01191 line = lline;
01192 sscanf( line, "%d %d %d",&yyyy,&mm,&dd);
01193
01194
01195 Jdate_read = julday( mm,dd,yyyy,0,0,0.0 );
01196 if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
01197 printf (" \n***Error\nStructure TP data (%s) init date > simulation start date\n", modelFileName); exit (-1);
01198 }
01199 if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
01200 {
01201 line = Scip( line, TAB );
01202 line = Scip( line, TAB );
01203 for ( j = 0; j < numTPser; j++ ) {
01204 line = Scip( line, TAB );
01205 if (*line == TAB) Hist_data[IDhist[j]].arrayP[((int)(lincnt))] = 0;
01206 else sscanf( line, "%f", &Hist_data[IDhist[j]].arrayP[((int)(lincnt))] );
01207 Hist_data[IDhist[j]].arrayP[((int)(lincnt))] *= 1.0E-6;
01208
01209 }
01210
01211 lincnt++;
01212 }
01213
01214 }
01215
01216 if (lincnt == 0 ) {
01217 printf (" \n***Error\nTP Time series data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
01218 }
01219 else if (lincnt != PORnumday ) {
01220 printf (" \n***Error\nTP Time series data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
01221 }
01222 sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
01223
01224 fclose (F_struct_TP);
01225 }
01226
01227
01228
01229
01234 void ReadStructTS_PtSer( char* filename )
01235 {
01236 char lline[2001], *line;
01237 char scenario[20];
01238 int i, j, k;
01239 int lincnt;
01240
01241 char structName[20];
01242 int yyyy, mm, dd;
01243 double Jdate_read;
01244
01245 usrErr0 ( "\tControl structures' Tracer concen. time series...");
01246
01247 { if(H_OPSYS==UNIX)
01248 sprintf( modelFileName, "%s/%s/Data/%s.struct_TS", ModelPath, ProjName, filename );
01249 else sprintf( modelFileName, "%s%s:Data:%s.struct_TS", ModelPath, ProjName, filename );
01250 }
01251
01252
01253 if ( ( F_struct_TS = fopen( modelFileName, "r" ) ) == NULL )
01254 {
01255 printf( "Can't open %s file!\n", modelFileName ) ;
01256 exit(-1) ;
01257 }
01258
01259 fgets(lline, 2000, F_struct_TS);
01260 line=lline;
01261 sscanf(line, "%s",&scenario);
01262 if (strcmp(scenario,SimAlt) != 0) {
01263 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);
01264 exit(-1);
01265 }
01266 fgets(lline, 2000, F_struct_TS);
01267 line=lline;
01268 line = Scip( line, TAB );
01269 line = Scip( line, TAB );
01270 line = Scip( line, TAB );
01271
01272
01273 for ( j = 0; j < numTSser; j++ ) {
01274 IDhist[j] = -1;
01275 strcpy(structName,"NULL");
01276
01277 sscanf ( line,"%s\t",&structName);
01278 for ( i = 0; i < num_struct_hist+disAggNum; i++ ) {
01279 if (strcmp(Hist_data[i].S_nam , structName) == 0) {
01280 IDhist[j] = i;
01281 line = Scip( line, TAB );
01282 }
01283 }
01284 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);
01285 exit(-1);
01286 }
01287 }
01288
01289 lincnt = 0;
01290
01291 while ( fgets(lline, 2000, F_struct_TS) != NULL && !feof( F_struct_TS ) )
01292 {
01293 line = lline;
01294 sscanf( line, "%d %d %d",&yyyy,&mm,&dd);
01295
01296
01297 Jdate_read = julday( mm,dd,yyyy,0,0,0.0 );
01298 if ((lincnt == 0) && (Jdate_read > Jdate_init) ) {
01299 printf (" \n***Error\nStructure Tracer data (%s) init date > simulation start date\n", modelFileName); exit (-1);
01300 }
01301 if ((Jdate_read >= Jdate_init) && (Jdate_read <= Jdate_end))
01302 {
01303 line = Scip( line, TAB );
01304 line = Scip( line, TAB );
01305 for ( j = 0; j < numTSser; j++ ) {
01306 line = Scip( line, TAB );
01307 if (*line == TAB) Hist_data[IDhist[j]].arrayS[((int)(lincnt))] = 0;
01308 else sscanf( line, "%f", &Hist_data[IDhist[j]].arrayS[((int)(lincnt))] );
01309 Hist_data[IDhist[j]].arrayS[((int)(lincnt))] *= 1.0;
01310
01311 }
01312
01313 lincnt++;
01314 }
01315
01316 }
01317
01318 if (lincnt == 0 ) {
01319 printf (" \n***Error\nTracer Time series data (%s) has date problem (0 records read)\n", modelFileName); exit (-1);
01320 }
01321 else if (lincnt != PORnumday ) {
01322 printf (" \n***Error\nTracer Time series data (%s) has problem (%d records read, %d expected)\n", modelFileName, lincnt, PORnumday); exit (-1);
01323 }
01324 sprintf (msgStr, "done, found %d records.",lincnt); usrErr(msgStr);
01325
01326 fclose (F_struct_TS);
01327 }
01328
01329
01330
01331
01332
01342 struct Schedule *Read_schedule ( char *sched_name, char *filename, float BASE_DATUM )
01343 {
01344 char ss[601], *s, scenario[30], modnam[30];
01345 int i, count = 1, sched_found=0;
01346 struct Schedule *Sce_Graph;
01347
01348 { if(H_OPSYS==UNIX)
01349 sprintf( modelFileName, "%s/%s/Data/%s.graph", ModelPath, ProjName, filename );
01350 else sprintf( modelFileName, "%s%s:Data:%s.graph", ModelPath, ProjName, filename );
01351 }
01352
01353
01354 if ( ( schedFile = fopen( modelFileName, "r" ) ) == NULL )
01355 {
01356 printf( "Can't open %s file!\n", modelFileName ) ;
01357 exit(-1) ;
01358 }
01359
01360 fgets( ss, 600, schedFile );
01361 fgets( ss, 600, schedFile );
01362 fgets( ss, 600, schedFile );
01363 fgets( ss, 600, schedFile ); sscanf( ss,"%s %s", &modnam, &scenario );
01364 if (strcmp(scenario,SimAlt) != 0) {
01365 sprintf(msgStr, "The scenario (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
01366 scenario, modelFileName, SimAlt); usrErr(msgStr);
01367 exit(-1);
01368 }
01369 if (strcmp(modnam,modelName) != 0) {
01370 sprintf(msgStr, "The model name (%s) found in the %s file doesn't match the one (%s) you asked for in Driver.parm!",
01371 modnam, modelFileName, modelName); usrErr(msgStr);
01372 exit(-1);
01373 }
01374
01375
01376 while ( fgets( ss, 600, schedFile ) != NULL && !feof( schedFile ) )
01377 { if ( !Wrdcmp ( ss, sched_name ) ) continue;
01378 sched_found = True;
01379
01380 fgets( ss, 600, schedFile );
01381 s = ss;
01382
01383 while ( *s++ != EOL ) if ( *s == '(' ) count++;
01384
01385
01386 if ( (Sce_Graph = (struct Schedule *) malloc( (size_t) sizeof( struct Schedule ))) == NULL )
01387 {
01388 printf( "Failed to allocate memory for Managed flow Schedule graph in (%s)\n ", sched_name ) ;
01389 exit( -2 ) ;
01390 };
01391
01392
01393 if ( (Sce_Graph->graph_points =
01394 (struct Points *) malloc( (size_t) sizeof( struct Points ) * count)) == NULL )
01395 {
01396 printf( "Failed to allocate memory for Managed flow Schedule graph in (%s)\n ", sched_name ) ;
01397 exit( -2 ) ;
01398 };
01399
01400 Sce_Graph->num_point = count;
01401
01402 s = ss;
01403 for ( i = 0; i < count && *s!=EOS; i++ )
01404 {
01405 s = Scip( s,'(');
01406 sscanf( s, "%f", &Sce_Graph->graph_points[i].time );
01407 s = Scip ( s, ',');
01408 sscanf( s, "%f", &Sce_Graph->graph_points[i].value );
01409 Sce_Graph->graph_points[i].value += BASE_DATUM;
01410 }
01411 }
01412
01413 if (!sched_found) {
01414 sprintf(msgStr, "The schedule name (%s) was not found in the %s schedule file! (Was named in the CanalData.struct file.)",
01415 sched_name, modelFileName); usrErr(msgStr);
01416 exit(-1);
01417 }
01418
01419
01420 sprintf(msgStr, "Read target schedule graph: %s, which contained %d time points.",strcpy(Sce_Graph->schedName, sched_name), Sce_Graph->num_point);
01421 WriteMsg(msgStr,1);
01422 NumScheds = NumScheds+1;
01423
01424
01425 fclose (schedFile);
01426 return Sce_Graph;
01427 }
01428
01429
01441 void Channel_configure ( float *ElevMap, struct Chan *channel_first )
01442 { int i, j;
01443 int cellLoc, cellLoc_i, cellLoc_j;
01444 int HV, k, L_range, R_range;
01445 float T_length;
01446 float C_length;
01447 float distTot;
01448 float A1, B1, A2, B2;
01449 float L, p_middle, m, i45, x, y, x1, y1, length, RoIL, RoIR;
01450 float cell_step = 1.0;
01451
01452 float a0, b0, a1, b1;
01453 int c_num, z;
01454 struct Cells *cell;
01455 struct Chan *channel;
01456 struct Chan_segm *Segm;
01457
01458 int numChek;
01459
01460
01461 int *marked;
01462 marked = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"marked");
01463 init_pvar(marked,NULL,'d',0);
01464
01465
01466 sprintf( modelFileName, "%s/%s/Output/Debug/CanalCells_interaction.txt", OutputPath, ProjName );
01467 if ( ( CanalCellInteractDebug = fopen( modelFileName, "w" ) ) == NULL )
01468 {
01469 printf( "Can't open %s file!\n", modelFileName ) ;
01470 exit(-1) ;
01471 }
01472
01473
01474 if ( (Chan_list =
01475 (struct Chan **) malloc( (size_t) sizeof( struct Chan *) * num_chan)) == NULL )
01476 {
01477 printf( "Failed to allocate memory for next canal (%d)\n ", channel->number ) ;
01478 exit( -2 ) ;
01479 };
01480
01481
01482 if ( (cell = (struct Cells *) malloc( (size_t) sizeof( struct Cells ))) == NULL )
01483 {
01484 printf( "Failed to allocate memory for cell structure\n " ) ;
01485 exit( -2 ) ;
01486 }
01487
01488
01489 channel = channel_first;
01490
01491
01492 while (channel != NULL)
01493 {
01494 Segm = channel->segments;
01495 c_num = channel->number;
01496 z = channel->levee;
01497 RoIL = channel->roil;
01498 RoIR = channel->roir;
01499
01500 Chan_list[c_num - CHo] = channel;
01501
01502 channel->cells = cell;
01503 distTot = 0;
01504 C_length = 0;
01505 C_Mark = 0;
01506
01507
01508 if (channel->width < 0) goto SKIPCHAN;
01509
01510 fprintf ( CanalCellInteractDebug, "N = %d levee = %d\n", c_num, z );
01511
01512
01513 cellLoc_i=(int)(Segm->x0+1);
01514 cellLoc_j=(int)(Segm->y0+1);
01515 cellLoc = T(cellLoc_i,cellLoc_j );
01516 Chan_list[c_num - CHo]->elev_start = ElevMap[cellLoc];
01517 if ( !ON_MAP[cellLoc] ) {
01518 printf( "Starting elevation of Reach %d is out of the active domain. Fix the reach location.\n",c_num);
01519 exit(-2);
01520 }
01521 fprintf ( CanalCellInteractDebug, " %d, %d, StartElev = %f\n", cellLoc_i,cellLoc_j, Chan_list[c_num - CHo]->elev_start );
01522
01523
01524 while (Segm != NULL)
01525 {
01526 a0 = Segm->x0;
01527 b0 = Segm->y0;
01528 a1 = Segm->x1;
01529 b1 = Segm->y1;
01530
01531 x1 = 0; y1 = 0;
01532
01533 T_length = sqrt( (b1 - b0)*(b1 - b0) + (a1 - a0)*(a1 - a0) );
01534 C_length += T_length;
01535
01536 if ( a1 != a0 && b1 != b0 )
01537 {
01538 A1 = (b1 - b0)/(a1 - a0);
01539 B1 = b0 - A1*a0;
01540 A2 = 1/A1; B2 = - B1/A1;
01541 }
01542
01543 x = a0; y = b0;
01544 i = Round (x); j = Round (y);
01545
01546
01547 while ( x1 != a1 || y1 != b1 )
01548 {
01549 if ( a1 == a0 )
01550 { x1 = a0;
01551 y1 = (j == Round (b1) ) ? b1 : j + STEP(b1 - b0)*cell_step;
01552 length = y1 - y;
01553 }
01554
01555 else if ( b1 == b0 )
01556 { x1 = (i == Round (a1) ) ? a1 : i + STEP(a1 - a0)*cell_step;
01557 y1 = b0;
01558 length = x1 - x;
01559 }
01560
01561 else if ( i == Round (a1) && j == Round (b1) )
01562 { x1 = a1; y1 = b1;
01563 length = T_length * (y1 - y)/(b1 - b0);
01564 C_Mark = 0;
01565 }
01566
01567 else
01568 { y1 = A1*(i + STEP(a1 - a0)*cell_step) + B1; length = T_length * (y1 - y)/(b1 - b0);
01569 x1 = A2*(j + STEP(b1 - b0)*cell_step) + B2; L = T_length * (x1 - x)/(a1 - a0);
01570
01571 if ( length < L )
01572 x1 = i + STEP(a1 - a0)*cell_step;
01573 else
01574 { y1 = j + STEP(b1 - b0)*cell_step; length = L;
01575 }
01576 }
01577
01578 if ( Abs (a1 - a0) < Abs (b1 - b0) )
01579
01580 { p_middle = (x1 + x)/2;
01581 m = i + cell_step/2;
01582 HV = 1;
01583 }
01584 else
01585 { p_middle = (y1 + y)/2;
01586 m = j + cell_step/2;
01587 HV = 0;
01588 }
01589
01590 distTot = distTot + length*celWid;
01591
01592
01593
01594 if ( Round (a1) == Round (a0) )
01595
01596 { if (b0 < b1)
01597 {
01598
01599 if ( p_middle > m )
01600 { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01601 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01602 }
01603 else
01604 { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i - 1, j, c_num, marked, distTot );
01605 L_range = (int) RoIL; R_range = (int)(RoIR - 1);
01606 }
01607 for ( k = 0; k < L_range; k++ )
01608 cell = MarkCell ( i - k - 1, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01609 for ( k = 0; k < R_range; k++ )
01610 cell = MarkCell ( i + k + 1, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01611 }
01612 else
01613 {
01614
01615 if ( p_middle > m )
01616 { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01617 L_range = (int)RoIL; R_range = (int)(RoIR - 1);
01618 }
01619 else
01620 { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i - 1, j, c_num, marked, distTot );
01621 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01622 }
01623 for ( k = 0; k < L_range; k++ )
01624 cell = MarkCell ( i + k + 1, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01625 for ( k = 0; k < R_range; k++ )
01626 cell = MarkCell ( i - k - 1, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01627 }
01628 }
01629 else if ( Round (b1) == Round (b0) )
01630 { if (a0 < a1)
01631 {
01632
01633 if ( p_middle > m )
01634 { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01635 L_range = (int)RoIL; R_range = (int)(RoIR -1);
01636 }
01637 else
01638 { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j - 1, c_num, marked, distTot );
01639 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01640 }
01641 for ( k = 0; k < L_range; k++ )
01642 cell = MarkCell ( i, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01643 for ( k = 0; k < R_range; k++ )
01644 cell = MarkCell ( i, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01645 }
01646 else
01647 {
01648
01649 if ( p_middle > m )
01650 { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01651 L_range = (int)(RoIL - 1); R_range = (int)RoIR;
01652 }
01653 else
01654 { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j - 1, c_num, marked, distTot );
01655 L_range = (int)RoIL; R_range = (int)(RoIR - 1);
01656 }
01657 for ( k = 0; k < L_range; k++ )
01658 cell = MarkCell ( i, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01659 for ( k = 0; k < R_range; k++ )
01660 cell = MarkCell ( i, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01661 }
01662 }
01663 else
01664
01665 { if ( A1>0 )
01666 i45 = ( HV ) ? p_middle - m : m - p_middle;
01667 else
01668 i45 = p_middle - m;
01669
01670
01671 if ( b0 < b1 )
01672 { if ( a0 < a1 )
01673 { if ( i45 > 0 )
01674 { cell = MarkCell ( i, j, LEFT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01675 cell = MarkCell ( i + 1, j - 1, RIGHT, length, cell, SOUTH, z, i + 1, j - 1, c_num, marked, distTot );
01676 }
01677 else
01678 { cell = MarkCell ( i, j, RIGHT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01679 cell = MarkCell ( i - 1, j + 1, LEFT, length, cell, EAST, z, i - 1, j + 1, c_num, marked, distTot );
01680 }
01681 for ( k = 0; k < RoIR-1; k++ )
01682 cell = MarkCell ( i + k + 1, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01683 for ( k = 0; k < RoIL-1; k++ )
01684 cell = MarkCell ( i - k - 1, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01685 }
01686 else
01687 { if ( i45 > 0 )
01688 { cell = MarkCell ( i, j, LEFT, length, cell, NONE, z, i, j, c_num, marked, distTot );
01689 cell = MarkCell ( i + 1, j + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01690 }
01691 else
01692 { cell = MarkCell ( i, j, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01693 cell = MarkCell ( i - 1, j - 1, LEFT, length, cell, NONE, z, i - 1, j - 1, c_num, marked, distTot );
01694 }
01695 for ( k = 0; k < RoIR-1; k++ )
01696 cell = MarkCell ( i + k + 1, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01697 for ( k = 0; k < RoIL-1; k++ )
01698 cell = MarkCell ( i - k - 1, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01699 }
01700 }
01701
01702 else
01703 { if ( a0 > a1 )
01704 { if ( i45 < 0 )
01705 { cell = MarkCell ( i, j, LEFT, length, cell, SOUTH, z, i, j, c_num, marked, distTot );
01706 cell = MarkCell ( i - 1, j + 1, RIGHT, length, cell, EAST, z, i - 1, j + 1, c_num, marked, distTot );
01707 }
01708 else
01709 { cell = MarkCell ( i, j, RIGHT, length, cell, EAST, z, i, j, c_num, marked, distTot );
01710 cell = MarkCell ( i + 1, j - 1, LEFT, length, cell, SOUTH, z, i + 1, j - 1, c_num, marked, distTot );
01711 }
01712
01713 for ( k = 0; k < RoIL-1; k++ )
01714 cell = MarkCell ( i + k + 1, j - k - 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01715 for ( k = 0; k < RoIR-1; k++ )
01716 cell = MarkCell ( i - k - 1, j + k + 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01717 }
01718 else
01719 { if ( i45 < 0 )
01720 { cell = MarkCell ( i, j, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01721 cell = MarkCell ( i - 1, j - 1, RIGHT, length, cell, NONE, z, i - 1, j - 1, c_num, marked, distTot );
01722 }
01723 else
01724 { cell = MarkCell ( i, j, RIGHT, length, cell, NONE, z, i, j, c_num, marked, distTot );
01725 cell = MarkCell ( i + 1, j + 1, LEFT, length, cell, ALL, z, i + 1, j + 1, c_num, marked, distTot );
01726 }
01727
01728 for ( k = 0; k < (int)(RoIL-1); k++ )
01729 cell = MarkCell ( i + k + 1, j + k + 1, LEFT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01730 for ( k = 0; k < (int)(RoIR-1); k++ )
01731 cell = MarkCell ( i - k - 1, j - k - 1, RIGHT, length, cell, ALL, z, i, j, c_num, marked, distTot );
01732 }
01733 }
01734 }
01735
01736 num_cell++;
01737
01738
01739
01740
01741 x = x1; y = y1;
01742 i = Round (x); if ( a1 < a0 && i == x ) i--;
01743
01744 j = Round (y); if ( b1 < b0 && j == y ) j--;
01745 C_Mark = 1;
01746
01747
01748 }
01749
01750 Segm = Segm->next_segm;
01751 }
01752
01753
01754 cellLoc_i=(int)(a1+1);
01755 cellLoc_j=(int)(b1+1);
01756 cellLoc = T(cellLoc_i,cellLoc_j );
01757 Chan_list[c_num - CHo]->elev_end = ElevMap[cellLoc];
01758
01759 Chan_list[c_num - CHo]->num_of_cells = num_cell;
01760 num_cell = 0;
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770 Chan_list[c_num - CHo]->length = C_length*celWid;
01771 Chan_list[c_num - CHo]->area = Chan_list[c_num - CHo]->width*Chan_list[c_num - CHo]->length;
01772 Chan_list[c_num - CHo]->minVol = Chan_list[c_num - CHo]->area*MINDEPTH;
01773 Chan_list[c_num - CHo]->seg_len = Chan_list[c_num - CHo]->length/Chan_list[c_num - CHo]->num_of_cells;
01774 Chan_list[c_num - CHo]->SW_flow_coef = sqrt(Chan_list[c_num - CHo]->seg_len) * sec_per_day * C_F;
01775 Chan_list[c_num - CHo]->SPG_flow_coef = Chan_list[c_num - CHo]->cond * GP_calibGWat;
01776 Chan_list[c_num - CHo]->elev_drop = Chan_list[c_num - CHo]->elev_start - Chan_list[c_num - CHo]->elev_end;
01777 Chan_list[c_num - CHo]->elev_slope = Chan_list[c_num - CHo]->elev_drop / Chan_list[c_num - CHo]->length ;
01778
01779 if ( !ON_MAP[cellLoc] ) {
01780 printf( "Ending elevation of Reach %d is out of the active domain. Fix the reach location.\n",c_num);
01781 exit(-2);
01782 }
01783 fprintf ( CanalCellInteractDebug, " %d, %d, EndingElev = %f, ElevSlope= %f, ReachLength= %f\n", cellLoc_i,cellLoc_j,
01784 Chan_list[c_num - CHo]->elev_end, Chan_list[c_num - CHo]->elev_slope, Chan_list[c_num - CHo]->length );
01785
01786 getCanalElev(c_num - CHo);
01787
01788 SKIPCHAN:
01789 channel = channel->next_in_list;
01790
01791 cell_last->next_cell = NULL;
01792
01793 }
01794
01795 fclose ( CanalCellInteractDebug );
01796
01797 return;
01798 }
01799
01800
01812 void getCanalElev (int chan_n)
01813 {
01814 int i, celcount=0;
01815 Chan_list[chan_n]->elev_avg = 0.0 ;
01816
01817 cell = Chan_list[chan_n]->cells;
01818 i = T(cell->x, cell->y);
01819
01820 while ( cell != NULL )
01821 {
01822 celcount++;
01823 cell->reachElev = Chan_list[chan_n]->elev_start - cell->reachDistDown * Chan_list[chan_n]->elev_slope;
01824 Chan_list[chan_n]->elev_avg = Chan_list[chan_n]->elev_avg + cell->reachElev;
01825 fprintf ( CanalCellInteractDebug, " Reach= %d, Cell (%d, %d), ReachElev= %f\n", chan_n+CHo, cell->x, cell->y, cell->reachElev );
01826 cell = cell->next_cell;
01827 }
01828 if (celcount == 0) { printf(" Error in getting cells along canal for elevation avg\n"); Exit(-1); }
01829 Chan_list[chan_n]->elev_avg = Chan_list[chan_n]->elev_avg/celcount;
01830 }
01831
01832
01833
01858 struct Cells *MarkCell (int x, int y, int c_ind, float length, struct Cells *cell, int direction, int levee,
01859 int xm, int ym, int c_num, int *marked, float distTot )
01860 { struct Cells *cell_next;
01861 int N, cellLoc;
01862
01863 if ( x > s0 || y > s1 )
01864 { printf ( "Error in definition of canal in cell row=%d col=%d ", x, y );
01865 exit (-2);
01866 }
01867
01868 N = T(xm+1, ym+1);
01869 if ( !ON_MAP[N] ) return cell;
01870
01871 if ( levee && direction != ALL )
01872 ON_MAP[N] = ((ON_MAP[N]-1) & (direction+100)) + 1;
01873
01874
01875
01876
01877 if ( !C_Mark ) return cell;
01878
01879
01880
01881
01882 cellLoc = T(x+1,y+1);
01883 if (marked[cellLoc] == c_num ) return cell;
01884 marked[cellLoc] = c_num;
01885
01886
01887
01888
01889 cell->x = x+1;
01890 cell->y = y+1;
01891 cell->ind = c_ind;
01892 cell->length = length*celWid;
01893 cell->reachDistDown = distTot;
01894
01895 fprintf ( CanalCellInteractDebug, " %d, %d, l/r = %d, distance= %7.1f\n", cell->x, cell->y, c_ind, cell->reachDistDown );
01896
01897
01898 if ( (cell_next = (struct Cells *) malloc((size_t) sizeof( struct Cells ))) == NULL )
01899 {
01900 printf( "Failed to allocate memory for cell structure\n " ) ;
01901 exit( -2 ) ;
01902 }
01903 cell->next_cell = cell_next;
01904 cell_last = cell;
01905
01906 return cell_next;
01907 }
01908
01909
01935 void FluxChannel ( int chan_n, float *SWH, float *ElevMap, float *MC,
01936 float *GWH, float *poros, float *GWcond,
01937 double *NA, double *PA, double *SA, double *GNA, double *GPA, double *GSA,
01938 float *Unsat, float *sp_yield)
01939 {
01940
01972 int converged = 0;
01973 float error = 1.0;
01974 float prev_error = 1.0;
01975 int iter = 0;
01976 int max_iter = CHANNEL_MAX_ITER;
01977 float init_factor = 0.5;
01978 float factor = init_factor;
01979
01980 int i;
01981 int attempt=0;
01982 double CanWatDep = Chan_list[chan_n]->wat_depth;
01983 double T_flux_S, T_flux_G, T_flux_L;
01984 double fluxO, fluxCk;
01985 double fluxS, fluxG, fluxL;
01986 double SW_coef, SPG_coef, GW_coef;
01987 double N_fl, P_fl, S_fl;
01988
01989 double I_Length;
01990 double seg_area;
01991 double MIN_VOL;
01992 double GW_head;
01993 double h_GWflow, h_SPflow;
01994 double dh;
01995 double Qin, Qout, Qout1;
01996 double SWater;
01997 double tot_head;
01998 double CH_vol;
01999 double CH_vol_R;
02000 double CH_bottElev;
02001 double H_rad_ch, H_rad_cell;
02002 double can_cell_M, can_cell_A;
02003 float errorConv;
02004
02005
02006 int equilib, revers=1;
02007 float H_pot_don, H_pot_rec;
02008
02009 float fluxTOunsat_don, fluxHoriz;
02010 float UnsatZ_don, UnsatCap_don, UnsatPot_don, Sat_vs_unsat;
02011 float UnsatZ_rec, UnsatCap_rec, UnsatPot_rec;
02012 float sfTOsat_don, unsatTOsat_don, unsatTOsat_rec, satTOsf_rec, horizTOsat_rec;
02013
02014 float sedwat_don, sedwat_rec;
02015 float sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec;
02016 float sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don;
02017 float sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz;
02018
02019
02020
02021
02022 Qin = 0; Qout = 0;
02023 structs = struct_first;
02024 while ( structs != NULL )
02025 {
02026
02027 if ( structs->canal_fr == chan_n+CHo )
02028 {
02029 if ( structs->flow > 0 ) Qout += structs->flow;
02030 else Qin -= structs->flow;
02031 }
02032 if ( structs->canal_to == chan_n+CHo )
02033 {
02034 if ( structs->flow > 0 ) Qin += structs->flow;
02035 else Qout -= structs->flow;
02036 }
02037
02038 structs = structs->next_in_list;
02039 }
02040
02041 CH_vol_R = Chan_list[chan_n]->area * CanWatDep;
02042 MIN_VOL = Chan_list[chan_n]->minVol;
02043 SPG_coef = Chan_list[chan_n]->SPG_flow_coef;
02044
02045 I_Length = Chan_list[chan_n]->seg_len;
02046 seg_area = I_Length * Chan_list[chan_n]->width;
02047 can_cell_M = Chan_list[chan_n]->area*CELL_SIZE;
02048 can_cell_A = Chan_list[chan_n]->area+CELL_SIZE;
02049
02050
02051
02052
02053 do
02054 {
02055 if ( Abs( error) < F_ERROR )
02056 {
02057 converged = 1;
02058
02059 errorConv=error;
02060 factor = 0;
02061 }
02062 else if ( iter >= (max_iter-1) )
02063 {
02064 if (attempt == 4) {
02065 converged = 1;
02066 sprintf(msgStr,"Day %6.1f: ERROR - couldn't converge on Canal %d, HALTING the iterative relaxation after 4 tries, error = %f.",
02067 SimTime.TIME, Chan_list[chan_n]->number, error);
02068 WriteMsg( msgStr,True ); dynERRORnum++;
02069 }
02070 else {
02071 if (debug > 3) { sprintf(msgStr,"Day %6.1f: Warning - couldn't converge on Canal %d, redoing the iterative relaxation, error = %f.",
02072 SimTime.TIME, Chan_list[chan_n]->number, error);
02073 WriteMsg( msgStr,True ); }
02074
02075 attempt++;
02076 iter = 0;
02077 if (attempt == 1) { factor = init_factor*2.0; max_iter = CHANNEL_MAX_ITER * 2; }
02078 else if (attempt == 2) { factor = init_factor*0.5; max_iter = CHANNEL_MAX_ITER * 2; }
02079 else if (attempt == 3) { factor = init_factor*4.0; max_iter = CHANNEL_MAX_ITER * 4; }
02080 else if (attempt == 4) { factor = init_factor*0.05; max_iter = CHANNEL_MAX_ITER * 4; }
02081 }
02082 }
02083
02084 else
02085 {
02086
02087 if (error * prev_error < 0 && iter != 1)
02088 factor = - sgn(error) * Abs( factor ) * 0.5;
02089
02090
02091 if (iter == 1 && error > 0) factor = - factor;
02092
02093 }
02094
02095 if (iter != 0 && !converged )
02096 {
02097 CanWatDep += factor;
02098 prev_error = error;
02099 if ( CanWatDep < MINDEPTH )
02100 {
02101 if (debug > 3) { sprintf(msgStr,"Day %6.1f: Warning - estimate is below minimum depth, revising estimate. (Canal %d, depth %f)",
02102 SimTime.TIME, Chan_list[chan_n]->number,CanWatDep);
02103 WriteMsg( msgStr,True ); }
02104
02105 CanWatDep = MINDEPTH;
02106
02107 factor = -factor*0.1;
02108 }
02109 }
02110
02111 T_flux_G = T_flux_S = T_flux_L = 0.0;
02112 CH_vol = Chan_list[chan_n]->area * CanWatDep;
02113
02114
02115 cell = Chan_list[chan_n]->cells;
02116
02117 while ( cell != NULL )
02118 {
02119 i = T(cell->x, cell->y);
02120 if ( !ON_MAP[i] )
02121
02122 { cell = cell->next_cell; continue;
02123 }
02124 if (basn[i] == 0) {
02125 sprintf(msgStr,"ERROR - Cell (%d, %d) is out-of-system.", cell->x, cell->y);
02126 WriteMsg( msgStr,True ); dynERRORnum++; }
02127
02128 fluxS = 0.; fluxG = 0.; fluxL = 0.;
02129
02130
02131
02132
02133 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] );
02134 GW_coef = GWcond[i] * poros[HAB[i]];
02135 SWater = SWH[i];
02136
02137 if (debug>0 && (SWater < -GP_MinCheck) ) {
02138 sprintf(msgStr,"Day %6.1f: capacityERR - neg sfwater along canal %d, candepth %f, cell (%d, %d), celldepth %f",
02139 SimTime.TIME,Chan_list[chan_n]->number,CanWatDep, cell->x, cell->y, SWater);
02140 WriteMsg( msgStr,True ); dynERRORnum++; }
02141
02142 GW_head = Min( GWH[i]/poros[HAB[i]], ElevMap[i]);
02143 tot_head = ramp(SWater) + GW_head;
02144 CH_bottElev = cell->reachElev - Chan_list[chan_n]->depth;
02145
02146 dh = ( CH_bottElev + CanWatDep ) - tot_head;
02147
02148
02149 H_rad_ch = Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
02150 SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);
02151
02152 H_rad_cell = Max ( ( seg_area * ramp(CanWatDep - Chan_list[chan_n]->depth) +
02153 SWater * (CELL_SIZE-seg_area) ) / CELL_SIZE, 0.0);
02154
02155
02156 if (dh > 0.0) {
02157 h_GWflow = Min(Chan_list[chan_n]->depth, CanWatDep );
02158 h_SPflow = Max(CH_bottElev + CanWatDep - ElevMap[i], 0.0);
02159 }
02160 else if (dh < 0.0) {
02161 h_GWflow = Max(GW_head-CH_bottElev, 0.0);
02162 h_SPflow = Max(tot_head-ElevMap[i], 0.0);
02163 }
02164 else {
02165 h_GWflow = 0.0;
02166 h_SPflow = 0.0;
02167 }
02168
02169
02170
02171 switch( Chan_list[chan_n]->levee )
02172 {
02173
02174 case 2:
02175
02176 fluxL = (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02177 fluxG = (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0);
02178 break;
02179 case 0:
02180
02181 if ( dh > 0 ) {
02182
02183 fluxS = f_Manning ( dh, H_rad_ch, SW_coef);
02184 }
02185 else if (SWater>GP_DetentZ) {
02186 fluxS = f_Manning ( dh, H_rad_cell, SW_coef) ;
02187 if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02188 fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02189 }
02190
02191
02192 fluxG = (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0);
02193 break;
02194 case 1:
02195 if ( cell->ind < 0 )
02196 {
02197 if ( dh > 0 ) {
02198
02199 fluxS = f_Manning ( dh, H_rad_ch, SW_coef);
02200 }
02201 else if (SWater>GP_DetentZ) {
02202 fluxS = f_Manning ( dh, H_rad_cell, SW_coef) ;
02203 if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02204 fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02205 }
02206
02207 }
02208 else
02209
02210 fluxL = (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02211
02212
02213 fluxG = (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0);
02214 break;
02215 case -1:
02216 if ( cell->ind > 0 )
02217 {
02218 if ( dh > 0 ) {
02219
02220 fluxS = f_Manning( dh, H_rad_ch, SW_coef);
02221 }
02222 else if (SWater>GP_DetentZ) {
02223 fluxS = f_Manning ( dh, H_rad_cell, SW_coef);
02224 if (-fluxS > (SWater-GP_DetentZ) *CELL_SIZE )
02225 fluxS = -(SWater-GP_DetentZ)*CELL_SIZE;
02226 }
02227 }
02228 else
02229
02230 fluxL = (h_SPflow > 0.0) ? (f_Ground ( dh, h_SPflow, SPG_coef, I_Length) ) : (0.0);
02231
02232
02233 fluxG = (h_GWflow > 0.0) ? (f_Ground ( dh, h_GWflow, GW_coef, I_Length) ) : (0.0);
02234 break;
02235 default:
02236 printf ("Error in levee location: Channel N %d", chan_n+CHo );
02237 exit (-1);
02238 break;
02239 }
02240
02241
02242
02243
02244
02245
02246
02247 if (fluxS>0 && ( (fluxCk = tot_head + fluxS/CELL_SIZE - (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) ) > 0.0 ) ) {
02248 fluxS = (can_cell_M*(CH_bottElev + CanWatDep)/can_cell_A - can_cell_M*tot_head/can_cell_A );
02249 fluxCk = tot_head + fluxS/CELL_SIZE - (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area);
02250 if (fluxCk > F_ERROR && debug>3) {
02251 sprintf(msgStr,"Day %6.1f: Warning - excessive head in cell after flux from Canal %d (%f m diff)",
02252 SimTime.TIME, Chan_list[chan_n]->number,fluxCk);
02253 WriteMsg( msgStr,True );
02254 }
02255 }
02256
02257
02258 if ( fluxS<0 && ( (fluxCk = (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) - (tot_head + fluxS/CELL_SIZE) ) > 0.0 ) ) {
02259 fluxS = (can_cell_M*(CH_bottElev + CanWatDep)/can_cell_A - can_cell_M*tot_head/can_cell_A );
02260 fluxCk = (CH_bottElev + CanWatDep - fluxS/Chan_list[chan_n]->area) - (tot_head + fluxS/CELL_SIZE);
02261 if (fluxCk > F_ERROR && debug>3) {
02262 sprintf(msgStr,"Day %6.1f: Warning - excessive head in Canal %d after flux from cell (%f m diff)",
02263 SimTime.TIME, Chan_list[chan_n]->number,fluxCk);
02264 WriteMsg( msgStr,True );
02265 }
02266 }
02267
02268
02269 if ( fluxS > 0 || fluxG > 0 || fluxL > 0)
02270 {
02271 fluxO = (fluxS > 0) ? fluxS : 0;
02272 fluxO += (fluxG > 0) ? fluxG : 0;
02273 fluxO += (fluxL > 0) ? fluxL : 0;
02274
02275
02276 if ( (fluxCk = CH_vol - MIN_VOL - fluxS - fluxG - fluxL) < 0.0 )
02277 { fluxCk = (fluxO>0) ? ( ( fluxO + fluxCk )/fluxO) : (0.0);
02278 if (fluxS>0) fluxS *= fluxCk;
02279 if (fluxG>0) fluxG *= fluxCk;
02280 if (fluxL>0) fluxL *= fluxCk;
02281 if ( (fluxCk = (CH_vol - MIN_VOL - fluxS - fluxG - fluxL)/Chan_list[chan_n]->area) < -GP_MinCheck ) {
02282 sprintf(msgStr,"Day %6.1f: capacityERR - excessive draining from Canal %d (%f m below min)", SimTime.TIME,Chan_list[chan_n]->number,fluxCk);
02283 WriteMsg( msgStr,True ); dynERRORnum++; }
02284 }
02285 }
02286
02287
02288 if (debug>99 && fluxG != 0.0) do {
02289
02290
02291
02292
02293
02294
02295 fluxHoriz = fluxG;
02296
02297
02298
02299 UnsatZ_rec = (fluxG>0) ? ( ElevMap[i] - GWH[i] / poros[HAB[i]] ) : (0.0) ;
02300 if (debug>2 && (UnsatZ_rec < -0.001 ) ) {
02301 sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in recipient cell (%d,%d) in canal-cell gw fluxes!",
02302 SimTime.TIME, UnsatZ_rec, cell->x, cell->y ); WriteMsg( msgStr,True ); dynERRORnum++; }
02303
02304 UnsatCap_rec = (fluxG>0) ? (UnsatZ_rec * poros[HAB[i]]) : (0.0);
02305
02306 UnsatPot_rec = (fluxG>0) ? (UnsatCap_rec - Unsat[i] ) : (0.0);
02307
02308
02309 horizTOsat_rec = (fluxG>0) ? (fluxHoriz) : (0.0);
02310 satTOsf_rec = (fluxG>0) ? (Max(fluxHoriz - UnsatPot_rec, 0.0) ) : (0.0);
02311 if (fluxG>0) unsatTOsat_rec = (UnsatZ_rec > 0.0) ?
02312
02313 ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[i]] ) / UnsatZ_rec * Unsat[i] ) :
02314 (0.0);
02315 else
02316 unsatTOsat_rec = 0.0;
02317
02318
02319 if (fluxG>0) H_pot_rec = (GWH[i] + horizTOsat_rec + unsatTOsat_rec - satTOsf_rec)
02320 / poros[HAB[i]] + (SWater + satTOsf_rec);
02321
02322
02323 equilib = 1;
02324
02325 } while (!equilib);
02326
02327
02328
02329
02330 CH_vol -= ( fluxS + fluxG + fluxL);
02331
02332 if ( converged )
02333 {
02334
02335
02336
02337 if ( fluxS != 0 ) {
02338 if ( fluxS > 0 )
02339 {
02340
02341 N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxS) : (0.0);
02342 P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxS) : (0.0);
02343 S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxS) : (0.0);
02344
02345 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02346 if ( (Chan_list[chan_n]->family != basn_list[basn[i]]->family) && (debug > 0) ) {
02347 sprintf(msgStr,"Day %6.1f: ERROR - no basin-basin canal-cell overland flow! (Reach %d, basn %d, with cell %d,%d, basin %d)",
02348 SimTime.TIME,Chan_list[chan_n]->number, Chan_list[chan_n]->basin, cell->x, cell->y, basn[i]);
02349 WriteMsg( msgStr,True ); dynERRORnum++;
02350 }
02351 else {
02352
02353 if ( !Chan_list[chan_n]->parent ) {
02354
02355 VOL_OUT_OVL[Chan_list[chan_n]->basin] += fluxS;
02356 P_OUT_OVL[Chan_list[chan_n]->basin] += P_fl;
02357 S_OUT_OVL[Chan_list[chan_n]->basin] += S_fl;
02358 }
02359 if ( !basn_list[basn[i]]->parent ) {
02360
02361 VOL_IN_OVL[basn[i]] += fluxS;
02362 P_IN_OVL[basn[i]] += P_fl;
02363 S_IN_OVL[basn[i]] += S_fl;
02364 }
02365 }
02366 }
02367 }
02368 else
02369 {
02370
02371 N_fl = ( SWH[i] > 0 ) ? NA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
02372 P_fl = ( SWH[i] > 0 ) ? PA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
02373 S_fl = ( SWH[i] > 0 ) ? SA[i]*fluxS/CELL_SIZE/SWH[i] : 0.;
02374 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02375
02376 if ( (Chan_list[chan_n]->family != basn_list[basn[i]]->family) && (debug > 0) ) {
02377 sprintf(msgStr,"Day %6.1f: ERROR - no basin-basin canal-cell overland flow! (Reach %d, basn %d, with cell %d,%d, basin %d)",
02378 SimTime.TIME,Chan_list[chan_n]->number, Chan_list[chan_n]->basin, cell->x, cell->y, basn[i]);
02379 WriteMsg( msgStr,True ); dynERRORnum++;
02380 }
02381 else {
02382
02383 if ( !basn_list[basn[i]]->parent ) {
02384
02385 VOL_OUT_OVL[basn[i]] -= fluxS;
02386 P_OUT_OVL[basn[i]] -= P_fl;
02387 S_OUT_OVL[basn[i]] -= S_fl;
02388 }
02389 if ( !Chan_list[chan_n]->parent ) {
02390
02391 VOL_IN_OVL[Chan_list[chan_n]->basin] -= fluxS;
02392 P_IN_OVL[Chan_list[chan_n]->basin] -= P_fl;
02393 S_IN_OVL[Chan_list[chan_n]->basin] -= S_fl;
02394 }
02395 }
02396 }
02397
02398
02399 }
02400
02401 NA[i] += N_fl;
02402 PA[i] += P_fl;
02403 SA[i] += S_fl;
02404 Chan_list[chan_n]->N -= N_fl;
02405 Chan_list[chan_n]->P -= P_fl;
02406 Chan_list[chan_n]->S -= S_fl;
02407 }
02408
02409
02410 if ( fluxG != 0 ) {
02411 if ( fluxG > 0 )
02412 {
02413 N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxG) : (0.0);
02414 P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxG) : (0.0);
02415 S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxG) : (0.0);
02416
02417 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02418 if (Chan_list[chan_n]->family !=
02419 basn_list[basn[i]]->family ) {
02420 if ( !Chan_list[chan_n]->parent ) {
02421
02422 VOL_OUT_GW[Chan_list[chan_n]->family] += fluxG;
02423 P_OUT_GW[Chan_list[chan_n]->family] += P_fl;
02424 S_OUT_GW[Chan_list[chan_n]->family] += S_fl;
02425 }
02426 if ( !basn_list[basn[i]]->parent ) {
02427
02428 VOL_IN_GW[basn_list[basn[i]]->family] += fluxG;
02429 P_IN_GW[basn_list[basn[i]]->family] += P_fl;
02430 S_IN_GW[basn_list[basn[i]]->family] += S_fl;
02431 }
02432 VOL_OUT_GW[Chan_list[chan_n]->basin] += fluxG;
02433 VOL_IN_GW[basn[i]] += fluxG;
02434 P_OUT_GW[Chan_list[chan_n]->basin] += P_fl;
02435 P_IN_GW[basn[i]] += P_fl;
02436 S_OUT_GW[Chan_list[chan_n]->basin] += S_fl;
02437 S_IN_GW[basn[i]] += S_fl;
02438 }
02439 else {
02440
02441 if ( !Chan_list[chan_n]->parent ) {
02442
02443 VOL_OUT_GW[Chan_list[chan_n]->basin] += fluxG;
02444 P_OUT_GW[Chan_list[chan_n]->basin] += P_fl;
02445 S_OUT_GW[Chan_list[chan_n]->basin] += S_fl;
02446 }
02447 if ( !basn_list[basn[i]]->parent ) {
02448
02449 VOL_IN_GW[basn[i]] += fluxG;
02450 P_IN_GW[basn[i]] += P_fl;
02451 S_IN_GW[basn[i]] += S_fl;
02452 }
02453 }
02454 }
02455 }
02456
02457
02458 else
02459 {
02460 N_fl = ( GWH[i] > 0 ) ? GNA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
02461 P_fl = ( GWH[i] > 0 ) ? GPA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
02462 S_fl = ( GWH[i] > 0 ) ? GSA[i]*fluxG/CELL_SIZE/GWH[i] : 0.;
02463 if (Chan_list[chan_n]->basin != basn[i] && debug > -999) {
02464 if (Chan_list[chan_n]->family !=
02465 basn_list[basn[i]]->family ) {
02466 if ( !basn_list[basn[i]]->parent ) {
02467
02468 VOL_OUT_GW[basn_list[basn[i]]->family] -= fluxG;
02469 P_OUT_GW[basn_list[basn[i]]->family] -= P_fl;
02470 S_OUT_GW[basn_list[basn[i]]->family] -= S_fl;
02471 }
02472 if ( !Chan_list[chan_n]->parent ) {
02473
02474 VOL_IN_GW[Chan_list[chan_n]->family] -= fluxG;
02475 P_IN_GW[Chan_list[chan_n]->family] -= P_fl;
02476 S_IN_GW[Chan_list[chan_n]->family] -= S_fl;
02477 }
02478 VOL_IN_GW[Chan_list[chan_n]->basin] -= fluxG;
02479 VOL_OUT_GW[basn[i]] -= fluxG;
02480 P_IN_GW[Chan_list[chan_n]->basin] -= P_fl;
02481 P_OUT_GW[basn[i]] -= P_fl;
02482 S_IN_GW[Chan_list[chan_n]->basin] -= S_fl;
02483 S_OUT_GW[basn[i]] -= S_fl;
02484 }
02485 else {
02486
02487 if ( !basn_list[basn[i]]->parent ) {
02488
02489 VOL_OUT_GW[basn[i]] -= fluxG;
02490 P_OUT_GW[basn[i]] -= P_fl;
02491 S_OUT_GW[basn[i]] -= S_fl;
02492 }
02493 if ( !Chan_list[chan_n]->parent ) {
02494
02495 VOL_IN_GW[Chan_list[chan_n]->basin] -= fluxG;
02496 P_IN_GW[Chan_list[chan_n]->basin] -= P_fl;
02497 S_IN_GW[Chan_list[chan_n]->basin] -= S_fl;
02498 }
02499 }
02500 }
02501 }
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514 GNA[i] += N_fl;
02515 GPA[i] += P_fl;
02516 GSA[i] += S_fl;
02517 Chan_list[chan_n]->N -= N_fl;
02518 Chan_list[chan_n]->P -= P_fl;
02519 Chan_list[chan_n]->S -= S_fl;
02520
02521
02522
02523
02524
02525 }
02526
02527
02528 if ( fluxL != 0 ) {
02529 if ( fluxL > 0 )
02530 {
02531 N_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->N/CH_vol_R*fluxL) : (0.0);
02532 P_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->P/CH_vol_R*fluxL) : (0.0);
02533 S_fl = (CH_vol_R > 0) ? (Chan_list[chan_n]->S/CH_vol_R*fluxL) : (0.0);
02534 if (Chan_list[chan_n]->basin != basn[i] && debug > -999 ) {
02535 if (Chan_list[chan_n]->family !=
02536 basn_list[basn[i]]->family ) {
02537 if ( !Chan_list[chan_n]->parent ) {
02538
02539 VOL_OUT_SPG[Chan_list[chan_n]->family] += fluxL;
02540 P_OUT_SPG[Chan_list[chan_n]->family] += P_fl;
02541 S_OUT_SPG[Chan_list[chan_n]->family] += S_fl;
02542 }
02543 if ( !basn_list[basn[i]]->parent ) {
02544
02545 VOL_IN_SPG[basn_list[basn[i]]->family] += fluxL;
02546 P_IN_SPG[basn_list[basn[i]]->family] += P_fl;
02547 S_IN_SPG[basn_list[basn[i]]->family] += S_fl;
02548 }
02549 VOL_OUT_SPG[Chan_list[chan_n]->basin] += fluxL;
02550 VOL_IN_SPG[basn[i]] += fluxL;
02551 P_OUT_SPG[Chan_list[chan_n]->basin] += P_fl;
02552 P_IN_SPG[basn[i]] += P_fl;
02553 S_OUT_SPG[Chan_list[chan_n]->basin] += S_fl;
02554 S_IN_SPG[basn[i]] += S_fl;
02555 }
02556 else {
02557
02558 if ( !Chan_list[chan_n]->parent ) {
02559
02560 VOL_OUT_SPG[Chan_list[chan_n]->basin] += fluxL;
02561 P_OUT_SPG[Chan_list[chan_n]->basin] += P_fl;
02562 S_OUT_SPG[Chan_list[chan_n]->basin] += S_fl;
02563 }
02564 if ( !basn_list[basn[i]]->parent ) {
02565
02566 VOL_IN_SPG[basn[i]] += fluxL;
02567 P_IN_SPG[basn[i]] += P_fl;
02568 S_IN_SPG[basn[i]] += S_fl;
02569 }
02570 }
02571 }
02572 if ( tot_head > GW_head ) {
02573 NA[i] += N_fl;
02574 PA[i] += P_fl;
02575 SA[i] += S_fl;
02576 }
02577 else {
02578 GNA[i] += N_fl;
02579 GPA[i] += P_fl;
02580 GSA[i] += S_fl;
02581 }
02582 }
02583 else if (fluxL < 0)
02584 {
02585 if ( tot_head > GW_head ) {
02586 N_fl = ( SWH[i] > 0 ) ? NA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
02587 P_fl = ( SWH[i] > 0 ) ? PA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
02588 S_fl = ( SWH[i] > 0 ) ? SA[i]*fluxL/CELL_SIZE/SWH[i] : 0.;
02589 NA[i] += N_fl;
02590 PA[i] += P_fl;
02591 SA[i] += S_fl;
02592 }
02593 else {
02594 N_fl = ( GWH[i] > 0 ) ? GNA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
02595 P_fl = ( GWH[i] > 0 ) ? GPA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
02596 S_fl = ( GWH[i] > 0 ) ? GSA[i]*fluxL/CELL_SIZE/GWH[i] : 0.;
02597 GNA[i] += N_fl;
02598 GPA[i] += P_fl;
02599 GSA[i] += S_fl;
02600 }
02601 if (Chan_list[chan_n]->basin != basn[i] && debug > -999) {
02602 if (Chan_list[chan_n]->family !=
02603 basn_list[basn[i]]->family ) {
02604 if ( !basn_list[basn[i]]->parent ) {
02605
02606 VOL_OUT_SPG[basn_list[basn[i]]->family] -= fluxL;
02607 P_OUT_SPG[basn_list[basn[i]]->family] -= P_fl;
02608 S_OUT_SPG[basn_list[basn[i]]->family] -= S_fl;
02609 }
02610 if ( !Chan_list[chan_n]->parent ) {
02611
02612 VOL_IN_SPG[Chan_list[chan_n]->family] -= fluxL;
02613 P_IN_SPG[Chan_list[chan_n]->family] -= P_fl;
02614 S_IN_SPG[Chan_list[chan_n]->family] -= S_fl;
02615 }
02616 VOL_IN_SPG[Chan_list[chan_n]->basin] -= fluxL;
02617 VOL_OUT_SPG[basn[i]] -= fluxL;
02618 P_IN_SPG[Chan_list[chan_n]->basin] -= P_fl;
02619 P_OUT_SPG[basn[i]] -= P_fl;
02620 S_IN_SPG[Chan_list[chan_n]->basin] -= S_fl;
02621 S_OUT_SPG[basn[i]] -= S_fl;
02622 }
02623 else {
02624
02625 if ( !basn_list[basn[i]]->parent ) {
02626
02627 VOL_OUT_SPG[basn[i]] -= fluxL;
02628 P_OUT_SPG[basn[i]] -= P_fl;
02629 S_OUT_SPG[basn[i]] -= S_fl;
02630 }
02631 if ( !Chan_list[chan_n]->parent ) {
02632
02633 VOL_IN_SPG[Chan_list[chan_n]->basin] -= fluxL;
02634 P_IN_SPG[Chan_list[chan_n]->basin] -= P_fl;
02635 S_IN_SPG[Chan_list[chan_n]->basin] -= S_fl;
02636 }
02637 }
02638 }
02639 }
02640
02641 Chan_list[chan_n]->N -= N_fl;
02642 Chan_list[chan_n]->P -= P_fl;
02643 Chan_list[chan_n]->S -= S_fl;
02644 }
02645
02646
02647
02648 SWH[i] += (fluxS)/CELL_SIZE;
02649 GWH[i] += (fluxG)/CELL_SIZE;
02650 if ( tot_head > GW_head ) SWH[i] += fluxL/CELL_SIZE;
02651 else GWH[i] += fluxL/CELL_SIZE;
02652
02653
02654 }
02655
02656 T_flux_S += fluxS;
02657 T_flux_G += fluxG;
02658 T_flux_L += fluxL;
02659
02660 cell = cell->next_cell;
02661 }
02662
02663
02664 error = (CanWatDep - Chan_list[chan_n]->wat_depth) - (Qin - Qout - T_flux_S - T_flux_G - T_flux_L)/ Chan_list[chan_n]->area;
02665 if (converged && error != errorConv ) {
02666
02667 if (fabs(error-errorConv)>F_ERROR ) {
02668 sprintf(msgStr,"Day %6.1f: ERROR - converg error != to error after cell-canal fluxes (Chan %d, est= %5.3f, last depth= %5.3f m)",
02669 SimTime.TIME, Chan_list[chan_n]->number, CanWatDep, Chan_list[chan_n]->wat_depth);
02670 WriteMsg( msgStr,True ); dynERRORnum++;
02671 }
02672
02673 if (CanWatDep>Chan_list[chan_n]->depth*3.0 || CanWatDep < 0.0) {
02674 CanWatDep = MINDEPTH;
02675 sprintf(msgStr," Reset Chan %d depth to minimum allowed (%5.3f m)",
02676 Chan_list[chan_n]->number, MINDEPTH);
02677 WriteMsg( msgStr,True ); }
02678
02679
02680
02681
02682 if (debug>4) {
02683 printf("\nDay %6.1f: ERROR - converg error != to error after cell-canal fluxes (Chan %d, est= %5.3f m)\n",
02684 SimTime.TIME, Chan_list[chan_n]->number, CanWatDep);
02685 exit (-1);
02686
02687 }
02688
02689 }
02690
02691 ++iter;
02692
02693 } while (!converged );
02694
02695 if ( factor != 0 && Qout > 0)
02696 {
02697 sprintf(msgStr,"Day %6.1f: ERROR - hit bottom of Canal %d, stopped the iterative relaxation", SimTime.TIME, Chan_list[chan_n]->number);
02698 WriteMsg( msgStr,True );
02699 dynERRORnum++;
02700
02701
02702
02703 Qout1 = (Chan_list[chan_n]->wat_depth - CanWatDep)*Chan_list[chan_n]->area + Qin - T_flux_S - T_flux_G - T_flux_L;
02704 structs = struct_first;
02705 while ( structs != NULL )
02706 {
02707 if ( structs->canal_fr == chan_n+CHo ) structs->flow *= Qout1/Qout;
02708
02709 structs = structs->next_in_list;
02710 }
02711 Qout = Qout1;
02712 }
02713
02714 Chan_list[chan_n]->wat_depth = CanWatDep;
02715 CH_vol = CanWatDep * Chan_list[chan_n]->area;
02716
02717 if( debug > 3 )
02718 {
02719
02720 Chan_list[chan_n]->P_con = (CH_vol > 0) ? (Chan_list[chan_n]->P/CH_vol*1000.0) : 0.0;
02721 Chan_list[chan_n]->S_con = (CH_vol > 0) ? (Chan_list[chan_n]->S/CH_vol) : 0.0;
02722 fprintf( canDebugFile,
02723 "%4d %8.3f %8.4f %8.4f %12.1f %12.1f %12.1f %12.1f %12.1f %5d %9.1f %7.5f \n",
02724 Chan_list[chan_n]->number, Chan_list[chan_n]->wat_depth, Chan_list[chan_n]->P_con,
02725 Chan_list[chan_n]->S_con, T_flux_S, T_flux_G, T_flux_L, Qin, Qout, iter, Chan_list[chan_n]->area, error );
02726
02727 }
02728
02729 return;
02730 }
02731
02732
02741 float f_Manning( float delta, float Water, float SW_coef )
02742 { float a_delta;
02743
02744 a_delta = Abs (delta);
02745 return ( sgn( delta ) * SW_coef * pow(Water, GP_mannDepthPow) * sqrt(a_delta) * canstep);
02746
02747 }
02748
02749
02759 float f_Ground ( float dh, float height, float GW_coef, float I_Length )
02760 {
02761 float w;
02762
02763
02764 w = dh * I_Length * GW_coef / (0.5*celWid) * height * canstep;
02765
02766 return ( w);
02767
02768 }
02769
02770
02811 void Flows_in_Structures ( float *SWH, float *GW, float *poros, float *Elev, double *NA, double *PA, double *SA )
02812 {
02816 struct Structure *structs, *struct_hist_start;
02817 int flow_code;
02818
02819 structs = struct_first;
02820
02821 while ( structs != NULL )
02822 {
02823 if ( structs->flag < 0 || (structs->flag >1) ) {
02824 structs->flow = 0.0;
02825 structs->conc_P = 0.0;
02826 structs->conc_N = 0.0;
02827 structs->conc_S = 0.0;
02828 }
02829 else {
02830
02833 if ( structs->flagTWsched == 1 )
02834 {
02835 structs->TW_stage = GetGraph ( structs->TW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 0.0) );
02836
02837
02838 if (structs->flagTWtide == 1) structs->TW_stage = structs->TW_stage + (SimTime.TIME/365) * GP_SLRise ;
02839 }
02840
02841 if ( structs->flagHWsched == 1 )
02842 {
02843 structs->HW_stage = GetGraph ( structs->HW_graph, ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 0.0) );
02844
02845
02846 if (structs->flagHWtide == 1) structs->HW_stage = structs->HW_stage + (SimTime.TIME/365) * GP_SLRise ;
02847 }
02848
02849 if (strcmp(structs->SrcDest,"CanCan") == 0) flow_code = 0;
02850 if (strcmp(structs->SrcDest,"CelCel") == 0) flow_code = 1;
02851 if (strcmp(structs->SrcDest,"CanCel") == 0) flow_code = 2;
02852 if (strcmp(structs->SrcDest,"CelCan") == 0) flow_code = 3;
02853
02854
02855 switch(flow_code) {
02856 case 0:
02857 if (structs->flag == 1) flowData_CanCan(structs, struct_hist_start, SWH, NA, PA, SA);
02858 else flowCalc_CanCan(structs, SWH, GW, poros, Elev, NA, PA, SA);
02859 break;
02860 case 1:
02861 if (structs->flag == 1) flowData_CelCel(structs, SWH, NA, PA, SA);
02862 else flowCalc_CelCel(structs, SWH, Elev, NA, PA, SA);
02863 break;
02864 case 2:
02865 if (structs->flag == 1) flowData_CanCel(structs, struct_hist_start, SWH, NA, PA, SA);
02866 else flowCalc_CanCel(structs, SWH, GW, poros, Elev, NA, PA, SA);
02867 break;
02868 case 3:
02869 if (structs->flag == 1) flowData_CelCan(structs, SWH, NA, PA, SA);
02870 else flowCalc_CelCan(structs, SWH, GW, poros, Elev, NA, PA, SA);
02871 break;
02872 default: {
02873 printf ("\nError in source-destination data for structure N %s: impossible condition! \n", structs->S_nam );
02874 exit (-1);
02875 }
02876 }
02877 }
02878
02879
02880 structs->SumFlow += structs->flow;
02881 structs->Sum_P += structs->conc_P * structs->flow;
02882
02883 structs->Sum_S += structs->conc_S * structs->flow;
02884
02885 structs->multiOut = 0;
02886
02887
02888 if ( printchan ) {
02889 fprintf( WstructOutFile, "%7.0f\t", (structs->SumFlow)/(2446.848) );
02890 fprintf( WstructOutFile_P, "%7.4f\t",
02891 (structs->SumFlow > 0.0) ? (structs->Sum_P/structs->SumFlow*1000.0) : (0.0) );
02892
02893
02894 fprintf( WstructOutFile_S, "%7.4f\t",
02895 (structs->SumFlow > 0.0) ? (structs->Sum_S/structs->SumFlow) : (0.0) );
02896
02897 structs->SumFlow = 0.0;
02898 structs->Sum_P = 0.0;
02899
02900 structs->Sum_S = 0.0;
02901 }
02902
02903
02904 structs = structs->next_in_list;
02905 }
02906
02907 if (printchan) {
02908 fflush( WstructOutFile );
02909 fflush( WstructOutFile_P );
02910
02911 fflush( WstructOutFile_S );
02912 }
02913
02914 return;
02915
02916 }
02917
02918
02919
02927 void flowData_CanCan( struct Structure *structs, struct Structure *struct_hist_start, float *SWH, double *NA, double *PA, double *SA)
02928 {
02929 int cell_to, cell_fr, can_to, can_fr;
02930 int cell_struct;
02931 double CH_vol;
02932 double Nconc, Pconc, Sconc;
02933 double N_fl, P_fl, S_fl;
02934 double chan_over;
02935 double ChanHistOut;
02936
02937 can_fr = structs->canal_fr - CHo;
02938 can_to = structs->canal_to - CHo;
02939 cell_struct = T(structs->str_cell_i,structs->str_cell_j);
02940
02941
02942 if (structs->multiOut == 1) return;
02943
02944
02945 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
02946 if (structs->flow == 0.0) {
02947 structs->conc_P = 0.0;
02948 structs->conc_N = 0.0;
02949 structs->conc_S = 0.0;
02950 return;
02951 }
02952 else if (structs->flow < 0.0) {
02953 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
02954 SimTime.TIME, Chan_list[can_to]->number );
02955 WriteMsg( msgStr,True );
02956 dynERRORnum++;
02957 return;
02958 }
02959
02960
02961
02962
02963
02964 struct_hist_start = structs;
02965 ChanHistOut = 0.0;
02966
02967 while ( structs != NULL )
02968 {
02969 if ( structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
02970 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
02971 ChanHistOut += structs->flow;
02972 }
02973
02974 structs = structs->next_in_list;
02975 }
02976
02977 structs = struct_hist_start;
02978
02979
02980
02981
02982
02983
02984
02985 CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth;
02986
02987 Pconc = (structs->TPser == -1) ?
02988 ( (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0) ) :
02989 ( (structs->TPser == 0) ? (structs->TP) : ( Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
02990 Nconc = (structs->TNser == -1) ?
02991 ( (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0) ) :
02992 ( (structs->TNser == 0) ? (structs->TN) : ( Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
02993 Sconc = (structs->TSser == -1) ?
02994 ( (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0) ) :
02995 ( (structs->TSser == 0) ? (structs->TS) : ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
02996
02997
02998
02999 CH_vol = Chan_list[can_fr]->area*(ramp(Chan_list[can_fr]->wat_depth-MINDEPTH) );
03000 chan_over = (ChanHistOut > 0.0 ) ? ( CH_vol/ChanHistOut ) : 1.0;
03001
03002 do
03003 {
03004 if (structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
03005
03006 structs->multiOut = 1;
03007
03008
03009 if ( chan_over < 1.0 ) {
03010 if( debug > 0 ) {
03011 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",
03012 SimTime.TIME, structs->S_nam, ChanHistOut, Chan_list[can_fr]->number,
03013 CH_vol, (ChanHistOut - CH_vol), Chan_list[can_fr]->wat_depth );
03014
03015 }
03016
03017 structs->flow = (structs->flow) * chan_over ;
03018 }
03019
03020
03021 Chan_list[can_fr]->sumHistOut += structs->flow;
03022
03023
03024 can_to = structs->canal_to - CHo;
03025 cell_to = T(structs->cell_i_to,structs->cell_j_to);
03026
03027 structs->conc_P = Pconc;
03028 structs->conc_N = Nconc;
03029 structs->conc_S = (structs->TSser == -1) ?
03030 (Sconc) :
03031 ( (structs->TSser == 0) ? (structs->TS) :
03032 ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
03033
03034 P_fl = structs->flow * structs->conc_P;
03035 N_fl = structs->flow * structs->conc_N;
03036 S_fl = structs->flow * structs->conc_S;
03037
03038 Chan_list[can_fr]->N -= N_fl ;
03039 Chan_list[can_fr]->P -= P_fl ;
03040 Chan_list[can_fr]->S -= (structs->TSser == -1) ? (S_fl) : (0);
03041
03042
03043
03044
03045
03046 if (can_to > 0) {
03047
03048
03049 Chan_list[can_to]->sumHistIn += structs->flow;
03050
03051 Chan_list[can_to]->N += N_fl ;
03052 Chan_list[can_to]->P += P_fl ;
03053 Chan_list[can_to]->S += S_fl ;
03054
03055
03056
03057 if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > -999 ) {
03058 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03059 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow;
03060 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03061 P_IN_STR[Chan_list[can_to]->basin] += P_fl;
03062 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0);
03063 S_IN_STR[Chan_list [can_to]->basin] += S_fl;
03064 }
03065
03066
03067
03068 }
03069 else if (cell_to > 0) {
03070
03071
03072 if (structs->cell_i_to > 2 ) {
03073 SWH[cell_to] += structs->flow/CELL_SIZE;
03074 PA[cell_to] += P_fl;
03075 NA[cell_to] += N_fl;
03076 SA[cell_to] += S_fl;
03077
03078 if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
03079 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03080 VOL_IN_STR[basn[cell_to]] += structs->flow;
03081 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03082 P_IN_STR[basn[cell_to]] += P_fl;
03083 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0);
03084 S_IN_STR[basn[cell_to]] += S_fl;
03085
03086
03087
03088
03089
03090
03091
03092
03093 }
03094
03095
03096
03097 }
03098 else if (structs->cell_i_to == 2 && debug > -999 ) {
03099 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03100 VOL_OUT_STR[0] += structs->flow;
03101 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03102 P_OUT_STR[0] += P_fl;
03103 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0);
03104 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0);
03105 }
03106 }
03107
03108 if (structs->flow == 0.0) {
03109 structs->conc_P = 0.0;
03110 structs->conc_N = 0.0;
03111 structs->conc_S = 0.0;
03112 }
03113 else if (structs->flow < 0.0) {
03114 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from another cell/chan into chan %d",
03115 SimTime.TIME,Chan_list[can_fr]->number );
03116 WriteMsg( msgStr,True );
03117 dynERRORnum++;
03118 }
03119 }
03120
03121 structs = structs->next_in_list;
03122 } while (structs != NULL);
03123
03124 structs = struct_hist_start;
03125
03126
03127
03128 }
03129
03130
03131
03156 void flowCalc_CanCan( struct Structure *structs, float *SWH, float *GW,
03157 float *poros, float *Elev, double *NA, double *PA, double *SA)
03158 {
03159 int can_to, can_fr;
03160 int cell_struct;
03161 float HeadH, HeadT, HeadH2, HeadT2;
03162 double CH_vol;
03163 double N_fl, P_fl, S_fl;
03164 float HeadH_drop, HeadT_drop;
03165
03166 float HWcheck, TWcheck, HWdiff, TWdiff;
03167 int struct_open;
03168
03169 can_fr = structs->canal_fr - CHo;
03170 can_to = structs->canal_to - CHo;
03171 cell_struct = T(structs->str_cell_i,structs->str_cell_j);
03172
03173
03174 if ( structs->w_coef == 0.0) {
03175 HeadH_drop = 0.5 * Abs(Chan_list[can_fr]->elev_drop);
03176 HeadT_drop = 0.5 * Abs(Chan_list[can_to]->elev_drop);
03177
03178
03179 HeadH = HeadH_drop +
03180 Elev[cell_struct] - Chan_list[can_fr]->depth + Chan_list[can_fr]->wat_depth
03181 + (Chan_list[can_fr]->sumHistIn - Chan_list[can_fr]->sumHistOut)/Chan_list[can_fr]->area
03182 + (Chan_list[can_fr]->sumRuleIn - Chan_list[can_fr]->sumRuleOut)/Chan_list[can_fr]->area;
03183
03184 HeadT = -HeadT_drop +
03185 Elev [cell_struct] - Chan_list[can_to]->depth + Chan_list[can_to]->wat_depth
03186 + (Chan_list[can_to]->sumHistIn - Chan_list[can_to]->sumHistOut)/Chan_list[can_to]->area
03187 + (Chan_list[can_to]->sumRuleIn - Chan_list[can_to]->sumRuleOut)/Chan_list[can_to]->area;
03188
03189 if ( HeadH>HeadT ) {
03190 structs->flow =
03191 Chan_list[can_fr]->area*Chan_list[can_to]->area/
03192 (Chan_list[can_fr]->area+Chan_list[can_to]->area) * (HeadH-HeadT) ;
03193 if (debug >3) {
03194 HeadH2 = HeadH - structs->flow/Chan_list[can_fr]->area;
03195 HeadT2 = HeadT + structs->flow/Chan_list[can_to]->area;
03196 sprintf( msgStr, "Day %6.1f: virtual flow from chan %d (%f m) to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
03197 SimTime.TIME,Chan_list[can_fr]->number,HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth,
03198 Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
03199 WriteMsg( msgStr,True );
03200 if (Abs(HeadH2 - HeadT2) > GP_MinCheck)
03201 { sprintf( msgStr, "Day %6.1f: ERROR virtual flow from chan %d to chan %d, New headH-headT=%f, Orig HeadH-HeadT=%f",
03202 SimTime.TIME,Chan_list[can_fr]->number,Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
03203 WriteMsg( msgStr,True ); dynERRORnum++; }
03204
03205 }
03206
03207 CH_vol = Chan_list[can_fr]->area* ramp((HeadH-HeadH_drop)-Elev[cell_struct]+Chan_list[can_fr]->depth-MINDEPTH);
03208 if (structs->flow > CH_vol) {
03209 structs->flow = CH_vol;
03210 HeadH2 = HeadH - structs->flow/Chan_list[can_fr]->area;
03211 HeadT2 = HeadT + structs->flow/Chan_list[can_to]->area;
03212 if (debug>3 && Abs(HeadH2 - HeadT2) > GP_MinCheck) {
03213 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",
03214 SimTime.TIME,Chan_list[can_fr]->number,CH_vol/Chan_list[can_fr]->area,Chan_list[can_to]->number,(HeadH2-HeadT2),(HeadH-HeadT) );
03215 WriteMsg( msgStr,True ); }
03216 }
03217
03218
03219 CH_vol = CH_vol + MINDEPTH*Chan_list[can_fr]->area;
03220
03221 structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
03222 structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
03223 structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
03224 }
03225
03226 else structs->flow = 0.0;
03227
03228 }
03229
03230
03231
03232 else
03233 {
03234 if ( (structs->flagTWsched != 1) || (structs->flagHWsched != 1) ) {
03235 sprintf (msgStr,"\nYou must assign both Headwater and Tailwater schedules for the rule-based canal-to-canal structure %s!\n", structs->S_nam);
03236 usrErr( msgStr );
03237 exit(-1);
03238 }
03239 if ( (structs->cell_HW == 0) || (structs->cell_TW == 0) ) {
03240 sprintf (msgStr,"\nYou must assign cells to check for both the Headwater and Tailwater schedules for the rule-based canal-to-canal structure %s!\n", structs->S_nam);
03241 usrErr( msgStr );
03242 exit(-1);
03243 }
03244
03245 struct_open = 0;
03246
03247
03248
03249
03250 HWcheck = GW[structs->cell_HW]/poros[HAB[structs->cell_HW]] + SWH[structs->cell_HW] ;
03251 TWcheck = GW[structs->cell_TW]/poros[HAB[structs->cell_TW]] + SWH[structs->cell_TW] ;
03252
03253
03254 struct_open = ( (structs->TW_stage > TWcheck) && (structs->HW_stage < HWcheck) ) ? (1) : (0);
03255
03256 if ( struct_open == 1)
03257 {
03258 HWdiff = HWcheck - structs->HW_stage;
03259 TWdiff = structs->TW_stage - TWcheck;
03260 CH_vol = Chan_list[can_fr]->area * ramp(Chan_list[can_fr]->wat_depth-MINDEPTH);
03261 structs->flow = Min(structs->w_coef * ( HWdiff+TWdiff ) * canstep, CH_vol );
03262 }
03263 else structs->flow = 0.0;
03264
03265
03266 }
03267
03268
03269 if (structs->flow == 0.0) {
03270 structs->conc_P = 0.0;
03271 structs->conc_N = 0.0;
03272 structs->conc_S = 0.0;
03273 return;
03274 }
03275 else if (structs->flow < 0.0 && structs->w_coef != 0) {
03276 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG RULE-BASED FLOW from chan %d",
03277 SimTime.TIME,Chan_list[can_to]->number);
03278 WriteMsg( msgStr,True );
03279 return;
03280 }
03281
03282
03283 CH_vol = CH_vol + MINDEPTH*Chan_list[can_fr]->area;
03284
03285 structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
03286 structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
03287 structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
03288
03289
03290
03291
03292 Chan_list[can_fr]->sumRuleOut += structs->flow;
03293 Chan_list[can_to]->sumRuleIn += structs->flow;
03294
03295
03296 P_fl = structs->flow * structs->conc_P;
03297 N_fl = structs->flow * structs->conc_N;
03298 S_fl = structs->flow * structs->conc_S;
03299
03300 Chan_list[can_fr]->N -= N_fl ;
03301 Chan_list[can_fr]->P -= P_fl ;
03302 Chan_list[can_fr]->S -= S_fl ;
03303
03304 Chan_list[can_to]->N += N_fl ;
03305 Chan_list[can_to]->P += P_fl ;
03306 Chan_list[can_to]->S += S_fl ;
03307
03308
03309 if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > 0 ) {
03310 sprintf( msgStr, "Day %6.1f: ERROR - no hydrologic basin-basin flows in canal-to-canal structure, canals %d - %d",
03311 SimTime.TIME, Chan_list[can_to]->number, Chan_list[can_fr]->number);
03312 WriteMsg( msgStr,True ); dynERRORnum++;
03313 }
03314
03315 }
03316
03317
03318
03325 void flowData_CelCel(struct Structure *structs, float *SWH, double *NA, double *PA, double *SA)
03326 {
03327 int cell_to, cell_fr;
03328
03329 double cell_vol;
03330 double N_fl, P_fl, S_fl;
03331
03332 cell_to = T(structs->cell_i_to,structs->cell_j_to);
03333 cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);
03334
03335
03336
03337 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03338 if (structs->flow == 0.0) {
03339 structs->conc_P = 0.0;
03340 structs->conc_N = 0.0;
03341 structs->conc_S = 0.0;
03342 return;
03343 }
03344 else if (structs->flow < 0.0) {
03345 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from cell (%d,%d)",
03346 SimTime.TIME, structs->cell_i_to, structs->cell_j_to);
03347 WriteMsg( msgStr,True );
03348 return;
03349 }
03350
03351 if (structs->cell_i_fr != 2 )
03352 {
03353
03354 cell_vol = SWH[cell_fr] * CELL_SIZE;
03355 if (structs->flow > cell_vol) {
03356 if( debug > 0 ) {
03357 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",
03358 SimTime.TIME, structs->S_nam, structs->flow, structs->cell_i_fr, structs->cell_j_fr, cell_vol, (structs->flow - cell_vol) );
03359
03360 }
03361
03362 structs->flow = cell_vol;
03363 }
03364 }
03365
03366
03367
03368 if (structs->cell_i_fr == 2) {
03369 structs->conc_P = (structs->TPser == 1) ?
03370 Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
03371 ((structs->TPser == 0) ? (structs->TP) : (0.0) );
03372 structs->conc_N = (structs->TNser == 1) ?
03373 Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
03374 ((structs->TNser == 0) ? (structs->TN) : (0.0) );
03375 structs->conc_S = (structs->TSser == 1) ?
03376 Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
03377 ((structs->TSser == 0) ? (structs->TS) : (0.0) );
03378 }
03379 else
03380 {
03381 structs->conc_P = (structs->TPser == -1) ?
03382 ( (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03383 ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
03384 structs->conc_N = ( structs->TNser == -1) ?
03385 ( (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03386 ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
03387 structs->conc_S = ( structs->TSser == -1) ?
03388 ( (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03389 ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
03390 }
03391
03392
03393
03394 P_fl = structs->flow * structs->conc_P;
03395 N_fl = structs->flow * structs->conc_N;
03396 S_fl = structs->flow * structs->conc_S;
03397
03398
03399 if (structs->cell_i_to == 2) {
03400 if (debug > -999) {
03401 VOL_OUT_STR[basn[cell_fr]] += structs->flow;
03402 VOL_OUT_STR[0] += structs->flow;
03403 P_OUT_STR[basn[cell_fr]] += P_fl;
03404 P_OUT_STR[0] += P_fl;
03405 S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0);
03406 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0);
03407 }
03408 PA[cell_fr] -= P_fl;
03409 NA[cell_fr] -= N_fl;
03410 SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03411
03412
03413
03414
03415
03416 SWH[cell_fr] -= structs->flow/CELL_SIZE;
03417 }
03418 else if (structs->cell_i_fr == 2) {
03419 if (debug > -999) {
03420 VOL_IN_STR[basn[cell_to]] += structs->flow;
03421 VOL_IN_STR[0] += structs->flow;
03422 P_IN_STR[basn[cell_to]] += P_fl;
03423 P_IN_STR[0] += P_fl;
03424 S_IN_STR[basn[cell_to]] += S_fl;
03425 S_IN_STR[0] += S_fl;
03426 }
03427 PA[cell_to] += P_fl;
03428 NA[cell_to] += N_fl;
03429 SA[cell_to] += S_fl;
03430
03431 SWH[cell_to] += structs->flow/CELL_SIZE;
03432 }
03433 else {
03434 PA[cell_fr] -= P_fl;
03435 NA[cell_fr] -= N_fl;
03436 SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03437
03438
03439
03440
03441
03442 PA[cell_to] += P_fl;
03443 NA[cell_to] += N_fl;
03444 SA[cell_to] += S_fl;
03445
03446 SWH[cell_fr] -= structs->flow/CELL_SIZE;
03447 SWH[cell_to] += structs->flow/CELL_SIZE;
03448 if (basn[cell_to] != basn[cell_fr] && debug > -999 ) {
03449 VOL_OUT_STR[basn[cell_fr]] += structs->flow;
03450 VOL_IN_STR[basn[cell_to]] += structs->flow;
03451 P_OUT_STR[basn[cell_fr]] += P_fl;
03452 P_IN_STR[basn[cell_to]] += P_fl;
03453 S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0);
03454 S_IN_STR[basn[cell_to]] += S_fl;
03455
03456
03457 if ( (basn_list[basn[cell_fr]]->family == basn_list[basn[cell_to]]->family) && (debug > 0) ) {
03458 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",
03459 SimTime.TIME, structs->cell_i_fr,structs->cell_j_fr,structs->cell_i_to,structs->cell_j_to );
03460 WriteMsg( msgStr,True ); dynERRORnum++;
03461 }
03462 }
03463
03464 }
03465 return;
03466 }
03467
03468
03481 void flowCalc_CelCel(struct Structure *structs, float *SWH, float *Elev, double *NA, double *PA, double *SA)
03482 {
03483 int cell_to, cell_fr;
03484
03485 double cell_vol;
03486 double N_fl, P_fl, S_fl;
03487 float FFlux;
03488 extern float *SF_WT_VEL_mag;
03489
03490 cell_to = T(structs->cell_i_to,structs->cell_j_to);
03491 cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);
03492
03493
03494 if ( structs->w_coef != 0.0) {
03495 printf ("\nSorry - no rule-based cell-cell flows in this model version!\n"); exit(-1);
03496 }
03497 else {
03498
03499
03500
03501
03502 FFlux = Flux_SWcells(structs->cell_i_fr,structs->cell_j_fr,
03503 structs->cell_i_to,structs->cell_j_to,SWH,Elev,MCopen);
03504
03505
03506 structs->flow = FFlux*CELL_SIZE;
03507 structs->conc_P = (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ;
03508
03509 structs->conc_S = (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ;
03510
03511 if (FFlux != 0.0) Flux_SWstuff ( structs->cell_i_fr,structs->cell_j_fr,
03512 structs->cell_i_to,structs->cell_j_to,FFlux,SWH,SA,NA,PA,SF_WT_VEL_mag);
03513
03514
03515 }
03516 }
03517
03518
03519
03520
03528 void flowData_CanCel(struct Structure *structs, struct Structure *struct_hist_start, float *SWH, double *NA, double *PA, double *SA)
03529 {
03530 int cell_to, can_to, can_fr;
03531 int cell_struct;
03532 double CH_vol;
03533 double Nconc, Pconc, Sconc;
03534 double N_fl, P_fl, S_fl;
03535 double chan_over;
03536 double ChanHistOut;
03537
03538 can_fr = structs->canal_fr - CHo;
03539 cell_to = T(structs->cell_i_to,structs->cell_j_to);
03540 cell_struct = T(structs->str_cell_i,structs->str_cell_j);
03541
03542
03543
03544 if (structs->multiOut == 1) return;
03545
03546 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03547
03548 if (structs->flow == 0.0) {
03549 structs->conc_P = 0.0;
03550 structs->conc_N = 0.0;
03551 structs->conc_S = 0.0;
03552 return;
03553 }
03554 else if (structs->flow < 0.0) {
03555 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from cell (%d,%d)",
03556 SimTime.TIME, structs->cell_i_to, structs->cell_j_to );
03557 WriteMsg( msgStr,True ); dynERRORnum++;
03558 return;
03559 }
03560
03561
03562
03563
03564 struct_hist_start = structs;
03565 ChanHistOut = 0.0;
03566
03567 while ( structs != NULL )
03568 {
03569 if ( structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
03570 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03571 ChanHistOut += structs->flow;
03572 }
03573
03574 structs = structs->next_in_list;
03575 }
03576
03577 structs = struct_hist_start;
03578
03579
03580
03581
03582
03583 CH_vol = Chan_list[can_fr]->area*Chan_list[can_fr]->wat_depth;
03584
03585 Pconc = (structs->TPser == -1) ?
03586 ( (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0) ) :
03587 ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) ) ;
03588 Nconc = (structs->TNser == -1) ?
03589 ( (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0) ) :
03590 ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) ) ;
03591 Sconc = (structs->TSser == -1) ?
03592 ( (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0) ) :
03593 ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) ) ;
03594
03595
03596 CH_vol = Chan_list[can_fr]->area*(ramp(Chan_list[can_fr]->wat_depth-MINDEPTH) );
03597 chan_over = (ChanHistOut > 0.0 ) ? ( CH_vol/ChanHistOut ) : 1.0;
03598
03599 do
03600 {
03601 if ( structs->flag == 1 && struct_hist_start->canal_fr == structs->canal_fr ) {
03602
03603 structs->multiOut = 1;
03604
03605
03606 if ( chan_over < 1.0 ) {
03607 if( debug > 0 ) {
03608 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",
03609 SimTime.TIME, structs->S_nam, ChanHistOut, Chan_list[can_fr]->number,
03610 CH_vol, (ChanHistOut - CH_vol), Chan_list[can_fr]->wat_depth );
03611
03612 }
03613
03614 structs->flow = (structs->flow) * chan_over ;
03615 }
03616
03617
03618 Chan_list[can_fr]->sumHistOut += structs->flow;
03619
03620
03621 can_to = structs->canal_to - CHo;
03622 cell_to = T(structs->cell_i_to,structs->cell_j_to);
03623
03624 structs->conc_P = Pconc;
03625 structs->conc_N = Nconc;
03626 structs->conc_S = (structs->TSser == -1) ?
03627 (Sconc) :
03628 ( (structs->TSser == 0) ? (structs->TS) :
03629 ( Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
03630
03631 P_fl = structs->flow * structs->conc_P;
03632 N_fl = structs->flow * structs->conc_N;
03633 S_fl = structs->flow * structs->conc_S;
03634
03635 Chan_list[can_fr]->N -= N_fl ;
03636 Chan_list[can_fr]->P -= P_fl ;
03637 Chan_list[can_fr]->S -= (structs->TSser == -1) ? (S_fl) : (0);
03638
03639
03640
03641
03642
03643 if (can_to > 0) {
03644
03645
03646 Chan_list[can_to]->sumHistIn += structs->flow;
03647
03648 Chan_list[can_to]->N += N_fl ;
03649 Chan_list[can_to]->P += P_fl ;
03650 Chan_list[can_to]->S += S_fl ;
03651
03652 if (Chan_list[can_fr]->basin != Chan_list[can_to]->basin && debug > -999 ) {
03653 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03654 VOL_IN_STR[Chan_list [can_to]->basin] += structs->flow;
03655 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03656 P_IN_STR[Chan_list [can_to]->basin] += P_fl;
03657 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0);
03658 S_IN_STR[Chan_list [can_to]->basin] += S_fl;
03659 }
03660 }
03661 else if (cell_to > 0) {
03662
03663
03664 if (structs->cell_i_to > 2 ) {
03665 SWH[cell_to] += structs->flow/CELL_SIZE;
03666 PA[cell_to] += P_fl;
03667 NA[cell_to] += N_fl;
03668 SA[cell_to] += S_fl;
03669
03670 if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
03671 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03672 VOL_IN_STR[basn[cell_to]] += structs->flow;
03673 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03674 P_IN_STR[basn[cell_to]] += P_fl;
03675 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0);
03676 S_IN_STR[basn[cell_to]] += S_fl;
03677
03678
03679
03680
03681
03682
03683
03684 }
03685 }
03686 else if (structs->cell_i_to == 2 && debug > -999 ) {
03687 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03688 VOL_OUT_STR[0] += structs->flow;
03689 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03690 P_OUT_STR[0] += P_fl;
03691 S_OUT_STR[Chan_list[can_fr]->basin] += (structs->TSser == -1) ? (S_fl) : (0);
03692 S_OUT_STR[0] += (structs->TSser == -1) ? (S_fl) : (0);
03693 }
03694 }
03695
03696 if (structs->flow == 0.0) {
03697 structs->conc_P = 0.0;
03698 structs->conc_N = 0.0;
03699 structs->conc_S = 0.0;
03700 }
03701 else if (structs->flow < 0.0) {
03702 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from another chan/cell into chan %d",
03703 SimTime.TIME, Chan_list[can_fr]->number );
03704 WriteMsg( msgStr,True ); dynERRORnum++;
03705 }
03706 }
03707
03708 structs = structs->next_in_list;
03709 } while (structs != NULL);
03710
03711 structs = struct_hist_start;
03712
03713
03714 return;
03715
03716 }
03717
03718
03719
03747 void flowCalc_CanCel(struct Structure *structs, float *SWH, float *GW,
03748 float *poros, float *Elev, double *NA, double *PA, double *SA)
03749 {
03750 int cell_to, can_fr;
03751 int cell_struct;
03752 float HeadH, HeadT;
03753 double CH_vol;
03754 double N_fl, P_fl, S_fl;
03755 float flow_rev;
03756
03757
03758 extern float* boundcond_depth;
03759
03760 float HWcheck, TWcheck, Hdiff;
03761 int struct_open;
03762
03763
03764 can_fr = structs->canal_fr - CHo;
03765 cell_to = T(structs->cell_i_to,structs->cell_j_to);
03766 cell_struct = T(structs->str_cell_i,structs->str_cell_j);
03767
03768 if (structs->cell_i_to != 2 ) {
03769 sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow is restricted to outflow from the model domain.",
03770 structs->S_nam );
03771 usrErr( msgStr ); exit(-1);
03772 }
03773 if (structs->flagTWsched == 1 && structs->flagHWsched == 1 ) {
03774 sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow can not have both HeadWater and TailWater Schedules.",
03775 structs->S_nam );
03776 usrErr( msgStr); exit(-1);
03777 }
03778 if ( structs->flagHWsched != 1 && structs->flagTWsched != 1 && structs->cell_i_TW == 0 ) {
03779 sprintf( msgStr, "ERROR in flow configuration for %s structure : without a head or tail water schedule, this rule-based flow must have a tailwater-check-cell defined (to acquire data from the boundary-condition model).",
03780 structs->S_nam );
03781 usrErr( msgStr ); exit(-1);
03782 }
03783
03784
03785
03786 HeadH = Chan_list[can_fr]->elev_avg - Chan_list[can_fr]->depth + Chan_list[can_fr]->wat_depth
03787 + (Chan_list[can_fr]->sumHistIn - Chan_list[can_fr]->sumHistOut)/Chan_list[can_fr]->area
03788 + (Chan_list[can_fr]->sumRuleIn - Chan_list[can_fr]->sumRuleOut)/Chan_list[can_fr]->area;
03789
03790 HeadT = (structs->flagTWsched == 1) ? (structs->TW_stage) : ( boundcond_depth[T(structs->cell_i_TW,structs->cell_j_TW)] + Elev[cell_struct]);
03791
03792 Hdiff = HeadH - HeadT;
03793
03794 CH_vol = Chan_list[can_fr]->area* ramp(HeadH-Chan_list[can_fr]->elev_avg+Chan_list[can_fr]->depth-MINDEPTH);
03795
03796 struct_open = 1;
03797
03798 if (structs->flagHWsched == 1 && structs->cell_HW != 0) {
03799
03800
03801
03802 HWcheck = GW[structs->cell_HW]/poros[HAB[structs->cell_HW]] + SWH[structs->cell_HW] ;
03803 struct_open = (structs->HW_stage < HWcheck) ? (1) : (0);
03804 Hdiff = HWcheck - structs->HW_stage ;
03805 }
03806
03807
03808 if ( (Hdiff > 0.0) && struct_open == 1)
03809 {
03810 structs->flow = Min(structs->w_coef * ( Hdiff ) * canstep, CH_vol );
03811
03812
03813
03814 flow_rev = ( structs->cell_i_to != 2 ) ?
03815 (HeadH - structs->flow/Chan_list[can_fr]->area) - (HeadT + structs->flow/CELL_SIZE) :
03816 (0.0);
03817
03818
03819 if (flow_rev < 0.0) { structs->flow =
03820 ramp( Chan_list[can_fr]->area*CELL_SIZE*HeadH/(Chan_list[can_fr]->area+CELL_SIZE) -
03821 Chan_list[can_fr]->area*CELL_SIZE*HeadT/(Chan_list[can_fr]->area+CELL_SIZE) );
03822 }
03823 }
03824 else structs->flow = 0.0;
03825
03826
03827
03828
03829 if (structs->flow == 0.0) {
03830 structs->conc_P = 0.0;
03831 structs->conc_N = 0.0;
03832 structs->conc_S = 0.0;
03833 return;
03834 }
03835 else if (structs->flow < 0) {
03836 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG FLOW from cell (%d,%d)",
03837 SimTime.TIME, structs->cell_i_to, structs->cell_j_to );
03838 WriteMsg( msgStr,True ); dynERRORnum++;
03839 return;
03840 }
03841
03842
03843
03844 CH_vol = Chan_list[can_fr]->area* ramp(HeadH-Elev[cell_struct]+Chan_list[can_fr]->depth);
03845
03846 structs->conc_P = (CH_vol>0) ? (Chan_list[can_fr]->P/CH_vol) : (0.0);
03847 structs->conc_N = (CH_vol>0) ? (Chan_list[can_fr]->N/CH_vol) : (0.0);
03848 structs->conc_S = (CH_vol>0) ? (Chan_list[can_fr]->S/CH_vol) : (0.0);
03849
03850 P_fl = structs->flow * structs->conc_P;
03851 N_fl = structs->flow * structs->conc_N;
03852 S_fl = structs->flow * structs->conc_S;
03853
03854
03855
03856 Chan_list[can_fr]->sumRuleOut += structs->flow;
03857
03858
03859 if (structs->cell_i_to > 2 ) {
03860 SWH[cell_to] += structs->flow/CELL_SIZE;
03861 PA[cell_to] += P_fl;
03862 NA[cell_to] += N_fl;
03863 SA[cell_to] += S_fl;
03864
03865 if (Chan_list[can_fr]->basin != basn[cell_to] && debug > -999 ) {
03866 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03867 VOL_IN_STR[basn[cell_to]] += structs->flow;
03868 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03869 P_IN_STR[basn[cell_to]] += P_fl;
03870 S_OUT_STR[Chan_list[can_fr]->basin] += S_fl;
03871 S_IN_STR[basn[cell_to]] += S_fl;
03872
03873
03874
03875
03876
03877
03878
03879 }
03880 }
03881 else if (structs->cell_i_to == 2 && debug > -999 ) {
03882 VOL_OUT_STR[Chan_list[can_fr]->basin] += structs->flow;
03883 VOL_OUT_STR[0] += structs->flow;
03884 P_OUT_STR[Chan_list[can_fr]->basin] += P_fl;
03885 P_OUT_STR[0] += P_fl;
03886 S_OUT_STR[Chan_list[can_fr]->basin] += S_fl;
03887 S_OUT_STR[0] += S_fl;
03888 }
03889
03890
03891
03892 Chan_list[can_fr]->N -= N_fl ;
03893 Chan_list[can_fr]->P -= P_fl ;
03894 Chan_list[can_fr]->S -= S_fl ;
03895 }
03896
03897
03898
03899
03907 void flowData_CelCan( struct Structure *structs, float *SWH, double *NA, double *PA, double *SA)
03908 {
03909 int cell_fr, can_to;
03910 int cell_struct;
03911 double cell_vol;
03912 double N_fl, P_fl, S_fl;
03913
03914 can_to = structs->canal_to - CHo;
03915 cell_fr = T(structs->cell_i_fr,structs->cell_j_fr);
03916 cell_struct = T(structs->str_cell_i,structs->str_cell_j);
03917
03918
03919
03920 if (structs->multiOut == 1) {
03921 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
03922 SimTime.TIME, Chan_list[can_to]->number );
03923 WriteMsg( msgStr,True ); dynERRORnum++;
03924 return;
03925 }
03926
03927 structs->flow = Hist_data[structs->histID-1].arrayWat[((int)(SimTime.TIME))]*canstep;
03928
03929 if (structs->flow == 0.0) {
03930 structs->conc_P = 0.0;
03931 structs->conc_N = 0.0;
03932 structs->conc_S = 0.0;
03933 return;
03934 }
03935 else if (structs->flow <0) {
03936 sprintf( msgStr, "Day %6.1f: ERROR IN ASSUMPTION OF NO NEG HIST FLOW from chan %d",
03937 SimTime.TIME, Chan_list[can_to]->number );
03938 WriteMsg( msgStr,True ); dynERRORnum++;
03939 return;
03940 }
03941
03942 if (structs->cell_i_fr != 2 )
03943 {
03944
03945 cell_vol = SWH[cell_fr] * CELL_SIZE;
03946 if (structs->flow > cell_vol) {
03947 if( debug > 0 ) {
03948 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",
03949 SimTime.TIME, structs->S_nam, structs->flow, structs->cell_i_fr, structs->cell_j_fr, cell_vol, (structs->flow - cell_vol) );
03950
03951 }
03952 structs->flow = cell_vol;
03953 }}
03954
03955
03956 if (structs->cell_i_fr == 2) {
03957 structs->conc_P = (structs->TPser == 1) ?
03958 Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
03959 ((structs->TPser == 0) ? (structs->TP) : (0.0) );
03960 structs->conc_N = (structs->TNser == 1) ?
03961 Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
03962 ((structs->TNser == 0) ? (structs->TN) : (0.0) );
03963 structs->conc_S = (structs->TSser == 1) ?
03964 Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
03965 ((structs->TSser == 0) ? (structs->TS) : (0.0) );
03966 }
03967 else
03968 {
03969 structs->conc_P = (structs->TPser == -1) ?
03970 ( (SWH[cell_fr] > 0) ? PA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03971 ( (structs->TPser == 0) ? (structs->TP) : (Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] ) );
03972 structs->conc_N = ( structs->TNser == -1) ?
03973 ( (SWH[cell_fr] > 0) ? NA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03974 ( (structs->TNser == 0) ? (structs->TN) : (Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] ) );
03975 structs->conc_S = ( structs->TSser == -1) ?
03976 ( (SWH[cell_fr] > 0) ? SA[cell_fr]/CELL_SIZE/SWH[cell_fr] : 0.0 ) :
03977 ( (structs->TSser == 0) ? (structs->TS) : (Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] ) );
03978 }
03979
03980
03981 P_fl = structs->flow * structs->conc_P;
03982 N_fl = structs->flow * structs->conc_N;
03983 S_fl = structs->flow * structs->conc_S;
03984 if (structs->cell_i_fr == 2 ) {
03985 if (debug > -999) {
03986 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow;
03987 VOL_IN_STR[0] += structs->flow;
03988 P_IN_STR[Chan_list[can_to]->basin] += P_fl;
03989 P_IN_STR[0] += P_fl;
03990 S_IN_STR[Chan_list[can_to]->basin] += S_fl;
03991 S_IN_STR[0] += S_fl;
03992 } }
03993 else {
03994 SWH[cell_fr] -= structs->flow/CELL_SIZE;
03995 PA[cell_fr] -= P_fl;
03996 NA[cell_fr] -= N_fl;
03997 SA[cell_fr] -= (structs->TSser == -1) ? (S_fl) : (0);
03998
03999
04000
04001
04002
04003 if (basn[cell_fr] != Chan_list[can_to]->basin && debug > -999 ) {
04004 VOL_OUT_STR[basn[cell_fr]] += structs->flow;
04005 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow;
04006 P_OUT_STR[basn[cell_fr]] += P_fl;
04007 P_IN_STR[Chan_list[can_to]->basin] += P_fl;
04008 S_OUT_STR[basn[cell_fr]] += (structs->TSser == -1) ? (S_fl) : (0);
04009 S_IN_STR[Chan_list[can_to]->basin] += S_fl;
04010
04011
04012 if ( (Chan_list[can_to]->family == basn_list[basn[cell_fr]]->family) && (debug > 0) ) {
04013 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",
04014 SimTime.TIME, structs->cell_i_fr, structs->cell_j_fr, Chan_list[can_to]->number );
04015 WriteMsg( msgStr,True ); dynERRORnum++;
04016 }
04017 }
04018 }
04019
04020
04021
04022 Chan_list[can_to]->sumHistIn += structs->flow;
04023
04024
04025 Chan_list[can_to]->N += N_fl ;
04026 Chan_list[can_to]->P += P_fl ;
04027 Chan_list[can_to]->S += S_fl ;
04028 return;
04029
04030 }
04031
04032
04033
04061 void flowCalc_CelCan( struct Structure *structs, float *SWH, float *GW,
04062 float *poros, float *Elev, double *NA, double *PA, double *SA )
04063 {
04064 int cell_fr, can_to;
04065 int cell_struct;
04066 float HeadH, HeadT;
04067 double N_fl, P_fl, S_fl;
04068
04069
04070 extern float* boundcond_depth;
04071
04072 float HWcheck, TWcheck, Hdiff;
04073 int struct_open;
04074
04075
04076 can_to = structs->canal_to - CHo;
04077
04078 if (structs->cell_i_fr != 2 ) {
04079 sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow is restricted to inflow from outside the model domain.",
04080 structs->S_nam );
04081 usrErr( msgStr); exit(-1);
04082 }
04083 if (structs->flagTWsched == 1 && structs->flagHWsched == 1 ) {
04084 sprintf( msgStr, "ERROR in flow configuration for %s structure : this rule-based flow can not have both HeadWater and TailWater Schedules.",
04085 structs->S_nam );
04086 usrErr( msgStr); exit(-1);
04087 }
04088 if ( structs->flagTWsched != 1 && structs->flagHWsched != 1 && structs->cell_i_HW == 0 ) {
04089 sprintf( msgStr, "ERROR in flow configuration for %s structure : without a head or tail water schedule, this rule-based flow must have a headwater-check-cell defined (to acquire data from the boundary-condition model).",
04090 structs->S_nam );
04091 usrErr( msgStr ); exit(-1);
04092 }
04093
04094 HeadH = (structs->flagHWsched == 1) ? (structs->HW_stage) : ( boundcond_depth[T(structs->cell_i_HW,structs->cell_j_HW)] + Elev[cell_struct]);
04095
04096
04097 HeadT = Chan_list[can_to]->elev_avg - Chan_list[can_to]->depth + Chan_list[can_to]->wat_depth
04098 + (Chan_list[can_to]->sumHistIn - Chan_list[can_to]->sumHistOut)/Chan_list[can_to]->area
04099 + (Chan_list[can_to]->sumRuleIn - Chan_list[can_to]->sumRuleOut)/Chan_list[can_to]->area;
04100
04101 Hdiff = HeadH - HeadT;
04102
04103 struct_open = 1;
04104
04105 if (structs->flagTWsched == 1 && structs->cell_TW != 0) {
04106
04107
04108
04109 TWcheck = GW[structs->cell_TW]/poros[HAB[structs->cell_TW]] + SWH[structs->cell_TW] ;
04110 struct_open = (structs->TW_stage > TWcheck) ? (1) : (0);
04111 Hdiff = structs->TW_stage - TWcheck;
04112 }
04113
04114
04115 if ( struct_open == 1)
04116 {
04117 structs->flow = Max(structs->w_coef * ( Hdiff ) * canstep, 0 );
04118
04119 }
04120 else structs->flow = 0.0;
04121
04122
04123 if (structs->flow == 0.0) {
04124 structs->conc_P = 0.0;
04125 structs->conc_N = 0.0;
04126 structs->conc_S = 0.0;
04127 return;
04128 }
04129
04130
04131
04132
04133 if (structs->cell_i_fr == 2) {
04134 structs->conc_P = (structs->TPser == 1) ?
04135 Hist_data[structs->histID-1].arrayP[((int)(SimTime.TIME))] :
04136 ((structs->TPser == 0) ? (structs->TP) : (0.0) );
04137 structs->conc_N = (structs->TNser == 1) ?
04138 Hist_data[structs->histID-1].arrayN[((int)(SimTime.TIME))] :
04139 ((structs->TNser == 0) ? (structs->TN) : (0.0) );
04140 structs->conc_S = (structs->TSser == 1) ?
04141 Hist_data[structs->histID-1].arrayS[((int)(SimTime.TIME))] :
04142 ((structs->TSser == 0) ? (structs->TS) : (0.0) );
04143 }
04144 P_fl = structs->flow * structs->conc_P;
04145 N_fl = structs->flow * structs->conc_N;
04146 S_fl = structs->flow * structs->conc_S;
04147
04148
04149
04150 Chan_list[can_to]->sumRuleIn += structs->flow;
04151
04152
04153
04154
04155 if (structs->cell_i_fr == 2 ) {
04156 VOL_IN_STR[Chan_list[can_to]->basin] += structs->flow;
04157 VOL_IN_STR[0] += structs->flow;
04158 P_IN_STR[Chan_list[can_to]->basin] += P_fl;
04159 P_IN_STR[0] += P_fl;
04160 S_IN_STR[Chan_list[can_to]->basin] += S_fl;
04161 S_IN_STR[0] += S_fl;
04162 }
04163 else {
04164 sprintf( msgStr, "Day %6.1f: ERROR in flow to chan %d: cannot have calculated cell->canal flows except from outside of model domain. ",
04165 SimTime.TIME, Chan_list[can_to]->number);
04166 WriteMsg( msgStr,True ); dynERRORnum++;
04167 return;
04168 }
04169
04170
04171
04172 Chan_list[can_to]->N += N_fl ;
04173 Chan_list[can_to]->P += P_fl ;
04174 Chan_list[can_to]->S += S_fl ;
04175
04176 }
04177
04178
04179
04180
04181
04190 float GetGraph (struct Schedule *graph, float x)
04191 {
04192 float s;
04193 int ig = 0, Np;
04194
04195 Np = graph->num_point;
04196
04197 while(1)
04198 {
04199 if (x <= graph->graph_points[ig].time)
04200 { if ( ig > 0 ) ig--;
04201 else return ( graph->graph_points[0].value );
04202 }
04203 else if (x > graph->graph_points[ig].time && x <= graph->graph_points[ig+1].time)
04204 {
04205 s = ( graph->graph_points[ig+1].value - graph->graph_points[ig].value )/
04206 ( graph->graph_points[ig+1].time - graph->graph_points[ig].time );
04207 return ( s * ( x - graph->graph_points[ig].time ) + graph->graph_points[ig].value );
04208 }
04209 else if (x > graph->graph_points[ig+1].time)
04210 { if ( ig < Np ) ig++;
04211 else return ( graph->graph_points[Np].value );
04212 }
04213 }
04214 }
04215
04216
04217
04223 float UTM2kmy ( double UTM )
04224 { return ( Abs(UTM - UTM_EOrigin)/celWid );
04225 }
04226
04232 float UTM2kmx ( double UTM )
04233 { return ( Abs(UTM - UTM_NOrigin)/celWid );
04234 }
04235
04236
04243 int Wrdcmp ( char *s, char *t )
04244 {
04245 for ( ; *s == *t; s++, t++ )
04246 if ( isspace (*(s+1)) || *(s+1) == EOS ) return 1;
04247
04248 return 0;
04249
04250 }
04251
04252