Fluxes.c

Go to the documentation of this file.

/* double precision for all constituent vars */

/* General NOTES on revisions to this source file:
        April 00 VERSION 2.1 - output available for application for water quality performance
        July 02 VERSIOIN 2.1a - first widely-available public release of code -
              - made a start at adding dynamic stage boundary conditions - INCOMPLETE, not functional
              - added some more code documentation, etc.

        July 2004 VERSION 2.2.1 - added dynamic stage boundary conditions
              - in globals (celldyn1) module of UnitMod.c, have read daily, coarse-grid boundary stage data
              - in this version, we have not turned on that new boundary condition code, still using v2.1 
                estimated boundary conditions
        Aug 2004 VERSION 2.3.0 - added dynamic boundary conditions 
                - used new spatial time series data input for dynamic boundary conditions
        Oct 2004 VERSION 2.3.1 - internal release 
                - checked that results reasonable, not fully checked
        Nov 2004 VERSION 2.3.2 - documentation upgrade
        Jan 2005 v2.4.0 - encoded the AND dispersion
        Oct 2006 v2.6.0: added more boundary condition functionality
                - added model-experiment parameters (ModExperimParms_NOM, boundary condition experiments for non-gridIO data)
            - omission found in surface water flows with external cells (e.g., boundary of Big Cypress in regional ELM), 
                where the external cell was expected to have stage (output from SFWMM), but it was provided (SFWMM) data on depth.
                (Stage was correctly derived for the groundwater flows with external cells, which is where the high-volume, 
                critical flows occur). 
        Nov 2006 v2.6.0: enhanced documentation for release
        Oct 2007 v2.7.0: added functionality of calc'ing surface water flow velocities
                - new variable of surface water flow velocities
                - fixed bug (identified in Oct '06) in auto-scaling of dispersion among model grids of diff sizes
                - fixed bug (identified in Oct '06) in boundary condition surface water flows
                

*/

#include "fluxes.h" 


/*************************************************************************************************/
void Flux_SWater(int it, float *SURFACE_WAT,float *SED_ELEV,float *HYD_MANNINGS_N, 
                 double *STUF1, double *STUF2, double *STUF3, float *SfWat_vel_day)
{ 
int ix, iy;
  float FFlux;


  
  /* check the donor and recipients cells for a) on-map, b) the cell attribute that allows sfwater
     boundary flow from the system and c) the attribute that indicates levee presence:
     the levee attribute of 1 (bitwise) allows flow to east;
     attribute of 2 (bitwise) allows flow to south (levee atts calc'd by WatMgmt.c) 
   */

  /* as always, x is row, y is column! */
  for(ix=1; ix<=s0; ix++) 
  {
    if (it%2)   /* alternate loop directions every other hyd_iter (it) */
    {
      for(iy=1; iy<=s1; iy++)  /* loop from west to east */
      {
        if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy+1)] && (int)(ON_MAP[T(ix,iy)]-1) & 1  ) || 
              BCondFlow[T(ix,iy+1)] == 3 ||  BCondFlow[T(ix,iy)] == 3  )
        {
          FFlux = Flux_SWcells(ix,iy,ix,iy+1,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N); 
          /* FFlux units = m */
          if (FFlux != 0.0) 
            Flux_SWstuff ( ix,iy,ix,iy+1,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3,SfWat_vel_day);

        } /* endof if */
      } /* end of for */
    }  /* end of if */
     

    else 
    { 
      for(iy=s1; iy>=1; iy--)   /* loop from east to west */
        if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy-1)] && (int)(ON_MAP[T(ix,iy-1)]-1) & 1 ) || 
              BCondFlow[T(ix,iy-1)] == 3 || BCondFlow[T(ix,iy)] == 3  )
        {
          FFlux = Flux_SWcells(ix,iy-1,ix,iy,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N); 
                                /* FFlux units = m */
          if (FFlux != 0.0) 
            Flux_SWstuff ( ix,iy-1,ix,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3,SfWat_vel_day);

        }  /* end of if */
    }  /* end of else */
  }
        
  for(iy=1; iy<=s1; iy++) 
  {
    if (it%2)   /* alternate loop directions every other hyd_iter (it) */
    {
      for(ix=1; ix<=s0; ix++)  /* loop from north to south */
        if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix+1,iy)]  && (int)(ON_MAP[T(ix,iy)]-1) & 2 ) ||
              BCondFlow[T(ix+1,iy)] == 3 || BCondFlow[T(ix,iy)] == 3  )
        { 
          FFlux = Flux_SWcells(ix,iy,ix+1,iy,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N); 
          /* FFlux units = m */
          if (FFlux != 0.0) 
            Flux_SWstuff ( ix,iy,ix+1,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3,SfWat_vel_day);
        }
    } 


    else 
    { 
      for(ix=s0; ix>=1; ix--)  /* loop from south to north */
        if( ( ON_MAP[T(ix,iy)] && ON_MAP[T(ix-1,iy)] && (int)(ON_MAP[T(ix-1,iy)]-1) & 2 ) ||
              BCondFlow[T(ix-1,iy)] == 3 || BCondFlow[T(ix,iy)] == 3  )
        {
             
          FFlux = Flux_SWcells(ix-1,iy,ix,iy,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N); 
          /* FFlux units = m */
          if (FFlux != 0.0) 
            Flux_SWstuff ( ix-1,iy,ix,iy,FFlux,SURFACE_WAT,STUF1,STUF2,STUF3,SfWat_vel_day);

        } /* end of if */
    } /* end of else */
  } /* end of for */

} /* end of function */




