Fluxes.c File Reference

Dynamic equations: Horizontal solutions of spatial raster fluxes. More...

#include "fluxes.h"

Include dependency graph for Fluxes.c:

Go to the source code of this file.

Functions
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)
	Surface water horizontal flux.
float	Flux_SWcells (int i0, int i1, int j0, int j1, float *SWater, float *Elevation, float *MC)
	Surface water flux calculations.
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)
	Flux of surface water constituents.
double	Veloc_Calc (float flux, float depth_i, float depth_j, float tim_step)
	Calculate horizontal flow velocities of surface water.
float	Disp_Calc (float flux, float depth_i, float depth_j, float tim_step, double sf_wt_veloc)
	Calculate dispersion for constituent fluxes.
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)
	Groundwater fluxing routine.
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)
	Ground water flux eqns, incl. integration of groundwater with surface water.

---

Detailed Description

Dynamic equations: Horizontal solutions of spatial raster fluxes.

This set of modules contains the cell-cell surface&ground water flux functions, and the integration of surface, unsat, and saturated water in the vertical dimension. The hydrologic and constituent budgets for basins/IRs are calculated for the fluxes.

The major functions in the sequence in which they are found in this file:
Flux_SWater (alternating direction explicit (ADE) calling function - surface water)
Flux_SWcells (calculate potential and actual fluxes among cells - surface water)
Flux_SWstuff (flux constituents associated with water fluxes - surface water)
Flux_GWater (alternating direction explicit (ADE) calling function - ground water)
Flux_GWcells (calculate potential and actual fluxes among cells - ground water; includes vertical surface-ground water integration, constituent fluxes)

Note: documented with Doxygen, which expects specific syntax within special comments.

The Everglades Landscape Model (ELM).
last updated: Mar 2008

Definition in file Fluxes.c.

---

Function Documentation

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
	)

Surface water horizontal flux.

This is the alternating direction function that first fluxes water in the E/W direction and then in the N/S direction. It sweeps from left to right, then goes back from right to left, then goes from top to bottom and finally from bottom to top. This alternates every hyd_iter.

Surface water model boundary exchanges may only occur across cells designated with designated attribute in the boundary condition map file.

The surface water variable is actually updated in the Flux_SWstuff function

Parameters:

	it	Horizontal iteration number
	SURFACE_WAT	SURFACE_WAT
	SED_ELEV	SED_ELEV
	HYD_MANNINGS_N	HYD_MANNINGS_N
	STUF1	SALT_SURF_WT
	STUF2	DINdummy
	STUF3	TP_SF_WT
	SfWat_vel_day	SF_WT_VEL_mag

Remarks:

Surface water flow - levee interaction rules:
ON_MAP value: if running managed canal network, ON_MAP is modified in WatMgmt.c, beyond just in/out (true/false) of model active domain, to these values (expressed in integer format):

• 0 Not in model active domain
• 1 Allow flow in no direction
• 2 Allow flow to east<->west
• 3 Allow flow to south<->north
• 4 Allow flow in all directions

ON_MAP values in cases of cells interacting with water control structure interactions:

• 101 water control structure interaction (allow no flow)
• 102 water control structure interaction (allow flow to east<->west)
• 103 water control structure interaction (allow flow to south<->north)
• 104 water control structure interaction (allow all flow)

• 201 water control structure interaction (allow no flow)
• 202 water control structure interaction (allow flow to east<->west)
• 203 water control structure interaction (allow flow to south<->north)
• 204 water control structure interaction (allow all flow) etc.
  

