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

Fluxes.c

Go to the documentation of this file.
00001 
00023 /* double precision for all constituent vars */
00024 
00025 /* General NOTES on revisions to this source file:
00026         April 00 VERSION 2.1 - output available for application for water quality performance
00027         July 02 VERSIOIN 2.1a - first widely-available public release of code -
00028               - made a start at adding dynamic stage boundary conditions - INCOMPLETE, not functional
00029               - added some more code documentation, etc.
00030 
00031         July 2004 VERSION 2.2.1 - added dynamic stage boundary conditions
00032               - in globals (celldyn1) module of UnitMod.c, have read daily, coarse-grid boundary stage data
00033               - in this version, we have not turned on that new boundary condition code, still using v2.1 
00034                 estimated boundary conditions
00035         Aug 2004 VERSION 2.3.0 - added dynamic boundary conditions 
00036                 - used new spatial time series data input for dynamic boundary conditions
00037         Oct 2004 VERSION 2.3.1 - internal release 
00038                 - checked that results reasonable, not fully checked
00039         Nov 2004 VERSION 2.3.2 - documentation upgrade
00040         Jan 2005 v2.4.0 - encoded the AND dispersion
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   /* check the donor and recipients cells for a) on-map, b) the cell attribute that allows sfwater
00121      boundary flow from the system and c) the attribute that indicates levee presence:
00122      the levee attribute of 1 (bitwise) allows flow to east;
00123      attribute of 2 (bitwise) allows flow to south (levee atts calc'd by WatMgmt.c) 
00124    */
00125 
00126   /* as always, x is row, y is column! */
00127   for(ix=1; ix<=s0; ix++) 
00128   {
00129     if (it%2)   /* alternate loop directions every other hyd_iter (it) */
00130     {
00131       for(iy=1; iy<=s1; iy++)  /* loop from west to east */
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           /* FFlux units = m */
00138           if (FFlux != 0.0) 
00139             Flux_SWstuff ( ix,iy,ix,iy+1,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00140 
00141         } /* endof if */
00142       } /* end of for */
00143     }  /* end of if */
00144      
00145 
00146     else 
00147     { 
00148       for(iy=s1; iy>=1; iy--)   /* loop from east to west */
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                                 /* FFlux units = m */
00154           if (FFlux != 0.0) 
00155             Flux_SWstuff ( ix,iy-1,ix,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00156 
00157         }  /* end of if */
00158     }  /* end of else */
00159   }
00160         
00161   for(iy=1; iy<=s1; iy++) 
00162   {
00163     if (it%2)   /* alternate loop directions every other hyd_iter (it) */
00164     {
00165       for(ix=1; ix<=s0; ix++)  /* loop from north to south */
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           /* FFlux units = m */
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--)  /* loop from south to north */
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           /* FFlux units = m */
00186           if (FFlux != 0.0) 
00187             Flux_SWstuff ( ix-1,iy,ix,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3);
00188 
00189         } /* end of if */
00190     } /* end of else */
00191   } /* end of for */
00192 
00193 } /* end of function */
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   /* If an on-map cell is marked 3, we are at a model boundary allowing surface water exchange */     
00236   if (!ON_MAP[cellLoci] && BCondFlow[cellLocj] == 3 ) 
00237   {
00238     /* the off-map cell given head 5 cm less than donor */
00239     /* Hi = Elevation[cellLocj] + Max(SWater[cellLocj]-0.05,0.0) ;  */ /* v2.3 not using this Hi */ 
00240 
00241     /* new dynamic boundary condition stage */
00242     Hi = boundcond_depth[cellLoci]; /* ?v2.2 not using this Hi */ 
00243 
00244     MC_cells = MC[cellLocj];       /* the mannings n is not avg, but the value of onmap boundary cell */
00245   }
00246 
00247   else 
00248     Hi = SWater[cellLoci] + Elevation[cellLoci];
00249 
00250   if(BCondFlow[cellLoci] == 3 && !ON_MAP[cellLocj] )   
00251   {  
00252     /* the off-map cell given head 5 cm less than donor  */
00253       /* Hj = Elevation[cellLoci] + Max(SWater[cellLoci]-0.05,0.0); */   /* v2.3 not using this Hj */ 
00254 
00255     Hj = boundcond_depth[cellLocj]; /* ?v2.2 not using this Hj */
00256 
00257     MC_cells = MC[cellLoci];      /* the mannings n is not avg, but the value of onmap boundary cell */
00258   }
00259 
00260   else 
00261     Hj = SWater[cellLocj] + Elevation[cellLocj];
00262 
00263   dh = Hi - Hj;         /* dh is "from --> to" */
00264   adh = Abs (dh);
00265                 
00266   if (dh > 0) 
00267   {
00268     if(SWater[cellLoci] < GP_DetentZ) 
00269       return 0.0; 
00270     /* step_Cell is constant calc'd in Generic_Driver.c at initialization ( m^(-1.5) * sec )*/
00271     Flux = (MC_cells != 0) ? 
00272         (pow(adh,GP_mannHeadPow) / MC_cells  * pow(SWater[cellLoci],GP_mannDepthPow)*step_Cell) : (0.0);
00273 
00274     /* ensure adequate volume avail */
00275     Flux =  ( Flux > ramp(SWater[cellLoci] - GP_DetentZ) ) ? (ramp(SWater[cellLoci] - GP_DetentZ)) : (Flux);
00276 
00277     /* check to ensure no flip/flops associated with depth */
00278     if ( ( Hi - Flux ) < ( Hj + Flux ) )        
00279       Flux = Min ( dh/2.0, ramp(SWater[cellLoci] - GP_DetentZ) );
00280 
00281   } /* end of dh > 0 */
00282 
00283   else
00284   {     
00285     if (SWater[cellLocj] < GP_DetentZ) 
00286       return 0.0; 
00287     /* step_Cell is constant calc'd in Generic_Driver.c at initialization ( m^(-1.5) * sec )*/
00288     /* Flux is negative in this case */
00289     Flux = (MC_cells != 0) ? 
00290            ( - pow(adh,GP_mannHeadPow) / MC_cells  * pow(SWater[cellLocj],GP_mannDepthPow)*step_Cell) : (0.0);
00291 
00292     /* ensure adequate volume avail */
00293     Flux =  ( -Flux > ramp(SWater[cellLocj] - GP_DetentZ) ) ? (-ramp(SWater[cellLocj] - GP_DetentZ)) : (Flux);
00294 
00295     /* check to ensure no flip/flops associated with depth */           
00296     if ( ( Hi - Flux ) > ( Hj + Flux ) )  
00297       Flux = - Min ( adh/2.0, ramp(SWater[cellLocj] - GP_DetentZ) );
00298   }
00299         
00300   return (Flux);        /* returns height 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; /* numerical dispersion (m2/d) */
00329     float veloc;    /* velocity (m/d) */
00330     float velocAdj;  /* velocity adjusted for numerical dispersion (m/d) */
00331     float FluxAdj; /* Flux adjusted for numerical dispersioin */
00332     double fl_prop_i, fl_prop_j; /* proportion of surface water constituents that should be fluxed with water (depends on magnitude of dispersion) */
00333     extern float dispParm_scaled; 
00334 
00335     cel_i = T(i0,i1);
00336     cel_j = T(j0,j1);
00337     
00338 /*    veloc =  Abs(Flux) * celWid/( (Flux>0.0)?(SURFACE_WAT[cel_i]):(SURFACE_WAT[cel_j])) / (sfstep) ;  
00339     disp_num = 0.5 * veloc * (celWid - veloc * sfstep)  ; 
00340     velocAdj = (veloc * celWid - disp_num)/celWid; 
00341     FluxAdj = dispParm_scaled * velocAdj * sfstep * ( (Flux>0.0)?(SURFACE_WAT[cel_i]):(SURFACE_WAT[cel_j]))/celWid;
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); /* pos Flux */
00346     fl_prop_j = (SURFACE_WAT[cel_j]>0.0) ? (Min(Flux+FluxAdj,0.0)/SURFACE_WAT[cel_j]) : (0.0); /* neg Flux */
00347         /*  if (i0==41 && i1==12 && veloc>100.0) printf ("\nveloc=%f m/d, disp_num=%f m2/d, FF=%f m; FFadj=%f m \n",veloc, disp_num, Flux, FluxAdj); */
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;  /* add the masses of constituents */
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; /* now update the surfwater depths */
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 /* mass balance sums */
00401  if (basn[cel_i] != basn[cel_j])  { 
00402 
00403         /* first do the normal case where all cells are on-map */
00404      if ( ON_MAP[cel_j] && ON_MAP[cel_i] ) { 
00405          if ( Flux > 0  ) { /* positive fluxes out of basn[cel_i] */
00406              if (basn_list[basn[cel_i]]->family !=  
00407                  basn_list[basn[cel_j]]->family ) {        /* if the flow is not within the family... */
00408                  if  ( !basn_list[basn[cel_i]]->parent  ) { /* and if the donor or recipient is a child... */
00409                      /* then find out about the flow for the family's sake */
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                      /* then find out about the flow for the family's sake */
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                      /* now sum the parents' flows */
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 {    /* if it's flow within a family, just keep
00430                           track of what the children do among themselves */            
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) { /* check for basin misconfiguration (allowable basin-basin flows) */
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 { /* negative fluxes out of basn[cel_j] */
00456              if (basn_list[basn[cel_i]]->family !=  
00457                  basn_list[basn[cel_j]]->family ) {        /* if the flow is not within the family... */
00458                  if  ( !basn_list[basn[cel_j]]->parent  ) { /* and if the donor or recipient is a child... */
00459                      /* then find out about the flow for the family's sake */
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                      /* then find out about the flow for the family's sake */
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                      /* now sum the parents' flows */
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 {    /* if it's flow within a family, just keep
00480                           track of what the children do among themselves */            
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) { /* check for basin misconfiguration (allowable basin-basin flows) */
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]) /* so now the j,j cell is off-map, recipient if pos flow */
00508          if (Flux > 0) { 
00509              if ( !basn_list[basn[cel_i]]->parent ) { /* child's play */
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              /* parents' play */
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 { /* negative flows */
00523              if ( !basn_list[basn[cel_i]]->parent ) {/* child's play */
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              /* parents' play */
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]) /* so now the i,i cell is off-map, donor if pos flow */
00537          if (Flux > 0) { 
00538              if ( !basn_list[basn[cel_j]]->parent ) {/* child's play */
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              /* parents' play */
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 { /*negative flows */
00552              if ( !basn_list[basn[cel_j]]->parent ) {/* child's play */
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              /* parents' play */
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; /* numerical dispersion (m2/d) */
00584     float veloc;    /* velocity (m/d) */
00585     float veloc_adj;  /* velocity adjusted for numerical dispersion (m/d) */
00586     float flux_adj; /* Flux adjusted for numerical dispersioin */
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; /* ix is row, iy is col! */
00639 /* we only check the donor cell for on-map, allowing losses from the system */
00640  for(ix=1; ix<=s0; ix++) 
00641  {
00642      if (it%2) {  /* alternate loop directions every other iteration (it = int(hyd_iter-1)/2 ) */
00643      for(iy=1; iy<=s1; iy++)  /* loop from west to east */
00644              /* allow boundary flow if donor cell is marked 4 or 9 */
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--)   /* loop from east to west */
00655              /* allow boundary flow if donor cell is marked 4 or 9 */
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) {  /* alternate loop directions every other iteration (it = int(hyd_iter-1)/2 ) */
00669      for(ix=1; ix<=s0; ix++)  /* loop from north to south */
00670              /* allow boundary flow if donor cell is marked 4 or 9 */
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--)  /* loop from south to north */
00681              /* allow boundary flow if donor cell is marked 4 or 9 */
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] ) { /* for an external i0,i1 cell */
00815         poros[HAB[cell_i]] = poros[HAB[cell_j]]; /* other variables = donor cell or zero */
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; /* no water in any potential unsat zone */
00820 
00821             /* outflow of groundwater from boundary cell*/
00822         if (BCondFlow[cell_j] == 4) { 
00823             if (boundcond_depth[cell_i] > 0.0) /* boundcond_depth=relative water depth; v2.3tmp instead of: boundcond_depth[cell_i] > elev[cell_i]  */
00824             {
00825                 SatWat[cell_i] =  elev[cell_i] * poros[HAB[cell_i]];    
00826                 SfWat[cell_i] = boundcond_depth[cell_i]; /* boundcond_depth=relative water depth; v2.3tmp instead of: boundcond_depth[cell_i] - elev[cell_i]  */
00827             }
00828         else
00829             {
00830                 SatWat[cell_i] =  (elev[cell_i] + boundcond_depth[cell_i]) * poros[HAB[cell_i]];  /* boundcond_depth=relative water depth; v2.3tmp instead of: 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] ) { /* for an external j0,j1 cell */
00840         poros[HAB[cell_j]] = poros[HAB[cell_i]]; /* other variables = donor cell or zero */
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; /* no water in any potential unsat zone */
00845             
00846             /* outflow of groundwater from boundary cell*/
00847         if (BCondFlow[cell_i] == 4) { 
00848             if (boundcond_depth[cell_j] > 0.0) /* boundcond_depth=relative water depth; v2.3tmp instead of: boundcond_depth[cell_j] > elev[cell_j]  */
00849             {
00850                 SatWat[cell_j] =  elev[cell_j] * poros[HAB[cell_j]];    
00851                 SfWat[cell_j] = boundcond_depth[cell_j]; /* boundcond_depth=relative water depth; v2.3tmp instead of: boundcond_depth[cell_j] - elev[cell_j]  */
00852             }
00853         else
00854             {
00855                 SatWat[cell_j] =  (elev[cell_j] + boundcond_depth[cell_j]) * poros[HAB[cell_j]];  /* boundcond_depth=relative water depth; v2.3tmp instead of: 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; /* average hyd conductivity of the two cells */
00865     
00866 /* hydraulic head difference, positive flux from cell_i to cell_j */
00867     dh = tot_head_i - tot_head_j; 
00868 
00869  
00870         /* determine donor and recipient cells based on head difference,
00871            and provide the sign of the flux;
00872            The surface water detention depth also used here as threshold
00873            head difference for fluxing (currently global, 1 cm)*/
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; /* no flux (or surface-groundwater integration) under minimal head diffs */
00888     
00889 
00890 /* Potential horizontal flux eqn (Darcy's eqn simplified to square cells).
00891    This is the maximum height (m) of water vol to flux under fully saturated condition.
00892    Note that the Min check is probably not important under current implementations, as SatWat includes water
00893    down to the base datum, which is currently 6 meters below sea level (actually, 1929 NGVD) */
00894     Flux = Min(Abs(dh) * AvgRate * SatWat[cell_don] / CELL_SIZE * gwstep , SatWat[cell_don]);
00895  
00896 /**** this do-while routine (1) integrates the surface, saturated, and unsaturated water,
00897   and (2) checks to ensure the heads do not reverse in a time step due to large fluxes.
00898   If heads do reverse, the total Flux is decremented until there is no reversal */
00899 /****/
00900     do {
00901         /* The total potential flux is apportioned to (1) the horizontal component
00902            that fluxes to an adjacent cell and (2) the vertical component that
00903            remains in the donor cell after the horizontal outflow from a donor cell.
00904            Thus, an unsaturated zone is created, or increased in size, with loss of
00905            saturated water from the donor cell; this lateral gravitational flow leaves behind
00906            the field capacity moisture in an unsat zone. (If donor-cell surface water is present,
00907            it potentially will replace the unsaturated soil capacity within the same time
00908            step in this routine).*/
00909         fluxTOunsat_don = Flux/poros[HAB[cell_don]] * (poros[HAB[cell_don]] - sp_yield[HAB[cell_don]])  ;
00910         fluxHoriz = Flux - fluxTOunsat_don;
00911 
00912 /**** given the total potential groundwater Flux, get the donor cell's
00913       new **post-flux** capacities */
00914         UnsatZ_don  =   elev[cell_don] - (SatWat[cell_don]-fluxHoriz /*Flux*/) / poros[HAB[cell_don]] ; /* unsat zone depth */
00915     
00916     /* ?v2.2 this check was against "-0.001", increased here to "-0.01", as it was only showing several mm capacityERR
00917         due to canal interactions (i.e., only when watmgmt turned on), and we haven't explicitly added sf-ground 
00918         integration to grid cells in canal code */
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]]; /* maximum pore space capacity (UnsatCap)
00924                                                               in the depth (UnsatZ) of new unsaturated zone */
00925         UnsatPot_don  = UnsatCap_don - (Unsat[cell_don]+fluxTOunsat_don); /* (height of) the volume of pore space
00926                                                                              (soil "removed") that is unoccupied in the unsat zone */
00927             /* determining the pathway of flow (to sat vs. unsat) of surface water depending on depth
00928                of an unsat zone relative to the surface water.  With a relatively deep unsat zone, this
00929                downflow tends to zero (infiltration occurs within the vertical hydrology module of UnitMod.c) */ 
00930         Sat_vs_unsat  = 1/Exp(100.0*Max((SfWat[cell_don]-UnsatZ_don),0.0)); 
00931 
00932     /* sf-unsat-sat fluxes */
00933         if (SfWat[cell_don] > 0.0) { /* surface water downflow is assumed to be as fast
00934                                         as horizontal groundwater outflows */
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            /* with downflow of surface water into an unsat zone, the proportion of that height
00940               that is made into saturated storage is allocated to the sat storage variable */
00941             if (sfTOsat_don >=  UnsatPot_don) {
00942                     sfTOsat_don = UnsatPot_don ; /* downflow constrained to the avail soil capacity */
00943                     unsatTOsat_don = Unsat[cell_don]; /* allocate unsat to sat in this case */
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); /*  allocate to saturated storage whatever proportion of
00949                                unsat zone that is now saturated by sfwat downflow */
00950             }
00951         }
00952         else { /* w/o surface water, these vertical fluxes are zero */
00953             sfTOsat_don = unsatTOsat_don = 0.0;
00954         }
00955 /**** potential new head in donor cell will be ... */
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 /**** recipient cell's **pre-flux** capacities */
00961         UnsatZ_rec  =   elev[cell_rec] - SatWat[cell_rec] / poros[HAB[cell_rec]] ; /* unsat zone depth */
00962     
00963     /* ?v2.2 this check was against "-0.001", increased here to "-0.01", as it was only showing several mm capacityERR
00964         due to canal interactions (i.e., only when watmgmt turned on), and we haven't explicitly added sf-ground 
00965         integration to grid cells in canal code */
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]]; /* maximum pore space capacity (UnsatCap)
00971                                                               in the depth (UnsatZ) of new unsaturated zone */
00972         UnsatPot_rec  = UnsatCap_rec - Unsat[cell_rec]; /* (height of) the volume of pore space (soil "removed")
00973                                                            that is unoccupied in the unsat zone */
00974      /* sf-unsat-sat fluxes */
00975         horizTOsat_rec = fluxHoriz; /* lateral inflow to soil into saturated storage */
00976         satTOsf_rec = Max(fluxHoriz - UnsatPot_rec, 0.0); /* upwelling beyond soil capacity */
00977         unsatTOsat_rec = (UnsatZ_rec > 0.0) ? /* incorporation of unsat moisture into sat storage with
00978                                                  rising water table due to horiz inflow */
00979             ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[cell_rec]] ) / UnsatZ_rec * Unsat[cell_rec] ) :
00980             (0.0);
00981 /**** potential new head in recipient cell will be ... */
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 /**** check for a head reversal - if so, reduce flux to achieve equilibrium (donor >= recip) */
00986         if ( (Flux != 0.0) && ((H_pot_rec - H_pot_don) > GP_MinCheck) ) {
00987             revers++;
00988             Flux *= 0.90;
00989             equilib = 0; 
00990                /* if the unsat water storage causes an unfixable (very rare(/never?)) head reversal due to the
00991                    sfwat and gwater integration alone, set the Flux to zero, then allow the vertical
00992                    sf/ground integration to be done with no horiz flux and be done with it */
00993             if (revers>1) { /* v2.3 temp, was >4 */
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); /* end of routine for integrating groundwater-surface water and preventing head reversals */
01004         
01005 /**********
01006   finished calc'ing the water fluxes
01007   *************/
01008 /* calc the flux of the constituents between sed/soil across cells in horiz dir,
01009  but don't update the state vars until the conc's for vertical flows have been calc'd */
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 /* pass along the constituents between sed/soil and surface water in vertical dir;
01016    for "simplicity", this does not account for the phosphorus active-zone conc. gradient,
01017    but assumes the homogeneity of conc within the entire soil column */
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 /* update the groundwater constituents due to horiz flows */
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 /* update the hydro state variables */
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); /* (horizTOsat_rec + satTOsf_rec) = fluxHoriz */
01059         
01060 /*******BUDGET
01061 **************/
01062 /*  sum the flow of groundwater and constituents (nutrients, salinity, etc) among basins, with budget calcs.*/
01063     if ( basn[cell_i] != basn[cell_j] ) { 
01064 
01065         fluxHoriz = sign*fluxHoriz*CELL_SIZE; /* signed water flux volume (m^3) */
01066         sed_stuf1fl_horz = sign*sed_stuf1fl_horz; /* signed constituent flux (kg) */
01067         sed_stuf3fl_horz = sign*sed_stuf3fl_horz; /* signed constituent flux (kg) */
01068 
01069         /* first do the normal case where all cells are on-map */
01070         if ( ON_MAP[cell_j] && ON_MAP[cell_i] ) { 
01071             if ( fluxHoriz > 0  ) { /* fluxes out of basn[cell_i] */
01072              if (basn_list[basn[cell_i]]->family !=  
01073                  basn_list[basn[cell_j]]->family ) {        /* if the flow is not within the family... */
01074                  if  ( !basn_list[basn[cell_i]]->parent  ) { /* and if the donor or recipient is a child... */
01075                      /* then find out about the flow for the family's sake */
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                      /* then find out about the flow for the family's sake */
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                      /* now sum the parents' flows */
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 {    /* if it's flow within a family, just keep
01097                           track of what the children do among themselves */            
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 { /* negative fluxes out of basn[cell_j] */
01112              if (basn_list[basn[cell_i]]->family !=  
01113                  basn_list[basn[cell_j]]->family ) {        /* if the flow is not within the family... */
01114                  if  ( !basn_list[basn[cell_j]]->parent  ) { /* and if the donor or recipient is a child... */
01115                      /* then find out about the flow for the family's sake */
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                      /* then find out about the flow for the family's sake */
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                      /* now sum the parents' flows */
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 {    /* if it's flow within a family, just keep
01136                           track of what the children do among themselves */            
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]) {/* so now the j,j cell is off-map, recipient if pos flow */
01153             if (fluxHoriz > 0) { 
01154              if ( !basn_list[basn[cell_i]]->parent ) { /* child's play */
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              /* parents' play */
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 {/* negative flows */
01168              if ( !basn_list[basn[cell_i]]->parent ) {/* child's play */
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              /* parents' play */
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]) { /* so now the i,i cell is off-map, donor if pos flow */
01184             if (fluxHoriz > 0) { 
01185              if ( !basn_list[basn[cell_j]]->parent ) {/* child's play */
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              /* parents' play */
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 {/*negative flows */
01199              if ( !basn_list[basn[cell_j]]->parent ) {/* child's play */
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              /* parents' play */
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 } /* end Flux_GWcells */
01219 
01220   
01221 
01222 
01223 /***************************************************************************************/

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