/************************************************************************************************/
float Flux_SWcells(int i0,int i1,int j0,int j1, float *SWater,float *Elevation,float *MC)
{
  float dh, adh; 
  float MC_cells;
  float Hi, Hj, stage_i, stage_j; 
  float Flux = 0.; 
  int cellLoci = T(i0,i1);
  int cellLocj = T(j0,j1);
  float garb, garb2;
  float GP_BC_SfWat_addOut; /* v2.6 added new calculated var, used only in model experiments */

   
  MC_cells = (MC[cellLoci] + MC[cellLocj] )/2.;

  /* If an on-map cell is marked 3, we are at a model boundary allowing surface water exchange */     
  if (!ON_MAP[cellLoci] && BCondFlow[cellLocj] == 3 ) 
  {
    /* v2.6, ONLY in model experiments (isModExperim), calculate synthetic boundary stages  */
        /* assume equivalent land surface elevation in cell external to model, 
           and induce a constant stage diff */
           /* v2.6 parameters found in special "ModExperimParms_NOM" file */
    if (isModExperim) { 
        if (j0 < GP_BC_InputRow) { 
                SWater[cellLoci] = SWater[cellLocj] + GP_BC_SfWat_add;
        } 
        else { 
                GP_BC_SfWat_addOut = ( SWater[cellLocj] > GP_BC_SfWat_targ ) ? 
                        ( (GP_BC_SfWat_addHi+GP_BC_SfWat_add) * ( pow( (SWater[cellLocj] / (GP_BC_SfWat_targ+0.001)) ,2.0) ) ) : 
                        (GP_BC_SfWat_add);
                SWater[cellLoci] = Max( SWater[cellLocj] - GP_BC_SfWat_addOut, 0.0);
        }
        stage_i = Elevation[cellLocj] + SWater[cellLoci];   
    }
    else
    {
        /* v2.6 - error in v2.5 found here, where Hi did not have Elevation added (boundcond_depth was originally intended to be stage, w/ stage in name).  
        Regional ELM v2.5 did not have significant overland flow across domain boundaries in the two regions where that operated (northern BCNP, Model Lands north of C-111, used to be no-flow areas for overland boundary flows). */ 
        stage_i = Elevation[cellLocj] + Max( Max(boundcond_depth[cellLoci],0.0) + GP_BC_SfWat_add,0.0);  /* v2.6 GP_BC_SfWat_add is non-zero only when performing model experiments (ModExperimParms parm file exists) */
        /*stage_i = boundcond_depth[cellLoci]; */ /* v2.5 */
    }

    Hi =  stage_i ;  /* whether model experiment or standard sim, get the value of Hi for overland flux calcs */

    MC_cells = MC[cellLocj];       /* the mannings n is not avg, but the value of onmap boundary cell */
  }

  else 
    Hi = SWater[cellLoci] + Elevation[cellLoci];

  if(BCondFlow[cellLoci] == 3 && !ON_MAP[cellLocj] )   
  {  
    /* v2.6, ONLY in model experiments (isModExperim), calculate synthetic boundary stages  */
        /* assume equivalent land surface elevation in cell external to model, 
           and induce a constant stage diff */
           /* v2.6 parameters found in special "ModExperimParms_NOM" file */
    if (isModExperim) { 
        if (i0 < GP_BC_InputRow) {
                SWater[cellLocj] = SWater[cellLoci] + GP_BC_SfWat_add;
        } 
        else { 
                GP_BC_SfWat_addOut = ( SWater[cellLoci] > GP_BC_SfWat_targ ) ? 
                        ( (GP_BC_SfWat_addHi+GP_BC_SfWat_add) * ( pow( (SWater[cellLoci] / (GP_BC_SfWat_targ+0.001)) ,2.0) ) ) : 
                        (GP_BC_SfWat_add);
                SWater[cellLocj] = Max(SWater[cellLoci] - GP_BC_SfWat_addOut, 0.0);
        }
        stage_j = Elevation[cellLoci] + SWater[cellLocj];    
    }
    else
    {
        /* v2.6 - error in v2.5 found here, where Hi did not have Elevation added (boundcond_depth was originally intended to be stage, w/ stage in name).  
        Regional ELM v2.5 did not have significant overland flow across domain boundaries in the two regions where that operated (northern BCNP, Model Lands north of C-111, used to be no-flow areas for overland boundary flows). */ 
        stage_j = Elevation[cellLoci] + Max( Max(boundcond_depth[cellLocj],0.0) + GP_BC_SfWat_add,0.0);  /* v2.6 GP_BC_SfWat_add is non-zero only when performing model experiments (ModExperimParms parm file exists) */ 
        /*stage_j = boundcond_depth[cellLocj];*/ /* v2.5 */
    }

    Hj = stage_j;  /* whether model experiment or standard sim, get the value of Hj for overland flux calcs */

    MC_cells = MC[cellLoci];      /* the mannings n is not avg, but the value of onmap boundary cell */
  }

  else 
    Hj = SWater[cellLocj] + Elevation[cellLocj];

  dh = Hi - Hj;         /* dh is "from --> to" */
  adh = Abs (dh);
                
  if (dh > 0) 
  {
    if(SWater[cellLoci] < GP_DetentZ) 
      return 0.0; 
    /* step_Cell is constant calc'd in Generic_Driver.c at initialization ( m^(-1.5) * sec )*/
    Flux = (MC_cells != 0) ? 
        (pow(adh,GP_mannHeadPow) / MC_cells  * pow(SWater[cellLoci],GP_mannDepthPow)*step_Cell) : (0.0);
//garb=Flux;
     /* ensure adequate volume avail */
    Flux =  ( Flux > ramp(SWater[cellLoci] - GP_DetentZ) ) ? (ramp(SWater[cellLoci] - GP_DetentZ)) : (Flux);

    /* check to ensure no flip/flops associated with depth */
    if ( ( Hi - Flux ) < ( Hj + Flux ) ) Flux = Min ( 0.5*dh, ramp(SWater[cellLoci] - GP_DetentZ) ); /* v2.8.4 was dh/2.0 */
//if (Flux > 0) printf("\nsfwt= %f; N= %f; dh= %f; Flux1= %f; Flux= %f; ",SWater[cellLoci],MC_cells, dh, garb, Flux); 

 } /* end of dh > 0 */

  else
  {     
    if (SWater[cellLocj] < GP_DetentZ) 
      return 0.0; 
    /* step_Cell is constant calc'd in Generic_Driver.c at initialization ( m^(-1.5) * sec )*/
    /* Flux is negative in this case */
    Flux = (MC_cells != 0) ? 
           ( - pow(adh,GP_mannHeadPow) / MC_cells  * pow(SWater[cellLocj],GP_mannDepthPow)*step_Cell) : (0.0);

    /* ensure adequate volume avail */
    Flux =  ( -Flux > ramp(SWater[cellLocj] - GP_DetentZ) ) ? (-ramp(SWater[cellLocj] - GP_DetentZ)) : (Flux);

    /* check to ensure no flip/flops associated with depth */           
    if ( ( Hi - Flux ) > ( Hj + Flux ) ) Flux = - Min ( 0.5*adh, ramp(SWater[cellLocj] - GP_DetentZ) ); /* v2.8.4 was adh/2.0 */
// if (Flux < 0) printf("\nsfwt= %f; N= %f; dh= %f; Flux= %f; ",SWater[cellLoci],MC_cells, dh, Flux); 

  }
        
  return (Flux);        /* returns height flux */
}