The results of basic bitwise operations shown here (for those of us who don't routinely work in this mode):

• (ON_MAP) (ON_MAP-1) (ON_MAP-1)&1 (ON_MAP-1)&2
• \000 -1 1 2
• \001 0 0 0
• \002 1 1 0
• \003 2 0 2
• \004 3 1 2

Boundary condition flow allowances (read from BoundCond.BIN input map, values all within "ON_MAP" model domain):

• Value of (integer) BCondFlow:
• 1 Allow no flows external to model domain
• 3 Allow surface- and ground- water flows to/from external boundary cells
• 4 Allow groundwater (but not surface) flows to/from external boundary cells
• (9) Archaic, will remove (allowing groundwater flows, with static external stage conditions)

Variables local to function

ix, iy Model grid cell row, column, respectively
FFlux Water flux (m) between grid cells

Definition at line 121 of file Fluxes.c.

References BCondFlow, Flux_SWcells(), Flux_SWstuff(), ON_MAP, s0, s1, and T.

Referenced by horizFlow().

{ 
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 */

Here is the call graph for this function:

float Flux_SWcells	(	int	i0,
		int	i1,
		int	j0,
		int	j1,
		float *	SWater,
		float *	Elevation,
		float *	MC
	)

Surface water flux calculations.

Application of Manning's eqn to calculate flux between two adjacent cells (i0,i1) and (j0,j1) Returns height flux. Flux is positive if flow is from i to j. Checks for available volume, and that flow cannot make the head in recepient cell higher than in the donor one.

Parameters:

	i0	Row of cell from (positive flow)
	i1	Column of cell from (positive flow)
	j0	Row of cell to (positive flow)
	j1	Column of cell to (positive flow)
	SWater	SURFACE_WAT
	Elevation	SED_ELEV
	MC	HYD_MANNINGS_N

Returns:

Flux Water flux (m)

Variables local to function

dh, adh The head difference, and absolute value of head difference (m)
MC_cells The mean manning's coefficient for a pair of cells or a cell
Hi, Hj The stage height in the i'th and j'th cells, respectively (m)
cellLoci, cellLocj The array location of the i'th and j'th cells, respectively

Definition at line 231 of file Fluxes.c.

References Abs, BCondFlow, boundcond_depth, GP_BC_InputRow, GP_BC_SfWat_add, GP_BC_SfWat_addHi, GP_BC_SfWat_targ, GP_DetentZ, GP_mannDepthPow, GP_mannHeadPow, isModExperim, Max, Min, ON_MAP, ramp, step_Cell, and T.

Referenced by flowCalc_CelCel(), and Flux_SWater().

{
  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
	)

Flux of surface water constituents.

The constituents (nutrients, salinity, etc), are passed among cells, and surface water variable is updated, along with making budget calcs. These are units of constituent mass being fluxed from i0,i1 to j0,j1 Flux & SURFACE_WAT are in units of height (m)

Parameters:

	i0	Row of cell from (positive flow)
	i1	Column of cell from (positive flow)
	j0	Row of cell to (positive flow)
	j1	Column of cell to (positive flow)
	Flux	Water flux (m) between grid cells
	SURFACE_WAT	SURFACE_WAT
	STUF1	SALT_SURF_WT
	STUF2	DINdummy
	STUF3	TP_SF_WT
	SfWat_vel_day	SF_WT_VEL_mag

Definition at line 383 of file Fluxes.c.

References basins, basn, basn_list, CELL_SIZE, debug, Disp_Calc(), dynERRORnum, basndef::family, basndef::FLok, GP_BC_Pconc, GP_BC_Sconc, GP_MinCheck, Max, Min, msgStr, basndef::numFLok, ON_MAP, P_IN_OVL, P_OUT_OVL, S_IN_OVL, S_OUT_OVL, sfstep, SimTime, T, simTime::TIME, True, Veloc_Calc(), VOL_IN_OVL, VOL_OUT_OVL, WatMgmtOn, and WriteMsg().

Referenced by flowCalc_CelCel(), and Flux_SWater().

{
    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;
        
}

Here is the call graph for this function:

double Veloc_Calc	(	float	flux,
		float	depth_i,
		float	depth_j,
		float	tim_step
	)

Calculate horizontal flow velocities of surface water.

Parameters:

	flux	Water flux (m)
	depth_i	Water depth (m) in i'th location
	depth_j	Water depth (m) in j'th location
	tim_step	Time step of calculation (d)

Returns:

sf_wt_veloc Return the horizontal velocity of advected water (m/d)

Definition at line 649 of file Fluxes.c.

References Abs, and celWid.

Referenced by Flux_SWstuff().

{
        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
	)

Calculate dispersion for constituent fluxes.

Parameters:

	flux	Water flux (m)
	depth_i	Water depth (m) in i'th location
	depth_j	Water depth (m) in j'th location
	tim_step	Time step of calculation (d)
	sf_wt_veloc	orizontal velocity of advected water (m/d)

Returns:

flux_adj Return the adjusted flux (height, m), adjusted (increased/decreased) due to constituent dispersion

Definition at line 670 of file Fluxes.c.

References celWid, GP_dispLenRef, and GP_dispParm.

Referenced by Flux_SWstuff().

{
    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
	)

Groundwater fluxing routine.

This is the alternating direction function that first fluxes water in the E/W direction and then in the N/S direction. It sweeps from left to right, then goes back from right to left, then goes from top to bottom and finally from bottom to top. This alternates every iteration. Water is fluxed with mass conservation, allowing outflow from (no inflow to) the model system.

Parameters:

	it	Horizontal iteration number
	SatWat	SAT_WATER
	Unsat	UNSAT_WATER
	SfWat	SURFACE_WAT
	rate	HYD_RCCONDUCT
	poros	HP_HYD_POROSITY
	sp_yield	HP_HYD_SPEC_YIELD
	elev	SED_ELEV
	gwSTUF1	SALT_SED_WT
	gwSTUF2	DINdummy
	gwSTUF3	TP_SED_WT
	swSTUF1	SALT_SURF_WT
	swSTUF2	DINdummy
	swSTUF3	TP_SF_WT

Remarks:

Groundwater fluxing (this function) is called every other horizontal hyd_iter (it = int(hyd_iter-1)/2 ), or half as frequently as the surface water flow calculations.

Boundary condition flow allowances (read from BoundCond.BIN input map, values all within "ON_MAP" model domain): Value of (integer) BCondFlow:

• 1 Allow no flows external to model domain
• 3 Allow surface- and ground- water flows to/from external boundary cells
• 4 Allow groundwater (but not surface) flows to/from external boundary cells
• (9) Archaic, will remove (allowing groundwater flows, with static external stage conditions)

Definition at line 721 of file Fluxes.c.

References BCondFlow, Flux_GWcells(), ON_MAP, s0, s1, and T.

Referenced by horizFlow().

{ 
  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); 
         }
     } 
 }

}

