00001
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "fluxes.h"
00044
00045
00046
00107 void Flux_SWater(int it, float *SURFACE_WAT,float *SED_ELEV,float *HYD_MANNINGS_N,
00108 double *STUF1, double *STUF2, double *STUF3)
00109 {
00115 int ix, iy;
00116 float FFlux;
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 for(ix=1; ix<=s0; ix++)
00128 {
00129 if (it%2)
00130 {
00131 for(iy=1; iy<=s1; iy++)
00132 {
00133 if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy+1)] && (int)(ON_MAP[T(ix,iy)]-1) & 1 ) ||
00134 BCondFlow[T(ix,iy+1)] == 3 || BCondFlow[T(ix,iy)] == 3 )
00135 {
00136 FFlux = Flux_SWcells(ix,iy,ix,iy+1,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N);
00137
00138 if (FFlux != 0.0)
00139 Flux_SWstuff ( ix,iy,ix,iy+1,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00140
00141 }
00142 }
00143 }
00144
00145
00146 else
00147 {
00148 for(iy=s1; iy>=1; iy--)
00149 if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy-1)] && (int)(ON_MAP[T(ix,iy-1)]-1) & 1 ) ||
00150 BCondFlow[T(ix,iy-1)] == 3 || BCondFlow[T(ix,iy)] == 3 )
00151 {
00152 FFlux = Flux_SWcells(ix,iy-1,ix,iy,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N);
00153
00154 if (FFlux != 0.0)
00155 Flux_SWstuff ( ix,iy-1,ix,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00156
00157 }
00158 }
00159 }
00160
00161 for(iy=1; iy<=s1; iy++)
00162 {
00163 if (it%2)
00164 {
00165 for(ix=1; ix<=s0; ix++)
00166 if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix+1,iy)] && (int)(ON_MAP[T(ix,iy)]-1) & 2 ) ||
00167 BCondFlow[T(ix+1,iy)] == 3 || BCondFlow[T(ix,iy)] == 3 )
00168 {
00169 FFlux = Flux_SWcells(ix,iy,ix+1,iy,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N);
00170
00171 if (FFlux != 0.0)
00172 Flux_SWstuff ( ix,iy,ix+1,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00173 }
00174 }
00175
00176
00177 else
00178 {
00179 for(ix=s0; ix>=1; ix--)
00180 if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix-1,iy)] && (int)(ON_MAP[T(ix-1,iy)]-1) & 2 ) ||
00181 BCondFlow[T(ix-1,iy)] == 3 || BCondFlow[T(ix,iy)] == 3 )
00182 {
00183
00184 FFlux = Flux_SWcells(ix-1,iy,ix,iy,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N);
00185
00186 if (FFlux != 0.0)
00187 Flux_SWstuff ( ix-1,iy,ix,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00188
00189 }
00190 }
00191 }
00192
00193 }
00194
00195
00196
00197
00198
00217 float Flux_SWcells(int i0,int i1,int j0,int j1, float *SWater,float *Elevation,float *MC)
00218 {
00225 float dh, adh;
00226 float MC_cells;
00227 float Hi, Hj;
00228 float Flux = 0.;
00229 int cellLoci = T(i0,i1);
00230 int cellLocj = T(j0,j1);
00231
00232
00233 MC_cells = (MC[cellLoci] + MC[T(j0,j1)])/2.;
00234
00235
00236 if (!ON_MAP[cellLoci] && BCondFlow[cellLocj] == 3 )
00237 {
00238
00239
00240
00241
00242 Hi = boundcond_depth[cellLoci];
00243
00244 MC_cells = MC[cellLocj];
00245 }
00246
00247 else
00248 Hi = SWater[cellLoci] + Elevation[cellLoci];
00249
00250 if(BCondFlow[cellLoci] == 3 && !ON_MAP[cellLocj] )
00251 {
00252
00253
00254
00255 Hj = boundcond_depth[cellLocj];
00256
00257 MC_cells = MC[cellLoci];
00258 }
00259
00260 else
00261 Hj = SWater[cellLocj] + Elevation[cellLocj];
00262
00263 dh = Hi - Hj;
00264 adh = Abs (dh);
00265
00266 if (dh > 0)
00267 {
00268 if(SWater[cellLoci] < GP_DetentZ)
00269 return 0.0;
00270
00271 Flux = (MC_cells != 0) ?
00272 (pow(adh,GP_mannHeadPow) / MC_cells * pow(SWater[cellLoci],GP_mannDepthPow)*step_Cell) : (0.0);
00273
00274
00275 Flux = ( Flux > ramp(SWater[cellLoci] - GP_DetentZ) ) ? (ramp(SWater[cellLoci] - GP_DetentZ)) : (Flux);
00276
00277
00278 if ( ( Hi - Flux ) < ( Hj + Flux ) )
00279 Flux = Min ( dh/2.0, ramp(SWater[cellLoci] - GP_DetentZ) );
00280
00281 }
00282
00283 else
00284 {
00285 if (SWater[cellLocj] < GP_DetentZ)
00286 return 0.0;
00287
00288
00289 Flux = (MC_cells != 0) ?
00290 ( - pow(adh,GP_mannHeadPow) / MC_cells * pow(SWater[cellLocj],GP_mannDepthPow)*step_Cell) : (0.0);
00291
00292
00293 Flux = ( -Flux > ramp(SWater[cellLocj] - GP_DetentZ) ) ? (-ramp(SWater[cellLocj] - GP_DetentZ)) : (Flux);
00294
00295
00296 if ( ( Hi - Flux ) > ( Hj + Flux ) )
00297 Flux = - Min ( adh/2.0, ramp(SWater[cellLocj] - GP_DetentZ) );
00298 }
00299
00300 return (Flux);
00301 }
00302
00303
00304
00324 void Flux_SWstuff(int i0,int i1,int j0,int j1, float Flux, float *SURFACE_WAT, double *STUF1, double *STUF2, double *STUF3)
00325 {
00326 float m1=0.0, m2=0.0, m3=0.0;
00327 int ii, flo_chek, cel_i, cel_j;
00328 float disp_num;
00329 float veloc;
00330 float velocAdj;
00331 float FluxAdj;
00332 double fl_prop_i, fl_prop_j;
00333 extern float dispParm_scaled;
00334
00335 cel_i = T(i0,i1);
00336 cel_j = T(j0,j1);
00337
00338
00339
00340
00341
00342
00343 FluxAdj = Disp_Calc(Flux, SURFACE_WAT[cel_i], SURFACE_WAT[cel_j], sfstep);
00344
00345 fl_prop_i = (SURFACE_WAT[cel_i]>0.0) ? (Max(Flux-FluxAdj,0.0)/SURFACE_WAT[cel_i]) : (0.0);
00346 fl_prop_j = (SURFACE_WAT[cel_j]>0.0) ? (Min(Flux+FluxAdj,0.0)/SURFACE_WAT[cel_j]) : (0.0);
00347
00348 fl_prop_i = Min(fl_prop_i, 1.0);
00349 fl_prop_j = Min(fl_prop_j, 1.0);
00350
00351 if (Flux >0.0) {
00352 m1 = STUF1[cel_i]*fl_prop_i;
00353 m2 = STUF2[cel_i]*fl_prop_i;
00354 m3 = STUF3[cel_i]*fl_prop_i;
00355 }
00356 else {
00357 m1 = STUF1[cel_j]*fl_prop_j;
00358 m2 = STUF2[cel_j]*fl_prop_j;
00359 m3 = STUF3[cel_j]*fl_prop_j;
00360 }
00361
00362 STUF1[cel_j] += m1;
00363 STUF2[cel_j] += m2;
00364 STUF3[cel_j] += m3;
00365 STUF1[cel_i] -= m1;
00366 STUF2[cel_i] -= m2;
00367 STUF3[cel_i] -= m3;
00368 SURFACE_WAT[cel_j] += Flux;
00369 SURFACE_WAT[cel_i] -= Flux;
00370
00371 if (debug > 2) {
00372 if (STUF3[cel_j] < -GP_MinCheck) {
00373 sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water P (%f kg) in cell (%d,%d) after SWfluxes!",
00374 SimTime.TIME, STUF3[cel_j], j0,j1 );
00375 WriteMsg( msgStr,True ); dynERRORnum++;}
00376 if (STUF3[cel_i] < -GP_MinCheck) {
00377 sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water P (%f kg) in cell (%d,%d) after SWfluxes!",
00378 SimTime.TIME, STUF3[cel_i], i0,i1 );
00379 WriteMsg( msgStr,True ); dynERRORnum++; }
00380 if (STUF1[cel_j] < -GP_MinCheck) {
00381 sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water S (%f kg) in cell (%d,%d) after SWfluxes!",
00382 SimTime.TIME, STUF1[cel_j], j0,j1 );
00383 WriteMsg( msgStr,True ); dynERRORnum++; }
00384 if (STUF1[cel_i] < -GP_MinCheck) {
00385 sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water S (%f kg) in cell (%d,%d) after SWfluxes!",
00386 SimTime.TIME, STUF1[cel_i], i0,i1 );
00387 WriteMsg( msgStr,True ); dynERRORnum++; }
00388 if (SURFACE_WAT[cel_j] < -GP_MinCheck) {
00389 sprintf(msgStr,"Day %6.1f: capacityERR - negative surface water (%f m) in cell (%d,%d)!",
00390 SimTime.TIME, SURFACE_WAT[cel_j], j0,j1 ) ;
00391 WriteMsg( msgStr,True ); dynERRORnum++;}
00392
00393 if (SURFACE_WAT[cel_i] < -GP_MinCheck) {
00394 sprintf(msgStr,"Day %6.1f: capacityERR - negative surface water (%f m) in cell (%d,%d)!",
00395 SimTime.TIME, SURFACE_WAT[cel_i], i0,i1 ) ;
00396 WriteMsg( msgStr,True ); dynERRORnum++; }
00397 }
00398
00399
00400
00401 if (basn[cel_i] != basn[cel_j]) {
00402
00403
00404 if ( ON_MAP[cel_j] && ON_MAP[cel_i] ) {
00405 if ( Flux > 0 ) {
00406 if (basn_list[basn[cel_i]]->family !=
00407 basn_list[basn[cel_j]]->family ) {
00408 if ( !basn_list[basn[cel_i]]->parent ) {
00409
00410 VOL_OUT_OVL[basn_list[basn[cel_i]]->family] += Flux*CELL_SIZE;
00411 P_OUT_OVL[basn_list[basn[cel_i]]->family] += m3;
00412 S_OUT_OVL[basn_list[basn[cel_i]]->family] += m1;
00413 }
00414 if ( !basn_list[basn[cel_j]]->parent ) {
00415
00416 VOL_IN_OVL[basn_list[basn[cel_j]]->family] += Flux*CELL_SIZE;
00417 P_IN_OVL[basn_list[basn[cel_j]]->family] += m3;
00418 S_IN_OVL[basn_list[basn[cel_j]]->family] += m1;
00419 }
00420
00421 VOL_OUT_OVL[basn[cel_i]] += Flux*CELL_SIZE;
00422 P_OUT_OVL[basn[cel_i]] += m3;
00423 S_OUT_OVL[basn[cel_i]] += m1;
00424 VOL_IN_OVL[basn[cel_j]] += Flux*CELL_SIZE;
00425 P_IN_OVL[basn[cel_j]] += m3;
00426 S_IN_OVL[basn[cel_j]] += m1;
00427
00428 }
00429 else {
00430
00431 if ( !basn_list[basn[cel_i]]->parent ) {
00432 VOL_OUT_OVL[basn[cel_i]] += Flux*CELL_SIZE;
00433 P_OUT_OVL[basn[cel_i]] += m3;
00434 S_OUT_OVL[basn[cel_i]] += m1;
00435 }
00436 if ( !basn_list[basn[cel_j]]->parent ) {
00437 VOL_IN_OVL[basn[cel_j]] += Flux*CELL_SIZE;
00438 P_IN_OVL[basn[cel_j]] += m3;
00439 S_IN_OVL[basn[cel_j]] += m1;
00440 }
00441 }
00442
00443 if (debug > 0 && WatMgmtOn) {
00444 basins = basn_list[basn[cel_i]];
00445 flo_chek = 0;
00446 for (ii=0; ii<basins->numFLok; ii++) { if (basn[cel_j] == basins->FLok[ii] ) flo_chek = 1; }
00447 if (!flo_chek) {
00448 sprintf(msgStr,"Day %5.3f: ERROR - no (pos) SW flow from cell (%d,%d) of basin %d into cell (%d,%d) of basin %d!",
00449 SimTime.TIME, i0,i1,basn[cel_i], j0,j1, basn[cel_j]);
00450 WriteMsg( msgStr,True ); dynERRORnum++; }
00451 }
00452
00453
00454 }
00455 else {
00456 if (basn_list[basn[cel_i]]->family !=
00457 basn_list[basn[cel_j]]->family ) {
00458 if ( !basn_list[basn[cel_j]]->parent ) {
00459
00460 VOL_OUT_OVL[basn_list[basn[cel_j]]->family] -= Flux*CELL_SIZE;
00461 P_OUT_OVL[basn_list[basn[cel_j]]->family] -= m3;
00462 S_OUT_OVL[basn_list[basn[cel_j]]->family] -= m1;
00463 }
00464 if ( !basn_list[basn[cel_i]]->parent ) {
00465
00466 VOL_IN_OVL[basn_list[basn[cel_i]]->family] -= Flux*CELL_SIZE;
00467 P_IN_OVL[basn_list[basn[cel_i]]->family] -= m3;
00468 S_IN_OVL[basn_list[basn[cel_i]]->family] -= m1;
00469 }
00470
00471 VOL_OUT_OVL[basn[cel_j]] -= Flux*CELL_SIZE;
00472 P_OUT_OVL[basn[cel_j]] -= m3;
00473 S_OUT_OVL[basn[cel_j]] -= m1;
00474 VOL_IN_OVL[basn[cel_i]] -= Flux*CELL_SIZE;
00475 P_IN_OVL[basn[cel_i]] -= m3;
00476 S_IN_OVL[basn[cel_i]] -= m1;
00477
00478 }
00479 else {
00480
00481 if ( !basn_list[basn[cel_j]]->parent ) {
00482 VOL_OUT_OVL[basn[cel_j]] -= Flux*CELL_SIZE;
00483 P_OUT_OVL[basn[cel_j]] -= m3;
00484 S_OUT_OVL[basn[cel_j]] -= m1;
00485 }
00486 if ( !basn_list[basn[cel_i]]->parent ) {
00487 VOL_IN_OVL[basn[cel_i]] -= Flux*CELL_SIZE;
00488 P_IN_OVL[basn[cel_i]] -= m3;
00489 S_IN_OVL[basn[cel_i]] -= m1;
00490 }
00491 }
00492
00493
00494 if (debug > 0 && WatMgmtOn) {
00495 basins = basn_list[basn[cel_j]];
00496 flo_chek = 0;
00497 for (ii=0; ii<basins->numFLok; ii++) { if (basn[cel_i] == basins->FLok[ii] ) flo_chek = 1; }
00498 if (!flo_chek) {
00499 sprintf(msgStr,"Day %5.3f: ERROR - no (neg) SW flow from cell (%d,%d) of basin %d into cell (%d,%d) of basin %d!",
00500 SimTime.TIME, i0,i1, basn[cel_j], j0,j1, basn[cel_i]);
00501 WriteMsg( msgStr,True ); dynERRORnum++; }
00502 }
00503
00504
00505 }
00506 }
00507 else if ( !ON_MAP[cel_j])
00508 if (Flux > 0) {
00509 if ( !basn_list[basn[cel_i]]->parent ) {
00510 VOL_OUT_OVL[basn_list[basn[cel_i]]->family] += Flux*CELL_SIZE;
00511 P_OUT_OVL[basn_list[basn[cel_i]]->family] += m3;
00512 S_OUT_OVL[basn_list[basn[cel_i]]->family] += m1;
00513 }
00514
00515 VOL_OUT_OVL[basn[cel_i]] += Flux*CELL_SIZE;
00516 VOL_OUT_OVL[0]+= Flux*CELL_SIZE;
00517 P_OUT_OVL[basn[cel_i]] += m3;
00518 P_OUT_OVL[0]+= m3;
00519 S_OUT_OVL[basn[cel_i]] += m1;
00520 S_OUT_OVL[0]+= m1;
00521 }
00522 else {
00523 if ( !basn_list[basn[cel_i]]->parent ) {
00524 VOL_IN_OVL[basn_list[basn[cel_i]]->family] -= Flux*CELL_SIZE;
00525 P_IN_OVL[basn_list[basn[cel_i]]->family] -= m3;
00526 S_IN_OVL[basn_list[basn[cel_i]]->family] -= m1;
00527 }
00528
00529 VOL_IN_OVL[basn[cel_i]] -= Flux*CELL_SIZE;
00530 VOL_IN_OVL[0] -= Flux*CELL_SIZE;
00531 P_IN_OVL[basn[cel_i]] -= m3;
00532 P_IN_OVL[0] -= m3;
00533 S_IN_OVL[basn[cel_i]] -= m1;
00534 S_IN_OVL[0] -= m1;
00535 }
00536 else if ( !ON_MAP[cel_i])
00537 if (Flux > 0) {
00538 if ( !basn_list[basn[cel_j]]->parent ) {
00539 VOL_IN_OVL[basn_list[basn[cel_j]]->family] += Flux*CELL_SIZE;
00540 P_IN_OVL[basn_list[basn[cel_j]]->family] += m3;
00541 S_IN_OVL[basn_list[basn[cel_j]]->family] += m1;
00542 }
00543
00544 VOL_IN_OVL[basn[cel_j]] += Flux*CELL_SIZE;
00545 VOL_IN_OVL[0] += Flux*CELL_SIZE;
00546 P_IN_OVL[basn[cel_j]] += m3;
00547 P_IN_OVL[0] += m3;
00548 S_IN_OVL[basn[cel_j]] += m1;
00549 S_IN_OVL[0] += m1;
00550 }
00551 else {
00552 if ( !basn_list[basn[cel_j]]->parent ) {
00553 VOL_OUT_OVL[basn_list[basn[cel_j]]->family] -= Flux*CELL_SIZE;
00554 P_OUT_OVL[basn_list[basn[cel_j]]->family] -= m3;
00555 S_OUT_OVL[basn_list[basn[cel_j]]->family] -= m1;
00556 }
00557
00558 VOL_OUT_OVL[basn[cel_j]] -= Flux*CELL_SIZE;
00559 VOL_OUT_OVL[0] -= Flux*CELL_SIZE;
00560 P_OUT_OVL[basn[cel_j]] -= m3;
00561 P_OUT_OVL[0] -= m3;
00562 S_OUT_OVL[basn[cel_j]] -= m1;
00563 S_OUT_OVL[0] -= m1;
00564 }
00565 }
00566
00567 return;
00568
00569 }
00570
00571
00581 float Disp_Calc(float flux, float depth_i, float depth_j, float tim_step)
00582 {
00583 float disp_num;
00584 float veloc;
00585 float veloc_adj;
00586 float flux_adj;
00587 extern float dispParm_scaled;
00588
00589 veloc = Abs(flux) * celWid/( (flux>0.0) ? (depth_i) : (depth_j) ) / (tim_step) ;
00590 disp_num = 0.5 * veloc * (celWid - veloc * tim_step) ;
00591 veloc_adj = (veloc * celWid - disp_num)/celWid;
00592 flux_adj = dispParm_scaled * veloc_adj * tim_step * ( (flux>0.0) ? (depth_i) : (depth_j) )/celWid;
00593 return (flux_adj);
00594 }
00595
00596
00597
00632 void Flux_GWater(int it, float *SatWat, float *Unsat, float *SfWat,
00633 float *rate, float *poros, float *sp_yield, float *elev,
00634 double *gwSTUF1, double *gwSTUF2, double *gwSTUF3,
00635 double *swSTUF1, double *swSTUF2, double *swSTUF3)
00636
00637 {
00638 int ix, iy;
00639
00640 for(ix=1; ix<=s0; ix++)
00641 {
00642 if (it%2) {
00643 for(iy=1; iy<=s1; iy++)
00644
00645 if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy+1)]) ||
00646 ((BCondFlow[T(ix,iy+1)]+1) % 5 == 0) || ((BCondFlow[T(ix,iy)]+1) % 5 == 0) )
00647 {
00648 Flux_GWcells(ix,iy,ix,iy+1,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
00649 gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3);
00650 }
00651 }
00652
00653 else {
00654 for(iy=s1; iy>=1; iy--)
00655
00656 if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy-1)]) ||
00657 ((BCondFlow[T(ix,iy-1)]+1) % 5 == 0) || ((BCondFlow[T(ix,iy)]+1) % 5 == 0) )
00658 {
00659
00660 Flux_GWcells(ix,iy,ix,iy-1,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
00661 gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3);
00662 }
00663 }
00664 }
00665
00666 for(iy=1; iy<=s1; iy++)
00667 {
00668 if (it%2) {
00669 for(ix=1; ix<=s0; ix++)
00670
00671 if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix+1,iy)]) ||
00672 ((BCondFlow[T(ix+1,iy)]+1) % 5 == 0) || ((BCondFlow[T(ix,iy)]+1) % 5 == 0) )
00673 {
00674
00675 Flux_GWcells(ix,iy,ix+1,iy,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
00676 gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3);
00677 }
00678 }
00679 else {
00680 for(ix=s0; ix>=1; ix--)
00681
00682 if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix-1,iy)]) ||
00683 ((BCondFlow[T(ix-1,iy)]+1) % 5 == 0) || ((BCondFlow[T(ix,iy)]+1) % 5 == 0) )
00684 {
00685
00686 Flux_GWcells(ix,iy,ix-1,iy,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
00687 gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3);
00688 }
00689 }
00690 }
00691
00692 }
00693
00694
00726 void Flux_GWcells( int i0, int i1, int j0, int j1, float *SatWat, float *Unsat, float *SfWat,
00727 float *rate, float *poros, float *sp_yield, float *elev,
00728 double *gwSTUF1, double *gwSTUF2, double *gwSTUF3,
00729 double *swSTUF1, double *swSTUF2, double *swSTUF3)
00730 {
00731
00766 int cell_don, cell_rec, cell_i, cell_j;
00767 int sign, equilib, revers=1;
00768 float GW_head_i,GW_head_j,tot_head_i,tot_head_j;
00769 float dh;
00770 float AvgRate;
00771
00772 float H_pot_don, H_pot_rec;
00773
00774 double Flux;
00775 double fluxTOunsat_don;
00776 double fluxHoriz;
00777
00778 float Sat_vs_unsat;
00779
00780 float UnsatZ_don, UnsatCap_don, UnsatPot_don;
00781 float UnsatZ_rec, UnsatCap_rec, UnsatPot_rec;
00782 double sfTOsat_don, unsatTOsat_don, unsatTOsat_rec, satTOsf_rec, horizTOsat_rec;
00783
00784 float sedwat_don, sedwat_rec;
00785 double sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don;
00786 double sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz;
00787 double sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec;
00788
00789
00790 cell_i = T(i0,i1);
00791 cell_j = T(j0,j1);
00792
00793 if (ON_MAP[cell_i] ) {
00794 GW_head_i = SatWat[cell_i] / poros[HAB[cell_i]];
00795 tot_head_i = GW_head_i + SfWat[cell_i];
00796 }
00797
00798 if (ON_MAP[cell_j] ) {
00799 GW_head_j = SatWat[cell_j] / poros[HAB[cell_j]];
00800 tot_head_j = GW_head_j + SfWat[cell_j];
00801 }
00802
00814 if (!ON_MAP[cell_i] ) {
00815 poros[HAB[cell_i]] = poros[HAB[cell_j]];
00816 sp_yield[HAB[cell_i]] = sp_yield[HAB[cell_j]];
00817 elev[cell_i] = elev[cell_j];
00818 rate[cell_i] = rate[cell_j];
00819 Unsat[cell_i] = 0.0;
00820
00821
00822 if (BCondFlow[cell_j] == 4) {
00823 if (boundcond_depth[cell_i] > 0.0)
00824 {
00825 SatWat[cell_i] = elev[cell_i] * poros[HAB[cell_i]];
00826 SfWat[cell_i] = boundcond_depth[cell_i];
00827 }
00828 else
00829 {
00830 SatWat[cell_i] = (elev[cell_i] + boundcond_depth[cell_i]) * poros[HAB[cell_i]];
00831 SfWat[cell_i] = 0.0;
00832 }
00833 }
00834
00835 GW_head_i = SatWat[cell_i] / poros[HAB[cell_i]];
00836 tot_head_i = GW_head_i + SfWat[cell_i];
00837 }
00838
00839 if (!ON_MAP[cell_j] ) {
00840 poros[HAB[cell_j]] = poros[HAB[cell_i]];
00841 sp_yield[HAB[cell_j]] = sp_yield[HAB[cell_i]];
00842 elev[cell_j] = elev[cell_i];
00843 rate[cell_j] = rate[cell_i];
00844 Unsat[cell_j] = 0.0;
00845
00846
00847 if (BCondFlow[cell_i] == 4) {
00848 if (boundcond_depth[cell_j] > 0.0)
00849 {
00850 SatWat[cell_j] = elev[cell_j] * poros[HAB[cell_j]];
00851 SfWat[cell_j] = boundcond_depth[cell_j];
00852 }
00853 else
00854 {
00855 SatWat[cell_j] = (elev[cell_j] + boundcond_depth[cell_j]) * poros[HAB[cell_j]];
00856 SfWat[cell_j] = 0.0;
00857 }
00858 }
00859
00860 GW_head_j = SatWat[cell_j] / poros[HAB[cell_j]];
00861 tot_head_j = GW_head_j + SfWat[cell_j];
00862 }
00863
00864 AvgRate = ( rate[cell_i] + rate[cell_j] )/2.0;
00865
00866
00867 dh = tot_head_i - tot_head_j;
00868
00869
00870
00871
00872
00873
00874 if (dh > GP_DetentZ)
00875 {
00876 cell_don=cell_i;
00877 cell_rec=cell_j;
00878 sign = 1;
00879 }
00880 else if (dh < -GP_DetentZ)
00881 {
00882 cell_don=cell_j;
00883 cell_rec=cell_i;
00884 sign = -1;
00885 }
00886 else
00887 return;
00888
00889
00890
00891
00892
00893
00894 Flux = Min(Abs(dh) * AvgRate * SatWat[cell_don] / CELL_SIZE * gwstep , SatWat[cell_don]);
00895
00896
00897
00898
00899
00900 do {
00901
00902
00903
00904
00905
00906
00907
00908
00909 fluxTOunsat_don = Flux/poros[HAB[cell_don]] * (poros[HAB[cell_don]] - sp_yield[HAB[cell_don]]) ;
00910 fluxHoriz = Flux - fluxTOunsat_don;
00911
00912
00913
00914 UnsatZ_don = elev[cell_don] - (SatWat[cell_don]-fluxHoriz ) / poros[HAB[cell_don]] ;
00915
00916
00917
00918
00919 if (debug>3 && (UnsatZ_don < -0.01 ) ) {
00920 sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in donor cell (%d,%d) in GWfluxes!",
00921 SimTime.TIME, UnsatZ_don, i0,i1 ); WriteMsg( msgStr,True ); dynERRORnum++; }
00922
00923 UnsatCap_don = UnsatZ_don * poros[HAB[cell_don]];
00924
00925 UnsatPot_don = UnsatCap_don - (Unsat[cell_don]+fluxTOunsat_don);
00926
00927
00928
00929
00930 Sat_vs_unsat = 1/Exp(100.0*Max((SfWat[cell_don]-UnsatZ_don),0.0));
00931
00932
00933 if (SfWat[cell_don] > 0.0) {
00934
00935 sfTOsat_don =
00936 ( (1.0-Sat_vs_unsat)*UnsatPot_don>SfWat[cell_don] ) ?
00937 ( SfWat[cell_don] ) :
00938 ( (1.0-Sat_vs_unsat)*UnsatPot_don);
00939
00940
00941 if (sfTOsat_don >= UnsatPot_don) {
00942 sfTOsat_don = UnsatPot_don ;
00943 unsatTOsat_don = Unsat[cell_don];
00944 }
00945 else {
00946 unsatTOsat_don = (UnsatZ_don > 0.0) ?
00947 ( (sfTOsat_don/poros[HAB[cell_don]] ) / UnsatZ_don * Unsat[cell_don] ) :
00948 (0.0);
00949
00950 }
00951 }
00952 else {
00953 sfTOsat_don = unsatTOsat_don = 0.0;
00954 }
00955
00956 H_pot_don = (SatWat[cell_don] - fluxTOunsat_don - fluxHoriz + sfTOsat_don + unsatTOsat_don )
00957 / poros[HAB[cell_don]] +(SfWat[cell_don] - sfTOsat_don) ;
00958
00959
00960
00961 UnsatZ_rec = elev[cell_rec] - SatWat[cell_rec] / poros[HAB[cell_rec]] ;
00962
00963
00964
00965
00966 if (debug>3 && (UnsatZ_rec < -0.01 ) ) {
00967 sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in recipient cell (%d,%d) in GWfluxes!",
00968 SimTime.TIME, UnsatZ_rec, i0,i1 ); WriteMsg( msgStr,True ); dynERRORnum++; }
00969
00970 UnsatCap_rec = UnsatZ_rec * poros[HAB[cell_rec]];
00971
00972 UnsatPot_rec = UnsatCap_rec - Unsat[cell_rec];
00973
00974
00975 horizTOsat_rec = fluxHoriz;
00976 satTOsf_rec = Max(fluxHoriz - UnsatPot_rec, 0.0);
00977 unsatTOsat_rec = (UnsatZ_rec > 0.0) ?
00978
00979 ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[cell_rec]] ) / UnsatZ_rec * Unsat[cell_rec] ) :
00980 (0.0);
00981
00982 H_pot_rec = (SatWat[cell_rec] + horizTOsat_rec + unsatTOsat_rec - satTOsf_rec)
00983 / poros[HAB[cell_rec]] + (SfWat[cell_rec] + satTOsf_rec) ;
00984
00985
00986 if ( (Flux != 0.0) && ((H_pot_rec - H_pot_don) > GP_MinCheck) ) {
00987 revers++;
00988 Flux *= 0.90;
00989 equilib = 0;
00990
00991
00992
00993 if (revers>1) {
00994 if (debug>2) { sprintf(msgStr,"Day %6.1f: FYI - head reversal (%f m) due to surf/ground integration, associated with cell (%d,%d). Total flux was to be %f. Fixed with 0 flux. ",
00995 SimTime.TIME, (H_pot_don - H_pot_rec), i0,i1,Flux*sign ); WriteMsg( msgStr,True );}
00996 Flux = 0.0;
00997 }
00998 }
00999 else {
01000 equilib = 1;
01001 }
01002
01003 } while (!equilib);
01004
01005
01006
01007
01008
01009
01010 sedwat_don = SatWat[cell_don] + Unsat[cell_don];
01011 sed_stuf1fl_horz = (sedwat_don > 0.0) ? (gwSTUF1[cell_don] / sedwat_don*fluxHoriz) : (0.0);
01012 sed_stuf2fl_horz = (sedwat_don > 0.0) ? (gwSTUF2[cell_don] / sedwat_don*fluxHoriz) : (0.0);
01013 sed_stuf3fl_horz = (sedwat_don > 0.0) ? (gwSTUF3[cell_don] / sedwat_don*fluxHoriz) : (0.0);
01014
01015
01016
01017
01018 if (sfTOsat_don > 0.0) {
01019 sf_stuf1fl_don = (SfWat[cell_don] > 0.0) ? (swSTUF1[cell_don] / SfWat[cell_don]*sfTOsat_don) : (0.0);
01020 sf_stuf2fl_don = (SfWat[cell_don] > 0.0) ? (swSTUF2[cell_don] / SfWat[cell_don]*sfTOsat_don) : (0.0);
01021 sf_stuf3fl_don = (SfWat[cell_don] > 0.0) ? (swSTUF3[cell_don] / SfWat[cell_don]*sfTOsat_don) : (0.0);
01022 gwSTUF1[cell_don] += sf_stuf1fl_don;
01023 gwSTUF2[cell_don] += sf_stuf2fl_don;
01024 gwSTUF3[cell_don] += sf_stuf3fl_don;
01025 swSTUF1[cell_don] -= sf_stuf1fl_don;
01026 swSTUF2[cell_don] -= sf_stuf2fl_don;
01027 swSTUF3[cell_don] -= sf_stuf3fl_don;
01028
01029 }
01030 if (satTOsf_rec > 0.0) {
01031 sedwat_rec = SatWat[cell_rec] + Unsat[cell_rec];
01032 sed_stuf1fl_rec = (sedwat_rec > 0.0) ? (gwSTUF1[cell_rec] / sedwat_rec*satTOsf_rec) : (0.0);
01033 sed_stuf2fl_rec = (sedwat_rec > 0.0) ? (gwSTUF2[cell_rec] / sedwat_rec*satTOsf_rec) : (0.0);
01034 sed_stuf3fl_rec = (sedwat_rec > 0.0) ? (gwSTUF3[cell_rec] / sedwat_rec*satTOsf_rec) : (0.0);
01035 gwSTUF1[cell_rec] -= sed_stuf1fl_rec;
01036 gwSTUF2[cell_rec] -= sed_stuf2fl_rec;
01037 gwSTUF3[cell_rec] -= sed_stuf3fl_rec;
01038 swSTUF1[cell_rec] += sed_stuf1fl_rec;
01039 swSTUF2[cell_rec] += sed_stuf2fl_rec;
01040 swSTUF3[cell_rec] += sed_stuf3fl_rec;
01041 }
01042
01043
01044 gwSTUF1[cell_don] -= sed_stuf1fl_horz;
01045 gwSTUF2[cell_don] -= sed_stuf2fl_horz;
01046 gwSTUF3[cell_don] -= sed_stuf3fl_horz;
01047 gwSTUF1[cell_rec] += sed_stuf1fl_horz;
01048 gwSTUF2[cell_rec] += sed_stuf2fl_horz;
01049 gwSTUF3[cell_rec] += sed_stuf3fl_horz;
01050
01051
01052
01053 SfWat[cell_don] += (-sfTOsat_don);
01054 Unsat[cell_don] += ( fluxTOunsat_don - unsatTOsat_don) ;
01055 SatWat[cell_don] += (sfTOsat_don + unsatTOsat_don - fluxTOunsat_don - fluxHoriz);
01056 SfWat[cell_rec] += ( satTOsf_rec);
01057 Unsat[cell_rec] += (-unsatTOsat_rec);
01058 SatWat[cell_rec] += (horizTOsat_rec + unsatTOsat_rec - satTOsf_rec);
01059
01060
01061
01062
01063 if ( basn[cell_i] != basn[cell_j] ) {
01064
01065 fluxHoriz = sign*fluxHoriz*CELL_SIZE;
01066 sed_stuf1fl_horz = sign*sed_stuf1fl_horz;
01067 sed_stuf3fl_horz = sign*sed_stuf3fl_horz;
01068
01069
01070 if ( ON_MAP[cell_j] && ON_MAP[cell_i] ) {
01071 if ( fluxHoriz > 0 ) {
01072 if (basn_list[basn[cell_i]]->family !=
01073 basn_list[basn[cell_j]]->family ) {
01074 if ( !basn_list[basn[cell_i]]->parent ) {
01075
01076 VOL_OUT_GW[basn_list[basn[cell_i]]->family] += fluxHoriz;
01077 P_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf3fl_horz;
01078 S_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf1fl_horz;
01079 }
01080 if ( !basn_list[basn[cell_j]]->parent ) {
01081
01082 VOL_IN_GW[basn_list[basn[cell_j]]->family] += fluxHoriz;
01083 P_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf3fl_horz;
01084 S_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf1fl_horz;
01085 }
01086
01087 VOL_OUT_GW[basn[cell_i]] += fluxHoriz;
01088 VOL_IN_GW[basn[cell_j]] += fluxHoriz;
01089 P_OUT_GW[basn[cell_i]] += sed_stuf3fl_horz;
01090 P_IN_GW[basn[cell_j]] += sed_stuf3fl_horz;
01091 S_OUT_GW[basn[cell_i]] += sed_stuf1fl_horz;
01092 S_IN_GW[basn[cell_j]] += sed_stuf1fl_horz;
01093 }
01094
01095
01096 else {
01097
01098 if ( !basn_list[basn[cell_i]]->parent ) {
01099 VOL_OUT_GW[basn[cell_i]] += fluxHoriz;
01100 P_OUT_GW[basn[cell_i]] += sed_stuf3fl_horz;
01101 S_OUT_GW[basn[cell_i]] += sed_stuf1fl_horz;
01102 }
01103 if ( !basn_list[basn[cell_j]]->parent ) {
01104 VOL_IN_GW[basn[cell_j]] += fluxHoriz;
01105 P_IN_GW[basn[cell_j]] += sed_stuf3fl_horz;
01106 S_IN_GW[basn[cell_j]] += sed_stuf1fl_horz;
01107 }
01108 }
01109
01110 }
01111 else {
01112 if (basn_list[basn[cell_i]]->family !=
01113 basn_list[basn[cell_j]]->family ) {
01114 if ( !basn_list[basn[cell_j]]->parent ) {
01115
01116 VOL_OUT_GW[basn_list[basn[cell_j]]->family] -= fluxHoriz;
01117 P_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf3fl_horz;
01118 S_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf1fl_horz;
01119 }
01120 if ( !basn_list[basn[cell_i]]->parent ) {
01121
01122 VOL_IN_GW[basn_list[basn[cell_i]]->family] -= fluxHoriz;
01123 P_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf3fl_horz;
01124 S_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf1fl_horz;
01125 }
01126
01127 VOL_OUT_GW[basn[cell_j]] -= fluxHoriz;
01128 VOL_IN_GW[basn[cell_i]] -= fluxHoriz;
01129 P_OUT_GW[basn[cell_j]] -= sed_stuf3fl_horz;
01130 P_IN_GW[basn[cell_i]] -= sed_stuf3fl_horz;
01131 S_OUT_GW[basn[cell_j]] -= sed_stuf1fl_horz;
01132 S_IN_GW[basn[cell_i]] -= sed_stuf1fl_horz;
01133
01134 }
01135 else {
01136
01137 if ( !basn_list[basn[cell_j]]->parent ) {
01138 VOL_OUT_GW[basn[cell_j]] -= fluxHoriz;
01139 P_OUT_GW[basn[cell_j]] -= sed_stuf3fl_horz;
01140 S_OUT_GW[basn[cell_j]] -= sed_stuf1fl_horz;
01141 }
01142 if ( !basn_list[basn[cell_i]]->parent ) {
01143 VOL_IN_GW[basn[cell_i]] -= fluxHoriz;
01144 P_IN_GW[basn[cell_i]] -= sed_stuf3fl_horz;
01145 S_IN_GW[basn[cell_i]] -= sed_stuf1fl_horz;
01146 }
01147 }
01148 }
01149
01150
01151 }
01152 else if ( !ON_MAP[cell_j]) {
01153 if (fluxHoriz > 0) {
01154 if ( !basn_list[basn[cell_i]]->parent ) {
01155 VOL_OUT_GW[basn_list[basn[cell_i]]->family] += fluxHoriz;
01156 P_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf3fl_horz;
01157 S_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf1fl_horz;
01158 }
01159
01160 VOL_OUT_GW[basn[cell_i]] += fluxHoriz;
01161 VOL_OUT_GW[0] += fluxHoriz;
01162 P_OUT_GW[basn[cell_i]] += sed_stuf3fl_horz;
01163 P_OUT_GW[0] += sed_stuf3fl_horz;
01164 S_OUT_GW[basn[cell_i]] += sed_stuf1fl_horz;
01165 S_OUT_GW[0] += sed_stuf1fl_horz;
01166 }
01167 else {
01168 if ( !basn_list[basn[cell_i]]->parent ) {
01169 VOL_IN_GW[basn_list[basn[cell_i]]->family] -= fluxHoriz;
01170 P_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf3fl_horz;
01171 S_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf1fl_horz;
01172 }
01173
01174 VOL_IN_GW[basn[cell_i]] -= fluxHoriz;
01175 VOL_IN_GW[0] -= fluxHoriz;
01176 P_IN_GW[basn[cell_i]] -= sed_stuf3fl_horz;
01177 P_IN_GW[0] -= sed_stuf3fl_horz;
01178 S_IN_GW[basn[cell_i]] -= sed_stuf1fl_horz;
01179 S_IN_GW[0] -= sed_stuf1fl_horz;
01180 }
01181 }
01182
01183 else if ( !ON_MAP[cell_i]) {
01184 if (fluxHoriz > 0) {
01185 if ( !basn_list[basn[cell_j]]->parent ) {
01186 VOL_IN_GW[basn_list[basn[cell_j]]->family] += fluxHoriz;
01187 P_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf3fl_horz;
01188 S_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf1fl_horz;
01189 }
01190
01191 VOL_IN_GW[basn[cell_j]] += fluxHoriz;
01192 VOL_IN_GW[0] += fluxHoriz;
01193 P_IN_GW[basn[cell_j]] += sed_stuf3fl_horz;
01194 P_IN_GW[0] += sed_stuf3fl_horz;
01195 S_IN_GW[basn[cell_j]] += sed_stuf1fl_horz;
01196 S_IN_GW[0] += sed_stuf1fl_horz;
01197 }
01198 else {
01199 if ( !basn_list[basn[cell_j]]->parent ) {
01200 VOL_OUT_GW[basn_list[basn[cell_j]]->family] -= fluxHoriz;
01201 P_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf3fl_horz;
01202 S_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf1fl_horz;
01203 }
01204
01205 VOL_OUT_GW[basn[cell_j]] -= fluxHoriz;
01206 VOL_OUT_GW[0] -= fluxHoriz;
01207 P_OUT_GW[basn[cell_j]] -= sed_stuf3fl_horz;
01208 P_OUT_GW[0] -= sed_stuf3fl_horz;
01209 S_OUT_GW[basn[cell_j]] -= sed_stuf1fl_horz;
01210 S_OUT_GW[0] -= sed_stuf1fl_horz;
01211 }
01212 }
01213 }
01214
01215
01216 return ;
01217
01218 }
01219
01220
01221
01222
01223