/*************************************************************************************************/
void Flux_SWstuff(int i0,int i1,int j0,int j1, float Flux, float *SURFACE_WAT, double *STUF1, double *STUF2, double *STUF3, float *SfWat_vel_day)
{
    float m1=0.0, m2=0.0, m3=0.0;
    int  ii, flo_chek, cel_i, cel_j;
    float FluxAdj; /* Flux adjusted for numerical dispersioin */
    double fl_prop_i, fl_prop_j; /* proportion of surface water constituents that should be fluxed with water (depends on magnitude of dispersion) */
    double sf_wt_veloc=0.0; /* v2.7.0  horizontal velocity of advected water (m/d) */

    cel_i = T(i0,i1);
    cel_j = T(j0,j1);
    
    /*sf_wt_veloc =  Abs(Flux) * celWid/( (Flux>0.0) ? (SURFACE_WAT[cel_i]) : (SURFACE_WAT[cel_j]) ) / (sfstep) ; */ 
    sf_wt_veloc = Veloc_Calc(Flux, SURFACE_WAT[cel_i], SURFACE_WAT[cel_j], sfstep);
    
    FluxAdj = Disp_Calc(Flux, SURFACE_WAT[cel_i], SURFACE_WAT[cel_j], sfstep, sf_wt_veloc);
    /* FluxAdj is the Flux that was adjusted (increase/decrease) to accomodate the targeted level of constituent dispersion */
    
    if (Flux > 0.0) {
        fl_prop_i = (SURFACE_WAT[cel_i]>0.0) ? (FluxAdj/SURFACE_WAT[cel_i]) : (0.0); /* pos Flux */
        fl_prop_i = Min(fl_prop_i, 1.0);
     }
    else {
        fl_prop_j = (SURFACE_WAT[cel_j]>0.0) ? (-FluxAdj/SURFACE_WAT[cel_j]) : (0.0); /* neg Flux */
        fl_prop_j = Max(fl_prop_j, -1.0); /* v2.6 error found? (was Min(fl_prop_j, 1.0) )*/
    }
 
 if (Flux >0.0) {  /* the i,i cell is donor if pos flow */ /* v2.6 added the boundary flux conc */
     m1 = (ON_MAP[cel_i]) ? (STUF1[cel_i]*fl_prop_i) : (Flux * CELL_SIZE * GP_BC_Sconc); 
     m2 = STUF2[cel_i]*fl_prop_i; /* inactive */
     m3 = (ON_MAP[cel_i]) ? (STUF3[cel_i]*fl_prop_i) : (Flux * CELL_SIZE * GP_BC_Pconc);
 }
 else  {   /* the j,j cell is donor if neg flow */
     m1 = (ON_MAP[cel_j]) ? (STUF1[cel_j]*fl_prop_j) : (Flux * CELL_SIZE * GP_BC_Sconc);
     m2 = STUF2[cel_j]*fl_prop_j; /* inactive */
     m3 = (ON_MAP[cel_j]) ? (STUF3[cel_j]*fl_prop_j) : (Flux * CELL_SIZE * GP_BC_Pconc);
 }
        
 if (ON_MAP[cel_j])  {
        STUF1[cel_j] += m1;  /* add the masses of constituents */
        STUF2[cel_j] += m2;
        STUF3[cel_j] += m3;
        SURFACE_WAT[cel_j] += Flux; /* now update the surfwater depths */
 }
 if (ON_MAP[cel_i])  {
        STUF1[cel_i] -= m1;
        STUF2[cel_i] -= m2;
        STUF3[cel_i] -= m3;
        SURFACE_WAT[cel_i] -= Flux;
 }
 
 /* v2.7.0 sum of the velocities during the multiple iterations within a day (will be averaged in UnitMod.c) */
 /* v2.8.4 the sum of velocities included flows into and out of cell, which, during looping, may have (?) tended to negate the desired sum
    When conditional with "if (Flux >0.0)  " to only use fluxes in a single direction, no real diff, so have left it as it was */
 SfWat_vel_day[cel_i] += sf_wt_veloc; 

 if (debug > 2) {
     if (STUF3[cel_j] < -GP_MinCheck) {
         sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water P (%f kg) in cell (%d,%d) after SWfluxes!", 
                 SimTime.TIME, STUF3[cel_j], j0,j1 ); 
         WriteMsg( msgStr,True );  dynERRORnum++;}
     if (STUF3[cel_i] < -GP_MinCheck) {
         sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water P (%f kg) in cell (%d,%d) after SWfluxes!", 
                 SimTime.TIME, STUF3[cel_i], i0,i1 ); 
         WriteMsg( msgStr,True );  dynERRORnum++; }
     if (STUF1[cel_j] < -GP_MinCheck) {
         sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water S (%f kg) in cell (%d,%d) after SWfluxes!", 
                 SimTime.TIME, STUF1[cel_j], j0,j1 ); 
         WriteMsg( msgStr,True );  dynERRORnum++; }
     if (STUF1[cel_i] < -GP_MinCheck) {
         sprintf(msgStr,"Day %6.1f: capacityERR - neg surface water S (%f kg) in cell (%d,%d) after SWfluxes!", 
                 SimTime.TIME, STUF1[cel_i], i0,i1 ); 
         WriteMsg( msgStr,True );  dynERRORnum++; }
     if (SURFACE_WAT[cel_j] < -GP_MinCheck) {
         sprintf(msgStr,"Day %6.1f: capacityERR - negative surface water (%f m) in cell (%d,%d)!", 
                 SimTime.TIME, SURFACE_WAT[cel_j], j0,j1 ) ; 
         WriteMsg( msgStr,True );  dynERRORnum++;}

     if (SURFACE_WAT[cel_i] < -GP_MinCheck) {
         sprintf(msgStr,"Day %6.1f: capacityERR - negative surface water (%f m) in cell (%d,%d)!", 
                 SimTime.TIME, SURFACE_WAT[cel_i], i0,i1 ) ; 
         WriteMsg( msgStr,True );  dynERRORnum++; }
 }
        