Here is the call graph for this function:

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
	)

Ground water flux eqns, incl. integration of groundwater with surface water.

Application of Darcy's eqn to calculate flux between two adjacent cells (i0,i1) and (j0,j1) Flux is positive if flow is from i to j. Checks for available volume, and that flow cannot make the head in recepient cell higher than in the donor one.

Integrates the vertical exchange among surface, unsaturated, and saturated water; updates all hydro and constituent state variables. Calcs the budget sums for water and constituents.

Parameters:

	i0	Row of cell from (positive flow)
	i1	Column of cell from (positive flow)
	j0	Row of cell to (positive flow)
	j1	Column of cell to (positive flow)
	SatWat	SAT_WATER
	Unsat	UNSAT_WATER
	SfWat	SURFACE_WAT
	rate	HYD_RCCONDUCT
	poros	HP_HYD_POROSITY
	sp_yield	HP_HYD_SPEC_YIELD
	elev	SED_ELEV
	gwSTUF1	SALT_SED_WT
	gwSTUF2	DINdummy
	gwSTUF3	TP_SED_WT
	swSTUF1	SALT_SURF_WT
	swSTUF2	DINdummy
	swSTUF3	TP_SF_WT

Variables local to function

cell_don, cell_rec, cell_i, cell_j The array location of the donor, recipient, i'th, and j'th cells, respectively
sign, equilib, revers The sign of a result, equilibrium flag, and head-reversal flag, respectively
GW_head_i,GW_head_j,tot_head_i,tot_head_j (m) the GroundWater head in the i'th and j'th cells, and the total water head in the i'th and j'th cells, respectively
dh (m) the difference in stage heights between a pair of cells
AvgRate (m/d) the average hydraulic conductivity of the two cells