/* mass balance sums */
 if (basn[cel_i] != basn[cel_j])  { 

        /* first do the normal case where all cells are on-map */
     if ( ON_MAP[cel_j] && ON_MAP[cel_i] ) { 
         if ( Flux > 0  ) { /* positive fluxes out of basn[cel_i] */
             if (basn_list[basn[cel_i]]->family !=  
                 basn_list[basn[cel_j]]->family ) {        /* if the flow is not within the family... */
                 if  ( !basn_list[basn[cel_i]]->parent  ) { /* and if the donor or recipient is a child... */
                     /* then find out about the flow for the family's sake */
                     VOL_OUT_OVL[basn_list[basn[cel_i]]->family] += Flux*CELL_SIZE;
                     P_OUT_OVL[basn_list[basn[cel_i]]->family] += m3;
                     S_OUT_OVL[basn_list[basn[cel_i]]->family] += m1;
                 }
                 if ( !basn_list[basn[cel_j]]->parent ) {                 
                     /* then find out about the flow for the family's sake */
                     VOL_IN_OVL[basn_list[basn[cel_j]]->family] += Flux*CELL_SIZE;
                     P_IN_OVL[basn_list[basn[cel_j]]->family] += m3;
                     S_IN_OVL[basn_list[basn[cel_j]]->family] += m1;
                 }
                     /* now sum the parents' flows */
                 VOL_OUT_OVL[basn[cel_i]] += Flux*CELL_SIZE; 
                 P_OUT_OVL[basn[cel_i]] += m3; 
                 S_OUT_OVL[basn[cel_i]] += m1; 
                 VOL_IN_OVL[basn[cel_j]]  += Flux*CELL_SIZE; 
                 P_IN_OVL[basn[cel_j]]  += m3; 
                 S_IN_OVL[basn[cel_j]]  += m1; 
                 
             }
             else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                 if  ( !basn_list[basn[cel_i]]->parent  ) { 
                     VOL_OUT_OVL[basn[cel_i]] += Flux*CELL_SIZE; 
                     P_OUT_OVL[basn[cel_i]] += m3; 
                     S_OUT_OVL[basn[cel_i]] += m1; 
                 }
                 if ( !basn_list[basn[cel_j]]->parent ) {                 
                     VOL_IN_OVL[basn[cel_j]]  += Flux*CELL_SIZE; 
                     P_IN_OVL[basn[cel_j]]  += m3; 
                     S_IN_OVL[basn[cel_j]]  += m1; 
                 }
             }
                        
             if (debug > 0 && WatMgmtOn) { /* check for basin misconfiguration (allowable basin-basin flows) */
                 basins = basn_list[basn[cel_i]];
                 flo_chek = 0;
                 for (ii=0; ii<basins->numFLok; ii++) { if (basn[cel_j] == basins->FLok[ii] ) flo_chek = 1; }
                 if (!flo_chek) {
                     sprintf(msgStr,"Day %5.3f: ERROR - no (pos) SW flow from cell (%d,%d) of basin %d into cell (%d,%d) of basin %d!", 
                             SimTime.TIME, i0,i1,basn[cel_i], j0,j1, basn[cel_j]); 
                     WriteMsg( msgStr,True );  dynERRORnum++; }
             }
             

         }
         else { /* negative fluxes out of basn[cel_j] */
             if (basn_list[basn[cel_i]]->family !=  
                 basn_list[basn[cel_j]]->family ) {        /* if the flow is not within the family... */
                 if  ( !basn_list[basn[cel_j]]->parent  ) { /* and if the donor or recipient is a child... */
                     /* then find out about the flow for the family's sake */
                     VOL_OUT_OVL[basn_list[basn[cel_j]]->family] -= Flux*CELL_SIZE;
                     P_OUT_OVL[basn_list[basn[cel_j]]->family] -= m3;
                     S_OUT_OVL[basn_list[basn[cel_j]]->family] -= m1;
                 }
                 if ( !basn_list[basn[cel_i]]->parent ) {                 
                     /* then find out about the flow for the family's sake */
                     VOL_IN_OVL[basn_list[basn[cel_i]]->family] -= Flux*CELL_SIZE;
                     P_IN_OVL[basn_list[basn[cel_i]]->family] -= m3;
                     S_IN_OVL[basn_list[basn[cel_i]]->family] -= m1;
                 }
                     /* now sum the parents' flows */
                 VOL_OUT_OVL[basn[cel_j]] -= Flux*CELL_SIZE; 
                 P_OUT_OVL[basn[cel_j]] -= m3; 
                 S_OUT_OVL[basn[cel_j]] -= m1; 
                 VOL_IN_OVL[basn[cel_i]]  -= Flux*CELL_SIZE; 
                 P_IN_OVL[basn[cel_i]]  -= m3; 
                 S_IN_OVL[basn[cel_i]]  -= m1; 
                 
             }
             else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                 if  ( !basn_list[basn[cel_j]]->parent  ) { 
                     VOL_OUT_OVL[basn[cel_j]] -= Flux*CELL_SIZE; 
                     P_OUT_OVL[basn[cel_j]] -= m3; 
                     S_OUT_OVL[basn[cel_j]] -= m1; 
                 }
                 if ( !basn_list[basn[cel_i]]->parent ) {                 
                     VOL_IN_OVL[basn[cel_i]]  -= Flux*CELL_SIZE; 
                     P_IN_OVL[basn[cel_i]]  -= m3; 
                     S_IN_OVL[basn[cel_i]]  -= m1; 
                 }
             }


             if (debug > 0 && WatMgmtOn) { /* check for basin misconfiguration (allowable basin-basin flows) */
                 basins = basn_list[basn[cel_j]];
                 flo_chek = 0;
                 for (ii=0; ii<basins->numFLok; ii++) { if (basn[cel_i] == basins->FLok[ii] ) flo_chek = 1; }
                 if (!flo_chek) {
                     sprintf(msgStr,"Day %5.3f: ERROR - no (neg) SW flow from cell (%d,%d) of basin %d into cell (%d,%d) of basin %d!", 
                             SimTime.TIME, i0,i1, basn[cel_j], j0,j1, basn[cel_i]); 
                     WriteMsg( msgStr,True );  dynERRORnum++; }
             }
                 
             
         }
     }
     else  if ( !ON_MAP[cel_j]) /* so now the j,j cell is off-map, recipient if pos flow */
         if (Flux > 0) { 
             if ( !basn_list[basn[cel_i]]->parent ) { /* child's play */
                 VOL_OUT_OVL[basn_list[basn[cel_i]]->family] += Flux*CELL_SIZE;
                 P_OUT_OVL[basn_list[basn[cel_i]]->family] += m3;
                 S_OUT_OVL[basn_list[basn[cel_i]]->family] += m1;
             }
             /* parents' play */
             VOL_OUT_OVL[basn[cel_i]] += Flux*CELL_SIZE; 
             VOL_OUT_OVL[0]+= Flux*CELL_SIZE;
             P_OUT_OVL[basn[cel_i]] += m3; 
             P_OUT_OVL[0]+= m3;
             S_OUT_OVL[basn[cel_i]] += m1; 
             S_OUT_OVL[0]+= m1;
         }
         else { /* negative flows */
             if ( !basn_list[basn[cel_i]]->parent ) {/* child's play */
                 VOL_IN_OVL[basn_list[basn[cel_i]]->family] -= Flux*CELL_SIZE;
                 P_IN_OVL[basn_list[basn[cel_i]]->family] -= m3;
                 S_IN_OVL[basn_list[basn[cel_i]]->family] -= m1;
             }
             /* parents' play */
             VOL_IN_OVL[basn[cel_i]]  -= Flux*CELL_SIZE;
             VOL_IN_OVL[0]               -= Flux*CELL_SIZE;
             P_IN_OVL[basn[cel_i]]  -= m3;
             P_IN_OVL[0]               -= m3;
             S_IN_OVL[basn[cel_i]]  -= m1;
             S_IN_OVL[0]               -= m1;
         }
     else  if ( !ON_MAP[cel_i]) /* so now the i,i cell is off-map, donor if pos flow */
         if (Flux > 0) { 
             if ( !basn_list[basn[cel_j]]->parent ) {/* child's play */
                 VOL_IN_OVL[basn_list[basn[cel_j]]->family] += Flux*CELL_SIZE;
                 P_IN_OVL[basn_list[basn[cel_j]]->family] += m3;
                 S_IN_OVL[basn_list[basn[cel_j]]->family] += m1;
             }
             /* parents' play */
             VOL_IN_OVL[basn[cel_j]] += Flux*CELL_SIZE;
             VOL_IN_OVL[0]              += Flux*CELL_SIZE;
             P_IN_OVL[basn[cel_j]] += m3;
             P_IN_OVL[0]              += m3;
             S_IN_OVL[basn[cel_j]] += m1;
             S_IN_OVL[0]              += m1;
         }
         else { /*negative flows */
             if ( !basn_list[basn[cel_j]]->parent ) {/* child's play */
                 VOL_OUT_OVL[basn_list[basn[cel_j]]->family] -= Flux*CELL_SIZE;
                 P_OUT_OVL[basn_list[basn[cel_j]]->family] -= m3;
                 S_OUT_OVL[basn_list[basn[cel_j]]->family] -= m1;
             }
             /* parents' play */
             VOL_OUT_OVL[basn[cel_j]]  -= Flux*CELL_SIZE;
             VOL_OUT_OVL[0]               -= Flux*CELL_SIZE;
             P_OUT_OVL[basn[cel_j]]  -= m3;
             P_OUT_OVL[0]               -= m3;
             S_OUT_OVL[basn[cel_j]]  -= m1;
             S_OUT_OVL[0]               -= m1;
         }
 }
                
 return;
        
}

/*************************************************************************************************/
double Veloc_Calc(float flux, float depth_i, float depth_j, float tim_step)
{
        double sf_wt_veloc_mag;
        
        sf_wt_veloc_mag =  Abs(flux) * celWid/( (flux>0.0) ? (depth_i) : (depth_j) ) / (tim_step) ;
//if (flux > 0) printf("depth= %f; Flux= %f; Velo= %f",depth_i,flux,sf_wt_veloc_mag);
        return (sf_wt_veloc_mag);
}

/*************************************************************************************************/
float Disp_Calc(float flux, float depth_i, float depth_j, float tim_step, double sf_wt_veloc)
{
    float disp_num; /* numerical dispersion in current model grid (m2/d) */
    float disp_num_targ; /*  targeted numerical dispersion (associated w/ grid size of the targeted dispersion length parameter) (m2/d) */
    float veloc_adj;  /* horizontal velocity associated w/ constituent flow, adjusted for the targeted numerical dispersion (m/d) */
    float flux_adj; /* flux (height) associated w/ constituents, adjusted for numerical dispersioin */
    extern float GP_dispLenRef, GP_dispParm; /* from unitmod.h */

    disp_num = 0.5 * sf_wt_veloc * (celWid - sf_wt_veloc * tim_step)  ; 
    disp_num_targ = 0.5 * sf_wt_veloc * (GP_dispLenRef - sf_wt_veloc * tim_step)  ; 
    veloc_adj = GP_dispParm * (sf_wt_veloc * celWid - disp_num + disp_num_targ )/celWid; 
    flux_adj = veloc_adj * tim_step * ( (flux>0.0) ? (depth_i) : (depth_j) )/celWid;
    return (flux_adj); /* return the height-flux that was adjusted (decreased/increased) due to dispersion */
}


/*************************************************************************************************/
void Flux_GWater(int it, float *SatWat, float *Unsat, float *SfWat,
                  float *rate, float *poros, float *sp_yield, float *elev,
                  double *gwSTUF1, double *gwSTUF2, double *gwSTUF3,
                  double *swSTUF1, double *swSTUF2, double *swSTUF3)

{ 
  int ix, iy; /* ix is row, iy is col! */
/* we only check the donor cell for on-map, allowing losses from the system */
 for(ix=1; ix<=s0; ix++) 
 {
     if (it%2) {  /* alternate loop directions every other iteration (it = int(hyd_iter-1)/2 ) */
     for(iy=1; iy<=s1; iy++)  /* loop from west to east */
             /* allow boundary flow if donor cell is marked 3 or higher (v2.5 had specific markings of 4 or 9, no longer pertinent in v2.6+) */
             if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy+1)]) || 
                ((BCondFlow[T(ix,iy+1)]) > 2) || ((BCondFlow[T(ix,iy)]) > 2) ) 
         {
             Flux_GWcells(ix,iy,ix,iy+1,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
                          gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3); 
         }
     } 
             
     else { 
     for(iy=s1; iy>=1; iy--)   /* loop from east to west */
             /* allow boundary flow if donor cell is marked 3 or higher  (v2.5 had specific markings of 4 or 9, no longer pertinent in v2.6+) */
         if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix,iy-1)]) || 
                ((BCondFlow[T(ix,iy-1)]) > 2) || ((BCondFlow[T(ix,iy)]) > 2) ) 
         {
         
             Flux_GWcells(ix,iy,ix,iy-1,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
                          gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3); 
         }
     } 
 }
        
 for(iy=1; iy<=s1; iy++) 
 {
     if (it%2) {  /* alternate loop directions every other iteration (it = int(hyd_iter-1)/2 ) */
     for(ix=1; ix<=s0; ix++)  /* loop from north to south */
             /* allow boundary flow if donor cell is marked 3 or higher  (v2.5 had specific markings of 4 or 9, no longer pertinent in v2.6+) */
         if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix+1,iy)]) || 
                ((BCondFlow[T(ix+1,iy)]) > 2) || ((BCondFlow[T(ix,iy)]) > 2) ) 
         {
         
             Flux_GWcells(ix,iy,ix+1,iy,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
                          gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3); 
         }
         } 
     else { 
     for(ix=s0; ix>=1; ix--)  /* loop from south to north */
             /* allow boundary flow if donor cell is marked 3 or higher  (v2.5 had specific markings of 4 or 9, no longer pertinent in v2.6+) */
         if (( ON_MAP[T(ix,iy)] && ON_MAP[T(ix-1,iy)]) || 
                ((BCondFlow[T(ix-1,iy)]) > 2) || ((BCondFlow[T(ix,iy)]) > 2) ) 
         {
         
             Flux_GWcells(ix,iy,ix-1,iy,SatWat, Unsat, SfWat, rate, poros, sp_yield, elev,
                          gwSTUF1, gwSTUF2, gwSTUF3, swSTUF1, swSTUF2, swSTUF3); 
         }
     } 
 }

}

/************************************************************************************************/
void Flux_GWcells( int i0, int i1, int j0, int j1, float *SatWat, float *Unsat, float *SfWat, 
                    float *rate, float *poros, float *sp_yield, float *elev,
                    double *gwSTUF1, double *gwSTUF2, double *gwSTUF3,
                    double *swSTUF1, double *swSTUF2, double *swSTUF3)
{       
            
    int cell_don, cell_rec, cell_i, cell_j;
    int sign, equilib, revers=1;
    float GW_head_i,GW_head_j,tot_head_i,tot_head_j;
    float dh;
    float AvgRate;

    float H_pot_don, H_pot_rec;
    
    double Flux;
    double fluxTOunsat_don;
    double fluxHoriz; 

    float Sat_vs_unsat;

    double UnsatZ_don, UnsatCap_don, UnsatPot_don; 
    double UnsatZ_rec, UnsatCap_rec, UnsatPot_rec; 
    double sfTOsat_don, unsatTOsat_don, unsatTOsat_rec, satTOsf_rec, horizTOsat_rec; 
    
    double sedwat_don, sedwat_rec;
    double sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don; 
    double sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz; 
    double sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec; 

    
    cell_i = T(i0,i1); 
    cell_j = T(j0,j1);
    
    if (ON_MAP[cell_i] ) {
        GW_head_i = SatWat[cell_i] / poros[HAB[cell_i]]; 
        tot_head_i = GW_head_i + SfWat[cell_i];
    }
    
    if (ON_MAP[cell_j] ) {
        GW_head_j = SatWat[cell_j] / poros[HAB[cell_j]]; 
        tot_head_j = GW_head_j + SfWat[cell_j];
    }
    
    if (!ON_MAP[cell_i] ) { /* for an external i0,i1 cell */
        poros[HAB[cell_i]] = poros[HAB[cell_j]]; /* other variables = donor cell or zero */
        sp_yield[HAB[cell_i]] = sp_yield[HAB[cell_j]]; 
        elev[cell_i] = elev[cell_j];
        rate[cell_i] = rate[cell_j];
        Unsat[cell_i] = 0.0; /* no water in any potential unsat zone */

            /* outflow of groundwater from boundary cell*/
        if (BCondFlow[cell_j] > 2)  { /* v2.6  (v2.5 had specific markings of 4 or 9, no longer pertinent in v2.6+)  */
            if (boundcond_depth[cell_i] > 0.0) /* boundcond_depth= positive or negative water depth relative to land surface elev */
            {
                SatWat[cell_i] =  elev[cell_i] * poros[HAB[cell_i]];    
                SfWat[cell_i] = boundcond_depth[cell_i]; 
            }
        else
            {
                SatWat[cell_i] =  (elev[cell_i] + boundcond_depth[cell_i]) * poros[HAB[cell_i]];  
                SfWat[cell_i] = 0.0; 
            }
        }
    
    GW_head_i = SatWat[cell_i] / poros[HAB[cell_i]]; 
    tot_head_i = GW_head_i + SfWat[cell_i];
    
            /* v2.6 when we don't have data on boundary condition heads */
    if ( isModExperim ) { 
        if (j0 < GP_BC_InputRow) {
                tot_head_i = tot_head_j /* + GP_BC_SatWat_add */ ;
        } 
        else {
                tot_head_i = (SfWat[cell_j] > GP_BC_SatWat_add) ? (tot_head_j - GP_BC_SatWat_add) : (tot_head_j);
        }
    } /* end of isModExperim */
    
    }

    if (!ON_MAP[cell_j] ) { /* for an external j0,j1 cell */
        poros[HAB[cell_j]] = poros[HAB[cell_i]]; /* other variables = donor cell or zero */
        sp_yield[HAB[cell_j]] = sp_yield[HAB[cell_i]]; 
        elev[cell_j] = elev[cell_i];
        rate[cell_j] = rate[cell_i];
        Unsat[cell_j] = 0.0; /* no water in any potential unsat zone */
            
            /* outflow of groundwater from boundary cell*/
        if ( BCondFlow[cell_i] > 2 )  { /* v2.6  (v2.5 had specific markings of 4 or 9, no longer pertinent in v2.6+)  */
            if (boundcond_depth[cell_j] > 0.0) /* boundcond_depth=relative water depth; v2.3tmp instead of: boundcond_depth[cell_j] > elev[cell_j]  */
            {
                SatWat[cell_j] =  elev[cell_j] * poros[HAB[cell_j]];    
                SfWat[cell_j] = boundcond_depth[cell_j]; /* boundcond_depth=relative water depth; v2.3tmp instead of: boundcond_depth[cell_j] - elev[cell_j]  */
            }
        else
            {
                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]]  */  
                SfWat[cell_j] = 0.0; 
            }
        }
            
    GW_head_j = SatWat[cell_j] / poros[HAB[cell_j]]; 
    tot_head_j = GW_head_j + SfWat[cell_j];

            /* v2.6 when we don't have data on boundary condition heads */
    if ( isModExperim ) {
        if (j0 < GP_BC_InputRow) { 
                tot_head_j = tot_head_i /* + GP_BC_SatWat_add */;
        } 
        else { 
                tot_head_j = (SfWat[cell_i] > GP_BC_SatWat_add) ? (tot_head_i - GP_BC_SatWat_add) : (tot_head_i);
        }
    } /* end of isModExperim */
    
    }

    AvgRate = ( rate[cell_i] + rate[cell_j] )/2.0; /* average hyd conductivity of the two cells */
    
/* hydraulic head difference, positive flux from cell_i to cell_j */
    dh = tot_head_i - tot_head_j; 

 
        /* determine donor and recipient cells based on head difference,
           and provide the sign of the flux;
           The surface water detention depth also used here as threshold
           head difference for fluxing (currently global, 1 cm)*/
    if (dh > GP_DetentZ) 
    {
        cell_don=cell_i;
        cell_rec=cell_j;
        sign = 1; 
    }
    else if (dh < -GP_DetentZ) 
    {   
        cell_don=cell_j;
        cell_rec=cell_i;
        sign = -1;
    }
    else 
        return; /* no flux (or surface-groundwater integration) under minimal head diffs */
    

/* Potential horizontal flux eqn (Darcy's eqn simplified to square cells).
   This is the maximum height (m) of water vol to flux under fully saturated condition.
   Note that the Min check is probably not important under current implementations, as SatWat includes water
   down to the model's base datum, which is below the land elevation vertical datum (NGVD29 or NAVD88) 
   by the value of GP_DATUM_DISTANCE (has always been 6 m as of v2.8)  */
    Flux = Min(Abs(dh) * AvgRate * SatWat[cell_don] / CELL_SIZE * gwstep , SatWat[cell_don]);
 
/**** this do-while routine (1) integrates the surface, saturated, and unsaturated water,
  and (2) checks to ensure the heads do not reverse in a time step due to large fluxes.
  If heads do reverse, the total Flux is decremented until there is no reversal */
/****/
    do {
        /* The total potential flux is apportioned to (1) the horizontal component
           that fluxes to an adjacent cell and (2) the vertical component that
           remains in the donor cell after the horizontal outflow from a donor cell.
           Thus, an unsaturated zone is created, or increased in size, with loss of
           saturated water from the donor cell; this lateral gravitational flow leaves behind
           the field capacity moisture in an unsat zone. (If donor-cell surface water is present,
           it potentially will replace the unsaturated soil capacity within the same time
           step in this routine).*/
        fluxTOunsat_don = Flux/poros[HAB[cell_don]] * (poros[HAB[cell_don]] - sp_yield[HAB[cell_don]])  ;
        fluxHoriz = Flux - fluxTOunsat_don;

/**** given the total potential groundwater Flux, get the donor cell's
      new **post-flux** capacities */
        UnsatZ_don  =   elev[cell_don] - (SatWat[cell_don]-fluxHoriz /*Flux*/) / poros[HAB[cell_don]] ; /* unsat zone depth */
    
    /* ?v2.2 this check was against "-0.001", increased here to "-0.01", as it was only showing several mm capacityERR
        due to canal interactions (i.e., only when watmgmt turned on), and we haven't explicitly added sf-ground 
        integration to grid cells in canal code */
        if (debug>3 && (UnsatZ_don < -0.01 ) ) { 
            sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in donor cell (%d,%d) in GWfluxes!", 
                    SimTime.TIME, UnsatZ_don, i0,i1 ); WriteMsg( msgStr,True ); dynERRORnum++; }
            
        UnsatCap_don =  UnsatZ_don * poros[HAB[cell_don]]; /* maximum pore space capacity (UnsatCap)
                                                              in the depth (UnsatZ) of new unsaturated zone */
        UnsatPot_don  = UnsatCap_don - (Unsat[cell_don]+fluxTOunsat_don); /* (height of) the volume of pore space
                                                                             (soil "removed") that is unoccupied in the unsat zone */
            /* determining the pathway of flow (to sat vs. unsat) of surface water depending on depth
               of an unsat zone relative to the surface water.  With a relatively deep unsat zone, this
               downflow tends to zero (infiltration occurs within the vertical hydrology module of UnitMod.c) */ 
        Sat_vs_unsat  = 1/Exp(100.0*Max((SfWat[cell_don]-UnsatZ_don),0.0)); 

    /* sf-unsat-sat fluxes */
        if (SfWat[cell_don] > 0.0) { /* surface water downflow is assumed to be as fast
                                        as horizontal groundwater outflows */
            sfTOsat_don  = 
                ( (1.0-Sat_vs_unsat)*UnsatPot_don>SfWat[cell_don] ) ? 
                ( SfWat[cell_don] ) : 
                ( (1.0-Sat_vs_unsat)*UnsatPot_don); 
           /* with downflow of surface water into an unsat zone, the proportion of that height
              that is made into saturated storage is allocated to the sat storage variable */
            if (sfTOsat_don >=  UnsatPot_don) {
                    sfTOsat_don = UnsatPot_don ; /* downflow constrained to the avail soil capacity */
                    unsatTOsat_don = Unsat[cell_don]; /* allocate unsat to sat in this case */
                }
            else { 
                unsatTOsat_don = (UnsatZ_don > 0.0) ?
                    ( (sfTOsat_don/poros[HAB[cell_don]] ) / UnsatZ_don * Unsat[cell_don] ) :
                    (0.0); /*  allocate to saturated storage whatever proportion of
                               unsat zone that is now saturated by sfwat downflow */
            }
        }
        else { /* w/o surface water, these vertical fluxes are zero */
            sfTOsat_don = unsatTOsat_don = 0.0;
        }
/**** potential new head in donor cell will be ... */
        H_pot_don = (SatWat[cell_don] - fluxTOunsat_don - fluxHoriz + sfTOsat_don + unsatTOsat_don )
            / poros[HAB[cell_don]] +(SfWat[cell_don] - sfTOsat_don) ; 
                

/**** recipient cell's **pre-flux** capacities */
        UnsatZ_rec  =   elev[cell_rec] - SatWat[cell_rec] / poros[HAB[cell_rec]] ; /* unsat zone depth */
    
    /* ?v2.2 this check was against "-0.001", increased here to "-0.01", as it was only showing several mm capacityERR
        due to canal interactions (i.e., only when watmgmt turned on), and we haven't explicitly added sf-ground 
        integration to grid cells in canal code */
        if (debug>3 && (UnsatZ_rec < -0.01 ) ) {
            sprintf(msgStr,"Day %6.1f: capacityERR - neg unsat depth (%f m) in recipient cell (%d,%d) in GWfluxes!", 
                    SimTime.TIME, UnsatZ_rec, i0,i1 ); WriteMsg( msgStr,True );  dynERRORnum++; }
            
        UnsatCap_rec =  UnsatZ_rec * poros[HAB[cell_rec]]; /* maximum pore space capacity (UnsatCap)
                                                              in the depth (UnsatZ) of new unsaturated zone */
        UnsatPot_rec  = UnsatCap_rec - Unsat[cell_rec]; /* (height of) the volume of pore space (soil "removed")
                                                           that is unoccupied in the unsat zone */
     /* sf-unsat-sat fluxes */
        horizTOsat_rec = fluxHoriz; /* lateral inflow to soil into saturated storage */
        satTOsf_rec = Max(fluxHoriz - UnsatPot_rec, 0.0); /* upwelling beyond soil capacity */
        unsatTOsat_rec = (UnsatZ_rec > 0.0) ? /* incorporation of unsat moisture into sat storage with
                                                 rising water table due to horiz inflow */
            ( ((horizTOsat_rec-satTOsf_rec)/poros[HAB[cell_rec]] ) / UnsatZ_rec * Unsat[cell_rec] ) :
            (0.0);
/**** potential new head in recipient cell will be ... */
        H_pot_rec = (SatWat[cell_rec] + horizTOsat_rec + unsatTOsat_rec - satTOsf_rec)
            / poros[HAB[cell_rec]] + (SfWat[cell_rec] + satTOsf_rec) ; 

/**** check for a head reversal - if so, reduce flux to achieve equilibrium (donor >= recip) */
        if ( (Flux != 0.0) && ((H_pot_rec - H_pot_don) > GP_MinCheck) ) {
            revers++;
            Flux *= 0.90;
            equilib = 0; 
               /* if the unsat water storage causes an unfixable (very rare(/never?)) head reversal due to the
                   sfwat and gwater integration alone, set the Flux to zero, then allow the vertical
                   sf/ground integration to be done with no horiz flux and be done with it */
            if (revers>1) { /* v2.3 temp, was >4 */
                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. ", 
                                       SimTime.TIME, (H_pot_don - H_pot_rec),  i0,i1,Flux*sign ); WriteMsg( msgStr,True );}
                Flux =  0.0;
                }
        }
        else {
            equilib = 1;
        }
        
    } while  (!equilib); /* end of routine for integrating groundwater-surface water and preventing head reversals */
        
/**********
  finished calc'ing the water fluxes
  *************/
/* calc the flux of the constituents between sed/soil across cells in horiz dir,
 but don't update the state vars until the conc's for vertical flows have been calc'd */
    sedwat_don = SatWat[cell_don] + Unsat[cell_don];
    sed_stuf1fl_horz = (sedwat_don > 0.0) ? (gwSTUF1[cell_don] / sedwat_don*fluxHoriz) : (0.0);
    sed_stuf2fl_horz = (sedwat_don > 0.0) ? (gwSTUF2[cell_don] / sedwat_don*fluxHoriz) : (0.0);
    sed_stuf3fl_horz = (sedwat_don > 0.0) ? (gwSTUF3[cell_don] / sedwat_don*fluxHoriz) : (0.0);

/* pass along the constituents between sed/soil and surface water in vertical dir;
   for "simplicity", this does not account for the phosphorus active-zone conc. gradient,
   but assumes the homogeneity of conc within the entire soil column */
    if (sfTOsat_don > 0.0) {
        sf_stuf1fl_don = (SfWat[cell_don] > 0.0) ? (swSTUF1[cell_don] / SfWat[cell_don]*sfTOsat_don) : (0.0);
        sf_stuf2fl_don = (SfWat[cell_don] > 0.0) ? (swSTUF2[cell_don] / SfWat[cell_don]*sfTOsat_don) : (0.0);
        sf_stuf3fl_don = (SfWat[cell_don] > 0.0) ? (swSTUF3[cell_don] / SfWat[cell_don]*sfTOsat_don) : (0.0);
        gwSTUF1[cell_don] += sf_stuf1fl_don;
        gwSTUF2[cell_don] += sf_stuf2fl_don;
        gwSTUF3[cell_don] += sf_stuf3fl_don;
        swSTUF1[cell_don] -= sf_stuf1fl_don;
        swSTUF2[cell_don] -= sf_stuf2fl_don;
        swSTUF3[cell_don] -= sf_stuf3fl_don;
        
    }
    if (satTOsf_rec > 0.0) {
        sedwat_rec = SatWat[cell_rec] + Unsat[cell_rec];
        sed_stuf1fl_rec = (sedwat_rec > 0.0) ? (gwSTUF1[cell_rec] / sedwat_rec*satTOsf_rec) : (0.0);
        sed_stuf2fl_rec = (sedwat_rec > 0.0) ? (gwSTUF2[cell_rec] / sedwat_rec*satTOsf_rec) : (0.0);
        sed_stuf3fl_rec = (sedwat_rec > 0.0) ? (gwSTUF3[cell_rec] / sedwat_rec*satTOsf_rec) : (0.0);
        gwSTUF1[cell_rec] -= sed_stuf1fl_rec;
        gwSTUF2[cell_rec] -= sed_stuf2fl_rec;
        gwSTUF3[cell_rec] -= sed_stuf3fl_rec;
        swSTUF1[cell_rec] += sed_stuf1fl_rec;
        swSTUF2[cell_rec] += sed_stuf2fl_rec;
        swSTUF3[cell_rec] += sed_stuf3fl_rec;
    }
    
/* update the groundwater constituents due to horiz flows */
    gwSTUF1[cell_don] -= sed_stuf1fl_horz;
    gwSTUF2[cell_don] -= sed_stuf2fl_horz;
    gwSTUF3[cell_don] -= sed_stuf3fl_horz;
    gwSTUF1[cell_rec] += sed_stuf1fl_horz;
    gwSTUF2[cell_rec] += sed_stuf2fl_horz;
    gwSTUF3[cell_rec] += sed_stuf3fl_horz;
        

/* update the hydro state variables */
    SfWat[cell_don]  += (-sfTOsat_don);
    Unsat[cell_don]  += ( fluxTOunsat_don - unsatTOsat_don) ;
    SatWat[cell_don] += (sfTOsat_don + unsatTOsat_don - fluxTOunsat_don - fluxHoriz);
    SfWat[cell_rec]  += ( satTOsf_rec); 
    Unsat[cell_rec]  += (-unsatTOsat_rec);
    SatWat[cell_rec] += (horizTOsat_rec + unsatTOsat_rec - satTOsf_rec); /* (horizTOsat_rec + satTOsf_rec) = fluxHoriz */
        
/*******BUDGET
**************/
/*  sum the flow of groundwater and constituents (nutrients, salinity, etc) among basins, with budget calcs.*/
    if ( basn[cell_i] != basn[cell_j] ) { 

        fluxHoriz = sign*fluxHoriz*CELL_SIZE; /* signed water flux volume (m^3) */
        sed_stuf1fl_horz = sign*sed_stuf1fl_horz; /* signed constituent flux (kg) */
        sed_stuf3fl_horz = sign*sed_stuf3fl_horz; /* signed constituent flux (kg) */

        /* first do the normal case where all cells are on-map */
        if ( ON_MAP[cell_j] && ON_MAP[cell_i] ) { 
            if ( fluxHoriz > 0  ) { /* fluxes out of basn[cell_i] */
             if (basn_list[basn[cell_i]]->family !=  
                 basn_list[basn[cell_j]]->family ) {        /* if the flow is not within the family... */
                 if  ( !basn_list[basn[cell_i]]->parent  ) { /* and if the donor or recipient is a child... */
                     /* then find out about the flow for the family's sake */
                     VOL_OUT_GW[basn_list[basn[cell_i]]->family] += fluxHoriz;
                     P_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf3fl_horz;
                     S_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf1fl_horz;
                 }
                 if ( !basn_list[basn[cell_j]]->parent ) {                 
                     /* then find out about the flow for the family's sake */
                     VOL_IN_GW[basn_list[basn[cell_j]]->family] += fluxHoriz;
                     P_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf3fl_horz;
                     S_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf1fl_horz;
                 }
                     /* now sum the parents' flows */
                 VOL_OUT_GW[basn[cell_i]] += fluxHoriz;
                 VOL_IN_GW[basn[cell_j]]  += fluxHoriz;
                 P_OUT_GW[basn[cell_i]] += sed_stuf3fl_horz;
                 P_IN_GW[basn[cell_j]]  += sed_stuf3fl_horz;
                 S_OUT_GW[basn[cell_i]] += sed_stuf1fl_horz;
                 S_IN_GW[basn[cell_j]]  += sed_stuf1fl_horz;
             }
             
                        
             else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                 if  ( !basn_list[basn[cell_i]]->parent  ) { 
                     VOL_OUT_GW[basn[cell_i]] += fluxHoriz; 
                     P_OUT_GW[basn[cell_i]] += sed_stuf3fl_horz; 
                     S_OUT_GW[basn[cell_i]] += sed_stuf1fl_horz; 
                 }
                 if ( !basn_list[basn[cell_j]]->parent ) {                 
                     VOL_IN_GW[basn[cell_j]]  += fluxHoriz; 
                     P_IN_GW[basn[cell_j]]  += sed_stuf3fl_horz; 
                     S_IN_GW[basn[cell_j]]  += sed_stuf1fl_horz; 
                 }
             }

            }
            else { /* negative fluxes out of basn[cell_j] */
             if (basn_list[basn[cell_i]]->family !=  
                 basn_list[basn[cell_j]]->family ) {        /* if the flow is not within the family... */
                 if  ( !basn_list[basn[cell_j]]->parent  ) { /* and if the donor or recipient is a child... */
                     /* then find out about the flow for the family's sake */
                     VOL_OUT_GW[basn_list[basn[cell_j]]->family] -= fluxHoriz;
                     P_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf3fl_horz;
                     S_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf1fl_horz;
                 }
                 if ( !basn_list[basn[cell_i]]->parent ) {                 
                     /* then find out about the flow for the family's sake */
                     VOL_IN_GW[basn_list[basn[cell_i]]->family] -= fluxHoriz;
                     P_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf3fl_horz;
                     S_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf1fl_horz;
                 }
                     /* now sum the parents' flows */
                VOL_OUT_GW[basn[cell_j]] -= fluxHoriz;
                VOL_IN_GW[basn[cell_i]]  -= fluxHoriz;
                P_OUT_GW[basn[cell_j]] -= sed_stuf3fl_horz;
                P_IN_GW[basn[cell_i]]  -= sed_stuf3fl_horz;
                S_OUT_GW[basn[cell_j]] -= sed_stuf1fl_horz;
                S_IN_GW[basn[cell_i]]  -= sed_stuf1fl_horz;

            }
             else {    /* if it's flow within a family, just keep
                          track of what the children do among themselves */            
                 if  ( !basn_list[basn[cell_j]]->parent  ) { 
                     VOL_OUT_GW[basn[cell_j]] -= fluxHoriz; 
                     P_OUT_GW[basn[cell_j]] -= sed_stuf3fl_horz; 
                     S_OUT_GW[basn[cell_j]] -= sed_stuf1fl_horz; 
                 }
                 if ( !basn_list[basn[cell_i]]->parent ) {                 
                     VOL_IN_GW[basn[cell_i]]  -= fluxHoriz; 
                     P_IN_GW[basn[cell_i]]  -= sed_stuf3fl_horz; 
                     S_IN_GW[basn[cell_i]]  -= sed_stuf1fl_horz; 
                 }
             }
            }
            

        }
        else  if ( !ON_MAP[cell_j]) {/* so now the j,j cell is off-map, recipient if pos flow */
            if (fluxHoriz > 0) { 
             if ( !basn_list[basn[cell_i]]->parent ) { /* child's play */
                 VOL_OUT_GW[basn_list[basn[cell_i]]->family] += fluxHoriz;
                 P_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf3fl_horz;
                 S_OUT_GW[basn_list[basn[cell_i]]->family] += sed_stuf1fl_horz;
             }
             /* parents' play */
                VOL_OUT_GW[basn[cell_i]] += fluxHoriz;
                VOL_OUT_GW[0]              += fluxHoriz;
                P_OUT_GW[basn[cell_i]] += sed_stuf3fl_horz;
                P_OUT_GW[0]              += sed_stuf3fl_horz;
                S_OUT_GW[basn[cell_i]] += sed_stuf1fl_horz;
                S_OUT_GW[0]              += sed_stuf1fl_horz;
            }
            else {/* negative flows */
             if ( !basn_list[basn[cell_i]]->parent ) {/* child's play */
                 VOL_IN_GW[basn_list[basn[cell_i]]->family] -= fluxHoriz;
                 P_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf3fl_horz;
                 S_IN_GW[basn_list[basn[cell_i]]->family] -= sed_stuf1fl_horz;
             }
             /* parents' play */
                VOL_IN_GW[basn[cell_i]]  -= fluxHoriz;
                VOL_IN_GW[0]               -= fluxHoriz;
                P_IN_GW[basn[cell_i]]  -= sed_stuf3fl_horz;
                P_IN_GW[0]               -= sed_stuf3fl_horz;
                S_IN_GW[basn[cell_i]]  -= sed_stuf1fl_horz;
                S_IN_GW[0]               -= sed_stuf1fl_horz;
            }
        }
        
        else  if ( !ON_MAP[cell_i]) { /* so now the i,i cell is off-map, donor if pos flow */
            if (fluxHoriz > 0) { 
             if ( !basn_list[basn[cell_j]]->parent ) {/* child's play */
                 VOL_IN_GW[basn_list[basn[cell_j]]->family] += fluxHoriz;
                 P_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf3fl_horz;
                 S_IN_GW[basn_list[basn[cell_j]]->family] += sed_stuf1fl_horz;
             }
             /* parents' play */
                VOL_IN_GW[basn[cell_j]] += fluxHoriz;
                VOL_IN_GW[0]              += fluxHoriz;
                P_IN_GW[basn[cell_j]] += sed_stuf3fl_horz;
                P_IN_GW[0]              += sed_stuf3fl_horz;
                S_IN_GW[basn[cell_j]] += sed_stuf1fl_horz;
                S_IN_GW[0]              += sed_stuf1fl_horz;
            }
            else {/*negative flows */
             if ( !basn_list[basn[cell_j]]->parent ) {/* child's play */
                 VOL_OUT_GW[basn_list[basn[cell_j]]->family] -= fluxHoriz;
                 P_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf3fl_horz;
                 S_OUT_GW[basn_list[basn[cell_j]]->family] -= sed_stuf1fl_horz;
             }
             /* parents' play */
                VOL_OUT_GW[basn[cell_j]]  -= fluxHoriz;
                VOL_OUT_GW[0]               -= fluxHoriz;
                P_OUT_GW[basn[cell_j]]  -= sed_stuf3fl_horz;
                P_OUT_GW[0]               -= sed_stuf3fl_horz;
                S_OUT_GW[basn[cell_j]]  -= sed_stuf1fl_horz;
                S_OUT_GW[0]               -= sed_stuf1fl_horz;
            }
    }
}


    return ; 
        
} /* end Flux_GWcells */

  


/***************************************************************************************/