H_pot_don, H_pot_rec (m) the potential new head in the donor and in the recipient cells, respectively

Flux (m/d) The potential horizontal flux between two cells
fluxTOunsat_don (m/d) donor cell, field capacity volume (height) remaining in unsaturated zone associated with a horizontal flux
fluxHoriz (m/d) the actual water volume (height) that may flux horizontally (leaving field capacity in donor cell)
Sat_vs_unsat (dimless) same form of control function (0,1) used in UnitMod.c Hydrologic Module cell_dyn7 to determine relative pathway of flow from surface storage (into saturated vs. unsaturated)

UnsatZ_don (m) donor cell, new unsat zone depth after calculated groundwater flow
UnsatCap_don (m) donor cell, maximum pore space capacity in the depth of new unsaturated zone
UnsatPot_don (m) donor cell, (height of) the volume of pore space (soil "removed") that is unoccupied in the unsat zone

UnsatZ_rec (m) recipient cell, old unsat zone depth before calculated groundwater flow
UnsatCap_rec (m) recipient cell, maximum pore space capacity in the depth of old unsaturated zone
UnsatPot_rec (m) recipient cell, (height of) the volume of pore space (soil "removed") that is unoccupied in the unsat zone

sfTOsat_don (m/d) donor cell, surface to saturated flow
unsatTOsat_don (m/d) donor cell, unsaturated to saturated flow
unsatTOsat_rec (m/d) recipient cell, unsaturated to saturated flow
satTOsf_rec (m/d) recipient cell, upflow to surface beyond soil capacity
horizTOsat_rec (m/d) recipient cell, horizontal inflow to soil into saturated storage (== fluxHoriz)

sedwat_don, sedwat_rec (m) the sediment/soil water in the donor and recipient cells, respectively

sed_stuf1fl_horz, sed_stuf2fl_horz, sed_stuf3fl_horz (kg/d) the horizontal mass flux of stuff (salt, N (unused), and P respectively)
sf_stuf1fl_don, sf_stuf2fl_don, sf_stuf3fl_don (kg/d) donor cell, the vertical surface-to-saturated mass flux of stuff (salt, N (unused), and P respectively)
sed_stuf1fl_rec, sed_stuf2fl_rec, sed_stuf3fl_rec (kg/d) recipient cell, the vertical sedwat-to-surface mass flux of stuff (salt, N (unused), and P respectively)

Note:

The variable "boundcond_depth" was intended to be stage height data. It is not stage. The boundcond_depth = water depth (positive/negative) relative to local land surface elevation output from the SFWMM v5.4 (that model's "daily_stg_minus_lsel.bin"). For ELM v2.3 (and perhaps subsequent versions?), we use this variable to estimate external boundary condition stage in grid cells adjacent to the active (ON_MAP) ELM domain cells, using the land surface elevation from the ON_MAP ELM domain cells. This method/data-use is the simplest (i.e., no pre-processing of the SFWMM output data file), and it remains to be determined whether it is maintained or upgraded in subsequent versions (using a couple of simple math re-arrangements shown below). It is flagged in comments with "v2.3tmp"

Definition at line 815 of file Fluxes.c.

References Abs, basn, basn_list, BCondFlow, boundcond_depth, CELL_SIZE, debug, dynERRORnum, Exp, basndef::family, GP_BC_InputRow, GP_BC_SatWat_add, GP_DetentZ, GP_MinCheck, gwstep, HAB, isModExperim, Max, Min, msgStr, ON_MAP, P_IN_GW, P_OUT_GW, S_IN_GW, S_OUT_GW, SimTime, T, simTime::TIME, True, VOL_IN_GW, VOL_OUT_GW, and WriteMsg().

Referenced by Flux_GWater().

{       
            
    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 */

Here is the call graph for this function: