UnitMod.c

Go to the documentation of this file.

/*###############
To change among model grid scales of different resolution/extent: 
  no code changes are needed here or elsewhere
###############

 ***
 NOTE: The init_eqns module is not fully updated for properly
 turning off individual **interacting** *biological/chemical* modules at runtime.  
 If one *biological/chemical* module is turned off, 
 they all need to be turned off. 

 The following *biological/chemical* modules must be ON or OFF as a group:
 (cell_dyn2 + cell_dyn4 + cell_dyn8 + cell_dyn9 + cell_dyn12)
 cell_dyn13, the net settling rate module, can be turned on only when above module group is off
 ***

*/

  
/* General NOTES on revisions to this source file:
        april 2000 VERSION 2.1 - output available for application for water quality performance
        july 2002 VERSIOIN 2.1a - first widely-available public release of code -
              - misc modifications, such as revising method and types of
              parameters that are read in (reduced number of parms in "habitat-specific" dbase);
              - added some more code documentation, etc.
        July 2003 VERSION 2.2.0 - moved declarations into header file, other similar changes
              to help reviewers better understand organization of codes.
              - Incorporated new rainfall module (highly revised rain_inp.c) that generalizes the
              code to better accomodate implementations of all scale; NO DATA CHANGE - passes identical data
              from SFWMM data file to the ELM.
              - Included are some specific, minor code changes (that do not alter calculations)  
        July 2004 VERSION 2.2.1 - improved code documentation
******        - AN ESSENTIAL COMPANION TO THIS CODE DOCUMENTATION IS NOW THE "ModelOutlist_creator_version#.xls"
******            workbook (OpenOffice/Excel).  In this workbook, all of the variables
******            that are spatial arrays (and other types) are defined/described (including their units).
              - In order to avoid the "clutter" of comments interspersing virtually every equation,
              this workbook was used to fully define the use of the variables in the model.
              - The code here STILL contains comments that describe the variables that may be most 
              important or that are better understood by explicitly stating underlying assumptions. 
              - Included are several specific, minor code changes (still maintaining v2.1a calcs). 
        Aug 2004 VERSION 2.3.0 - added bathymetry map input (dynamic boundary conditions in other src files -> version upgrade)
                - added bathymetry map input
                - added new spatial time series data input for dynamic boundary conditions
                - version increment to 2.3 due to changes in other source files
        Oct 2004 VERSION 2.3.1 - internal release 
                - checked that results reasonable, not fully checked
        Nov 2004 VERSION 2.3.2 - documentation upgrade
                - removed standing detritus (was operative in v2.3.0 & earlier), consumer, and fire modules completely from code
                - removed extraneous comments, inoperative variables, inoperative parameters
                - added Doxygen tags in comments
        Apr 2005 v2.4.4: changed init maps: icMult removed, icMacBio added
        May 2006 v2.5.0: added some debug-level output for boundary condition stage
        Oct 2006 v2.6.0: added more boundary condition functionality
                - reinstated the ptsl (search for "v2.6 PTSL") on-the-fly spatial time series interpolators for a non-gridIO alternative for meteorological data
                - added model-experiment parameters (ModExperimParms_NOM, boundary condition experiments right now)
        Nov 2006 v2.6.0: added input map of atmospheric phosphorus deposition rate 
                - added a P deposition method of applying temporally-constant, spatially-varying 
                P depos rate (TP_AtmosDepos).  If the global parameter of TP conc in rain is >= 0.0, then
                the old of method of applying P conc to rainfall is used instead.
        Oct 2007 v2.7.0: added functionality of calc'ing surface water flow velocities
                - new variable SF_WT_VEL_mag, of net surface water flow velocity magnitude (all directions)
                - fixed bug in auto-scaling of dispersion among model grids of diff sizes
        Nov 2007 v2.7.1: added functionality of outputting BC model water depths
                - new variable BCmodel_sfwat of the BoundaryCondition model's water depth (currently the SFWMM), so
                that those depths may be output (for comparison to ELM), if desired
        Feb 2008 v2.8.1: added functionality of salt (chloride) inputs to domain from atmosphere
                - new variable SALT_Atm_Depos, calculated from rainfall and new global parm GP_CL_IN_RAIN, 
                or from new input map SALT_AtmosDepos of constant total atmospheric deposition
                - added this input to the salt budget (BudgStats.c)
        June 2008 v2.8.2: added new variables for output
                - new variable HydRelDepPosNeg (stage minus land surface elevation).  
                - new associated variable HydRelDepPosNegAvg added to BudgStats.c (average of HydRelDepPosNeg)
                - These are for reporting purposes only; data are pos/neg, and thus map output MUST be 4 byte floats
        Mar 2011 v2.8.4: added new parameter
                - new parameter GP_WQualMonitZ (threshold depth used to report 0.000 for salt and TP conc in surface water when below threshold depth,
                used to use the hydro detention depth).  
        July 2011 v2.8.4: altered soil initializations
                - new DEPOS_ORG_MAT initialization routines
                - new initial map for TP_SEDWT_CONCACT (instead of hab-specific parameter)
                - now (optionally) using hab-specific parameter of activezone depth of 10 cm (instead of 30 cm)
        Sept 2011 v2.8.5: added succession parm inputs for salinity, and algorithm for mac growth response 
                - changed macrophyte module to alter growth in response to salinity re. existing HP_MAC_SALIN_THRESH
                - see Success.c for succession module changes

*/

#include "unitmod.h"
#include "unitmod_vars.h"
#include "unitmod_habparms.h"
#include "unitmod_globparms.h"
#include "budgstats_statvars.h"  

int call_cell_dyn(int sector,int step)
 {
  int rv=0;

  switch(sector) {

    case 99: { stats(step); rv=1; } break;
    case 0: { horizFlow(step); rv=1; } break;
    case 1: { cell_dyn1(step); rv=1; } break;
    case 2: { cell_dyn2(step); rv=1; } break;
    case 4: { cell_dyn4(step); rv=1; } break;
    case 7: { cell_dyn7(step); rv=1; } break;
    case 8: { cell_dyn8(step); rv=1; } break;
    case 9: { cell_dyn9(step); rv=1; } break;
    case 10: { cell_dyn10(step); rv=1; } break;
    case 12: { cell_dyn12(step); rv=1; } break;
    case 13: { cell_dyn13(step); rv=1; } break;
            default:  printf("Warning, undefined sector number:%d\n",sector);
  }
  return rv;
}


/*###############
The following modules (cell_dyn*) are the dynamic calcs for the model.
They are called in the order determined by the call_cell_dyn function.
################*/


/*******/
/* Parameter definitions:
      global parameters in GlobalParms.xls OpenOffice/Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
******/

void cell_dyn1(int step)
{
   int fail = -1;
   int ix, iy, cellLoc, stat=1; 
   float fTemp;
   extern IsSubstituteModel; /* from generic_driver.h */
   
   /* v2.5 may be temporary? */
   FILE *boundaryFile; 
   char filename[40];

     /* SimTime.TIME is a (float) julian day counter since initialization (calc'd in main function, Generic_Driver.c source).  
         (SimTime.TIME includes fractions of days if the vertical DT<1.0, but it is
         unlikely that the vertical DT will deviate from 1 day). */

         if ( debug > 4 ) {
                if (SimTime.TIME>60)  printf("\nDebug level is very high - you're probably getting too much output data in %s/%s/Output/Debug/! Decrease your simulation length.\n", OutputPath, ProjName); 
                exit(-1); 
         }

/* v2.6 PTSL */
/* a choice of rainfall and pET input data routines */
    if (!isPTSL) { /* do the gridIO inputs */
        if (FMOD(SimTime.TIME,1.0) < 0.001) { /* daily re-mapping of coarse-grid input data */
                /* remap sfwmm (or other grid) rain data to ELM grid scale */
          stat=rain_data_wmm(boundcond_rain); 
          if(stat==fail)  
         { 
                  sprintf(msgStr,"Problem with rainfall data, time = %f.\n",SimTime.TIME); 
                         WriteMsg(msgStr,1); 
                         exit(fail);    
         }

     /* remap sfwmm (or other grid) potential ET data to ELM grid scale */
                  stat=evap_data_wmm(boundcond_ETp); 
                  if(stat==-1)
                 {
                   sprintf(msgStr,"Problem with ETp data, time = %f.\n",SimTime.TIME);
                   WriteMsg(msgStr,1);
                   exit(fail);
                          } 

     /* remap sfwmm (or other grid) stage/water depth data to ELM grid scale */
      stat=stage_data_wmm(boundcond_depth);  
     if(stat==-1)
     {
       sprintf(msgStr,"Problem with stage data, time = %f.\n",SimTime.TIME);
       WriteMsg(msgStr,1);
       exit(fail);
     } 
     

         for ( ix = 1; ix <= s0 ; ix++ ) 
            for ( iy = 1 ; iy <= s1 ; iy++ )
                {
                        cellLoc= T(ix,iy);

    /* v2.7.1 assign the external boundary condition model's depths to an ELM variable within the rectangular 
    ELM domain (i.e., on-map and off-map), for output only (this variable not used in any calcs) */
                        if (IsSubstituteModel == 0) BCmodel_sfwat[cellLoc] = (boundcond_depth[cellLoc] > 0.0) ? (boundcond_depth[cellLoc]) : (0.0) ;
                        
                        else { /* IsSubstituteModel is true */
                                if(ON_MAP[cellLoc])  {
       /* v2.8 - added feature to replace all ELM hydrologic calcs, substituting the Boundary Condition Model's
                        "stage_minus_landElevation" variable into ELM hydro variables.
                        This was done in order to post-process the SFWMM or NSM using the same procedures as ELM 
                        At initial setup, it was determined whether the user is requesting to only run cell_dyn1 (and stats, case #99) in UnitMod.c,
                        which would make the IsSubstituteModel flag True */

                                SURFACE_WAT[cellLoc] = (boundcond_depth[cellLoc] > 0.0) ? (boundcond_depth[cellLoc]) : (0.0) ;
                                SAT_WT_HEAD[cellLoc] = (boundcond_depth[cellLoc] <= 0.0) ? (SED_ELEV[cellLoc] + boundcond_depth[cellLoc]) : (SED_ELEV[cellLoc]) ;
                                HydTotHd[cellLoc]  = SAT_WT_HEAD[cellLoc]+SURFACE_WAT[cellLoc];
                                SAT_WATER[cellLoc] = SAT_WT_HEAD[cellLoc] * HP_HYD_POROSITY[HAB[cellLoc]];  /* this isn't needed for current purposes, but may be of interest later */
                                UNSAT_DEPTH[cellLoc] = Max(SED_ELEV[cellLoc]-SAT_WT_HEAD[cellLoc],0.0);
                                UNSAT_WATER[cellLoc] = UNSAT_DEPTH[cellLoc] * HP_HYD_POROSITY[HAB[cellLoc]] * 0.0; /* assume there is no moisture in pore spaces, thus never any unsat storage */
              /* v2.8.2 new variable, for reporting purposes only: positive/negative water depth relative to land elevation (stage minus land elevation) */
                        HydRelDepPosNeg[cellLoc] =  HydTotHd[cellLoc] - SED_ELEV[cellLoc];
                        }
                        }
                                
                        }   



    } /* end of mapping multi-grid, gridIO data */
    } /* end of !isPTSL */
    
    else { /* instead of gridIO spatial data, calculate on-the-fly spatial interpolations of rain and pET point time series data  */
                /* NOTE that this option does not use gridIO stage/depth boundary condition data, and thus should be used in conjunction
                with the ModExperimParms_NOM synthetic boundary condition data */
        fTemp = FMOD(SimTime.TIME,Timestep2*NPtTS2);
                PTSL_CreatePointMap(pTSList2,boundcond_rain,'f',(int)(fTemp/Timestep2),1.0);
        fTemp = FMOD(SimTime.TIME,Timestep3*NPtTS3);
                PTSL_CreatePointMap(pTSList3,boundcond_ETp,'f',(int)(fTemp/Timestep3),1.0);
    }


            
    /* the "julian" day counter within a 365 day "year"  */
    DAYJUL = ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 365.0);
    /* DAYLENGTH is not used */
    /* DAYLENGTH = AMPL*Sin((DAYJUL-79.0)*0.01721)+12.0; */ /* length of daylight (hours) */

                
/*  Nikolov & Zeller (1992) solar radiation algorithm 
   the algorithm and all parameters are published in the (Ecol. Mod., 61:149-168) manuscript,
   and the only modifiable parameters (in ELM database system) are local latitude and altitude */
    SOLDEC1 = 0.39785*Sin(4.868961+0.017203*DAYJUL   
                          +0.033446*Sin(6.224111+0.017202*DAYJUL));
    SOLCOSDEC = sqrt(1.0-SOLDEC1*SOLDEC1);
    SOLELEV_SINE = Sin(LATRAD)*SOLDEC1+Cos(LATRAD)*SOLCOSDEC;
    SOLALTCORR = (1.0-Exp(-0.014*(GP_ALTIT-274.0)/(SOLELEV_SINE*274.0)));
    SOLDEC = Arctan(SOLDEC1/sqrt(1.0-SOLDEC1*SOLDEC1));
    SOLRISSET_HA1 = -Tan(LATRAD)*Tan(SOLDEC);
    SOLRISSET_HA = ( (SOLRISSET_HA1==0.0) ) ? ( PI*0.5 ) : (   ( (SOLRISSET_HA1<0.0) ) ? 
                                                               ( PI+Arctan(sqrt(1.0-SOLRISSET_HA1*SOLRISSET_HA1)/SOLRISSET_HA1)  ) : 
                                                               (    Arctan(sqrt(1.0-SOLRISSET_HA1*SOLRISSET_HA1)/SOLRISSET_HA1)));
    SOLRADATMOS = 458.37*2.0*(1.0+0.033*Cos(360.0/365.0*PI/180.0*DAYJUL))
        * ( Cos(LATRAD)*Cos(SOLDEC)*Sin(SOLRISSET_HA) 
            + SOLRISSET_HA*180.0/(57.296*PI)*Sin(LATRAD)*Sin(SOLDEC));
            
  
        /* daily habitat switching */
    if ( (SimTime.TIME - (int)SimTime.TIME) < DT/2 && HabSwitchOn ) /* HabSwitchOn flag to invoke switching */
        for(ix=1; ix<=s0; ix++)  
            for(iy=1; iy<=s1; iy++)  
                if(ON_MAP[cellLoc= T(ix,iy)]) { /* TPtoSoil == kg/kg */
                    HAB [cellLoc] = HabSwitch (ix, iy, SURFACE_WAT, TPtoSOIL, (int*)FIREdummy, SAL_SF_WT, HAB); 
                } 
}



/*******/
/*
  NC_ALG[cellLoc] == carbon mass of the NonCalcareous (more appropriately, the eutrophic, or non-native) periphyton community (gC/m^2)
  C_ALG[cellLoc] == carbon mass of the Calcareous (more appropriately, the oligotrophic, or native) periphyton community (gC/m^2)

  NC_ALG_P[cellLoc] == phosphorus mass of NonCalcareous (more appropriately, the eutrophic, or non-native) periphyton community (gP/m^2)
  C_ALG_P[cellLoc] == phosphorus mass of Calcareous (more appropriately, the oligotrophic, or native) periphyton community (gP/m^2)
  
Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn2(int step)
 {
 int ix, iy, cellLoc; 
 float reduc, min_litTemp,I_ISat, Z_extinct; 
 double PO4Pconc, PO4P;
 float C_ALG_thresh_CF, mortPot;
 
  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) {

    if(ON_MAP[cellLoc= T(ix,iy)])  {

/* these thresholds need updating when a habitat type of a grid cell changes */
     ALG_REFUGE[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*GP_ALG_REF_MULT;
     ALG_SAT[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*0.9;

     NC_ALG_AVAIL_MORT[cellLoc]  = Max(NC_ALG[cellLoc]-ALG_REFUGE[cellLoc],0);
     C_ALG_AVAIL_MORT[cellLoc]  = Max(C_ALG[cellLoc]-ALG_REFUGE[cellLoc],0);
/* bio-avail P (PO4) is calc'd from TP, using pre-processed regression for predicting PO4 from TP */
/* assume that periphyton (microbial) alkaline phosphotase activity keeps PO4 at least 10% of TP conc */
     PO4Pconc = Max(TP_SFWT_CONC_MG[cellLoc]*GP_PO4toTP + GP_PO4toTPint,
                                0.10 * TP_SFWT_CONC_MG[cellLoc]); /* mg/L */

/* light, water, temperature controls apply to both calc and non-calc */
     ALG_LIGHT_EXTINCT[cellLoc]  = GP_alg_light_ext_coef; 
                     /* algal self-shading implicit in density-dependent constraint function later */
     ALG_INCID_LIGHT[cellLoc]  = SOLRADGRD[cellLoc]*Exp(-MAC_LAI[cellLoc]*GP_ALG_SHADE_FACTOR);
     Z_extinct = SURFACE_WAT[cellLoc]*ALG_LIGHT_EXTINCT[cellLoc];
     I_ISat = ALG_INCID_LIGHT[cellLoc]/GP_ALG_LIGHT_SAT;
                     /* averaged over whole water column (based on Steele '65) */
     ALG_LIGHT_CF[cellLoc]  = ( Z_extinct > 0.0 ) ? 
                     ( 2.718/Z_extinct * (Exp(-I_ISat * Exp(-Z_extinct)) - Exp(-I_ISat)) ) :
                     (I_ISat*Exp(1.0-I_ISat));
                     /* low-water growth constraint ready for something better based on data */
     ALG_WAT_CF[cellLoc]  = ( SURFACE_WAT[cellLoc]>0.0 ) ? ( 1.0 ) : ( 0.0);
     ALG_TEMP_CF[cellLoc]  = tempCF(0, 0.20, GP_ALG_TEMP_OPT, 5.0, 40.0, H2O_TEMP[cellLoc]);

      min_litTemp = Min(ALG_LIGHT_CF[cellLoc],ALG_TEMP_CF[cellLoc]);

/* the 2 communities have same form of growth response to avail phosphorus */
     NC_ALG_NUT_CF[cellLoc]  =
                     Exp(-GP_alg_uptake_coef * Max(GP_NC_ALG_KS_P-PO4Pconc, 0.0)/GP_NC_ALG_KS_P) ; /* mg/L */
     C_ALG_NUT_CF[cellLoc]  =
                     Exp(-GP_alg_uptake_coef * Max(GP_C_ALG_KS_P-PO4Pconc, 0.0)/GP_C_ALG_KS_P) ; /* mg/L */

         /* the form of the control function assumes that at very low
             P conc, the alkaline phosphotase activity of the microbial assemblage scavenges P, maintaining a minimum nutrient availability to community */
     NC_ALG_PROD_CF[cellLoc]  = Min(min_litTemp,ALG_WAT_CF[cellLoc])*Max(NC_ALG_NUT_CF[cellLoc], GP_alg_alkP_min);
     C_ALG_PROD_CF[cellLoc]   = Min(min_litTemp,ALG_WAT_CF[cellLoc])*Max(C_ALG_NUT_CF[cellLoc], GP_alg_alkP_min);

     NC_ALG_RESP_POT[cellLoc]  = 
         ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
         ( 0.0) :
         ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc]*NC_ALG_AVAIL_MORT[cellLoc] ); 
     C_ALG_RESP_POT[cellLoc]  = 
         ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
         ( 0.0) :
         ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc] *C_ALG_AVAIL_MORT[cellLoc] ); 

     NC_ALG_RESP[cellLoc]  =  
         ( NC_ALG_RESP_POT[cellLoc]*DT>NC_ALG_AVAIL_MORT[cellLoc] ) ? 
         ( NC_ALG_AVAIL_MORT[cellLoc]/DT ) : 
         ( NC_ALG_RESP_POT[cellLoc]);
     C_ALG_RESP[cellLoc]  =  
         ( C_ALG_RESP_POT[cellLoc]*DT>C_ALG_AVAIL_MORT[cellLoc] ) ? 
         ( C_ALG_AVAIL_MORT[cellLoc]/DT ) : 
         ( C_ALG_RESP_POT[cellLoc]);
                 
          /* this is the threshold control function that increases
              calcareous/native periph mortality (likely due to loss of calcareous sheath) as P conc. increases */
     C_ALG_thresh_CF = Min(exp(GP_alg_R_accel*Max( TP_SFWT_CONC_MG[cellLoc]-GP_C_ALG_threshTP,0.0)/GP_C_ALG_threshTP), 100.0);

     NC_ALG_MORT_POT[cellLoc]  = 
         ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
         ( NC_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT_DRY ) : 
         ( NC_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT);
     C_ALG_MORT_POT[cellLoc]  = 
         ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
         ( C_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT_DRY ) : 
         ( C_ALG_thresh_CF * C_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT);

     NC_ALG_MORT[cellLoc]  = 
         ( (NC_ALG_MORT_POT[cellLoc]+NC_ALG_RESP[cellLoc])*DT>NC_ALG_AVAIL_MORT[cellLoc] ) ? 
         ( (NC_ALG_AVAIL_MORT[cellLoc]-NC_ALG_RESP[cellLoc]*DT)/DT ) : 
         ( NC_ALG_MORT_POT[cellLoc]);
     C_ALG_MORT[cellLoc]  = 
         ( (C_ALG_MORT_POT[cellLoc]+C_ALG_RESP[cellLoc])*DT>C_ALG_AVAIL_MORT[cellLoc] ) ? 
         ( (C_ALG_AVAIL_MORT[cellLoc]-C_ALG_RESP[cellLoc]*DT)/DT ) : 
         ( C_ALG_MORT_POT[cellLoc]);
                 
/* gross production of the 2 communities (gC/m2/d, NOT kgC/m2/d) */
/* density constraint on both noncalc and calc, competition effect accentuated by calc algae */

     NC_ALG_GPP[cellLoc]  =  NC_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*NC_ALG[cellLoc]       
             *Max( (1.0-(GP_AlgComp*C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
     C_ALG_GPP[cellLoc]  =  C_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*C_ALG[cellLoc] 
       *Max( (1.0-(        C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
       
                     
/* P uptake is dependent on available P and is relative to a maximum P:C ratio for the tissue
   (g C/m^2/d * P:Cmax * dimless * dimless = gP/m2/d (NOT kg) )*/
     NC_ALG_GPP_P[cellLoc] = NC_ALG_GPP[cellLoc] *GP_ALG_PC * NC_ALG_NUT_CF[cellLoc]
                     * Max(1.0-NC_ALG_PC[cellLoc]/GP_ALG_PC, 0.0); 
     C_ALG_GPP_P[cellLoc] = C_ALG_GPP[cellLoc] * GP_ALG_PC * C_ALG_NUT_CF[cellLoc]
                     * Max(1.0-C_ALG_PC[cellLoc]/GP_ALG_PC, 0.0); 
                 
/* check for available P mass (the nutCF does not) */
     PO4P = ramp(Min(PO4Pconc * SFWT_VOL[cellLoc],  conv_kgTOg*TP_SF_WT[cellLoc]) ); /*g P available; v2.5 put in ramp (constrain non-neg) */
     reduc = ( (NC_ALG_GPP_P[cellLoc]+C_ALG_GPP_P[cellLoc]) > 0) ? 
                     (PO4P / ( (NC_ALG_GPP_P[cellLoc]+C_ALG_GPP_P[cellLoc])*CELL_SIZE*DT) ) :
                     (1.0);
    /* can have high conc, but low mass of P avail, in presence of high peri biomass and high demand */
    /* reduce the production proportionally if excess demand is found */
    if (reduc < 1.0) {
                     NC_ALG_GPP[cellLoc] *= reduc;   
                     NC_ALG_GPP_P[cellLoc] *= reduc;   
                     C_ALG_GPP[cellLoc] *= reduc; 
                     C_ALG_GPP_P[cellLoc] *= reduc; 
                 }

/* state variables calc'd (gC/m2, NOT kgC/m2) */
     NC_ALG[cellLoc] =  NC_ALG[cellLoc] 
         + (NC_ALG_GPP[cellLoc]
          - NC_ALG_RESP[cellLoc] - NC_ALG_MORT[cellLoc]) * DT;

     C_ALG[cellLoc] =  C_ALG[cellLoc] 
         + (C_ALG_GPP[cellLoc] 
         - C_ALG_RESP[cellLoc]  - C_ALG_MORT[cellLoc]) * DT;

/* carbon NPP not currently used elsewhere, only for output */
     NC_ALG_NPP[cellLoc]  = NC_ALG_GPP[cellLoc]-NC_ALG_RESP[cellLoc]; 
     C_ALG_NPP[cellLoc]  = C_ALG_GPP[cellLoc]-C_ALG_RESP[cellLoc];                 
/* carbon total biomass of both communities, not currently used elsewhere, only for output */
     ALG_TOT[cellLoc] = NC_ALG[cellLoc] + C_ALG[cellLoc];
                 
                 
/*  now calc the P  associated with the C fluxes (GPP_P already calc'd) */
     mortPot = (double) NC_ALG_MORT[cellLoc] * NC_ALG_PC[cellLoc];
     NC_ALG_MORT_P[cellLoc] = (mortPot*DT>NC_ALG_P[cellLoc]) ?
                     (NC_ALG_P[cellLoc]/DT) :
                     (mortPot);
     mortPot = (double) C_ALG_MORT[cellLoc] * C_ALG_PC[cellLoc];
     C_ALG_MORT_P[cellLoc] = (mortPot*DT>C_ALG_P[cellLoc]) ?
                     (C_ALG_P[cellLoc]/DT) :
                     (mortPot);
                 

/* state variables calc'd (gP/m2, NOT kgP/m2) */
     NC_ALG_P[cellLoc] =  NC_ALG_P[cellLoc] 
         + (NC_ALG_GPP_P[cellLoc] - NC_ALG_MORT_P[cellLoc]) * DT;

     C_ALG_P[cellLoc] =  C_ALG_P[cellLoc] 
         + (C_ALG_GPP_P[cellLoc] - C_ALG_MORT_P[cellLoc]) * DT;

     NC_ALG_PC[cellLoc] = (NC_ALG[cellLoc]>0.0) ?
                     (NC_ALG_P[cellLoc]/ NC_ALG[cellLoc]) :
                     ( GP_ALG_PC * 0.03); /* default to 3% of  max P:C */
     C_ALG_PC[cellLoc] = (C_ALG[cellLoc]>0.0) ?
                     (C_ALG_P[cellLoc]/ C_ALG[cellLoc]) :
                     ( GP_ALG_PC * 0.03 ); /* default to 3% of max P:C */
     NC_ALG_PCrep[cellLoc] = (float)NC_ALG_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */
     C_ALG_PCrep[cellLoc] = (float)C_ALG_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */

     if (debug > 0 && NC_ALG[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NC_ALG C biomass (%f g/m2) in cell (%d,%d)!", 
                              SimTime.TIME, NC_ALG[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
                 if (debug > 0 && C_ALG[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg C_ALG C biomass (%f g/m2) in cell (%d,%d)!", 
                              SimTime.TIME, C_ALG[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
                 if (debug > 0 && NC_ALG_P[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NC_ALG_P P biomass (%f g/m2) in cell (%d,%d)!", 
                              SimTime.TIME, NC_ALG_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
                 if (debug > 0 && C_ALG_P[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg C_ALG_P P biomass (%f g/m2) in cell (%d,%d)!", 
                              SimTime.TIME, C_ALG_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}

                     
     TP_SFWT_UPTAK[cellLoc]  = (NC_ALG_GPP_P[cellLoc]+C_ALG_GPP_P[cellLoc])
                     *0.001*CELL_SIZE; /* gP/m2 => kg P */
/* recalc P in surface water state variable (kg P) */
     TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] - TP_SFWT_UPTAK[cellLoc] * DT;
     TP_SFWT_CONC[cellLoc]  = 
         ( SFWT_VOL[cellLoc] > 0.0 ) ? 
         ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
         ( 0.0); /* used in P fluxes for mass balance */
     TP_SFWT_CONC_MG[cellLoc]  = 
         ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? 
         (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
         (0.0); /* (g/m3==mg/L) used for reporting and other modules to evaluate P conc when water is present */

          }
  }
  }
}




/*******/
/*     DEPOS_ORG_MAT[cellLoc] == mass of Deposited Organic Matter (DOM/AOM) of soil zone, w/o floc layer (kg OM/m^2)

     DOP[cellLoc] ==  mass of Deposited Organic Phosphorus of soil zone, w/o floc layer, w/o P sorbed (kg P/m^2)
  
Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn4(int step)
 {
int ix, iy, cellLoc;
 float TPsoil, TP_sedMin, TP_sfMin;
  
 

  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) {

    if(ON_MAP[cellLoc= T(ix,iy)])  {


/* inputs of organic matter (kg OM/m2)*/
     DOM_fr_nphBio[cellLoc] = nphbio_mort_OM[cellLoc];
     DOM_FR_FLOC[cellLoc]  =  FLOC_DEPO[cellLoc] ; 

/* losses of organic matter (kg OM/m2) */

     DOM_QUALITY_CF[cellLoc]  =
          Min(Exp(-GP_DOM_decomp_coef * Max(GP_DOM_DECOMP_POPT-TP_SEDWT_CONCACTMG[cellLoc], 0.0)
          /GP_DOM_DECOMP_POPT),1.0) ; /* mg/L */
          DOM_TEMP_CF[cellLoc] = tempCF(0, 0.20, GP_DOM_DECOMP_TOPT, 5.0, 40.0, H2O_TEMP[cellLoc]);
     DOM_SED_AEROB_Z[cellLoc]  = Min(Max(UNSAT_DEPTH[cellLoc],HP_DOM_AEROBTHIN[HAB[cellLoc]]),HP_DOM_MAXDEPTH[HAB[cellLoc]]);
     DOM_SED_ANAEROB_Z[cellLoc]  = HP_DOM_MAXDEPTH[HAB[cellLoc]]-DOM_SED_AEROB_Z[cellLoc];

                     /* GP_calibDecomp is an adjustable calib parm  */
     DOM_DECOMP_POT[cellLoc]  = (double) GP_calibDecomp*GP_DOM_RCDECOMP*DOM_QUALITY_CF[cellLoc]*DOM_TEMP_CF[cellLoc]*DEPOS_ORG_MAT[cellLoc]
         *(Min(DOM_SED_AEROB_Z[cellLoc]/HP_DOM_MAXDEPTH[HAB[cellLoc]],1.0)*soil_MOIST_CF[cellLoc]
         +GP_DOM_DECOMPRED*Min(DOM_SED_ANAEROB_Z[cellLoc]/HP_DOM_MAXDEPTH[HAB[cellLoc]],1.0) );
     DOM_DECOMP[cellLoc]  =  
         ( (DOM_DECOMP_POT[cellLoc])*DT>DEPOS_ORG_MAT[cellLoc] ) ? 
         ( (DEPOS_ORG_MAT[cellLoc])/DT ) : 
         ( DOM_DECOMP_POT[cellLoc]);
/* calc state var (kg OM/m2) */
     DEPOS_ORG_MAT[cellLoc] =  DEPOS_ORG_MAT[cellLoc] + 
         ( DOM_fr_nphBio[cellLoc] + DOM_FR_FLOC[cellLoc] 
         - DOM_DECOMP[cellLoc] ) * DT;

/* soil elevation */
     DOM_Z[cellLoc] = (double) DEPOS_ORG_MAT[cellLoc] / DOM_BD[cellLoc] ; /* (m) organic depth  */
     SED_ELEV[cellLoc]  = DOM_Z[cellLoc]+Inorg_Z[cellLoc]+SED_INACT_Z[cellLoc];  /* total land surface elevation, including model GP_DATUM_DISTANCE below zero of land elev datum (NGVD/NAVD) (m) */

/* P DOM stoich (kg P /m2) */
     DOP_nphBio[cellLoc] = nphbio_mort_P[cellLoc];    
     DOP_FLOC[cellLoc] = FlocP_DEPO[cellLoc]; 

     DOP_DECOMP[cellLoc] = (double) DOM_DECOMP[cellLoc] * DOM_P_OM[cellLoc]; 

/* calc state var of total "organic" P in soil (NOT including dissolved in pore water or sorbed) (kgP/m2) */
     DOP[cellLoc] =  DOP[cellLoc]  
         + ( DOP_nphBio[cellLoc]  + DOP_FLOC[cellLoc] 
         - DOP_DECOMP[cellLoc]) * DT; /* kgP/m2 */

/* now the P ratio */
     DOM_P_OM[cellLoc] = (DEPOS_ORG_MAT[cellLoc]>0.0) ? ( DOP[cellLoc] / DEPOS_ORG_MAT[cellLoc]) : (0.0); /* kgP/kgOM */
     TPsoil = DOP[cellLoc]*CELL_SIZE + TP_SORB[cellLoc]; /* kg TP in soil */
     TPtoSOIL[cellLoc] = ((DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc])>0.0) ?
         (  TPsoil / (DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc]) ) : (0.0); /* kgP/kgsoil */
     TPtoVOL[cellLoc] =  (CELL_SIZE * DOM_Z[cellLoc]>0.0) ?
         (TPsoil / (CELL_SIZE * DOM_Z[cellLoc]) ) :
         (0.0); /* kgP/m3 soil */
         
     TPtoSOIL_rep[cellLoc] = TPtoSOIL[cellLoc] * conv_kgTOmg; /* reporting purposes only (kg/kg->mg/kg) */
     TPtoVOL_rep[cellLoc] = TPtoVOL[cellLoc] * conv_kgTOg; /* reporting purposes only (kg/m3->g/m3 == ug/cm3) */

/* now the P gain in sed water with decomp
 a small proportion goes into surface water P (below) */
     TP_sedMin  =  (1.0 - HP_DOM_AEROBTHIN[HAB[cellLoc]]  / HP_DOM_MAXDEPTH[HAB[cellLoc]] )
          * DOP_DECOMP[cellLoc] * CELL_SIZE;
                 
/* calc P in sed water state variables (kg P) */
     TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] + TP_sedMin * DT;
             /* this is the active zone, where uptake, sorption, and mineralization take place */
     TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] + TP_sedMin * DT;

     TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
          (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
          (0.0);
     TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
          ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
          (TP_SED_CONC[cellLoc]); /* g/L */
     TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m3==mg/L */
              
/* now store the ratio of the conc in the active zone relative to total, prior to horiz fluxes
***** very simple constant, code in transition **** */              
     TP_Act_to_Tot[cellLoc] = 1.0 / HP_TP_CONC_GRAD[HAB[cellLoc]];
     TP_Act_to_TotRep[cellLoc] =  (float) TP_Act_to_Tot[cellLoc];
                 
/* now the P gain in surface water with decomp in the very thin upper layer of the soil */
                     /* if there is no surface water present, assume that this
                        relative contribution will be an additional sorbed component that
                        is introduced to surface water column immediately upon hydration
                        with surface water */
     TP_sfMin  =  HP_DOM_AEROBTHIN[HAB[cellLoc]] / HP_DOM_MAXDEPTH[HAB[cellLoc]]
           * DOP_DECOMP[cellLoc] * CELL_SIZE;
                 
/* calc P in surface water state variable (kg P) */
     TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] + TP_sfMin * DT;
     TP_SFWT_CONC[cellLoc]  = 
           ( SFWT_VOL[cellLoc] > 0.0 ) ? 
           ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
           ( 0.0); /* used in P fluxes for mass balance */
     TP_SFWT_CONC_MG[cellLoc]  = 
           ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? 
           (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
           (0.0); /* (g/m3==mg/L) used for reporting and other modules to evaluate P conc when water is present */

/* for reporting only: calc sum of all P storages in grid cells (budget calcs do same for Basin/Indicator-Regions) (g P /m^2) */
     P_SUM_CELL[cellLoc] = ( (C_ALG_P[cellLoc] + NC_ALG_P[cellLoc]) * 0.001 * CELL_SIZE + /* gP/m2 => kgP */
           (mac_nph_P[cellLoc] + mac_ph_P[cellLoc] )* CELL_SIZE + /* kgP/m2 => kgP */
           TP_SORB[cellLoc] +
           ( FlocP[cellLoc] + DOP[cellLoc]  ) * CELL_SIZE + /* kgP/m2 => kgP */
           TP_SED_WT[cellLoc] + TP_SF_WT[cellLoc] ) /* kgP */
           /CELL_SIZE * conv_kgTOg; /* kg P/m^2 => g P/m^2 */
                 
    }
  }
  }
}


/*******/
void horizFlow(int step)
{
  int it;
  int ix, iy; 
 
  /* v2.7.0 adding velocity output capabilities */
  for(ix=1; ix<=s0; ix++) {
    for(iy=1; iy<=s1; iy++) {
      if(ON_MAP[T(ix,iy)])  SF_WT_VEL_mag[T(ix,iy)]=0.0; 
    }
  }

  for ( it = 0; it < hyd_iter; it++ )
  {
    sprintf(msgStr,"\b\b\b%3d",it); 
    usrErr0(msgStr);
       
    /* horizontal raster-vector canal fluxes and water management functions
       Water Management switch set at runtime in Driver.parm
       this routine also integrates surface, unsaturated, and saturated exchanges
     */
    /* Nitrogen (DINdummy) argument is dummy argument placeholder */
        if (WatMgmtOn  ) {
      Run_Canal_Network(SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N,SAT_WATER,HP_HYD_POROSITY,
                        HYD_RCCONDUCT, DINdummy,TP_SF_WT,SALT_SURF_WT,DINdummy,TP_SED_WT,SALT_SED_WT,
                        UNSAT_WATER,HP_HYD_SPEC_YIELD ); 
    }
    

    /* Function for horiz fluxing of surface water, no exchange with sat/unsat water */
    /* if the order of the solute is changed, be sure to change the mass bal info in FluxSWstuff fcn */
    /* Nitrogen (DINdummy) argument is dummy argument placeholder */
    Flux_SWater(it,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N,SALT_SURF_WT,DINdummy,TP_SF_WT,SF_WT_VEL_mag);       
  
    /* Function for horiz fluxing of ground water and its vertical itegration with unsat and surface water 
       It is only called every other hyd_iter iteration, and passes the realized number of gwat iterations to the function.  
       If the order of the solutes is changed, be sure to change the mass bal info in FluxGWstuff fcn 
     */
    /* Nitrogen (DINdummy) argument is dummy argument placeholder */
    if ( it%2 ) 
      Flux_GWater((it-1)/2, SAT_WATER, UNSAT_WATER, SURFACE_WAT, HYD_RCCONDUCT, HP_HYD_POROSITY,
                   HP_HYD_SPEC_YIELD, SED_ELEV, SALT_SED_WT, DINdummy, TP_SED_WT, SALT_SURF_WT, DINdummy, TP_SF_WT);  

  }  /* end of hyd_iter */  
  
   /* v2.7.0 adding velocity output capabilities */
  for(ix=1; ix<=s0; ix++) {
    for(iy=1; iy<=s1; iy++) {
       SF_WT_VEL_mag[T(ix,iy)]= SF_WT_VEL_mag[T(ix,iy)] / hyd_iter / 2.0; /* v2.8.4 - no changes, but explored the /hyd_iter/2.0 etc */
     }
  }

}



/*******/
/*
  SURFACE_WAT[cellLoc] == water storage above the land surface elevation (meters)
  UNSAT_WATER[cellLoc] == water storage in the pore spaces of the unsaturated zone (meters)
  SAT_WATER[cellLoc] == water storage in the pore spaces of the saturated zone (meters)
  

Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
******/
void cell_dyn7(int step)
 {
    int ix, iy, cellLoc; 
    float SatWat_Root_CF, field_cap; 
    float mann_height, N_density, f_LAI_eff, sfwat_pr1;
    float cloudy=0.0;
    
/* the horizontal raster and vector fluxes are always called before this cell_dyn */
    for(ix=1; ix<=s0; ix++) {
        for(iy=1; iy<=s1; iy++) {
            if(ON_MAP[cellLoc= T(ix,iy)])  {

              if (debug > 3) { /* these are old, relatively coarse checks - surface-ground water integration occurs in Fluxes.c */
                  if (SAT_WT_HEAD[cellLoc] - 0.01 > SED_ELEV[cellLoc] ) {
                      sprintf(msgStr,"Day %6.1f. Warning - SAT_WT_HEAD exceeds elev in cell (%d,%d) by %f.", 
                              SimTime.TIME, ix,iy,(SAT_WT_HEAD[cellLoc] - SED_ELEV[cellLoc]) ); 
                      WriteMsg( msgStr,True ); }
                  if (SURFACE_WAT[cellLoc] > 0.2 && UNSAT_DEPTH[cellLoc] > 0.1 ) {
                      sprintf(msgStr,"Day: %6.1f.  Warning - large sfwat depth (%5.2f m) in presence of unsat= %5.2f m, %4.2f %% moist, in cell (%d,%d).", 
                              SimTime.TIME, SURFACE_WAT[cellLoc], UNSAT_DEPTH[cellLoc],UNSAT_MOIST_PRP[cellLoc], ix,iy ); 
                      WriteMsg( msgStr,True );  }
                  if (SAT_WATER[cellLoc] < -0.01) { /* this seems unnecessary but... */
                      sprintf(msgStr,"Day %6.1f: capacityERR - neg SAT_WATER (%f m) in cell (%d,%d) before cell_dyn7!", 
                              SimTime.TIME, SAT_WATER[cellLoc], ix,iy ); 
                      WriteMsg( msgStr,True );  dynERRORnum++; }
                  if (SURFACE_WAT[cellLoc] < -0.01) {
                      sprintf(msgStr,"Day %6.1f: capacityERR - neg SURFACE_WAT (%f m) in cell (%d,%d) before cell_dyn7!", 
                              SimTime.TIME, SURFACE_WAT[cellLoc], ix,iy ); 
                      WriteMsg( msgStr,True );  dynERRORnum++; }
              }

/*  note that rainfall during a time step is added to surface water storage and available */
/* for runoff (horizFlow) before the calc of infiltration & ET associated with that new input */
/* (infiltration/ET etc will be of avail water the next time step after a rainfall event and horiz flows) */
                  /* NSM/SFWMM rainfall input data, created in cell_dyn1, convert here from tenths of mm to meters */
              SF_WT_FROM_RAIN[cellLoc] = boundcond_rain[cellLoc]*0.0001;  /*  tenths of mm *0.0001 = m */

              /* solar radiation at altitude of 274m in atmosphere cal/cm2/d) */
              /* v2.2+, CLOUDY (cloudiness) spatiotemporal data "temporarily" unavailable, is constant in space and time, local var "cloudy" */
              SOLRAD274[cellLoc] = SOLRADATMOS*(SOLBETA-GP_SOLOMEGA* ( ( cloudy>0.0 ) ? ( cloudy ) : ( 0.0) ) ) -SOLALPHA;
              SOLRADGRD[cellLoc] = SOLRAD274[cellLoc]+((SOLRADATMOS+1.0)-SOLRAD274[cellLoc])*SOLALTCORR;
              H2O_TEMP[cellLoc] = AIR_TEMP[cellLoc]; /* v2.2+, temperature data "temporarily" unavailable, is constant in space and time */
              
/******** determine new unsat potentials */
              SAT_WT_HEAD[cellLoc]  = SAT_WATER[cellLoc]/HP_HYD_POROSITY[HAB[cellLoc]];
              UNSAT_DEPTH[cellLoc]  = SED_ELEV[cellLoc]-SAT_WT_HEAD[cellLoc];
              UNSAT_CAP[cellLoc]  =  UNSAT_DEPTH[cellLoc]*HP_HYD_POROSITY[HAB[cellLoc]];
          
              UNSAT_MOIST_PRP[cellLoc]  = 
                  ( UNSAT_CAP[cellLoc]>0.0 ) ? 
                  ( Min(UNSAT_WATER[cellLoc]/UNSAT_CAP[cellLoc],1.0) ) : 
                  ( 1.0); 
                  /* determining the pathway of flow of surface water depending on depth
                     of an unsat zone relative to the surface water  */ 
              SAT_VS_UNSAT[cellLoc]  = 1/Exp(100.0*Max((SURFACE_WAT[cellLoc]-UNSAT_DEPTH[cellLoc]),0.0)); 
     /* empirical data of a (0-1) control function, the proportion of maximum vertical water infiltration rate through soil (dependent var) as a function of soil moisture proportion (0-1)  (independent var) */
              UNSAT_HYD_COND_CF[cellLoc]  = graph7(0x0,UNSAT_MOIST_PRP[cellLoc] ); 
                     /* field capacity = porosity - specific yield; spec yield== proportion of total soil vol
                            that represents water that can be moved by gravity */
              field_cap = (HP_HYD_POROSITY[HAB[cellLoc]]-HP_HYD_SPEC_YIELD[HAB[cellLoc]]);
                  /* unsat_avail is proportion of water in pore space available for gravitational flow (above field capacity) */
                  /* e.g., when moisture prop in pore space <= field_cap/pore_space, no percolation */
                  /* using old moisture proportion (hasn't changed unless unsat zone was replaced by sat water) */
              UNSAT_AVAIL[cellLoc]  = Max(UNSAT_MOIST_PRP[cellLoc]
                                          -(field_cap)/HP_HYD_POROSITY[HAB[cellLoc]],0.0);
              UNSAT_WT_POT[cellLoc]  = Max(UNSAT_CAP[cellLoc]-UNSAT_WATER[cellLoc],0.0);

/******** now determine the potential total transpiration and evaporation  */
/* Potential ET is input data used in SFWMM v5.4 */
/* GP_calibET is an adjustable calibration parameter (close to 1.0, adjusted in global parameter input file)  */
             HYD_EVAP_CALC[cellLoc]  =  boundcond_ETp[cellLoc] * 0.0001*GP_calibET;  /*  tenths of mm *0.0001 = m */

    /*  Leaf Area Index (LAI) of emergent macrophytes: this effective LAI estimates leaf area index that is above ponded surface water  */
              LAI_eff[cellLoc] =  
                  (MAC_HEIGHT[cellLoc]>0.0) ? 
                  (Max(1.0 - SURFACE_WAT[cellLoc]/MAC_HEIGHT[cellLoc], 0.0)*MAC_LAI[cellLoc]) : 
                  (0.0)  ;
 
          /* control function (0-1) of relative magnitude of the effective Leaf Area Index  */
              f_LAI_eff = exp(-LAI_eff[cellLoc]); 
              
              
              HYD_TOT_POT_TRANSP[cellLoc]  = HYD_EVAP_CALC[cellLoc] * (1.0-f_LAI_eff); 

     /* 0-1 control function of saturated water available to roots - capillary draw when roots don't quite reach down to water table */
              SatWat_Root_CF =  Exp(-10.0* Max(UNSAT_DEPTH[cellLoc]-HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]],0.0) ); 
     /*  HYDrologic, control function (0-1) of proportion of WATer in upper soil profile that is AVAILable for plant uptake, including unsaturated storage withdrawal, and small capillary withdrawal from saturated storage, depending on relative depths */
              HYD_WATER_AVAIL[cellLoc]  = (UNSAT_DEPTH[cellLoc] > HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]] ) ? 
                  ( Max(UNSAT_MOIST_PRP[cellLoc], SatWat_Root_CF) ) :
                  ( 1.0 ); 
    /* empirical data of a (0-1) control function, the proportion of water available to plants (dependent var) as a function of proportion (0-1) of water available upper soil profile (independent var) (generally, simply 1:1 relationship) */
              MAC_WATER_AVAIL_CF[cellLoc]  = graph8(0x0,HYD_WATER_AVAIL[cellLoc]); 

/******** next calc the potential and actual flows */
/* unsatLoss(1) unsat to sat percolation */
  /*unsat to sat flow here only includes percolation (rising water table accomodated in update after horiz fluxes) */ 
              UNSAT_PERC[cellLoc]  = Min(HP_HYD_RCINFILT[HAB[cellLoc]]*UNSAT_HYD_COND_CF[cellLoc],UNSAT_AVAIL[cellLoc]*UNSAT_WATER[cellLoc]);
              UNSAT_TO_SAT_FL[cellLoc]  =  
                  ( (UNSAT_PERC[cellLoc])*DT > UNSAT_WATER[cellLoc] ) ? 
                  ( UNSAT_WATER[cellLoc]/DT ) : 
                  ( UNSAT_PERC[cellLoc]);
/* unsatLoss(2) unsat to atmosph transpiration */
              HYD_UNSAT_POT_TRANS[cellLoc]  = (UNSAT_DEPTH[cellLoc] > HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]] ) ?
                  ( HYD_TOT_POT_TRANSP[cellLoc]*MAC_WATER_AVAIL_CF[cellLoc] ) :
                  (0.0); /* no unsat transp if roots are in saturated zone */
              UNSAT_TRANSP[cellLoc]  = 
                  ( (HYD_UNSAT_POT_TRANS[cellLoc]+UNSAT_TO_SAT_FL[cellLoc])*DT>UNSAT_WATER[cellLoc] ) ? 
                  ( (UNSAT_WATER[cellLoc]-UNSAT_TO_SAT_FL[cellLoc]*DT)/DT ) : 
                  ( HYD_UNSAT_POT_TRANS[cellLoc]);

/* satLoss(1) sat to deep aquifer recharge **RATE parameter IS ALWAYS SET to ZERO  *****/
              SAT_WT_RECHG[cellLoc]  = 
                  ( GP_HYD_RCRECHG*HP_HYD_SPEC_YIELD[HAB[cellLoc]]/HP_HYD_POROSITY[HAB[cellLoc]]*DT>SAT_WATER[cellLoc] ) ? 
                  ( SAT_WATER[cellLoc]/DT ) : 
                  ( GP_HYD_RCRECHG*HP_HYD_SPEC_YIELD[HAB[cellLoc]]/HP_HYD_POROSITY[HAB[cellLoc]]); 
                 
/* sat to surf upflow  when gwater exceeds soil capacity due to lateral inflows
   accomodated in gwFluxes */

/* satLoss(2) sat to unsat with lowering water table due to recharge loss ONLY (horiz outflows accomodated in gwFluxes)
   (leaves field capacity amount in unsat zone)*/
              SAT_TO_UNSAT_FL[cellLoc]  =  
                  ( SAT_WT_RECHG[cellLoc]*field_cap/HP_HYD_POROSITY[HAB[cellLoc]]*DT > SAT_WATER[cellLoc] ) ? 
                  ( (SAT_WATER[cellLoc])/DT ) : 
                  ( SAT_WT_RECHG[cellLoc]*field_cap/HP_HYD_POROSITY[HAB[cellLoc]]) ;
/* satLoss(3) sat to atmosph */
              HYD_SAT_POT_TRANS[cellLoc]  = HYD_TOT_POT_TRANSP[cellLoc]*SatWat_Root_CF; 
              SAT_WT_TRANSP[cellLoc]  =  
                  ( (HYD_SAT_POT_TRANS[cellLoc]+SAT_TO_UNSAT_FL[cellLoc])*DT > SAT_WATER[cellLoc] ) ? 
                  ( (SAT_WATER[cellLoc]-SAT_TO_UNSAT_FL[cellLoc]*DT)/DT ) : 
                  ( HYD_SAT_POT_TRANS[cellLoc]);

/* sfwatLoss(1) sf to sat */
                     /* downflow to saturated assumed to occur in situations with small
                        unsat layer overlain by significant surface water (SAT_VS_UNSAT very small),
                        with immediate hydraulic connectivity betweent the two storages */
              SF_WT_TO_SAT_DOWNFLOW[cellLoc]  = 
                  ( (1.0-SAT_VS_UNSAT[cellLoc])*UNSAT_WT_POT[cellLoc]*DT>SURFACE_WAT[cellLoc] ) ? 
                  ( SURFACE_WAT[cellLoc]/DT ) : 
                  ( (1.0-SAT_VS_UNSAT[cellLoc])*UNSAT_WT_POT[cellLoc]); 
/* sfwatLoss(2) sf to unsat infiltration (sat_vs_unsat splits these losses to groundwater, but downflow to sat is given priority) */
              SF_WT_POT_INF[cellLoc]  = 
                  ( (SAT_VS_UNSAT[cellLoc]*HP_HYD_RCINFILT[HAB[cellLoc]]+SF_WT_TO_SAT_DOWNFLOW[cellLoc])*DT>SURFACE_WAT[cellLoc] ) ? 
                  ( (SURFACE_WAT[cellLoc]-SF_WT_TO_SAT_DOWNFLOW[cellLoc]*DT)/DT ) : 
                  ( SAT_VS_UNSAT[cellLoc]*HP_HYD_RCINFILT[HAB[cellLoc]]);
              SF_WT_INFILTRATION[cellLoc]  = 
                  ( SF_WT_POT_INF[cellLoc]*DT > (UNSAT_WT_POT[cellLoc]-SF_WT_TO_SAT_DOWNFLOW[cellLoc]*DT) ) ? 
                  ((UNSAT_WT_POT[cellLoc]-SF_WT_TO_SAT_DOWNFLOW[cellLoc]*DT)/DT ) : 
                  ( SF_WT_POT_INF[cellLoc]);
              sfwat_pr1 = SF_WT_INFILTRATION[cellLoc]+SF_WT_TO_SAT_DOWNFLOW[cellLoc];
/* sfwatLoss(3) sf to atmosph */
              SF_WT_EVAP[cellLoc]  =  
                  ( (f_LAI_eff*HYD_EVAP_CALC[cellLoc]+sfwat_pr1 )*DT>SURFACE_WAT[cellLoc] ) ? 
                  ( (SURFACE_WAT[cellLoc]-sfwat_pr1*DT)/DT ) : 
                  ( f_LAI_eff*HYD_EVAP_CALC[cellLoc]); 


/******** then update the state variable storages */

/* unsat loss priority:  percolation,  transpiration */
/* calc unsaturated storage state var (m) */
              UNSAT_WATER[cellLoc] = UNSAT_WATER[cellLoc] 
                  + (SAT_TO_UNSAT_FL[cellLoc] + SF_WT_INFILTRATION[cellLoc] 
                     - UNSAT_TO_SAT_FL[cellLoc] - UNSAT_TRANSP[cellLoc]) * DT;

/* sat loss priority:  recharge to deep aquifer, re-assign to unsat with lowered water table, transpiration */
/* calc saturated storage state var (m) */
              SAT_WATER[cellLoc] =  SAT_WATER[cellLoc] 
                  + (UNSAT_TO_SAT_FL[cellLoc] + SF_WT_TO_SAT_DOWNFLOW[cellLoc] 
                     - SAT_WT_TRANSP[cellLoc] - SAT_TO_UNSAT_FL[cellLoc] - SAT_WT_RECHG[cellLoc]) * DT;

/* sfwat loss priority: downflow to replace groundwater loss, infiltration to unsat, evaporation */
/* calc surface storage state var (m) */
              SURFACE_WAT[cellLoc] = SURFACE_WAT[cellLoc] 
                  + (SF_WT_FROM_RAIN[cellLoc]  
                     - SF_WT_EVAP[cellLoc] - SF_WT_INFILTRATION[cellLoc] - SF_WT_TO_SAT_DOWNFLOW[cellLoc]) * DT;

/******** lastly, update of head calcs, unsat capacity, moisture proportion, etc.
 in order to calc water in diff storages for solute concentrations */
              SAT_WT_HEAD[cellLoc]  = SAT_WATER[cellLoc]/HP_HYD_POROSITY[HAB[cellLoc]];
              UNSAT_DEPTH[cellLoc]  = Max(SED_ELEV[cellLoc]-SAT_WT_HEAD[cellLoc],0.0);
              UNSAT_CAP[cellLoc]  =  UNSAT_DEPTH[cellLoc]*HP_HYD_POROSITY[HAB[cellLoc]];

              UNSAT_MOIST_PRP[cellLoc]  = 
                  ( UNSAT_CAP[cellLoc]>0.0 ) ? 
                  ( Min(UNSAT_WATER[cellLoc]/UNSAT_CAP[cellLoc],1.0) ) : 
                  ( 1.0); 
              HYD_DOM_ACTWAT_VOL[cellLoc]  = ( Min(HP_DOM_MAXDEPTH[HAB[cellLoc]],UNSAT_DEPTH[cellLoc])*UNSAT_MOIST_PRP[cellLoc] +
                                               Max(HP_DOM_MAXDEPTH[HAB[cellLoc]]-UNSAT_DEPTH[cellLoc], 0.0)*HP_HYD_POROSITY[HAB[cellLoc]] )
                  *CELL_SIZE; 
                  /* flag for presence of small amount of water storage in the active zone must be present */ 
              HYD_DOM_ACTWAT_PRES[cellLoc]  = 
                  ( HYD_DOM_ACTWAT_VOL[cellLoc] > CELL_SIZE*GP_DetentZ ) ? 
                  ( 1.0 ) : ( 0.0); 
              HYD_SED_WAT_VOL[cellLoc]  = (SAT_WATER[cellLoc]+UNSAT_WATER[cellLoc])*CELL_SIZE;
              SFWT_VOL[cellLoc]  = SURFACE_WAT[cellLoc]*CELL_SIZE;

              HydTotHd[cellLoc]  = SAT_WT_HEAD[cellLoc]+SURFACE_WAT[cellLoc]; /* only used for reporting purposes */
              /* v2.8.2 new variable, for reporting purposes only: positive/negative water depth relative to land elevation (stage minus land elevation) */
              HydRelDepPosNeg[cellLoc] =  HydTotHd[cellLoc] - SED_ELEV[cellLoc];
              
                  /* at the square of xx% of the macrophyte's height, the manning's n
                     calc will indicate the macrophyte *starting* to bend over,
                     starting to offer increasingly less resistance */
              mann_height = (GP_mann_height_coef*MAC_HEIGHT[cellLoc])*(GP_mann_height_coef*MAC_HEIGHT[cellLoc]); 
              N_density = Max(HP_MAC_MAXROUGH[HAB[cellLoc]] * MAC_REL_BIOM[cellLoc], HP_MAC_MINROUGH[HAB[cellLoc]] );
                  /* manning's n for overland (horiz) flow */
              mann_height = Max(mann_height,0.01); /* ensure that even in absence of veg, there is  miniscule (1 cm in model grid cell) height related to some form of veg */   
              HYD_MANNINGS_N[cellLoc]  = Max(-Abs((N_density-HP_MAC_MINROUGH[HAB[cellLoc]])
                                                  *(pow(2.0,(1.0-SURFACE_WAT[cellLoc]/mann_height))-1.0) ) 
                                             + N_density,HP_MAC_MINROUGH[HAB[cellLoc]]);

                  /* sum of transpiration for output only */
              HYD_TRANSP[cellLoc]  = UNSAT_TRANSP[cellLoc]+SAT_WT_TRANSP[cellLoc];
              HYD_ET[cellLoc]  = HYD_TRANSP[cellLoc]+SF_WT_EVAP[cellLoc];

            }
  }
    }
 }


/*******/
/*   MAC_NOPH_BIOMAS[cellLoc] == carbon mass of live Non-Photosynthetic tissues of macrophytes (kgC/m^2)
   MAC_PH_BIOMAS[cellLoc] == carbon mass of live Photosynthetic tissues of macrophytes (kgC/m^2)

   mac_nph_P[cellLoc] == phosphorus mass of live Non-Photosynthetic tissues of macrophytes (kgP/m^2)
   mac_ph_P[cellLoc] == phosphorus mass of live Photosynthetic tissues of macrophytes (kgP/m^2)

   mac_nph_OM[cellLoc] == organic matter mass of live Non-Photosynthetic tissues of macrophytes (kgOM/m^2)
   mac_ph_OM[cellLoc] == organic matter mass of live Photosynthetic tissues of macrophytes (kgOM/m^2)

spatial (cell) location defines habitat ( arrayOf[HAB[cellLoc]] );
habitat type can switch when global dynamic function (in cell_dyn1) calls the succession function

Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn8(int step)
 {
int ix, iy, cellLoc; 
float reduc, NPP_P, min_litTemp, nphbio_ddep, phbio_ddep, MAC_PHtoNPH, MAC_PHtoNPH_Init;
float Sal_Mean;
#define salinLow 0.1 /* v2.8.5 minimal salinity value (mg/L=ppt), below which growth isn't really affected for any macrophyre species */


  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) {

    if(ON_MAP[cellLoc= T(ix,iy)])  { 
              
/* these thresholds need updating when a habitat type of a grid cell changes */
      MAC_MAX_BIO[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]+HP_PHBIO_MAX[HAB[cellLoc]];
      NPHBIO_REFUGE[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
      NPHBIO_SAT[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]*0.9;
      PHBIO_REFUGE[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
      PHBIO_SAT[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]]*0.9;
/* various control functions for production calc */
     MAC_LIGHT_CF[cellLoc]  = SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]
         *Exp(1.0-SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]);
     MAC_TEMP_CF[cellLoc]  = tempCF(0, 0.20, HP_MAC_TEMPOPT[HAB[cellLoc]], 5.0, 40.0, AIR_TEMP[cellLoc]);
     HP_MAC_WAT_TOLER[HAB[cellLoc]] = Max(HP_MAC_WAT_TOLER[HAB[cellLoc]],0.005); /* water tolerance is supposed to be non-zero; set to 5mm if user input a 0 */
     MAC_WATER_CF[cellLoc]  = Min(MAC_WATER_AVAIL_CF[cellLoc], 
         Max(1.0-Max( (SURFACE_WAT[cellLoc]-HP_MAC_WAT_TOLER[HAB[cellLoc]])
         /HP_MAC_WAT_TOLER[HAB[cellLoc]],0.0),0.0));
     MAC_NUT_CF[cellLoc]  = 
                     Exp(-GP_mac_uptake_coef * Max(HP_MAC_KSP[HAB[cellLoc]]-TP_SEDWT_CONCACTMG[cellLoc], 0.0)
                         /HP_MAC_KSP[HAB[cellLoc]]) ; /* mg/L */
/* v2.8.5, re-instituted use of salinity control on mac production */
/* NOTE the placeholder salinity CF here !!! */
     Sal_Mean  = (SAL_SED_WT[cellLoc]+SAL_SF_WT[cellLoc])/2.0; 
     MAC_SALT_CF[cellLoc]  =  (Sal_Mean > HP_MAC_SALIN_THRESH[HAB[cellLoc]]) ? (Max(1.0-0.1*Sal_Mean/HP_MAC_SALIN_THRESH[HAB[cellLoc]],0.0) ) : (1.0);
     Max( 1.0-Max(Sal_Mean-(HP_MAC_SALIN_THRESH[HAB[cellLoc]]+salinLow),0.0)/(HP_MAC_SALIN_THRESH[HAB[cellLoc]]+salinLow), 0.0)  ; 

     min_litTemp = Min(MAC_LIGHT_CF[cellLoc], MAC_TEMP_CF[cellLoc]);
     MAC_PROD_CF[cellLoc]  = Min(min_litTemp,MAC_WATER_CF[cellLoc])
          *MAC_NUT_CF[cellLoc] *MAC_SALT_CF[cellLoc];
/* net primary production, kg C/m2/d */
     PHBIO_NPP[cellLoc]  = HP_PHBIO_RCNPP[HAB[cellLoc]]*MAC_PROD_CF[cellLoc]*MAC_PH_BIOMAS[cellLoc]
         * Max( (1.0-MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc]), 0.0);
/* P uptake is dependent on available P and relative to a maximum P:C ratio for the tissue (kg C/m^2/d * P:Cmax * dimless = kgP/m2/d ) */
     NPP_P = PHBIO_NPP[cellLoc]  * HP_PHBIO_IC_PC[HAB[cellLoc]]  * Max(MAC_NUT_CF[cellLoc]*2.0,1.0)
                     * Max(1.0-mac_ph_PC[cellLoc]/HP_PHBIO_IC_PC[HAB[cellLoc]], 0.0);
/* check for available P mass that will be taken up from sed water in active zone (nutCF does not); v2.5 constrained TP_SED_WT_AZ to pos */
     reduc = (NPP_P > 0.0) ? 
                     (Max(TP_SED_WT_AZ[cellLoc],0.0) / ( NPP_P*CELL_SIZE*DT) ) : /* within-plant variable stoichiometry */
                     (1.0);
    /* reduce the production proportionally if excess demand is found */
                /* can have high conc, but low mass of P avail, in presence of high plant biomass and high demand */
     if (reduc < 1.0) {
                     PHBIO_NPP[cellLoc] *= reduc;
                     NPP_P *= reduc;
                 }
                 
/* losses from photobio */
     phbio_ddep = Max(1.0-Max( (PHBIO_SAT[cellLoc]-MAC_PH_BIOMAS[cellLoc]) 
                                           /(PHBIO_SAT[cellLoc]-PHBIO_REFUGE[cellLoc]),0.0),0.0);
     PHBIO_AVAIL[cellLoc]  = MAC_PH_BIOMAS[cellLoc]*phbio_ddep;
     MAC_PHtoNPH_Init = HP_PHBIO_MAX[HAB[cellLoc]] / HP_NPHBIO_MAX[HAB[cellLoc]] ; /*if habitat type changes */
     MAC_PHtoNPH = (MAC_NOPH_BIOMAS[cellLoc]>0.0) ? ( MAC_PH_BIOMAS[cellLoc] / MAC_NOPH_BIOMAS[cellLoc]) : (MAC_PHtoNPH_Init);
                 
     NPHBIO_TRANSLOC_POT[cellLoc]  = (MAC_PHtoNPH>MAC_PHtoNPH_Init) ?
                     (exp(HP_MAC_TRANSLOC_RC[HAB[cellLoc]] *(MAC_PHtoNPH-MAC_PHtoNPH_Init)) - 1.0) :
                     (0.0); 
     NPHBIO_TRANSLOC[cellLoc]  =  
         ( (NPHBIO_TRANSLOC_POT[cellLoc])*DT >PHBIO_AVAIL[cellLoc] ) ? 
         ( (PHBIO_AVAIL[cellLoc])/DT ) : 
         ( NPHBIO_TRANSLOC_POT[cellLoc]);

     PHBIO_MORT_POT[cellLoc]  = HP_PHBIO_RCMORT[HAB[cellLoc]] *PHBIO_AVAIL[cellLoc] 
         *(1.0 + (1.0-MAC_WATER_AVAIL_CF[cellLoc]) )/2.0;
     PHBIO_MORT[cellLoc]  =
                                ( (PHBIO_MORT_POT[cellLoc]+NPHBIO_TRANSLOC[cellLoc])*DT>PHBIO_AVAIL[cellLoc] ) ? 
         ( (PHBIO_AVAIL[cellLoc]-NPHBIO_TRANSLOC[cellLoc]*DT)/DT ) : 
         ( PHBIO_MORT_POT[cellLoc]);


/* losses from non-photobio  */
     nphbio_ddep = Max(1.0-Max((NPHBIO_SAT[cellLoc]-MAC_NOPH_BIOMAS[cellLoc])
                                          /(NPHBIO_SAT[cellLoc]-NPHBIO_REFUGE[cellLoc]),0.0),0.0);
     NPHBIO_AVAIL[cellLoc]  = MAC_NOPH_BIOMAS[cellLoc]*nphbio_ddep; 

                 PHBIO_TRANSLOC_POT[cellLoc]  = (MAC_PHtoNPH<MAC_PHtoNPH_Init) ?
                     (exp(HP_MAC_TRANSLOC_RC[HAB[cellLoc]] *(MAC_PHtoNPH_Init-MAC_PHtoNPH)) - 1.0) :
                     (0.0); 
     PHBIO_TRANSLOC[cellLoc]  =  
         ( (PHBIO_TRANSLOC_POT[cellLoc])*DT >NPHBIO_AVAIL[cellLoc] ) ? 
         ( (NPHBIO_AVAIL[cellLoc])/DT ) : 
         ( PHBIO_TRANSLOC_POT[cellLoc]);
     NPHBIO_MORT_POT[cellLoc]  = NPHBIO_AVAIL[cellLoc]*HP_PHBIO_RCMORT[HAB[cellLoc]]
                     * (1.0 + Max(1.0-MAC_PH_BIOMAS[cellLoc]/HP_PHBIO_MAX[HAB[cellLoc]],0.0) )/2.0; /* decreased mort w/ increasing photobio */
     NPHBIO_MORT[cellLoc]  =
                                ( (PHBIO_TRANSLOC[cellLoc]+NPHBIO_MORT_POT[cellLoc])*DT>NPHBIO_AVAIL[cellLoc] ) ? 
         ( (NPHBIO_AVAIL[cellLoc]-PHBIO_TRANSLOC[cellLoc]*DT)/DT ) : 
         ( NPHBIO_MORT_POT[cellLoc]);
                 

/* calc nonphotosynthetic biomass state var (kgC/m2) */
     MAC_NOPH_BIOMAS[cellLoc] =  MAC_NOPH_BIOMAS[cellLoc] 
                     + (NPHBIO_TRANSLOC[cellLoc] - NPHBIO_MORT[cellLoc] 
                        - PHBIO_TRANSLOC[cellLoc]  ) * DT;
                 
/* calc nonphotosynthetic biomass state var (kgC/m2) */
     MAC_PH_BIOMAS[cellLoc] = MAC_PH_BIOMAS[cellLoc] 
         + (PHBIO_TRANSLOC[cellLoc] + PHBIO_NPP[cellLoc] 
           - PHBIO_MORT[cellLoc] - NPHBIO_TRANSLOC[cellLoc]) * DT;

/* total biomass */
     MAC_TOT_BIOM[cellLoc]  = MAC_PH_BIOMAS[cellLoc]+MAC_NOPH_BIOMAS[cellLoc];
/* book-keeping calcs used in other modules */
     MAC_REL_BIOM[cellLoc]  = 
         ( MAC_TOT_BIOM[cellLoc] > 0.0 ) ? 
         MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc] : 
         0.0001;
     MAC_HEIGHT[cellLoc]  = pow(MAC_REL_BIOM[cellLoc],0.33)*HP_MAC_MAXHT[HAB[cellLoc]];
     MAC_LAI[cellLoc]  = MAC_REL_BIOM[cellLoc]*HP_MAC_MAXLAI[HAB[cellLoc]];
/* note that an "effective" LAI that accounts for water depth is calc'd in hydro module */

/*  now calc the P and organic matter associated with the C fluxes */
                 phbio_npp_P[cellLoc] = NPP_P;     /* within-plant variable stoichiometry */
                 phbio_npp_OM[cellLoc] = PHBIO_NPP[cellLoc] / HP_PHBIO_IC_CTOOM[HAB[cellLoc]]; /* habitat-specfic stoichiometry */

                 phbio_mort_P[cellLoc] = PHBIO_MORT[cellLoc] * mac_ph_PC[cellLoc];
                 phbio_mort_OM[cellLoc] = PHBIO_MORT[cellLoc] / mac_ph_CtoOM[cellLoc];

                 phbio_transl_P[cellLoc] = PHBIO_TRANSLOC[cellLoc] * mac_nph_PC[cellLoc];
                 phbio_transl_OM[cellLoc] = PHBIO_TRANSLOC[cellLoc] / mac_nph_CtoOM[cellLoc];

                 nphbio_transl_P[cellLoc] = NPHBIO_TRANSLOC[cellLoc] * mac_ph_PC[cellLoc];
                 nphbio_transl_OM[cellLoc] = NPHBIO_TRANSLOC[cellLoc] / mac_ph_CtoOM[cellLoc];
                 
                 nphbio_mort_P[cellLoc] = NPHBIO_MORT[cellLoc] * mac_nph_PC[cellLoc];
                 nphbio_mort_OM[cellLoc] = NPHBIO_MORT[cellLoc] / mac_nph_CtoOM[cellLoc];


/* state vars: now calc the P and OM associated with those C state vars */
                 mac_nph_P[cellLoc] = mac_nph_P[cellLoc]
                     + (nphbio_transl_P[cellLoc] - nphbio_mort_P[cellLoc]
                        - phbio_transl_P[cellLoc]  ) * DT;
                 mac_nph_PC[cellLoc] = ( (MAC_NOPH_BIOMAS[cellLoc] > 0.0) && (mac_nph_P[cellLoc] > 0.0)) ?
                     (mac_nph_P[cellLoc] / MAC_NOPH_BIOMAS[cellLoc]) : /* these second mass checks not needed */
                     0.3 * HP_NPHBIO_IC_PC[HAB[cellLoc]]; /* default to 0.3 of max for habitat */
                 mac_nph_PC_rep[cellLoc] = (float)mac_nph_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */
                 

                 mac_nph_OM[cellLoc] = mac_nph_OM[cellLoc]
                     + (nphbio_transl_OM[cellLoc] - nphbio_mort_OM[cellLoc]
                        - phbio_transl_OM[cellLoc] ) * DT;
                 mac_nph_CtoOM[cellLoc] = ( (mac_nph_OM[cellLoc] > 0.0) && (MAC_NOPH_BIOMAS[cellLoc]>0.0)) ?
                     (MAC_NOPH_BIOMAS[cellLoc] / mac_nph_OM[cellLoc]) :
                     HP_NPHBIO_IC_CTOOM[HAB[cellLoc]];

                 mac_ph_P[cellLoc] = mac_ph_P[cellLoc]
                     + (phbio_transl_P[cellLoc] + phbio_npp_P[cellLoc] - phbio_mort_P[cellLoc]
                        - nphbio_transl_P[cellLoc] ) * DT;
                 mac_ph_PC[cellLoc] = ( (MAC_PH_BIOMAS[cellLoc] > 0.0) && (mac_ph_P[cellLoc]>0.0)) ?
                     (mac_ph_P[cellLoc] / MAC_PH_BIOMAS[cellLoc]) :
                     0.3 * HP_PHBIO_IC_PC[HAB[cellLoc]]; /* default to 0.3 of max for habitat */
                 mac_ph_PC_rep[cellLoc] = (float)mac_ph_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */

                 mac_ph_OM[cellLoc] = mac_ph_OM[cellLoc]
                     + (phbio_transl_OM[cellLoc] + phbio_npp_OM[cellLoc] - phbio_mort_OM[cellLoc]
                        - nphbio_transl_OM[cellLoc]  ) * DT;
                 mac_ph_CtoOM[cellLoc] = ( (mac_ph_OM[cellLoc] > 0.0) && (MAC_PH_BIOMAS[cellLoc]>0.0)) ?
                     (MAC_PH_BIOMAS[cellLoc] / mac_ph_OM[cellLoc]) :
                     HP_PHBIO_IC_CTOOM[HAB[cellLoc]];

                 if (debug > 0 && MAC_NOPH_BIOMAS[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NPhoto C biomass (%f kg) in cell (%d,%d)!", 
                              SimTime.TIME, MAC_NOPH_BIOMAS[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
                 if (debug > 0 && MAC_PH_BIOMAS[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg Photo C biomass (%f kg) in cell (%d,%d)!", 
                              SimTime.TIME, MAC_PH_BIOMAS[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
                 if (debug > 0 && mac_nph_P[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NPhoto P biomass (%f kg) in cell (%d,%d)!", 
                              SimTime.TIME, mac_nph_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
                 if (debug > 0 && mac_ph_P[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg Photo P biomass (%f kg) in cell (%d,%d)!", 
                              SimTime.TIME, mac_ph_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}


/******** phosphorus removed from water here *************/
     TP_SEDWT_UPTAKE[cellLoc]  = phbio_npp_P[cellLoc]*CELL_SIZE; 
/* recalc P in sed water state variable (kg P) */
     TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] - (TP_SEDWT_UPTAKE[cellLoc]) * DT;
     TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
                     (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
                  (0.0);
                 
                      /* this is the active zone, where uptake, sorption, and mineralization take place */
                TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] - (TP_SEDWT_UPTAKE[cellLoc]) * DT;
                 TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
                     ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
                     (TP_SED_CONC[cellLoc]); /* g/L */
                 TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m3==mg/L */
              
          }
  }
  }
}


/*******/
/*
TP_SF_WT[cellLoc] == mass of P in surface water (kg P)
TP_SED_WT[cellLoc] == mass of P in sediment/soil (pore) water (kg P)
TP_SORB[cellLoc] == mass of P sorbed to soil (kg P)

Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn9(int step)
 {
 int ix, iy, cellLoc; 
 float TPpartic, TPsettlRat, TP_settl_pot;
 double PO4Pconc, nonPO4Pconc;
 
  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) {

    if(ON_MAP[cellLoc= T(ix,iy)])  {
/* calc concentrations after horiz fluxes */
              TP_SFWT_CONC[cellLoc]  = 
                  ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                  ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                  ( 0.0); /* used in P fluxes for mass balance */
             /* using regression for predicting PO4 from TP */
              PO4Pconc =  Max( TP_SFWT_CONC[cellLoc]*GP_PO4toTP + 0.001*GP_PO4toTPint,0.0);  /* g/l  (note that intercept may be <0)*/
      /* after spatial (horiz) fluxes, recalc the active zone P mass based on old active/total ratio */
              TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
                  (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
                  (0.0);
            /* this is the active zone, where uptake, sorption, and mineralization take place */
             TP_SED_WT_AZ[cellLoc] = TP_SED_CONC[cellLoc] * TP_Act_to_Tot[cellLoc] * HYD_DOM_ACTWAT_VOL[cellLoc];
              TP_SEDWT_CONCACT[cellLoc] =(HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
                  ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
                  ( TP_SED_CONC[cellLoc]);

/* inputs to surf  P (kg P)  */
              /* v2.6, using negative-parameter switch in v2.8:
                - either use rainfall-based P deposition (parameter >= 0.0), 
                or spatially varying, temporally constant, daily rate of total deposition (used if conc parameter is negative)*/
              /* TODO - investigate methods/data for wet and dry deposition */
              TP_Atm_Depos[cellLoc]  = (GP_TP_IN_RAIN < 0.0) ? (TP_AtmosDepos[cellLoc]) : (SF_WT_FROM_RAIN[cellLoc]*CELL_SIZE*GP_TP_IN_RAIN*0.001) ; /* P deposition, kgP/day */

/* calc various loss and/or vertical exchange potentials (kg P) */
              /* diffusive flux */
              TP_UPFLOW_POT[cellLoc]  = /*  advective upflow has been handled in surf-ground integration in fluxes.c  */
                  Max((TP_SEDWT_CONCACT[cellLoc]-PO4Pconc)
                        *GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0); /* document this diffusion constant multiplier (involving conversion from sec/day) */
              TP_UPFLOW[cellLoc]  = 
                  ( (TP_UPFLOW_POT[cellLoc])*DT>TP_SED_WT_AZ[cellLoc] ) ? 
                  ( (TP_SED_WT_AZ[cellLoc])/DT ) : 
                  ( TP_UPFLOW_POT[cellLoc]);
                  /* units = mgP/L */
              TP_K[cellLoc]  = Max(GP_TP_K_SLOPE*TP_SORBCONC[cellLoc]+GP_TP_K_INTER,0.0);
/*fix rate */
              TP_SORB_POT[cellLoc]  = 
                  ( HYD_DOM_ACTWAT_PRES[cellLoc]>0.0 ) ? 
                  ( (double) 0.001 
                    *(TP_K[cellLoc]*(pow(Max(TP_SEDWT_CONCACT[cellLoc],0.0),0.8) )
                      *0.001*(DEPOS_ORG_MAT[cellLoc]*CELL_SIZE+DIM[cellLoc])-TP_SORB[cellLoc] ) ) : 
                  ( 0.0);
              if (TP_SORB_POT[cellLoc]>0.0) {
                  TP_SORBTION[cellLoc]  =  
                      ( (TP_SORB_POT[cellLoc]+TP_UPFLOW[cellLoc])*DT>TP_SED_WT_AZ[cellLoc] ) ? 
                      ( (TP_SED_WT_AZ[cellLoc]-TP_UPFLOW[cellLoc]*DT)/DT ) : 
                      ( TP_SORB_POT[cellLoc]);
              }
              else { /* neg sorption, loss from sorb variable*/
                  TP_SORBTION[cellLoc]  =  
                      ( (-TP_SORB_POT[cellLoc])*DT>TP_SORB[cellLoc] ) ? 
                      ( (-TP_SORB[cellLoc])/DT ) : 
                      ( TP_SORB_POT[cellLoc]);
              }
              
              /* diffusive + advective flux */
              TP_DNFLOW_POT[cellLoc]  = (SF_WT_INFILTRATION[cellLoc]+SF_WT_TO_SAT_DOWNFLOW[cellLoc])
                  *CELL_SIZE*TP_SFWT_CONC[cellLoc]   
                  + Max((PO4Pconc-TP_SEDWT_CONCACT[cellLoc])
                        *GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0) ;
              TP_DNFLOW[cellLoc]  =  
                  ( ( TP_DNFLOW_POT[cellLoc])*DT > TP_SF_WT[cellLoc] ) ? 
                  ( ( TP_SF_WT[cellLoc])/DT ) : 
                  ( TP_DNFLOW_POT[cellLoc]);
/* calc P in sed water state variables (kg P) */
              TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] +
                  ( TP_DNFLOW[cellLoc] - TP_UPFLOW[cellLoc] - TP_SORBTION[cellLoc]) * DT;
             /* this is the active zone, where uptake, sorption, and mineralization take place */
              TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] +
                  (TP_DNFLOW[cellLoc] - TP_UPFLOW[cellLoc] - TP_SORBTION[cellLoc]) * DT;

/* calc P in surface water state variable (kg P) */
              TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] + 
                  (TP_UPFLOW[cellLoc] + TP_Atm_Depos[cellLoc] 
                   - TP_DNFLOW[cellLoc]) * DT;

/* calc P sorbed-to-soil state variable (kg P) */
              TP_SORB[cellLoc] = TP_SORB[cellLoc] + (TP_SORBTION[cellLoc]) * DT;

              TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
                  (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
                  (0.0); /* g/L */
              TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
                  ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
                  (TP_SED_CONC[cellLoc]); /* g/L */
              TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m^3==mg/L */

              TP_SORBCONC[cellLoc] = ((DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc])>0.0) ?
                  ( (double) TP_SORB[cellLoc]*conv_kgTOg / (DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc]) ) :
                  (0.0); /* gP/kgsoil */

              TP_SORBCONC_rep[cellLoc] = TP_SORBCONC[cellLoc] * conv_gTOmg; /* reporting purposes only (g/kg->mg/kg)*/

              TP_SFWT_CONC[cellLoc]  = 
                  ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                  ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                  ( 0.0); /* g/L used in P fluxes for mass balance */
              TP_SFWT_CONC_MG[cellLoc]  = 
                  ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? 
                  (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
                  (0.0); /* g/m^3==mg/L used for reporting and other modules to evaluate P conc when water is present */

/* the following is an empirical method to calculate settling of particulate P out of the water column
   into the sediments.  It uses the fixed ratio of PO4 to TP, but allows for a decreasing proportion of
   TP to be in (large) particulate form as TP concentration drops below a chosen threshold - the sum of
   the TP is considered to be dissolved plus large particulate plus small particulate (that cannot settle out) */
                  /* mass (kg) of particulate fraction of TP, available for settling to sediments */
                  /* using regression for predicting PO4 from TP */
              PO4Pconc =  Max(TP_SFWT_CONC_MG[cellLoc]*GP_PO4toTP + GP_PO4toTPint,0.0);  /* mg/l (note that intercept may be <0)*/
              nonPO4Pconc = Max(TP_SFWT_CONC_MG[cellLoc]-PO4Pconc,0.0); /* non-PO4 conc, mg/l */
              TPpartic = nonPO4Pconc * (1.0-exp(-nonPO4Pconc/GP_TPpart_thresh)) *0.001 * SFWT_VOL[cellLoc] ; /* kg particulate P */


                  /* settling rate, 1/d */
              TPsettlRat = ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ?
                  (GP_settlVel/SURFACE_WAT[cellLoc]) :
                  0.0;
              
                  /* potential settling of particulate TP (kg/d) */
              TP_settl_pot = TPsettlRat * TPpartic;
              TP_settl[cellLoc]  =  
                  ( ( TP_settl_pot)*DT > TPpartic ) ? 
                  ( (TPpartic)/DT ) : 
                  ( TP_settl_pot);
/* calc P in surface water state variable (kg P) */
              TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] - TP_settl[cellLoc] * DT;

/* various book-keeping calcs used in other modules */
/* conc surf and sed wat = kgP/m3==gP/L, another var calc for mg/L */
              TP_SFWT_CONC[cellLoc]  = 
                  ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                  ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                  ( 0.0); /* used in P fluxes for mass balance */
              TP_SFWT_CONC_MG[cellLoc]  = 
                  ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? 
                  (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
                  (0.0); /* g/m3==mg/L used for reporting and other modules to evaluate P conc when water is present */
              
              if (debug > 0 && TP_SF_WT[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg TP_SF_WT (%f kg) in cell (%d,%d)!", 
                                                                    SimTime.TIME, TP_SF_WT[cellLoc], ix,iy ); WriteMsg( msgStr,True ); usrErr( msgStr ); dynERRORnum++; }
              if (debug > 0 && TP_SED_WT_AZ[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg TP_SED_WT_AZ (%f kg) in cell (%d,%d)!", 
                                                                      SimTime.TIME, TP_SED_WT_AZ[cellLoc], ix,iy ); WriteMsg( msgStr,True ); usrErr( msgStr ); dynERRORnum++; }

    }
  }
  }
}

/*******/
/*
Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn13(int step)
 {
 int ix, iy, cellLoc; 
 float TPsettlRat, TP_settl_pot;
 
  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) { 

    if(ON_MAP[cellLoc= T(ix,iy)])  {
/* concentration of P in surface water used in P fluxes for mass transfer (kgP/m3==gP/L) */
              TP_SFWT_CONC[cellLoc]  = 
                  ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                  ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                  ( 0.0); 
/* inputs to surf  P (kg P)  */             
              /* see notes on changes, in cell_dyn9 module */
              TP_Atm_Depos[cellLoc]  = (GP_TP_IN_RAIN < 0.0) ? (TP_AtmosDepos[cellLoc]) : (SF_WT_FROM_RAIN[cellLoc]*CELL_SIZE*GP_TP_IN_RAIN*0.001) ; /* P deposition, kgP/day */

/* TP settling rate calculation (m/d) = 0 if water depth (m) less than depth threshold parm */
              if (SURFACE_WAT[cellLoc] > GP_WQMthresh ) {
                  TPsettlRat = WQMsettlVel[cellLoc];
                  TP_settlDays[cellLoc]++;
              }
              else
                  TPsettlRat = 0.0;
              
/* before we had to put in the day accumulator*/
/*               TPsettlRat = ( SURFACE_WAT[cellLoc] > GP_WQMthresh ) ?  */
/*                   (WQMsettlVel[cellLoc]) : 0.0; */
/* potential settling of particulate TP (m/d * m^2 * kg/m^3 = kg/d) */
              TP_settl_pot = TPsettlRat * CELL_SIZE * TP_SFWT_CONC[cellLoc];

/*  like EWQM, shows mass bal error when ->   TP_settl[cellLoc]  =   TP_settl_pot; */
              
/* constrain settling to be no more than kg P available in water column */
              TP_settl[cellLoc]  =   
                  ( ( TP_settl_pot)*DT > TP_SF_WT[cellLoc] ) ?  
                  ( (TP_SF_WT[cellLoc])/DT ) :  
                  ( TP_settl_pot); 
/* calc P in surface water state variable (kg P) */
              TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] +
                  ( TP_Atm_Depos[cellLoc] - TP_settl[cellLoc] ) * DT;

/* this was in EWQM!!! (mass balance error!):  if (TP_SF_WT[cellLoc]<0.0) TP_SF_WT[cellLoc]=0.0; */
              
/* concentration of P in surface water used in P fluxes for mass transfer (kgP/m3==gP/L) */
              TP_SFWT_CONC[cellLoc]  = 
                  ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                  ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                  ( 0.0); 
/* concentration used for reporting (e.g., in maps) when water is present. 
   Using a new (v2.8.4) depth parm also used for all
   concentration reporting thresholds
*/
              TP_SFWT_CONC_MG[cellLoc]  = 
                  ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ?  
                  (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
                  (0.0);  /* g/m3==mg/L */
    }
  }
  }
}



/*******/
/*
SALT_SURF_WT[cellLoc] == mass of salt dissolved in surface water (kg salt)
SALT_SED_WT[cellLoc] == mass of salt dissolved in sediment/soil (pore) water (kg salt)


Parameter definitions:
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn10(int step)
 {
int ix, iy, cellLoc;
 double SALT_SED_TO_SF_FLOW_pot;
 

  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) {

    if(ON_MAP[cellLoc= T(ix,iy)])  {
     SAL_SF_WT_mb[cellLoc]  = 
                     ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                     ( SALT_SURF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                     ( 0.0); /* used in salt fluxes for mass balance */
     SAL_SED_WT[cellLoc]  = 
                     ( HYD_SED_WAT_VOL[cellLoc]>0.0 ) ? 
                     ( SALT_SED_WT[cellLoc]/HYD_SED_WAT_VOL[cellLoc] ) : 
                     ( 0.0);
      /* v2.8:
         - either use rainfall-based salt (Cl) deposition (parameter >= 0.0), 
           or spatially varying, temporally constant, daily rate of total deposition (used if conc parameter is negative)*/
      /* TODO - investigate methods/data for wet and dry deposition */
              SALT_Atm_Depos[cellLoc]  = (GP_CL_IN_RAIN < 0.0) ? (SALT_AtmosDepos[cellLoc]) : (SF_WT_FROM_RAIN[cellLoc]*CELL_SIZE*GP_CL_IN_RAIN*0.001) ; /* chloride deposition, kgCl/day */

              /* diffusive + advective flux */ 
                 SALT_SFWAT_DOWNFL_POT[cellLoc]  = (SF_WT_INFILTRATION[cellLoc] + SF_WT_TO_SAT_DOWNFLOW[cellLoc])
                     * CELL_SIZE*SAL_SF_WT_mb[cellLoc]
                     + Max((SAL_SF_WT_mb[cellLoc]-SAL_SED_WT[cellLoc])
                           * GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0)  ; /* TODO: get cl diffusion parm; diffusion parms same as P */
     SALT_SFWAT_DOWNFL[cellLoc]  =  
                     ( SALT_SFWAT_DOWNFL_POT[cellLoc]*DT>SALT_SURF_WT[cellLoc] ) ? 
                     ( SALT_SURF_WT[cellLoc]/DT ) : 
                     ( SALT_SFWAT_DOWNFL_POT[cellLoc]);

              /* diffusive flux */  
                 SALT_SED_TO_SF_FLOW_pot  =  
                    /*  advective upflow has been handled in surf-ground integration in fluxes.c  */
                     Max((SAL_SED_WT[cellLoc]-SAL_SF_WT_mb[cellLoc])
                            *GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0)  ; /* TODO: get cl diffusion parm; diffusion parms same as P */
                 SALT_SED_TO_SF_FLOW[cellLoc]  =  
                     ( SALT_SED_TO_SF_FLOW_pot*DT>SALT_SED_WT[cellLoc] ) ? 
                     ( SALT_SED_WT[cellLoc]/DT ) : 
                     ( SALT_SED_TO_SF_FLOW_pot );
/* calc state vars (kg salt) */
     SALT_SED_WT[cellLoc] =  SALT_SED_WT[cellLoc]  
                     + (SALT_SFWAT_DOWNFL[cellLoc] - SALT_SED_TO_SF_FLOW[cellLoc]) * DT;

     SALT_SURF_WT[cellLoc] = SALT_SURF_WT[cellLoc] 
                     + (SALT_Atm_Depos[cellLoc] + SALT_SED_TO_SF_FLOW[cellLoc] - SALT_SFWAT_DOWNFL[cellLoc] ) * DT;

/* book-keeping concentration calcs, (kg/m3==g/L==ppt) used in other modules */
     SAL_SF_WT_mb[cellLoc]  = 
                     ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                     ( SALT_SURF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                     ( 0.0); /* used in salt fluxes for mass balance */
     SAL_SF_WT[cellLoc]  = 
                     ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? 
                     ( SAL_SF_WT_mb[cellLoc]  ) : 
                     ( 0.0); /* used for reporting and other modules to evaluate salinity when water is present */
     SAL_SED_WT[cellLoc]  = 
                     ( HYD_SED_WAT_VOL[cellLoc]>0.0 ) ? 
                     ( SALT_SED_WT[cellLoc]/HYD_SED_WAT_VOL[cellLoc] ) : 
                     ( 0.0);

    }
  }
  }
                  
}



/*
   FLOC[cellLoc] == mass of organic flocculent material at the interface between soil and surface water (kg OM/m^2) \n
   FlocP[cellLoc] == mass of phosphorus in the flocculent material at the interface between soil and surface water (kg P/m^2) \n

This routine originally was Suspended Organic Matter; was modified to instead
represent the organic matter in the flocculent sediment layer. \n

Parameter definitions: \n
      global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms) \n
      habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
*/
void cell_dyn12(int step)
 {
 int ix, iy, cellLoc; 
 float FlocP_DECOMP_pot, FlocP_DEPO_pot, FlocP_settl, Floc_settl;
 
  for(ix=1; ix<=s0; ix++) {
  for(iy=1; iy<=s1; iy++) {

    if(ON_MAP[cellLoc= T(ix,iy)])  {
/* inputs (kg OM / m2)  */
                  /* all periphyton mortality goes to floc */
              FLOC_FR_ALGAE[cellLoc]  = (double) (C_ALG_MORT[cellLoc]+NC_ALG_MORT[cellLoc])
                  /GP_ALG_C_TO_OM*0.001 ; 
                  /* all photobiomass mortality goes to floc */
              Floc_fr_phBio[cellLoc]  = phbio_mort_OM[cellLoc];

             /* all settling from water column goes to floc */
              FlocP_settl = TP_settl[cellLoc] / CELL_SIZE; /* kg P/m2; */
                  /* the particulate P settling is assumed a fixed Redfield P:OM ratio */
              Floc_settl =   FlocP_settl / GP_TP_P_OM; 
                 
          
/* outputs (kg OM / m2) */
              FLOC_DECOMP_QUAL_CF[cellLoc]  = /* use the avg conc of sed and sf water here */
                  Exp(-GP_DOM_decomp_coef * Max(GP_DOM_DECOMP_POPT-
                                 (TP_SFWT_CONC_MG[cellLoc]+TP_SEDWT_CONCACTMG[cellLoc])/2.0, 0.0)
                      /GP_DOM_DECOMP_POPT) ; /* mg/L */
/* decomp in surface water has higher rate than in soil,
 assuming this stock is of much higer substrate quality than the total soil layer */
/* GP_calibDecomp is an adjustable calib parm */
              soil_MOIST_CF[cellLoc]  =  (UNSAT_DEPTH[cellLoc]>HP_DOM_AEROBTHIN[HAB[cellLoc]]) ?
                     ( Max(UNSAT_MOIST_PRP[cellLoc],0.0) ) :
                     ( 1.0 );
              FLOC_DECOMP_POT[cellLoc]  = GP_calibDecomp * 10.0*GP_DOM_RCDECOMP*FLOC[cellLoc]
                   *DOM_TEMP_CF[cellLoc] *FLOC_DECOMP_QUAL_CF[cellLoc] * soil_MOIST_CF[cellLoc];
              FLOC_DECOMP[cellLoc]  = 
                  ( (FLOC_DECOMP_POT[cellLoc])*DT>FLOC[cellLoc] ) ? 
                  ( (FLOC[cellLoc])/DT ) : 
                  ( FLOC_DECOMP_POT[cellLoc]);

/* the incorporation of the floc layer into the "true" soil layer occurs
   at a rate dependent on the floc depth under flooded conditions, then constant rate under dry conditions */
              FLOC_DEPO_POT[cellLoc]  = 
                  ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
                  ( FLOC_Z[cellLoc]/GP_FlocMax * FLOC[cellLoc]*GP_Floc_rcSoil ) : 
                  ( FLOC[cellLoc]*GP_Floc_rcSoil);
              FLOC_DEPO[cellLoc]  = 
                  ( (FLOC_DEPO_POT[cellLoc]+FLOC_DECOMP[cellLoc])*DT>FLOC[cellLoc] ) ? 
                  ( (FLOC[cellLoc]-FLOC_DECOMP[cellLoc]*DT)/DT ) : 
                  ( FLOC_DEPO_POT[cellLoc]); 
/* calc the state var (kg OM / m2) */
              FLOC[cellLoc] =  FLOC[cellLoc] 
                  + ( Floc_settl + Floc_fr_phBio[cellLoc] + FLOC_FR_ALGAE[cellLoc]
                      - FLOC_DECOMP[cellLoc] - FLOC_DEPO[cellLoc] ) * DT;

/* the depth of floc is dependent on a fixed floc bulk density */
              FLOC_Z[cellLoc] = (double) FLOC[cellLoc] / GP_Floc_BD ;
                 

/* Floc phosphorus (kg P / m2)  */
              FlocP_FR_ALGAE[cellLoc]  = (double) (NC_ALG_MORT_P[cellLoc]
                                          + C_ALG_MORT_P[cellLoc])*0.001; /* kg P/m2 */
              FlocP_PhBio[cellLoc] = phbio_mort_P[cellLoc] ;    

              FlocP_DECOMP_pot =  FLOC_DECOMP[cellLoc] * FlocP_OM[cellLoc];
              FlocP_DECOMP[cellLoc]  = 
                  ( (FlocP_DECOMP_pot)*DT>FlocP[cellLoc] ) ? 
                  ( (FlocP[cellLoc])/DT ) : 
                  ( FlocP_DECOMP_pot); 
              FlocP_DEPO_pot =  FLOC_DEPO[cellLoc] * FlocP_OM[cellLoc];
              FlocP_DEPO[cellLoc]  = 
                  ( (FlocP_DEPO_pot+FlocP_DECOMP[cellLoc])*DT>FlocP[cellLoc] ) ? 
                  ( (FlocP[cellLoc]-FlocP_DECOMP[cellLoc]*DT)/DT ) : 
                  ( FlocP_DEPO_pot); 
              
/* calc the state var (kg P / m2) */
              FlocP[cellLoc] =  FlocP[cellLoc] 
                  + ( FlocP_settl + FlocP_PhBio[cellLoc] + FlocP_FR_ALGAE[cellLoc]
                      - FlocP_DECOMP[cellLoc] - FlocP_DEPO[cellLoc] ) * DT;

              FlocP_OM[cellLoc] = ( FLOC[cellLoc]>0.0) ? (FlocP[cellLoc]/FLOC[cellLoc]) : (0.0); /* kgP/kgOM */
              FlocP_OMrep[cellLoc] = (float) FlocP_OM[cellLoc] * conv_kgTOmg; /* mgP/kgOM, variable for output _rep-orting only */
              
              if (debug > 0 && FLOC[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg FLOC OM biomass (%f kg/m2) in cell (%d,%d)!", 
                                                                    SimTime.TIME, FLOC[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
              if (debug > 0 && FlocP[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg FLOC P biomass (%f kg/m2) in cell (%d,%d)!", 
                                                                      SimTime.TIME, FlocP[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}

/* now the P gain in sediment pore water with decomp - 90% goes to porewater, 10% to sfwat */
              TP_SED_MINER[cellLoc]  =  0.90 * FlocP_DECOMP[cellLoc] * CELL_SIZE ; 
/* calc P in sed water state variables (kg P) */
              TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] + 
                  (TP_SED_MINER[cellLoc]) * DT;
            /* this is the active zone, where uptake, sorption, and mineralization take place */
              TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] + 
                  (TP_SED_MINER[cellLoc]) * DT;

              TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
                  (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
                  (0.0);
               TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
                  ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
                  (TP_SED_CONC[cellLoc]);
              TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m3==mg/L */
              
              
/* now the P gain in surface water with decomp - 90% goes to porewater, 10% to sfwat */
              TP_SFWT_MINER[cellLoc]  = 0.10*FlocP_DECOMP[cellLoc] * CELL_SIZE ;  
/* calc P in surface water state variable (kg P) */
              TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] + 
                  (TP_SFWT_MINER[cellLoc]) * DT;
              TP_SFWT_CONC[cellLoc]  = 
                  ( SFWT_VOL[cellLoc] > 0.0 ) ? 
                  ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
                  ( 0.0); /* used in P fluxes for mass balance */
              TP_SFWT_CONC_MG[cellLoc]  = 
                  ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? 
                  (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
                  (0.0); /* g/m3==mg/L used for reporting and other modules to evaluate P conc when water is present */
    }
  }
  }
}


float tempCF(int form, float c1, float topt, float tmin, float tmax, float tempC) 
{
  if (form == 1) {
    /* Lassiter et al. 1975, where c1 = curvature parm */
    return (Exp(c1*(tempC-topt)) * pow(((tmax-tempC)/(tmax-topt)), (c1*(tmax-topt))) );
  }
  else {
    /* Jorgensen 1976 */
    return ( Exp(-2.3 * Abs((tempC-topt)/(topt-tmin))) );
  }
}

void init_static_data(void)
{
  usrErr0("Acquiring static spatial data..."); /* console message */
  
  read_map_file("ModArea",ON_MAP,'c',4.0,0.0);            /* defines model area, dimless unsigned char attributes, value 1 set to 4, 0=0 */
  read_map_file("BoundCond",BCondFlow,'d',1.0,0.0);       /* boundary condition flow restrictions, dimless integer attributes */
            /* ONLY when running the EWQM settling rate version of ELM */
  if (ESPmodeON) read_map_file("basinSetVel",WQMsettlVel,'f',0.0001,0.0);       /* basin-based settling rates (data in tenths mm/d, converted to m/d) */
  read_map_file("basins",basn,'d',1.0,0.0);               /* Basins/Indicator-Region map, dimless integer attributes  */

  alloc_mem_stats (); /* allocate memory for budget & stats arrays (need to have read the Basin/Indicator-Region map first) */
  BIRinit(); /* Set up the Basin/Indicator-Region (BIR) linkages/inheritances */
  BIRoutfiles(); /* Basin/Indicator-Region output files */
  
  usrErr("Done.");
} 


void init_dynam_data(void)
{
  usrErr0("Acquiring dynamic spatial data..."); /* console message */
  
  read_map_file("Elevation",ELEVATION,'f',0.01,0.0);      /*  positive elevation relative to zero of vertical datum; input data in cm, converted to m; (prior to ELMv2.8, datum was always NGVD 1929; some apps in v2.8+ are in NAVD88; datum is independent of source code ) */
  read_map_file("HAB",HAB,'c',1.0,0.0);                   /* habitat (veg classifications, dimless unsigned char attributes) */
  read_map_file("icMacBio",MAC_TOT_BIOM,'f',0.015,0.0);      /* initial total biomass of macrophytes (data in xyz, converted to kg C/m2) */
  read_map_file("icSfWt",SURFACE_WAT,'f',0.01,0.0);       /* initial ponded surface water (data in cm, converted to m) */
  read_map_file("icUnsat",UNSAT_DEPTH,'f',0.01,0.0);      /* initial depth of unsaturated zone (data in cm, converted to m) */
  read_map_file("soilTP",TPtoSOIL,'f',0.000001,0.0);  /* soil TP map (data in mgP/kgSoil, converted to kgP/kgSoil) */
  read_map_file("soilBD",BulkD,'f',1.0,0.0);         /* soil bulk dens map (data in mgSoil/cm3 == kgSoil/m3)  */
  read_map_file("soil_orgBD",DOM_BD,'f',1.0,0.0);    /* organic soil bulk dens map (data in mgOM/cm3 == kgOM/m3)  */
  read_map_file("soilTPpore",TP_SEDWT_CONCACT,'f',0.000001,0.0);    /* porewater TP concentration map (data in ug/L, converted to g/L == kg/m3)  */
  read_map_file("AtmosPdepos",TP_AtmosDepos,'f',1.0,0.0);    /* rate of atmospheric deposition of total phosphorus (data in mgP/m^2/yr, convert units later with eqn inits)  */

  /* 2 static mapps need to be here for re-initialing w/ a multiplier (sensitivity analysis) */ 
  read_map_file("HydrCond",HYD_RCCONDUCT,'f',1.0,0.0);   /* hydraulic conductivity (no conversion, data in m/d)  */
  read_map_file("Bathymetry",Bathymetry,'f',0.01,0.0);      /* positive bathymetry (depth) relative to zero of vertical datum used in land elevation surveys (i.e., positive depth below NGVD29 or NAVD88); all positive depth data in cm, converted to m) */

  usrErr("Done."); /* console message */
} 


void init_eqns(void)
 {
int ix,iy, cellLoc;
float tmp; /* used in checking nutrient availability for MichMent kinetics */
float min_litTemp; /* used to determine min of temper, light cf's for alg and macs */
float I_ISat, Z_extinct, PO4Pconc, PO4P;
float MAC_PHtoNPH, MAC_PHtoNPH_Init;
double soil_propOrg; 
#define satDensity 0.9 /* assign the relative proportion (0 - 1) of a population maximum to be the saturation density for a population */

  usrErr0("Initializing unit model equations..."); /* console message */
  
  SimTime.TIME = 0;
  DAYJUL = ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 365.0);
  LATRAD = ((int)(GP_LATDEG)+(GP_LATDEG-(int)(GP_LATDEG))*5.0/3.0)*PI/180.0;
  /* AMPL = Exp(7.42+0.045*LATRAD*180.0/PI)/3600.0; */
  /* DAYLENGTH = AMPL*Sin((DAYJUL-79.0)*0.01721)+12.0; */
  SOLALPHA = 32.9835-64.884*(1.0-1.3614*Cos(LATRAD));
  SOLDEC1 = 0.39785*Sin(4.868961+0.017203*DAYJUL   +0.033446*Sin(6.224111+0.017202*DAYJUL));
  SOLCOSDEC = sqrt(1.0-pow(SOLDEC1,2.0));
  SOLELEV_SINE = Sin(LATRAD)*SOLDEC1+Cos(LATRAD)*SOLCOSDEC;
  SOLALTCORR = (1.0-Exp(-0.014*(GP_ALTIT-274.0)/(SOLELEV_SINE*274.0)));
  SOLBETA = 0.715-0.3183*(1.0-1.3614*Cos(LATRAD));
  SOLDEC = Arctan(SOLDEC1/sqrt(1.0-pow(SOLDEC1,2.0)));
  SOLRISSET_HA1 = -Tan(LATRAD)*Tan(SOLDEC);
  SOLRISSET_HA = ( (SOLRISSET_HA1==0.0) ) ? ( PI*0.5 ) : (   ( (SOLRISSET_HA1<0.0) ) ? ( PI+Arctan(sqrt(1.0-pow(SOLRISSET_HA1,2.0))/SOLRISSET_HA1)  ) : (      Arctan(sqrt(1.0-pow(SOLRISSET_HA1,2.0))/SOLRISSET_HA1)));
  SOLRADATMOS = 458.37*2.0*(1.0+0.033*Cos(360.0/365.0*PI/180.0*DAYJUL))* (Cos(LATRAD)*Cos(SOLDEC)*Sin(SOLRISSET_HA) + SOLRISSET_HA*180.0/(57.296*PI)*Sin(LATRAD)*Sin(SOLDEC));
        
        for(ix=0; ix<=s0+1; ix++) {
            for(iy=0; iy<=s1+1; iy++) { 
                  cellLoc = T(ix,iy);
                  
                  AIR_TEMP[cellLoc] = 25.0;  /* spatial time series unavailable after 1995; globally constant in v2.2+ */
/* rainfall read from sfwmm data, remapped to ELM resolution */
                  boundcond_rain[cellLoc] =  SF_WT_FROM_RAIN[cellLoc] = boundcond_ETp[cellLoc] = boundcond_depth[cellLoc] = 0.0;
       
                                /* used to have cloudiness influence on GP_SOLOMEGA, now 0 */
                 SOLRAD274[cellLoc] = SOLRADATMOS*(SOLBETA-GP_SOLOMEGA*0.0 )-SOLALPHA;
                 SOLRADGRD[cellLoc] = SOLRAD274[cellLoc]+((SOLRADATMOS+1.0)-SOLRAD274[cellLoc])*SOLALTCORR;
                 H2O_TEMP[cellLoc] = AIR_TEMP[cellLoc];
         
                 ALG_REFUGE[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*GP_ALG_REF_MULT;
                 ALG_SAT[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*0.9;

              /* v2.3: with southern everglades topo, put bathymetry back into the model */
                 ELEVATION[cellLoc] = ELEVATION[cellLoc] * GP_IC_ELEV_MULT;
                 Bathymetry[cellLoc] = Bathymetry[cellLoc] * GP_IC_BATHY_MULT;
                 SED_ELEV[cellLoc] =  ELEVATION[cellLoc] - Bathymetry[cellLoc] + GP_DATUM_DISTANCE; 
                 SED_INACT_Z[cellLoc] = SED_ELEV[cellLoc]-HP_DOM_MAXDEPTH[HAB[cellLoc]]; 

                 MAC_MAX_BIO[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]+HP_PHBIO_MAX[HAB[cellLoc]];
                 NPHBIO_REFUGE[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
                 NPHBIO_SAT[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]] * satDensity;
                 PHBIO_REFUGE[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
                 PHBIO_SAT[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]] * satDensity;
            }
        }
        for(ix=1; ix<s0+1; ix++) {
            for(iy=1; iy<s1+1; iy++) { cellLoc = T(ix,iy);
/*initialization of major state variables */

              /* hydro */
                 HYD_RCCONDUCT[cellLoc] = HYD_RCCONDUCT[cellLoc] * GP_calibGWat;

/* map */        UNSAT_DEPTH[cellLoc] = Max(UNSAT_DEPTH[cellLoc] + GP_HYD_IC_UNSAT_ADD,0.0); /* m */
                 SAT_WT_HEAD[cellLoc] = SED_ELEV[cellLoc]- UNSAT_DEPTH[cellLoc]; /* m */
                 SAT_WATER[cellLoc] = SAT_WT_HEAD[cellLoc] * HP_HYD_POROSITY[HAB[cellLoc]]; /* m */

                 UNSAT_WATER[cellLoc] = HP_HYD_POROSITY[HAB[cellLoc]]*UNSAT_DEPTH[cellLoc] *GP_HYD_ICUNSATMOIST; /* m */
                 UNSAT_CAP[cellLoc] = UNSAT_DEPTH[cellLoc]*HP_HYD_POROSITY[HAB[cellLoc]]; /* m */
                 UNSAT_MOIST_PRP[cellLoc]  =  /* dimless proportion, 0-1 */
                     ( UNSAT_CAP[cellLoc]>0.0 ) ? 
                     ( Min(UNSAT_WATER[cellLoc]/UNSAT_CAP[cellLoc],1.0) ) : 
                     ( 1.0); 

/* map */        SURFACE_WAT[cellLoc] =  Max(SURFACE_WAT[cellLoc] + GP_HYD_IC_SFWAT_ADD, 0.0); /* m */
                 SFWT_VOL[cellLoc] = SURFACE_WAT[cellLoc]*CELL_SIZE;

                 HYD_DOM_ACTWAT_VOL[cellLoc] = Max(HP_DOM_MAXDEPTH[HAB[cellLoc]]-UNSAT_DEPTH[cellLoc]*
                                                   (1.0-UNSAT_MOIST_PRP[cellLoc]),0.0)*HP_HYD_POROSITY[HAB[cellLoc]]*CELL_SIZE;
                 HYD_SED_WAT_VOL[cellLoc] = (SAT_WATER[cellLoc]+UNSAT_WATER[cellLoc])*CELL_SIZE;
                 
              /* soil */
/* map */        DOM_BD[cellLoc] = DOM_BD[cellLoc] * GP_IC_DOM_BD_MULT; /* kg/m3 */
/* map */        BulkD[cellLoc] = BulkD[cellLoc] * GP_IC_BulkD_MULT;    /* kgOM/m3 */

                 soil_propOrg = (double) DOM_BD[cellLoc] / BulkD[cellLoc];
                 DIM[cellLoc] =   (1.0 - soil_propOrg) * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] * CELL_SIZE; /* kg inorganic matter */
                 DEPOS_ORG_MAT[cellLoc] = soil_propOrg * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] ; /* kg organic matter/m2 */
//                 DEPOS_ORG_MAT[cellLoc] = soil_propOrg * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] ; /* kg organic matter/m2 */
/* v2.8.4 new DEPOS_ORG_MAT init */

//                 Inorg_Z[cellLoc] = (1.0 - soil_propOrg) * HP_DOM_MAXDEPTH[HAB[cellLoc]]; /*  fixed inorganic depth (m) */
//                 DOM_Z[cellLoc] =   HP_DOM_MAXDEPTH[HAB[cellLoc]] - Inorg_Z[cellLoc]; /* m */
                 Inorg_Z[cellLoc] = DIM[cellLoc] / BulkD[cellLoc] * (1.0-soil_propOrg)  / CELL_SIZE ; /*  kgInorg / kgSoil/m3 * proprInorg / m2 == fixed inorganic depth (m) */
                 DOM_Z[cellLoc] =   HP_DOM_MAXDEPTH[HAB[cellLoc]] - Inorg_Z[cellLoc]; /* m */

// v2.8.4                 DEPOS_ORG_MAT[cellLoc] = DOM_BD[cellLoc]*DOM_Z[cellLoc]; /* (mgOM/cm3 ==> kgOM/m3) * m = kgOM/m2 */

                 DOM_SED_AEROB_Z[cellLoc] = Min(UNSAT_DEPTH[cellLoc]+HP_DOM_AEROBTHIN[HAB[cellLoc]],HP_DOM_MAXDEPTH[HAB[cellLoc]]); /* m */
                 DOM_SED_ANAEROB_Z[cellLoc] = HP_DOM_MAXDEPTH[HAB[cellLoc]]-DOM_SED_AEROB_Z[cellLoc]; /* m */

/* map */        TPtoSOIL[cellLoc] = TPtoSOIL[cellLoc] * GP_IC_TPtoSOIL_MULT; /* kgP/kgSoil */
                 DOP[cellLoc] =  (double) (1.0-GP_sorbToTP) * TPtoSOIL[cellLoc] * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] ; /* kgP/kg_soil * kg_soil/m3 * m == kgP/m2 */

                 TP_SORB[cellLoc] =   (double)  GP_sorbToTP * TPtoSOIL[cellLoc] * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] * CELL_SIZE; /* dimless * kgP/kg_soil * kg_soil/m3 * m * m^2 == kgP */

              /* floc layer of soil */
                 FLOC[cellLoc] = HP_FLOC_IC[HAB[cellLoc]]; /* kg OM/m2  */
                 FlocP[cellLoc] = FLOC[cellLoc]*HP_FLOC_IC_PC[HAB[cellLoc]]*HP_FLOC_IC_CTOOM[HAB[cellLoc]]; /* kg P /m2 */
                 FLOC_Z[cellLoc] = (double) FLOC[cellLoc] / GP_Floc_BD ; /* m */

              /* phosphorus */
                 TP_AtmosDepos[cellLoc]  = (double)TP_AtmosDepos[cellLoc] * conv_mgTOkg * CELL_SIZE / 365.0; /* P deposition, mgP/m2/yr * kg/mg * m2 / daysPerYear = kgP/day */
/* v2.4.4 */       TP_SFWT_CONC[cellLoc]  = GP_TP_ICSFWAT * conv_mgTOg; /* mg/L * g/mg  => g/L */
                 TP_SF_WT[cellLoc] =TP_SFWT_CONC[cellLoc] * SFWT_VOL[cellLoc]; /* kg P */
                 TP_SFWT_CONC_MG[cellLoc] = TP_SFWT_CONC[cellLoc] * conv_gTOmg; /* mg/L */
                      /* using regression for predicting PO4 from TP */
                 PO4Pconc = Max(TP_SFWT_CONC_MG[cellLoc]*GP_PO4toTP + GP_PO4toTPint, 0.10 * TP_SFWT_CONC_MG[cellLoc]); /* mg/L */
                 PO4P = PO4Pconc * SFWT_VOL[cellLoc] /conv_kgTOg; /*kg P available (from conc. in mg/l) */

/* v2.4.4 */ //      TP_SED_CONC[cellLoc] = GP_TP_ICSEDWAT * conv_mgTOg; /* mg/L * g/mg => g/L */
//                 TP_SED_WT[cellLoc] = TP_SED_CONC[cellLoc] * HYD_SED_WAT_VOL[cellLoc]; /* kg P */
//                      /* this is the active zone, where uptake, sorption, and mineralization take place */
//                  TP_Act_to_Tot[cellLoc] = 1.0 / HP_TP_CONC_GRAD[HAB[cellLoc]]; /* ratio of TP conc in active zone relative to total */
//                  TP_SED_WT_AZ[cellLoc] = TP_SED_CONC[cellLoc] * TP_Act_to_Tot[cellLoc] * HYD_DOM_ACTWAT_VOL[cellLoc]; /* kg P */
//                  TP_SEDWT_CONCACT[cellLoc] =(HYD_DOM_ACTWAT_VOL[cellLoc] > 0.0) ?
//                      ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
//                      ( TP_SED_CONC[cellLoc]); /* g/L */
//                  TP_SEDWT_CONCACTMG[cellLoc] = TP_SEDWT_CONCACT[cellLoc]*conv_gTOmg; /* mg/L */

/* v2.8.4 concentration in Active Zone (TP_SEDWT_CONCACT) is initialed by map (data here are g/L) */
                 TP_Act_to_Tot[cellLoc] = 1.0 / HP_TP_CONC_GRAD[HAB[cellLoc]]; /* ratio of TP conc in active zone relative to total */
                                 TP_SED_CONC[cellLoc] = (TP_Act_to_Tot[cellLoc]>0.0) ? (TP_SEDWT_CONCACT[cellLoc] / TP_Act_to_Tot[cellLoc]) : (0.0);
/* map */        TP_SEDWT_CONCACT[cellLoc] =(HYD_DOM_ACTWAT_VOL[cellLoc] > 0.0) ?
                     ( TP_SEDWT_CONCACT[cellLoc] ) :
                     ( 0.0 ); /* g/L */
                 TP_SEDWT_CONCACTMG[cellLoc] = TP_SEDWT_CONCACT[cellLoc]*conv_gTOmg; /* mg/L */
                 TP_SED_WT[cellLoc] = TP_SED_CONC[cellLoc] * HYD_SED_WAT_VOL[cellLoc]; /* kg P */
                 TP_SED_WT_AZ[cellLoc] = TP_SEDWT_CONCACT[cellLoc] * HYD_DOM_ACTWAT_VOL[cellLoc]; /* kg P */


              /* salt */
                 SALT_SED_WT[cellLoc] = HYD_SED_WAT_VOL[cellLoc]*HP_SALT_ICSEDWAT[HAB[cellLoc]];
                 SAL_SED_WT[cellLoc] = ( HYD_SED_WAT_VOL[cellLoc]>0.0 ) ? ( SALT_SED_WT[cellLoc]/HYD_SED_WAT_VOL[cellLoc] ) : ( 0.0);
                 SALT_SURF_WT[cellLoc] = SFWT_VOL[cellLoc]*HP_SALT_ICSFWAT[HAB[cellLoc]];
                 SAL_SF_WT[cellLoc] = ( SURFACE_WAT[cellLoc] > GP_WQualMonitZ ) ? ( SALT_SURF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : ( 0.0);

              /* periphyton */
/* 2.4.4 */       NC_ALG[cellLoc] = HP_ALG_MAX[HAB[cellLoc]] * GP_ALG_IC_MULT * GP_ALG_REF_MULT ; /* start w/ low, refuge-level of non-calcareous (eutrophic) periphyton, g C/m2 */
/* 2.4.4 */       C_ALG[cellLoc]  = HP_ALG_MAX[HAB[cellLoc]] * GP_ALG_IC_MULT * (1.0 - GP_ALG_REF_MULT); /* g C/m2 */
                 /* ic PC of periph in oligotrophic area is 3% of max P:C, varies across space via (0.1->1.0) map */
/* 2.4.4 */       NC_ALG_PC[cellLoc] = GP_ALG_PC; /* gP/ gC */
/* 2.4.4 */       C_ALG_PC[cellLoc]  = GP_ALG_PC; /* gP/ gC */

                 NC_ALG_P[cellLoc] = NC_ALG[cellLoc]*NC_ALG_PC[cellLoc];   /* g P/m2 */
                 C_ALG_P[cellLoc] = C_ALG[cellLoc]*C_ALG_PC[cellLoc];  /* g P/m2 */  

              /* macrophytes */
/* 2.4.4 */       MAC_PH_BIOMAS[cellLoc] = MAC_TOT_BIOM[cellLoc] * GP_MAC_IC_MULT * HP_PHBIO_MAX[HAB[cellLoc]]/MAC_MAX_BIO[cellLoc]; /* kg C/m2 */
                   /*  now calc the P and OM associated with that C */
/* 2.4.4 */       mac_ph_PC[cellLoc] = HP_PHBIO_IC_PC[HAB[cellLoc]]; 
                 mac_ph_P[cellLoc] = MAC_PH_BIOMAS[cellLoc] * mac_ph_PC[cellLoc]; /* kg P/m2 */
                 mac_ph_OM[cellLoc] = MAC_PH_BIOMAS[cellLoc] / HP_PHBIO_IC_CTOOM[HAB[cellLoc]];
                 mac_ph_CtoOM[cellLoc] = HP_PHBIO_IC_CTOOM[HAB[cellLoc]];
                 PHBIO_AVAIL[cellLoc] = MAC_PH_BIOMAS[cellLoc]*Max(1.0-Max((PHBIO_SAT[cellLoc]-MAC_PH_BIOMAS[cellLoc]) /(PHBIO_SAT[cellLoc]-PHBIO_REFUGE[cellLoc]),0.0),0.0);

/* 2.4.4 */       MAC_NOPH_BIOMAS[cellLoc] = MAC_TOT_BIOM[cellLoc] * GP_MAC_IC_MULT * HP_NPHBIO_MAX[HAB[cellLoc]]/MAC_MAX_BIO[cellLoc]; /* kg C/m2 */
                   /*  now calc the P and OM associated with that C */
/* 2.4.4 */       mac_nph_PC[cellLoc] = HP_NPHBIO_IC_PC[HAB[cellLoc]]; 
                 mac_nph_P[cellLoc] = MAC_NOPH_BIOMAS[cellLoc] * mac_nph_PC[cellLoc];  /* kg P/m2 */
                 mac_nph_OM[cellLoc] = MAC_NOPH_BIOMAS[cellLoc] / HP_NPHBIO_IC_CTOOM[HAB[cellLoc]];
                 mac_nph_CtoOM[cellLoc] = HP_NPHBIO_IC_CTOOM[HAB[cellLoc]];
                 NPHBIO_AVAIL[cellLoc] = MAC_NOPH_BIOMAS[cellLoc]*Max(1.0-Max((NPHBIO_SAT[cellLoc]-MAC_NOPH_BIOMAS[cellLoc])/(NPHBIO_SAT[cellLoc]-NPHBIO_REFUGE[cellLoc]),0.0),0.0);

                 MAC_REL_BIOM[cellLoc] = ( MAC_TOT_BIOM[cellLoc] > 0 ) ? MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc] : 0.0;
                 MAC_LAI[cellLoc] = MAC_REL_BIOM[cellLoc]*HP_MAC_MAXLAI[HAB[cellLoc]];
                 MAC_HEIGHT[cellLoc] = MAC_REL_BIOM[cellLoc]*HP_MAC_MAXHT[HAB[cellLoc]];
                 LAI_eff[cellLoc] =  (MAC_HEIGHT[cellLoc]>0.0) ? (Max(1.0 - SURFACE_WAT[cellLoc]/MAC_HEIGHT[cellLoc], 0.0)*MAC_LAI[cellLoc]) : (0.0)  ;                 

/* end of initialization of major state variables */
                 
/* *************************** */
/* These are the multiple calculations used if particular modules are turned off. \n
 NOTE: THIS section (of init_eqns() ) is not fully updated for properly
 turning off individual **interacting** *biological/chemical* modules.  
 If one *biological/chemical* module is turned off, 
 they all need to be turned off. (Note that cell_dyn's 3,5,6 should always be off). \n

 *** \n
 The following *biological/chemical* modules must be ON or OFF as a group:
 (cell_dyn2 + cell_dyn4 + cell_dyn8 + cell_dyn9  + cell_dyn12)
 cell_dyn13, the net settling rate module, can be turned on only when above module group is off
 *** \n
 
 In particular, the budget will
 not properly reflect actual dynamics if those 
 modules are not treated as a group.
*/
     NC_ALG_MORT_POT[cellLoc] = ( UNSAT_DEPTH[cellLoc]>0.05 ) ? ( NC_ALG[cellLoc]*GP_ALG_RC_MORT_DRY ) : ( NC_ALG[cellLoc]*GP_ALG_RC_MORT);
                     /* calcareous periphyton */
     C_ALG_MORT_POT[cellLoc] = ( UNSAT_DEPTH[cellLoc]>0.05 ) ? ( C_ALG[cellLoc]*GP_ALG_RC_MORT_DRY ) : ( C_ALG[cellLoc]*GP_ALG_RC_MORT);
                 ALG_TEMP_CF[cellLoc]  = tempCF(0, 0.20, GP_ALG_TEMP_OPT, 5.0, 40.0, H2O_TEMP[cellLoc]);
     NC_ALG_RESP_POT[cellLoc]  = 
         ( UNSAT_DEPTH[cellLoc]<0.05 ) ? 
         ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc]*NC_ALG[cellLoc] ) : 
         ( 0.0);
     NC_ALG_RESP[cellLoc]  =  
         ( NC_ALG_RESP_POT[cellLoc]*DT>NC_ALG[cellLoc] ) ? 
         ( NC_ALG[cellLoc]/DT ) : 
         ( NC_ALG_RESP_POT[cellLoc]);
                     /* calcareous periphyton */
     C_ALG_RESP_POT[cellLoc]  = 
         ( UNSAT_DEPTH[cellLoc]<0.05 ) ? 
         ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc]*C_ALG[cellLoc] ) : 
         ( 0.0);
     C_ALG_RESP[cellLoc]  =  
         ( C_ALG_RESP_POT[cellLoc]*DT>C_ALG[cellLoc] ) ? 
         ( C_ALG[cellLoc]/DT ) : 
         ( C_ALG_RESP_POT[cellLoc]);

     NC_ALG_AVAIL_MORT[cellLoc] = NC_ALG[cellLoc]-ALG_REFUGE[cellLoc];
     NC_ALG_MORT[cellLoc] = ( (NC_ALG_MORT_POT[cellLoc])*DT>NC_ALG_AVAIL_MORT[cellLoc] ) ? ( (NC_ALG_AVAIL_MORT[cellLoc])/DT ) : ( NC_ALG_MORT_POT[cellLoc]);
                     /* calcareous periphyton */
     C_ALG_AVAIL_MORT[cellLoc]  = C_ALG[cellLoc]-ALG_REFUGE[cellLoc];
     C_ALG_MORT[cellLoc]  = ( (C_ALG_MORT_POT[cellLoc])*DT>C_ALG_AVAIL_MORT[cellLoc] ) ? ( (C_ALG_AVAIL_MORT[cellLoc])/DT ) : ( C_ALG_MORT_POT[cellLoc]);

/* light, water, temperature controls apply to calc and non-calc periphyton */
     ALG_LIGHT_EXTINCT[cellLoc]  =  0.01; /* light extinction coef */
                     /* algal self-shading implicit in density-dependent constraint function later */
     ALG_INCID_LIGHT[cellLoc]  = SOLRADGRD[cellLoc]*Exp(-MAC_LAI[cellLoc]*GP_ALG_SHADE_FACTOR);
                 Z_extinct = SURFACE_WAT[cellLoc]*ALG_LIGHT_EXTINCT[cellLoc];
     I_ISat = ALG_INCID_LIGHT[cellLoc]/GP_ALG_LIGHT_SAT;
                     /*  averaged over whole water column (based on Steele '65) */
     ALG_LIGHT_CF[cellLoc]  = ( Z_extinct > 0.0 ) ? 
         ( 2.718/Z_extinct * (Exp(-I_ISat * Exp(-Z_extinct)) - Exp(-I_ISat)) ) :
                (I_ISat*Exp(1.0-I_ISat));
                     /*  low-water growth constraint ready for something better based on data */
     ALG_WAT_CF[cellLoc]  = ( SURFACE_WAT[cellLoc]>0.0 ) ? ( 1.0 ) : ( 0.0);
/* the 2 communities have same growth response to avail phosphorus - avail P is roughly calc'd from TP */
     NC_ALG_NUT_CF[cellLoc]  =  PO4Pconc/(PO4Pconc+GP_NC_ALG_KS_P) ;
     C_ALG_NUT_CF[cellLoc]  = PO4Pconc/(PO4Pconc+GP_C_ALG_KS_P); 
      min_litTemp = Min(ALG_LIGHT_CF[cellLoc],ALG_TEMP_CF[cellLoc]);
      NC_ALG_PROD_CF[cellLoc]  = Min(min_litTemp,ALG_WAT_CF[cellLoc])*NC_ALG_NUT_CF[cellLoc];
      C_ALG_PROD_CF[cellLoc]   = Min(min_litTemp,ALG_WAT_CF[cellLoc])*C_ALG_NUT_CF[cellLoc];
/* gross production of the 2 communities (gC/m2, NOT kgC/m2) */
                     /* density constraint contains both noncalc and calc, competition effect accentuated by calc algae */
                     /* used to increase calc growth by factor of 10 */
      NC_ALG_GPP[cellLoc]  =  NC_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*NC_ALG[cellLoc]       
             *Max( (1.0-(GP_AlgComp*C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
     C_ALG_GPP[cellLoc]  =  C_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*C_ALG[cellLoc] 
       *Max( (1.0-(    C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
/* check for available P mass (the MichMent kinetics nutCF do not) */
     tmp = ( (NC_ALG_GPP[cellLoc]+C_ALG_GPP[cellLoc]) > 0) ? 
         PO4P / ( (NC_ALG_GPP[cellLoc]+C_ALG_GPP[cellLoc]) 
         * 0.001*GP_ALG_PC*CELL_SIZE*DT) :
         1.0;
     if (tmp < 1.0) {
         NC_ALG_GPP[cellLoc] *= tmp;   
         C_ALG_GPP[cellLoc] *= tmp; 
    /* can have high conc, but low mass of P avail, in presence of high peri biomass and high demand */
    /* reduce the production proportionally if excess demand is found */
       }
/* total of calc and noncalc algae available and their total NPP */
     NC_ALG_NPP[cellLoc]  = NC_ALG_GPP[cellLoc]-NC_ALG_RESP[cellLoc]; 
     C_ALG_NPP[cellLoc]  = C_ALG_GPP[cellLoc]-C_ALG_RESP[cellLoc]; 


     DOM_QUALITY_CF[cellLoc]  = (Max(TP_SFWT_CONC_MG[cellLoc],TP_SEDWT_CONCACTMG[cellLoc]))
         /GP_DOM_DECOMP_POPT; /* mg/L */
     DOM_TEMP_CF[cellLoc] = Exp(0.20*(H2O_TEMP[cellLoc]-GP_DOM_DECOMP_TOPT))*pow(((40.0-H2O_TEMP[cellLoc])/(40.0-GP_DOM_DECOMP_TOPT)),(0.20*(40.0-GP_DOM_DECOMP_TOPT)));
     soil_MOIST_CF[cellLoc] = pow(Max(UNSAT_MOIST_PRP[cellLoc],0.0),0.75);
     DOM_DECOMP_POT[cellLoc] = GP_DOM_RCDECOMP*DOM_QUALITY_CF[cellLoc]*DOM_TEMP_CF[cellLoc]*DEPOS_ORG_MAT[cellLoc]*(Min(DOM_SED_AEROB_Z[cellLoc]/HP_DOM_MAXDEPTH[HAB[cellLoc]],1.0)*soil_MOIST_CF[cellLoc]+GP_DOM_DECOMPRED*Min(DOM_SED_ANAEROB_Z[cellLoc]/HP_DOM_MAXDEPTH[HAB[cellLoc]],1.0));
     DOM_DECOMP[cellLoc] =  ( (DOM_DECOMP_POT[cellLoc])*DT>DEPOS_ORG_MAT[cellLoc] ) ? ( (DEPOS_ORG_MAT[cellLoc])/DT ) : ( DOM_DECOMP_POT[cellLoc]);
/* added for P DOM stoich */
     DOP_DECOMP[cellLoc] = DOM_DECOMP[cellLoc] * DOM_P_OM[cellLoc]; 

     SAT_VS_UNSAT[cellLoc] = 1/Exp(100.0*Max((SURFACE_WAT[cellLoc]-UNSAT_DEPTH[cellLoc]),0.0));
     UNSAT_WT_POT[cellLoc] = Max(UNSAT_CAP[cellLoc]-UNSAT_WATER[cellLoc],0.0);
         SF_WT_TO_SAT_DOWNFLOW[cellLoc]  = ( (1.0-SAT_VS_UNSAT[cellLoc])*UNSAT_WT_POT[cellLoc]*DT>SURFACE_WAT[cellLoc] ) ? ( SURFACE_WAT[cellLoc]/DT ) : ( (1.0-SAT_VS_UNSAT[cellLoc])*UNSAT_WT_POT[cellLoc]); 
     SAT_WT_RECHG[cellLoc] = ( GP_HYD_RCRECHG*DT>SAT_WATER[cellLoc] ) ? ( SAT_WATER[cellLoc]/DT ) : ( GP_HYD_RCRECHG);
     SF_WT_POT_INF[cellLoc] = ( (SURFACE_WAT[cellLoc]<HP_HYD_RCINFILT[HAB[cellLoc]]*DT) ) ? ( SURFACE_WAT[cellLoc]/DT ) : ( HP_HYD_RCINFILT[HAB[cellLoc]]);
     SF_WT_INFILTRATION[cellLoc] = ( SF_WT_POT_INF[cellLoc]*SAT_VS_UNSAT[cellLoc]*DT>UNSAT_WT_POT[cellLoc] ) ? ( (UNSAT_WT_POT[cellLoc])/DT ) : ( SF_WT_POT_INF[cellLoc]*SAT_VS_UNSAT[cellLoc]);
     HYD_DOM_ACTWAT_PRES[cellLoc] = ( HYD_DOM_ACTWAT_VOL[cellLoc] > 0.01 ) ? ( 1.0 ) : ( 0.0);
     HYD_WATER_AVAIL[cellLoc] = Min(1.0, UNSAT_MOIST_PRP[cellLoc]+Exp(-10.0*Max(UNSAT_DEPTH[cellLoc]-HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]],0.0)));

     MAC_LIGHT_CF[cellLoc] = SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]*Exp(1.0-SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]);
     MAC_TEMP_CF[cellLoc] = Exp(0.20*(AIR_TEMP[cellLoc]-HP_MAC_TEMPOPT[HAB[cellLoc]])) *pow(((40.0-AIR_TEMP[cellLoc])/(40.0-HP_MAC_TEMPOPT[HAB[cellLoc]])),(0.20*(40.0-HP_MAC_TEMPOPT[HAB[cellLoc]])));
     MAC_WATER_AVAIL_CF[cellLoc] = graph8(0x0,HYD_WATER_AVAIL[cellLoc]);
     MAC_WATER_CF[cellLoc] = Min(MAC_WATER_AVAIL_CF[cellLoc], Max(1.0-Max((SURFACE_WAT[cellLoc]-HP_MAC_WAT_TOLER[HAB[cellLoc]])/HP_MAC_WAT_TOLER[HAB[cellLoc]],0.0),0.0));
     MAC_NUT_CF[cellLoc] =  TP_SEDWT_CONCACT[cellLoc]/(TP_SEDWT_CONCACT[cellLoc]+HP_MAC_KSP[HAB[cellLoc]]*0.001) ;

     MAC_SALT_CF[cellLoc] = ( HP_MAC_SALIN_THRESH[HAB[cellLoc]]>0.0 ) ? (  Max( 1.0-Max(SAL_SED_WT[cellLoc]-HP_MAC_SALIN_THRESH[HAB[cellLoc]],0.0)/HP_MAC_SALIN_THRESH[HAB[cellLoc]],0.0) ) : ( Max(1.0-SAL_SED_WT[cellLoc],0.0));
     min_litTemp = Min(MAC_LIGHT_CF[cellLoc], MAC_TEMP_CF[cellLoc]);
     MAC_PROD_CF[cellLoc]  = Min(min_litTemp,MAC_WATER_CF[cellLoc])
         *MAC_NUT_CF[cellLoc]*MAC_SALT_CF[cellLoc];
     PHBIO_NPP[cellLoc] = HP_PHBIO_RCNPP[HAB[cellLoc]]*MAC_PROD_CF[cellLoc]*MAC_PH_BIOMAS[cellLoc]* Max( (1.0-MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc]), 0.0);
/* check for available P mass that will be taken up (MichMent kinetics in nutCF does not) */
     tmp = (PHBIO_NPP[cellLoc] > 0) ? 
         TP_SED_WT[cellLoc] / ( PHBIO_NPP[cellLoc] * HP_NPHBIO_IC_PC[HAB[cellLoc]]*CELL_SIZE*DT) :
         1.0;
     if (tmp < 1.0) PHBIO_NPP[cellLoc] *= tmp;   
    /* reduce the production proportionally if excess demand is found */
    /* can have high conc, but low mass of P avail, in presence of high plant biomass and high demand */
/* now add the P and OM associated with that C */
     phbio_npp_P[cellLoc] = PHBIO_NPP[cellLoc] * HP_PHBIO_IC_PC[HAB[cellLoc]];     /* habitat-specfic stoichiometry */
     phbio_npp_OM[cellLoc] = PHBIO_NPP[cellLoc] / HP_PHBIO_IC_CTOOM[HAB[cellLoc]]; /* habitat-specfic stoichiometry */

/* init ("target") ph/nph ratio and new transloc algorithm */
     MAC_PHtoNPH_Init = HP_PHBIO_MAX[HAB[cellLoc]] / HP_NPHBIO_MAX[HAB[cellLoc]] ;
     MAC_PHtoNPH = MAC_PH_BIOMAS[cellLoc] / MAC_NOPH_BIOMAS[cellLoc];

     PHBIO_TRANSLOC_POT[cellLoc]  = 0.0; /* (MAC_PHtoNPH<MAC_PHtoNPH_Init) ? (exp(HP_MAC_TRANSLOC_RC[HAB[cellLoc]] *(MAC_PHtoNPH_Init-MAC_PHtoNPH)) - 1.0) : (0.0) */ 

     PHBIO_TRANSLOC[cellLoc] =  ( (PHBIO_TRANSLOC_POT[cellLoc])*DT >NPHBIO_AVAIL[cellLoc] ) ? ( (NPHBIO_AVAIL[cellLoc])/DT ) : ( PHBIO_TRANSLOC_POT[cellLoc]);
/*  now remove the P and OM associated with that C */
     phbio_transl_P[cellLoc] = PHBIO_TRANSLOC[cellLoc] * mac_nph_PC[cellLoc];
     phbio_transl_OM[cellLoc] = PHBIO_TRANSLOC[cellLoc] / mac_nph_CtoOM[cellLoc];
     NPHBIO_MORT_POT[cellLoc] = NPHBIO_AVAIL[cellLoc]*HP_PHBIO_RCMORT[HAB[cellLoc]]*(1.0-MAC_PH_BIOMAS[cellLoc]/HP_PHBIO_MAX[HAB[cellLoc]]);
     NPHBIO_MORT[cellLoc] =  ( (PHBIO_TRANSLOC[cellLoc]+NPHBIO_MORT_POT[cellLoc])*DT>NPHBIO_AVAIL[cellLoc] ) ? ( (NPHBIO_AVAIL[cellLoc]-PHBIO_TRANSLOC[cellLoc]*DT)/DT ) : ( NPHBIO_MORT_POT[cellLoc]);
     PHBIO_MORT_POT[cellLoc] = HP_PHBIO_RCMORT[HAB[cellLoc]] *PHBIO_AVAIL[cellLoc] *(1.0-MAC_WATER_AVAIL_CF[cellLoc]);
/* now remove the P and OM associated with that C */
     nphbio_mort_P[cellLoc] = NPHBIO_MORT[cellLoc] * mac_nph_PC[cellLoc];
     nphbio_mort_OM[cellLoc] = NPHBIO_MORT[cellLoc] / mac_nph_CtoOM[cellLoc];

     NPHBIO_TRANSLOC_POT[cellLoc] = 0.0; /* (MAC_PHtoNPH>MAC_PHtoNPH_Init) ? (exp(HP_MAC_TRANSLOC_RC[HAB[cellLoc]] *(MAC_PHtoNPH-MAC_PHtoNPH_Init)) - 1.0) : (0.0) */ 
     NPHBIO_TRANSLOC[cellLoc] = ( (NPHBIO_TRANSLOC_POT[cellLoc])*DT >MAC_PH_BIOMAS[cellLoc] ) ? ( (MAC_PH_BIOMAS[cellLoc])/DT ) : ( NPHBIO_TRANSLOC_POT[cellLoc]);
/*  now remove the P and OM associated with that C */
     nphbio_transl_P[cellLoc] = NPHBIO_TRANSLOC[cellLoc] * mac_ph_PC[cellLoc];
     nphbio_transl_OM[cellLoc] = NPHBIO_TRANSLOC[cellLoc] / mac_ph_CtoOM[cellLoc];
     PHBIO_MORT[cellLoc] = ( (PHBIO_MORT_POT[cellLoc]+NPHBIO_TRANSLOC[cellLoc])*DT>PHBIO_AVAIL[cellLoc] ) ? ( (PHBIO_AVAIL[cellLoc]-NPHBIO_TRANSLOC[cellLoc]*DT)/DT ) : ( PHBIO_MORT_POT[cellLoc]);
/*  now remove the P associated with that C */
     phbio_mort_P[cellLoc] = PHBIO_MORT[cellLoc] * mac_ph_PC[cellLoc];
     phbio_mort_OM[cellLoc] = PHBIO_MORT[cellLoc] / mac_ph_CtoOM[cellLoc];


     FLOC_DECOMP_QUAL_CF[cellLoc] = Min(TP_SFWT_CONC_MG[cellLoc]/GP_DOM_DECOMP_POPT,1.0) ;
     FLOC_DECOMP_POT[cellLoc] = GP_DOM_RCDECOMP*FLOC[cellLoc]*DOM_TEMP_CF[cellLoc] *FLOC_DECOMP_QUAL_CF[cellLoc];
     FLOC_DECOMP[cellLoc] = ( (FLOC_DECOMP_POT[cellLoc])*DT>FLOC[cellLoc] ) ? ( (FLOC[cellLoc])/DT ) : ( FLOC_DECOMP_POT[cellLoc]);
     FLOC_DEPO_POT[cellLoc] = (  SURFACE_WAT[cellLoc] > GP_DetentZ ) ? ( FLOC[cellLoc]*0.05 ) : ( FLOC[cellLoc]/DT);
     FLOC_DEPO[cellLoc] = ( (FLOC_DEPO_POT[cellLoc]+FLOC_DECOMP[cellLoc])*DT>FLOC[cellLoc] ) ? ( (FLOC[cellLoc]-FLOC_DECOMP[cellLoc]*DT)/DT ) : ( FLOC_DEPO_POT[cellLoc]);
 
     HYD_MANNINGS_N[cellLoc]  = Max(-Abs((HP_MAC_MAXROUGH[HAB[cellLoc]]-HP_MAC_MINROUGH[HAB[cellLoc]])*  (pow(2.0,(1.0-SURFACE_WAT[cellLoc]/MAC_HEIGHT[cellLoc]))-1.0)) + HP_MAC_MAXROUGH[HAB[cellLoc]],HP_MAC_MINROUGH[HAB[cellLoc]]);
  } /* spatial loop end */
  } /* spatial loop end */
  usrErr("Done.");

} /* end of init_eqns */


void init_canals(int irun) 
{
   if (irun == 1) {
      usrErr("Initializing Water Management...");
      Canal_Network_Init(GP_DATUM_DISTANCE,SED_ELEV); /* WatMgmt.c - initialize the canal network topology and data */
      usrErr("Done Water Management.");
   }
   else {
      reinitCanals();
   }

}

void init_succession(void) 
{
   HabSwitch_Init( ); 

}

void reinitBIR(void) 
{
   usrErr0("Re-initializing Basin/Indicator-Region info...");
   BIRstats_reset(); 
   BIRbudg_reset(); 
   Cell_reset_avg(); 
   Cell_reset_hydper(); 
   usrErr("Done.");
}

void reinitCanals(void) 
{
   usrErr0("Re-initializing Canal depths & constituent masses...");
   CanalReInit();  
   usrErr("Done.");
}


void gen_output(int step, ViewParm *view)
{
    #define numOutputs 50000
    static int iw[numOutputs];
    int oIndex;
    ViewParm   *vp;

    static outVar_struct tgen[] = {
      { (float**)&TP_settlDays, "TP_settlDays", 'f' },
      { (float**)&FLOC, "FLOC", 'f' },
      { (float**)&FLOC_DECOMP, "FLOC_DECOMP", 'f' },
      { (float**)&FLOC_DECOMP_POT, "FLOC_DECOMP_POT", 'f' },
      { (float**)&FLOC_DECOMP_QUAL_CF, "FLOC_DECOMP_QUAL_CF", 'f' },
      { (float**)&FLOC_DEPO, "FLOC_DEPO", 'f' },
      { (float**)&FLOC_DEPO_POT, "FLOC_DEPO_POT", 'f' },
      { (float**)&FLOC_FR_ALGAE, "FLOC_FR_ALGAE", 'f' },
      { (float**)&FLOC_Z, "FLOC_Z", 'f' },
      { (float**)&FlocP_OMrep, "FlocP_OMrep", 'f' },
      { (float**)&soil_MOIST_CF, "soil_MOIST_CF", 'f' },
      { (float**)&TP_SED_MINER, "TP_SED_MINER", 'f' },
      { (float**)&TP_SFWT_MINER, "TP_SFWT_MINER", 'f' },
      { (float**)&AIR_TEMP, "AIR_TEMP", 'f' },
      { (unsigned char**)&HAB, "HAB", 'c' },
      { (float**)&SOLRAD274, "SOLRAD274", 'f' },
      { (float**)&SOLRADGRD, "SOLRADGRD", 'f' },
      { (float**)&H2O_TEMP, "H2O_TEMP", 'f' },
      { (float**)&HYD_DOM_ACTWAT_PRES, "HYD_DOM_ACTWAT_PRES", 'f' },
      { (float**)&HYD_DOM_ACTWAT_VOL, "HYD_DOM_ACTWAT_VOL", 'f' },
      { (float**)&HYD_ET, "HYD_ET", 'f' },
      { (float**)&HYD_EVAP_CALC, "HYD_EVAP_CALC", 'f' },
      { (float**)&HYD_MANNINGS_N, "HYD_MANNINGS_N", 'f' },
      { (float**)&HYD_SAT_POT_TRANS, "HYD_SAT_POT_TRANS", 'f' },
      { (float**)&HYD_SED_WAT_VOL, "HYD_SED_WAT_VOL", 'f' },
      { (float**)&HYD_TOT_POT_TRANSP, "HYD_TOT_POT_TRANSP", 'f' },
      { (float**)&HYD_TRANSP, "HYD_TRANSP", 'f' },
      { (float**)&HYD_UNSAT_POT_TRANS, "HYD_UNSAT_POT_TRANS", 'f' },
      { (float**)&HYD_WATER_AVAIL, "HYD_WATER_AVAIL", 'f' },
      { (float**)&HydTotHd, "HydTotHd", 'f' },
      { (float**)&HydRelDepPosNeg, "HydRelDepPosNeg", 'f' },
      { (float**)&LAI_eff, "LAI_eff", 'f' },
      { (float**)&MAC_WATER_AVAIL_CF, "MAC_WATER_AVAIL_CF", 'f' },
      { (float**)&SAT_TO_UNSAT_FL, "SAT_TO_UNSAT_FL", 'f' },
      { (float**)&SAT_VS_UNSAT, "SAT_VS_UNSAT", 'f' },
      { (float**)&SAT_WATER, "SAT_WATER", 'f' },
      { (float**)&SAT_WT_HEAD, "SAT_WT_HEAD", 'f' },
      { (float**)&SAT_WT_RECHG, "SAT_WT_RECHG", 'f' },
      { (float**)&SAT_WT_TRANSP, "SAT_WT_TRANSP", 'f' },
      { (float**)&SF_WT_EVAP, "SF_WT_EVAP", 'f' },
      { (float**)&SF_WT_FROM_RAIN, "SF_WT_FROM_RAIN", 'f' },
      { (float**)&SF_WT_INFILTRATION, "SF_WT_INFILTRATION", 'f' },
      { (float**)&SF_WT_POT_INF, "SF_WT_POT_INF", 'f' },
      { (float**)&SF_WT_TO_SAT_DOWNFLOW, "SF_WT_TO_SAT_DOWNFLOW", 'f' },
      { (float**)&SF_WT_VEL_mag, "SF_WT_VEL_mag", 'f' },
      { (float**)&BCmodel_sfwat, "BCmodel_sfwat", 'f' },
      { (float**)&SFWT_VOL, "SFWT_VOL", 'f' },
      { (float**)&SURFACE_WAT, "SURFACE_WAT", 'f' },
      { (float**)&UNSAT_AVAIL, "UNSAT_AVAIL", 'f' },
      { (float**)&UNSAT_CAP, "UNSAT_CAP", 'f' },
      { (float**)&UNSAT_DEPTH, "UNSAT_DEPTH", 'f' },
      { (float**)&UNSAT_HYD_COND_CF, "UNSAT_HYD_COND_CF", 'f' },
      { (float**)&UNSAT_MOIST_PRP, "UNSAT_MOIST_PRP", 'f' },
      { (float**)&UNSAT_PERC, "UNSAT_PERC", 'f' },
      { (float**)&UNSAT_TO_SAT_FL, "UNSAT_TO_SAT_FL", 'f' },
      { (float**)&UNSAT_TRANSP, "UNSAT_TRANSP", 'f' },
      { (float**)&UNSAT_WATER, "UNSAT_WATER", 'f' },
      { (float**)&UNSAT_WT_POT, "UNSAT_WT_POT", 'f' },
      { (float**)&ELEVATION, "ELEVATION", 'f' },
      { (float**)&HYD_RCCONDUCT, "HYD_RCCONDUCT", 'f' },
      { (unsigned char**)&ON_MAP, "ON_MAP", 'c' },
      { (float**)&SED_INACT_Z, "SED_INACT_Z", 'f' },
      { (float**)&MAC_HEIGHT, "MAC_HEIGHT", 'f' },
      { (float**)&MAC_LAI, "MAC_LAI", 'f' },
      { (float**)&MAC_LIGHT_CF, "MAC_LIGHT_CF", 'f' },
      { (float**)&MAC_MAX_BIO, "MAC_MAX_BIO", 'f' },
      { (float**)&MAC_NOPH_BIOMAS, "MAC_NOPH_BIOMAS", 'f' },
      { (float**)&mac_nph_PC_rep, "mac_nph_PC_rep", 'f' },
      { (float**)&MAC_NUT_CF, "MAC_NUT_CF", 'f' },
      { (float**)&MAC_PH_BIOMAS, "MAC_PH_BIOMAS", 'f' },
      { (float**)&mac_ph_PC_rep, "mac_ph_PC_rep", 'f' },
      { (float**)&MAC_PROD_CF, "MAC_PROD_CF", 'f' },
      { (float**)&MAC_REL_BIOM, "MAC_REL_BIOM", 'f' },
      { (float**)&MAC_SALT_CF, "MAC_SALT_CF", 'f' },
      { (float**)&MAC_TEMP_CF, "MAC_TEMP_CF", 'f' },
      { (float**)&MAC_TOT_BIOM, "MAC_TOT_BIOM", 'f' },
      { (float**)&MAC_WATER_CF, "MAC_WATER_CF", 'f' },
      { (float**)&NPHBIO_AVAIL, "NPHBIO_AVAIL", 'f' },
      { (float**)&NPHBIO_MORT, "NPHBIO_MORT", 'f' },
      { (float**)&NPHBIO_MORT_POT, "NPHBIO_MORT_POT", 'f' },
      { (float**)&NPHBIO_REFUGE, "NPHBIO_REFUGE", 'f' },
      { (float**)&NPHBIO_SAT, "NPHBIO_SAT", 'f' },
      { (float**)&NPHBIO_TRANSLOC, "NPHBIO_TRANSLOC", 'f' },
      { (float**)&NPHBIO_TRANSLOC_POT, "NPHBIO_TRANSLOC_POT", 'f' },
      { (float**)&PHBIO_AVAIL, "PHBIO_AVAIL", 'f' },
      { (float**)&PHBIO_MORT, "PHBIO_MORT", 'f' },
      { (float**)&PHBIO_MORT_POT, "PHBIO_MORT_POT", 'f' },
      { (float**)&PHBIO_NPP, "PHBIO_NPP", 'f' },
      { (float**)&PHBIO_REFUGE, "PHBIO_REFUGE", 'f' },
      { (float**)&PHBIO_SAT, "PHBIO_SAT", 'f' },
      { (float**)&PHBIO_TRANSLOC, "PHBIO_TRANSLOC", 'f' },
      { (float**)&PHBIO_TRANSLOC_POT, "PHBIO_TRANSLOC_POT", 'f' },
      { (float**)&TP_SEDWT_UPTAKE, "TP_SEDWT_UPTAKE", 'f' },
      { (float**)&ALG_INCID_LIGHT, "ALG_INCID_LIGHT", 'f' },
      { (float**)&ALG_LIGHT_CF, "ALG_LIGHT_CF", 'f' },
      { (float**)&ALG_LIGHT_EXTINCT, "ALG_LIGHT_EXTINCT", 'f' },
      { (float**)&ALG_REFUGE, "ALG_REFUGE", 'f' },
      { (float**)&ALG_SAT, "ALG_SAT", 'f' },
      { (float**)&ALG_TEMP_CF, "ALG_TEMP_CF", 'f' },
      { (float**)&ALG_TOT, "ALG_TOT", 'f' },
      { (float**)&ALG_WAT_CF, "ALG_WAT_CF", 'f' },
      { (float**)&C_ALG, "C_ALG", 'f' },
      { (float**)&C_ALG_AVAIL_MORT, "C_ALG_AVAIL_MORT", 'f' },
      { (float**)&C_ALG_GPP, "C_ALG_GPP", 'f' },
      { (float**)&C_ALG_MORT, "C_ALG_MORT", 'f' },
      { (float**)&C_ALG_MORT_POT, "C_ALG_MORT_POT", 'f' },
      { (float**)&C_ALG_NPP, "C_ALG_NPP", 'f' },
      { (float**)&C_ALG_NUT_CF, "C_ALG_NUT_CF", 'f' },
      { (float**)&C_ALG_PROD_CF, "C_ALG_PROD_CF", 'f' },
      { (float**)&C_ALG_RESP, "C_ALG_RESP", 'f' },
      { (float**)&C_ALG_RESP_POT, "C_ALG_RESP_POT", 'f' },
      { (float**)&NC_ALG, "NC_ALG", 'f' },
      { (float**)&NC_ALG_AVAIL_MORT, "NC_ALG_AVAIL_MORT", 'f' },
      { (float**)&NC_ALG_GPP, "NC_ALG_GPP", 'f' },
      { (float**)&NC_ALG_MORT, "NC_ALG_MORT", 'f' },
      { (float**)&NC_ALG_MORT_POT, "NC_ALG_MORT_POT", 'f' },
      { (float**)&NC_ALG_NPP, "NC_ALG_NPP", 'f' },
      { (float**)&NC_ALG_NUT_CF, "NC_ALG_NUT_CF", 'f' },
      { (float**)&NC_ALG_PROD_CF, "NC_ALG_PROD_CF", 'f' },
      { (float**)&NC_ALG_RESP, "NC_ALG_RESP", 'f' },
      { (float**)&NC_ALG_RESP_POT, "NC_ALG_RESP_POT", 'f' },
      { (float**)&TP_SFWT_UPTAK, "TP_SFWT_UPTAK", 'f' },
      { (float**)&TP_Act_to_TotRep, "TP_Act_to_TotRep", 'f' },
      { (float**)&TP_DNFLOW, "TP_DNFLOW", 'f' },
      { (float**)&TP_DNFLOW_POT, "TP_DNFLOW_POT", 'f' },
      { (float**)&TP_Atm_Depos, "TP_Atm_Depos", 'f' },
      { (float**)&TP_K, "TP_K", 'f' },
      { (double**)&TP_SED_CONC, "TP_SED_CONC", 'd' },
      { (double**)&TP_SED_WT, "TP_SED_WT", 'd' },
      { (double**)&TP_SED_WT_AZ, "TP_SED_WT_AZ", 'd' },
      { (double**)&TP_SEDWT_CONCACT, "TP_SEDWT_CONCACT", 'd' },
      { (float**)&TP_SEDWT_CONCACTMG, "TP_SEDWT_CONCACTMG", 'f' },
      { (float**)&TP_settl, "TP_settl", 'f' },
      { (double**)&TP_SF_WT, "TP_SF_WT", 'd' },
      { (double**)&TP_SFWT_CONC, "TP_SFWT_CONC", 'd' },
      { (float**)&TP_SFWT_CONC_MG, "TP_SFWT_CONC_MG", 'f' },
      { (double**)&TP_SORB, "TP_SORB", 'd' },
      { (float**)&TP_SORB_POT, "TP_SORB_POT", 'f' },
      { (double**)&TP_SORBCONC, "TP_SORBCONC", 'd' },
      { (float**)&TP_SORBCONC_rep, "TP_SORBCONC_rep", 'f' },
      { (float**)&TP_SORBTION, "TP_SORBTION", 'f' },
      { (float**)&TP_UPFLOW, "TP_UPFLOW", 'f' },
      { (float**)&TP_UPFLOW_POT, "TP_UPFLOW_POT", 'f' },
      { (double**)&SAL_SED_WT, "SAL_SED_WT", 'd' },
      { (float**)&SAL_SF_WT, "SAL_SF_WT", 'f' },
      { (float**)&SALT_SED_TO_SF_FLOW, "SALT_SED_TO_SF_FLOW", 'f' },
      { (double**)&SALT_SED_WT, "SALT_SED_WT", 'd' },
      { (float**)&SALT_SFWAT_DOWNFL, "SALT_SFWAT_DOWNFL", 'f' },
      { (float**)&SALT_SFWAT_DOWNFL_POT, "SALT_SFWAT_DOWNFL_POT", 'f' },
      { (double**)&SALT_SURF_WT, "SALT_SURF_WT", 'd' },
      { (float**)&SALT_Atm_Depos, "SALT_Atm_Depos", 'f' },
      { (double**)&DEPOS_ORG_MAT, "DEPOS_ORG_MAT", 'd' },
      { (float**)&DOM_DECOMP, "DOM_DECOMP", 'f' },
      { (float**)&DOM_DECOMP_POT, "DOM_DECOMP_POT", 'f' },
      { (float**)&DOM_FR_FLOC, "DOM_FR_FLOC", 'f' },
      { (double**)&DOM_P_OM, "DOM_P_OM", 'd' },
      { (float**)&DOM_QUALITY_CF, "DOM_QUALITY_CF", 'f' },
      { (float**)&DOM_SED_AEROB_Z, "DOM_SED_AEROB_Z", 'f' },
      { (float**)&DOM_SED_ANAEROB_Z, "DOM_SED_ANAEROB_Z", 'f' },
      { (float**)&DOM_TEMP_CF, "DOM_TEMP_CF", 'f' },
      { (float**)&DOM_Z, "DOM_Z", 'f' },
      { (double**)&DOP, "DOP", 'd' },
      { (float**)&DOP_DECOMP, "DOP_DECOMP", 'f' },
      { (float**)&P_SUM_CELL, "P_SUM_CELL", 'f' },
      { (float**)&SED_ELEV, "SED_ELEV", 'f' },
      { (float**)&TPtoSOIL_rep, "TPtoSOIL_rep", 'f' },
      { (float**)&Floc_fr_phBioAvg, "Floc_fr_phBioAvg", 'f' },
      { (float**)&TPSfMinAvg, "TPSfMinAvg", 'f' },
      { (float**)&ETAvg, "ETAvg", 'f' },
      { (float**)&EvapAvg, "EvapAvg", 'f' },
      { (float**)&HydPerAnn, "HydPerAnn", 'f' },
      { (float**)&LAI_effAvg, "LAI_effAvg", 'f' },
      { (float**)&Manning_nAvg, "Manning_nAvg", 'f' },
      { (float**)&RainAvg, "RainAvg", 'f' },
      { (float**)&SfWatAvg, "SfWatAvg", 'f' },
      { (float**)&TotHeadAvg, "TotHeadAvg", 'f' },
      { (float**)&HydRelDepPosNegAvg, "HydRelDepPosNegAvg", 'f' },
      { (float**)&TranspAvg, "TranspAvg", 'f' },
      { (float**)&UnsatMoistAvg, "UnsatMoistAvg", 'f' },
      { (float**)&UnsatZavg, "UnsatZavg", 'f' },
      { (float**)&SF_WT_VEL_magAvg, "SF_WT_VEL_magAvg", 'f' },
      { (float**)&BCmodel_sfwatAvg, "BCmodel_sfwatAvg", 'f' },
      { (float**)&mac_nph_PCAvg, "mac_nph_PCAvg", 'f' },
      { (float**)&Mac_nphBioAvg, "Mac_nphBioAvg", 'f' },
      { (float**)&Mac_nphMortAvg, "Mac_nphMortAvg", 'f' },
      { (float**)&Mac_nppAvg, "Mac_nppAvg", 'f' },
      { (float**)&mac_ph_PCAvg, "mac_ph_PCAvg", 'f' },
      { (float**)&Mac_phBioAvg, "Mac_phBioAvg", 'f' },
      { (float**)&Mac_phMortAvg, "Mac_phMortAvg", 'f' },
      { (float**)&Mac_totBioAvg, "Mac_totBioAvg", 'f' },
      { (float**)&MacNutCfAvg, "MacNutCfAvg", 'f' },
      { (float**)&MacWatCfAvg, "MacWatCfAvg", 'f' },
      { (float**)&TPSedUptAvg, "TPSedUptAvg", 'f' },
      { (float**)&C_Peri_mortAvg, "C_Peri_mortAvg", 'f' },
      { (float**)&C_Peri_nppAvg, "C_Peri_nppAvg", 'f' },
      { (float**)&C_Peri_PCAvg, "C_Peri_PCAvg", 'f' },
      { (float**)&C_PeriAvg, "C_PeriAvg", 'f' },
      { (float**)&C_PeriNutCFAvg, "C_PeriNutCFAvg", 'f' },
      { (float**)&C_PeriRespAvg, "C_PeriRespAvg", 'f' },
      { (float**)&NC_Peri_mortAvg, "NC_Peri_mortAvg", 'f' },
      { (float**)&NC_Peri_nppAvg, "NC_Peri_nppAvg", 'f' },
      { (float**)&NC_Peri_PCAvg, "NC_Peri_PCAvg", 'f' },
      { (float**)&NC_PeriAvg, "NC_PeriAvg", 'f' },
      { (float**)&NC_PeriNutCFAvg, "NC_PeriNutCFAvg", 'f' },
      { (float**)&NC_PeriRespAvg, "NC_PeriRespAvg", 'f' },
      { (float**)&PeriAvg, "PeriAvg", 'f' },
      { (float**)&PeriLiteCFAvg, "PeriLiteCFAvg", 'f' },
      { (float**)&TPSfUptAvg, "TPSfUptAvg", 'f' },
      { (float**)&TP_settlAvg, "TP_settlAvg", 'f' },
      { (float**)&TPSedWatAvg, "TPSedWatAvg", 'f' },
      { (float**)&TPSfWatAvg, "TPSfWatAvg", 'f' },
      { (float**)&SaltSedAvg, "SaltSedAvg", 'f' },
      { (float**)&SaltSfAvg, "SaltSfAvg", 'f' },
      { (float**)&SedElevAvg, "SedElevAvg", 'f' },
      { (float**)&TPSedMinAvg, "TPSedMinAvg", 'f' },
      { (float**)&TPSorbAvg, "TPSorbAvg", 'f' },
      { (float**)&TPtoSOILAvg, "TPtoSOILAvg", 'f' },
      { (float**)&TPtoVOLAvg, "TPtoVOLAvg", 'f' },








      { NULL, NULL, '\0' },
  };
  
    outVar_struct *ptable;
    int  i;

    for (i = 0, ptable = tgen; ptable->pfvar != NULL; ptable++, i++)
    {
  vp = view + i;

/* TODO: develop flexible output of calendar-based and julian-based outsteps (jan 11, 2005) */
if (vp->step > 0)
                
            if (strcmp(ptable->pname,"HydPerAnn")!=0) { /* i.e., not the HydPerAnn variable */
                
                if  (step % vp->step == 0 && (vp->step != CalMonOut) ) {   /* standard julian-day outstep interval variables (note: the != CalMonOut needed for step=0) */ 
                    oIndex = iw[i]++;
                    write_output(oIndex, vp, *(ptable->pfvar), ptable->pname, ptable->ctype, step);
                }
                else if ( (avgPrint) && (vp->step == CalMonOut) ) { /* variables output at the special 1-calendar-month outstep interval */  
                    oIndex = iw[i]++;
                    write_output(oIndex, vp, *(ptable->pfvar), ptable->pname, ptable->ctype, step);
                }
            }
        
            else
                if (FMOD(DAYJUL, 273.0) ==0) { /* hydroperiod is printed at a special-case time (approximately Oct 1 every year) */
                    oIndex = iw[i]++;
                    write_output(oIndex, vp, *(ptable->pfvar), ptable->pname, ptable->ctype, step);
                }
        
     
    }

        /* after printing, zero the arrays holding averages or hydroperiods (v2.2.1note using avgPrint var)*/
    if (avgPrint) {
        Cell_reset_avg();
    }

    if (FMOD(DAYJUL, 273.0) ==0) {
        Cell_reset_hydper();
    }


} /* end of gen_output routine */


void ReadHabParms ( char* s_parm_name, int s_parm_relval )
{

 /* Periphyton/Algae */
  HP_ALG_MAX = get_hab_parm(s_parm_name, s_parm_relval, "HP_ALG_MAX"); 

/* Soil */
  HP_DOM_MAXDEPTH = get_hab_parm(s_parm_name, s_parm_relval,"HP_DOM_MAXDEPTH");  
  HP_DOM_AEROBTHIN = get_hab_parm(s_parm_name, s_parm_relval,"HP_DOM_AEROBTHIN"); 

/* Phosphorus */
  HP_TP_CONC_GRAD = get_hab_parm(s_parm_name, s_parm_relval,"HP_TP_CONC_GRAD"); 

/* Salt/tracer */
  HP_SALT_ICSEDWAT = get_hab_parm(s_parm_name, s_parm_relval,"HP_SALT_ICSEDWAT"); 
  HP_SALT_ICSFWAT = get_hab_parm(s_parm_name, s_parm_relval,"HP_SALT_ICSFWAT"); 
        
/* Macrophytes */
  HP_PHBIO_MAX = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_MAX"); 
  HP_NPHBIO_MAX = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_MAX"); 
  HP_MAC_MAXHT = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXHT"); 
  HP_NPHBIO_ROOTDEPTH = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_ROOTDEPTH"); 
  HP_MAC_MAXROUGH = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXROUGH"); 
  HP_MAC_MINROUGH = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MINROUGH"); 
  HP_MAC_MAXLAI = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXLAI"); 
  HP_MAC_MAXCANOPCOND = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXCANOPCOND"); /* unused in ELMv2.2, 2.3 */ 
  HP_MAC_CANOPDECOUP = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_CANOPDECOUP"); /* unused in ELMv2.2, 2.3 */ 
  HP_MAC_TEMPOPT = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_TEMPOPT"); 
  HP_MAC_LIGHTSAT = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_LIGHTSAT"); 
  HP_MAC_KSP = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_KSP"); 
  HP_PHBIO_RCNPP = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_RCNPP"); 
  HP_PHBIO_RCMORT = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_RCMORT"); 
  HP_MAC_WAT_TOLER = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_WAT_TOLER"); 
  HP_MAC_SALIN_THRESH = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_SALIN_THRESH"); /* unused in ELMv2.2, 2.3 */ 
  HP_PHBIO_IC_CTOOM = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_IC_CTOOM"); 
  HP_NPHBIO_IC_CTOOM = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_IC_CTOOM"); 
  HP_PHBIO_IC_PC = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_IC_PC"); 
  HP_NPHBIO_IC_PC = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_IC_PC"); 
  HP_MAC_TRANSLOC_RC = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_TRANSLOC_RC"); 

/* Hydrology */
  HP_HYD_RCINFILT = get_hab_parm(s_parm_name, s_parm_relval,"HP_HYD_RCINFILT"); 
  HP_HYD_SPEC_YIELD = get_hab_parm(s_parm_name, s_parm_relval,"HP_HYD_SPEC_YIELD");  
  HP_HYD_POROSITY = get_hab_parm(s_parm_name, s_parm_relval,"HP_HYD_POROSITY"); 

/* Floc */
  HP_FLOC_IC = get_hab_parm(s_parm_name, s_parm_relval,"HP_FLOC_IC"); 
  HP_FLOC_IC_CTOOM = get_hab_parm(s_parm_name, s_parm_relval,"HP_FLOC_IC_CTOOM"); 
  HP_FLOC_IC_PC = get_hab_parm(s_parm_name, s_parm_relval,"HP_FLOC_IC_PC"); 

/* Habitat succession */
  HP_SfDepthLo = get_hab_parm(s_parm_name, s_parm_relval,"HP_SfDepthLo"); 
  HP_SfDepthHi = get_hab_parm(s_parm_name, s_parm_relval,"HP_SfDepthHi"); 
  HP_SfDepthInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_SfDepthInt"); 
  HP_PhosLo = get_hab_parm(s_parm_name, s_parm_relval,"HP_PhosLo"); 
  HP_PhosHi = get_hab_parm(s_parm_name, s_parm_relval,"HP_PhosHi"); 
  HP_PhosInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_PhosInt"); 
  HP_FireInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_FireInt"); 
  HP_SalinLo = get_hab_parm(s_parm_name, s_parm_relval,"HP_SalinLo"); 
  HP_SalinHi = get_hab_parm(s_parm_name, s_parm_relval,"HP_SalinHi"); 
  HP_SalinInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_SalinInt"); 
  
}


void ReadGlobalParms( char* s_parm_name, int s_parm_relval )
 {
     
/* Geography, hydrology */
    GP_SOLOMEGA = get_global_parm(s_parm_name, s_parm_relval,"GP_SOLOMEGA"); 
    GP_ALTIT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALTIT"); 
    GP_LATDEG = get_global_parm(s_parm_name, s_parm_relval,"GP_LATDEG"); 
    GP_mannDepthPow = get_global_parm(s_parm_name, s_parm_relval,"GP_mannDepthPow"); 
    GP_mannHeadPow = get_global_parm(s_parm_name, s_parm_relval,"GP_mannHeadPow"); 
    GP_calibGWat = get_global_parm(s_parm_name, s_parm_relval,"GP_calibGWat"); 
    GP_IDW_pow = get_global_parm(s_parm_name, s_parm_relval,"GP_IDW_pow"); 
    GP_calibET = get_global_parm(s_parm_name, s_parm_relval,"GP_calibET"); 
    GP_DATUM_DISTANCE = get_global_parm(s_parm_name, s_parm_relval,"GP_DATUM_DISTANCE"); 
    GP_HYD_IC_SFWAT_ADD = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_IC_SFWAT_ADD"); 
    GP_HYD_IC_UNSAT_ADD = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_IC_UNSAT_ADD"); 
    GP_HYD_RCRECHG = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_RCRECHG"); 
    GP_HYD_ICUNSATMOIST = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_ICUNSATMOIST"); 
    GP_DetentZ = get_global_parm(s_parm_name, s_parm_relval,"GP_DetentZ"); 
    GP_WQualMonitZ = get_global_parm(s_parm_name, s_parm_relval,"GP_WQualMonitZ"); 
    GP_MinCheck = get_global_parm(s_parm_name, s_parm_relval,"GP_MinCheck"); 
    GP_SLRise = get_global_parm(s_parm_name, s_parm_relval,"GP_SLRise"); 
    GP_dispLenRef = get_global_parm(s_parm_name, s_parm_relval,"GP_dispLenRef"); 
    GP_dispParm = get_global_parm(s_parm_name, s_parm_relval,"GP_dispParm"); 
 
 /* Periphyton/Algae */
    GP_ALG_IC_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_IC_MULT"); 
    GP_alg_uptake_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_uptake_coef"); 
    GP_ALG_SHADE_FACTOR = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_SHADE_FACTOR"); 
    GP_algMortDepth = get_global_parm(s_parm_name, s_parm_relval,"GP_algMortDepth"); 
    GP_ALG_RC_MORT_DRY = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_MORT_DRY"); 
    GP_ALG_RC_MORT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_MORT"); 
    GP_ALG_RC_PROD = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_PROD"); 
    GP_ALG_RC_RESP = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_RESP"); 
    GP_alg_R_accel = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_R_accel"); 
    GP_AlgComp = get_global_parm(s_parm_name, s_parm_relval,"GP_AlgComp"); 
    GP_ALG_REF_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_REF_MULT"); 
    GP_NC_ALG_KS_P = get_global_parm(s_parm_name, s_parm_relval,"GP_NC_ALG_KS_P"); 
    GP_alg_alkP_min = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_alkP_min"); 
    GP_C_ALG_KS_P = get_global_parm(s_parm_name, s_parm_relval,"GP_C_ALG_KS_P"); 
    GP_ALG_TEMP_OPT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_TEMP_OPT"); 
    GP_C_ALG_threshTP = get_global_parm(s_parm_name, s_parm_relval,"GP_C_ALG_threshTP"); 
    GP_ALG_C_TO_OM = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_C_TO_OM"); 
    GP_alg_light_ext_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_light_ext_coef"); 
    GP_ALG_LIGHT_SAT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_LIGHT_SAT"); 
    GP_ALG_PC = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_PC"); 

/* Soil */
    GP_DOM_RCDECOMP = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_RCDECOMP"); 
    GP_DOM_DECOMPRED = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_DECOMPRED"); 
    GP_calibDecomp = get_global_parm(s_parm_name, s_parm_relval,"GP_calibDecomp"); 
    GP_DOM_decomp_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_decomp_coef"); 
    GP_DOM_DECOMP_POPT = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_DECOMP_POPT"); 
    GP_DOM_DECOMP_TOPT = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_DECOMP_TOPT"); 
    GP_sorbToTP = get_global_parm(s_parm_name, s_parm_relval,"GP_sorbToTP"); 
    GP_IC_ELEV_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_ELEV_MULT"); 
    GP_IC_BATHY_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_BATHY_MULT"); 
    GP_IC_DOM_BD_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_DOM_BD_MULT"); 
    GP_IC_BulkD_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_BulkD_MULT"); 
    GP_IC_TPtoSOIL_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_TPtoSOIL_MULT"); 
 
/* Macrophytes */
    GP_MAC_IC_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_MAC_IC_MULT"); 
    GP_MAC_REFUG_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_MAC_REFUG_MULT"); 
    GP_mac_uptake_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_mac_uptake_coef"); 
    GP_mann_height_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_mann_height_coef"); 

/* Floc */
    GP_Floc_BD = get_global_parm(s_parm_name, s_parm_relval,"GP_Floc_BD"); 
    GP_FlocMax = get_global_parm(s_parm_name, s_parm_relval,"GP_FlocMax"); 
    GP_TP_P_OM = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_P_OM"); 
    GP_Floc_rcSoil = get_global_parm(s_parm_name, s_parm_relval,"GP_Floc_rcSoil"); 
                
/* Phosphorus */
    GP_TP_DIFFCOEF = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_DIFFCOEF"); 
    GP_TP_K_INTER = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_K_INTER"); 
    GP_TP_K_SLOPE = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_K_SLOPE"); 
    GP_WQMthresh = get_global_parm(s_parm_name, s_parm_relval,"GP_WQMthresh"); 
    GP_PO4toTP = get_global_parm(s_parm_name, s_parm_relval,"GP_PO4toTP"); 
    GP_TP_IN_RAIN = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_IN_RAIN"); 
    GP_CL_IN_RAIN = get_global_parm(s_parm_name, s_parm_relval,"GP_CL_IN_RAIN"); 
    GP_PO4toTPint = get_global_parm(s_parm_name, s_parm_relval,"GP_PO4toTPint"); 
    GP_TP_ICSFWAT = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_ICSFWAT"); 
    GP_TP_ICSEDWAT = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_ICSEDWAT"); 
    GP_TPpart_thresh = get_global_parm(s_parm_name, s_parm_relval,"GP_TPpart_thresh"); 
    GP_TP_DIFFDEPTH = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_DIFFDEPTH"); 
    GP_settlVel = get_global_parm(s_parm_name, s_parm_relval,"GP_settlVel"); 
    
}

void ReadModExperimParms( char* s_parm_name, int s_parm_relval )
 {
    if (isModExperim) {
                 sprintf(msgStr,"####### Reading special model experiment parameters.  \nYou should NOT be using gridIO spatial time series data for rain, pET, and SFWMM-output depth #######..."); 
            WriteMsg(msgStr,1); 
            usrErr0(msgStr);
            
        GP_BC_SfWat_add = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_SfWat_add"); 
        GP_BC_SfWat_addHi = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_SfWat_addHi"); 
        GP_BC_SatWat_add = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_SatWat_add"); 
        GP_BC_SfWat_targ = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_SfWat_targ"); 
        GP_BC_InputRow = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_InputRow"); 
        GP_BC_Pconc = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_Pconc"); /* ug/L */
        GP_BC_Sconc = get_modexperim_parm(s_parm_name, s_parm_relval,"GP_BC_Sconc"); /* g/L */
        
        GP_BC_Pconc /= conv_kgTOmg; /* convert from ug/L to g/L (same as div by the kg->mg conversion) */

            sprintf(msgStr,"done."); 
            WriteMsg(msgStr,1); 
            usrErr(msgStr);
    }
    else {
        GP_BC_SfWat_add = 0.0; /* value of 0.0 has no effect on model results (adds nothing to boundary depth/stages) */
        GP_BC_SfWat_addHi = 0.0; /* value of 0.0 has no effect on model results (adds nothing to boundary depth/stages) */
        GP_BC_SatWat_add = 0.0; /* value of 0.0 has no effect on model results (adds nothing to boundary depth/stages) */
        GP_BC_SfWat_targ = 0.0; /* value of 0.0 has no effect on model results (only used in Oct06 specialized model experiment - code doesn't use this parm otherwise) */
        GP_BC_InputRow = 0; /* value of 0 has no effect on model results (only used in Oct06 specialized model experiment - code doesn't use this parm otherwise) */
        GP_BC_Pconc = 0.0;
        GP_BC_Sconc = 0.0;
    }
 } 

/* v2.6 PTSL */
void PtInterp_read() 
{
                /* rainfall and pET point data for models that do not use gridIO inputs, but use interpolated point time series */
            /* final arg in PTSL_ReadLists is the column number of of the data (not incl. the 3 date columns) */
        usrErr0 ("\n*** WARNING: You are not using gridIO spatial time series data for rain, pET, and SFWMM-output depth.  \nReading rainfall and pET point time series data... ");
                pTSList2 = PTSL_Init(50, 2);
                PTSL_ReadLists(pTSList2,"rain",TIndex2++,&Timestep2,&NPtTS2, 1); /* rainfall in data-col 1 (file's col 4) */
        
        pTSList3 = PTSL_Init(50, 2);
        PTSL_ReadLists(pTSList3,"pET",TIndex3++,&Timestep3,&NPtTS3, 1); /* pET data in data-col 1 (file's col 4) */

        usrErr("done.\n");

        if (99 == 66) { /* never-get-here: temporary v2.6 PTSL - Below (old code frags extracted from v2.1) are left here as template for other data (or re-implementing the ELM v2.1 pET calculations if data available) */
        pTSList4 = PTSL_Init(50, 2);
        PTSL_ReadLists(pTSList4,"pET",TIndex4++,&Timestep4,&NPtTS4, 2);  /* dew point temperature data can be in data-col 2 */
        pTSList5 = PTSL_Init(50, 2);
        PTSL_ReadLists(pTSList5,"pET",TIndex5++,&Timestep5,&NPtTS5, 3);  /* wind speed data can be in data-col 3 */
        pTSList6 = PTSL_Init(50, 2);
        PTSL_ReadLists(pTSList6,"pET",TIndex6++,&Timestep6,&NPtTS6, 4);  /* cloud cover data can be in data-col 4 */
        }
}


void get_map_dims() 
{
  read_map_dims("Elevation");
}

/* v2.8 a large number of variables are incorrect/inconsistent in initializing for double or float - corrected below, w/o comments */
void alloc_memory() 
{
  usrErr0("Allocating Memory...");  /* console message */
  
  ON_MAP = (unsigned char *) nalloc(sizeof(unsigned char)*(s0+2)*(s1+2),"ON_MAP");
  init_pvar(ON_MAP,NULL,'c',0.0);

  BCondFlow = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"BCondFlow");
  init_pvar(BCondFlow,NULL,'d',0.0);
  HAB = (unsigned char *) nalloc(sizeof(unsigned char)*(s0+2)*(s1+2),"HAB");
  init_pvar(HAB,NULL,'c',0.0);
  basn = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"basn"); /* Basin/Indicator-Region map */
  init_pvar(basn,NULL,'d',0.0);

  ALG_INCID_LIGHT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_INCID_LIGHT");
  init_pvar(ALG_INCID_LIGHT,NULL,'f',0.0);
  ALG_LIGHT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_LIGHT_CF");
  init_pvar(ALG_LIGHT_CF,NULL,'f',0.0);
  ALG_LIGHT_EXTINCT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_LIGHT_EXTINCT");
  init_pvar(ALG_LIGHT_EXTINCT,NULL,'f',0.0);
  ALG_WAT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_WAT_CF");
  init_pvar(ALG_WAT_CF,NULL,'f',0.0);
  ALG_TEMP_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_TEMP_CF");
  init_pvar(ALG_TEMP_CF,NULL,'f',0.0);
  ALG_REFUGE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_REFUGE");
  init_pvar(ALG_REFUGE,NULL,'f',0.0);
  ALG_SAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_SAT");
  init_pvar(ALG_SAT,NULL,'f',0.0);
  ALG_TOT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_TOT");
  init_pvar(ALG_TOT,NULL,'f',0.0);

  NC_ALG_AVAIL_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_AVAIL_MORT");
  init_pvar(NC_ALG_AVAIL_MORT,NULL,'f',0.0);
  NC_ALG_GPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_GPP");
  init_pvar(NC_ALG_GPP,NULL,'f',0.0);
  NC_ALG_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_MORT");
  init_pvar(NC_ALG_MORT,NULL,'f',0.0);
  NC_ALG_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_MORT_POT");
  init_pvar(NC_ALG_MORT_POT,NULL,'f',0.0);
  NC_ALG_NPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_NPP");
  init_pvar(NC_ALG_NPP,NULL,'f',0.0);
  NC_ALG_NUT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_NUT_CF");
  init_pvar(NC_ALG_NUT_CF,NULL,'f',0.0);
  NC_ALG_PROD_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_PROD_CF");
  init_pvar(NC_ALG_PROD_CF,NULL,'f',0.0);
  NC_ALG_RESP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_RESP");
  init_pvar(NC_ALG_RESP,NULL,'f',0.0);
  NC_ALG_RESP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_RESP_POT");
  init_pvar(NC_ALG_RESP_POT,NULL,'f',0.0);
  NC_ALG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG");
  init_pvar(NC_ALG,NULL,'f',0.0);
  C_ALG_AVAIL_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_AVAIL_MORT");
  init_pvar(C_ALG_AVAIL_MORT,NULL,'f',0.0);
  C_ALG_GPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_GPP");
  init_pvar(C_ALG_GPP,NULL,'f',0.0);
  C_ALG_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_MORT");
  init_pvar(C_ALG_MORT,NULL,'f',0.0);
  C_ALG_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_MORT_POT");
  init_pvar(C_ALG_MORT_POT,NULL,'f',0.0);
  C_ALG_NPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_NPP");
  init_pvar(C_ALG_NPP,NULL,'f',0.0);
  C_ALG_NUT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_NUT_CF");
  init_pvar(C_ALG_NUT_CF,NULL,'f',0.0);
  C_ALG_PROD_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_PROD_CF");
  init_pvar(C_ALG_PROD_CF,NULL,'f',0.0);
  C_ALG_RESP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_RESP");
  init_pvar(C_ALG_RESP,NULL,'f',0.0);
  C_ALG_RESP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_RESP_POT");
  init_pvar(C_ALG_RESP_POT,NULL,'f',0.0);
  C_ALG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG");
  init_pvar(C_ALG,NULL,'f',0.0);
  NC_ALG_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_P");
  init_pvar(NC_ALG_P,NULL,'f',0.0);
  C_ALG_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_P");
  init_pvar(C_ALG_P,NULL,'f',0.0);
  NC_ALG_GPP_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_GPP_P");
  init_pvar(NC_ALG_GPP_P,NULL,'f',0.0);
  C_ALG_GPP_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_GPP_P");
  init_pvar(C_ALG_GPP_P,NULL,'f',0.0);
  NC_ALG_MORT_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_MORT_P");
  init_pvar(NC_ALG_MORT_P,NULL,'f',0.0);
  C_ALG_MORT_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_MORT_P");
  init_pvar(C_ALG_MORT_P,NULL,'f',0.0);
  NC_ALG_PCrep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_PCrep");
  init_pvar(NC_ALG_PCrep,NULL,'f',0.0);
  C_ALG_PCrep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_PCrep");
  init_pvar(C_ALG_PCrep,NULL,'f',0.0);
  NC_ALG_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"NC_ALG_PC");
  init_pvar(NC_ALG_PC,NULL,'l',0.0);
  C_ALG_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"C_ALG_PC");
  init_pvar(C_ALG_PC,NULL,'l',0.0);

  DEPOS_ORG_MAT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DEPOS_ORG_MAT");
  init_pvar(DEPOS_ORG_MAT,NULL,'l',0.0);
        
  DOM_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_Z");
  init_pvar(DOM_Z,NULL,'f',0.0);
  DOM_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_DECOMP");
  init_pvar(DOM_DECOMP,NULL,'f',0.0);
  DOM_DECOMP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_DECOMP_POT");
  init_pvar(DOM_DECOMP_POT,NULL,'f',0.0);
  DOM_fr_nphBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_fr_nphBio");
  init_pvar(DOM_fr_nphBio,NULL,'f',0.0);
  
  Floc_fr_phBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"Floc_fr_phBio");
  init_pvar(Floc_fr_phBio,NULL,'f',0.0);
  DOM_FR_FLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_FR_FLOC");
  init_pvar(DOM_FR_FLOC,NULL,'f',0.0);
  soil_MOIST_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"soil_MOIST_CF");
  init_pvar(soil_MOIST_CF,NULL,'f',0.0);
  DOM_QUALITY_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_QUALITY_CF");
  init_pvar(DOM_QUALITY_CF,NULL,'f',0.0);
  DOM_SED_AEROB_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_SED_AEROB_Z");
  init_pvar(DOM_SED_AEROB_Z,NULL,'f',0.0);
  DOM_SED_ANAEROB_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_SED_ANAEROB_Z");
  init_pvar(DOM_SED_ANAEROB_Z,NULL,'f',0.0);
  DOM_TEMP_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_TEMP_CF");
  init_pvar(DOM_TEMP_CF,NULL,'f',0.0);
  ELEVATION = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ELEVATION");
  init_pvar(ELEVATION,NULL,'f',0.0);
  Bathymetry = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"Bathymetry");
  init_pvar(Bathymetry,NULL,'f',0.0);
  SED_ELEV = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SED_ELEV");
  init_pvar(SED_ELEV,NULL,'f',0.0);
  SED_INACT_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SED_INACT_Z");
  init_pvar(SED_INACT_Z,NULL,'f',0.0);
  DOP_FLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOP_FLOC");
  init_pvar(DOP_FLOC,NULL,'f',0.0);
  DOP_nphBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOP_nphBio");
  init_pvar(DOP_nphBio,NULL,'f',0.0);

  DOM_P_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DOM_P_OM");
  init_pvar(DOM_P_OM,NULL,'l',0.0);

  TPtoSOIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TPtoSOIL");
  init_pvar(TPtoSOIL,NULL,'f',0.0);
  TPtoSOIL_rep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TPtoSOIL_rep");
  init_pvar(TPtoSOIL_rep,NULL,'f',0.0);
  TPtoVOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TPtoVOL");
  init_pvar(TPtoVOL,NULL,'f',0.0);
  TPtoVOL_rep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TPtoVOL_rep");
  init_pvar(TPtoVOL_rep,NULL,'f',0.0);
  BulkD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"BulkD");
  init_pvar(BulkD,NULL,'f',0.0);
  DOM_BD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_BD");
  init_pvar(DOM_BD,NULL,'f',0.0);
  Inorg_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"Inorg_Z");
  init_pvar(Inorg_Z,NULL,'f',0.0);
  DIM = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DIM");
  init_pvar(DIM,NULL,'f',0.0);
  DOP_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOP_DECOMP");
  init_pvar(DOP_DECOMP,NULL,'f',0.0);

  DOP = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DOP");
  init_pvar(DOP,NULL,'l',0.0);

/* placeholder for fire */
  FIREdummy = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FIREdummy");
  init_pvar(FIREdummy,NULL,'f',0.0);
        
  HYD_DOM_ACTWAT_PRES = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_DOM_ACTWAT_PRES");
  init_pvar(HYD_DOM_ACTWAT_PRES,NULL,'f',0.0);
  HYD_DOM_ACTWAT_VOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_DOM_ACTWAT_VOL");
  init_pvar(HYD_DOM_ACTWAT_VOL,NULL,'f',0.0);
  HYD_EVAP_CALC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_EVAP_CALC");
  init_pvar(HYD_EVAP_CALC,NULL,'f',0.0);
  HYD_MANNINGS_N = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_MANNINGS_N");
  init_pvar(HYD_MANNINGS_N,NULL,'f',0.0);
  HYD_RCCONDUCT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_RCCONDUCT");
  init_pvar(HYD_RCCONDUCT,NULL,'f',0.0);
  HYD_SAT_POT_TRANS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_SAT_POT_TRANS");
  init_pvar(HYD_SAT_POT_TRANS,NULL,'f',0.0);
  HYD_SED_WAT_VOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_SED_WAT_VOL");
  init_pvar(HYD_SED_WAT_VOL,NULL,'f',0.0);
  HYD_TOT_POT_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_TOT_POT_TRANSP");
  init_pvar(HYD_TOT_POT_TRANSP,NULL,'f',0.0);
  HydTotHd = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HydTotHd");
  init_pvar(HydTotHd,NULL,'f',0.0);
  HydRelDepPosNeg = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HydRelDepPosNeg");
  init_pvar(HydRelDepPosNeg,NULL,'f',0.0);
  HYD_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_TRANSP");
  init_pvar(HYD_TRANSP,NULL,'f',0.0);
  HYD_ET = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_ET");
  init_pvar(HYD_ET,NULL,'f',0.0);
  HYD_UNSAT_POT_TRANS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_UNSAT_POT_TRANS");
  init_pvar(HYD_UNSAT_POT_TRANS,NULL,'f',0.0);
  HYD_WATER_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_WATER_AVAIL");
  init_pvar(HYD_WATER_AVAIL,NULL,'f',0.0);

/* sfwmm rainfall, stage, and pET mapped to elm grid */
  boundcond_rain = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"boundcond_rain");
  init_pvar(boundcond_rain,NULL,'f',0.0);
  boundcond_depth = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"boundcond_depth");
  init_pvar(boundcond_depth,NULL,'f',0.0);
  boundcond_ETp = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"boundcond_ETp");
  init_pvar(boundcond_ETp,NULL,'f',0.0);

  SAT_TO_UNSAT_FL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_TO_UNSAT_FL");
  init_pvar(SAT_TO_UNSAT_FL,NULL,'f',0.0);
  SAT_VS_UNSAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_VS_UNSAT");
  init_pvar(SAT_VS_UNSAT,NULL,'f',0.0);
  SAT_WATER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WATER");
  init_pvar(SAT_WATER,NULL,'f',0.0);
  SAT_WT_HEAD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WT_HEAD");
  init_pvar(SAT_WT_HEAD,NULL,'f',0.0);
  SAT_WT_RECHG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WT_RECHG");
  init_pvar(SAT_WT_RECHG,NULL,'f',0.0);
  SAT_WT_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WT_TRANSP");
  init_pvar(SAT_WT_TRANSP,NULL,'f',0.0);
  SF_WT_EVAP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_EVAP");
  init_pvar(SF_WT_EVAP,NULL,'f',0.0);
  SF_WT_FROM_RAIN = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_FROM_RAIN");
  init_pvar(SF_WT_FROM_RAIN,NULL,'f',0.0);
  SF_WT_INFILTRATION = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_INFILTRATION");
  init_pvar(SF_WT_INFILTRATION,NULL,'f',0.0);
  SF_WT_POT_INF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_POT_INF");
  init_pvar(SF_WT_POT_INF,NULL,'f',0.0);
  SF_WT_TO_SAT_DOWNFLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_TO_SAT_DOWNFLOW");
  init_pvar(SF_WT_TO_SAT_DOWNFLOW,NULL,'f',0.0);
  SFWT_VOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SFWT_VOL");
  init_pvar(SFWT_VOL,NULL,'f',0.0);
  SURFACE_WAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SURFACE_WAT");
  init_pvar(SURFACE_WAT,NULL,'f',0.0);
  SF_WT_VEL_mag = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_VEL_mag");
  init_pvar(SF_WT_VEL_mag,NULL,'f',0.0); 
  BCmodel_sfwat = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"BCmodel_sfwat");
  init_pvar(BCmodel_sfwat,NULL,'f',0.0); 

  UNSAT_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_AVAIL");
  init_pvar(UNSAT_AVAIL,NULL,'f',0.0);
  UNSAT_CAP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_CAP");
  init_pvar(UNSAT_CAP,NULL,'f',0.0);
  UNSAT_DEPTH = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_DEPTH");
  init_pvar(UNSAT_DEPTH,NULL,'f',0.0);
  UNSAT_HYD_COND_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_HYD_COND_CF");
  init_pvar(UNSAT_HYD_COND_CF,NULL,'f',0.0);
  UNSAT_MOIST_PRP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_MOIST_PRP");
  init_pvar(UNSAT_MOIST_PRP,NULL,'f',0.0);
  UNSAT_PERC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_PERC");
  init_pvar(UNSAT_PERC,NULL,'f',0.0);
  UNSAT_TO_SAT_FL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_TO_SAT_FL");
  init_pvar(UNSAT_TO_SAT_FL,NULL,'f',0.0);
  UNSAT_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_TRANSP");
  init_pvar(UNSAT_TRANSP,NULL,'f',0.0);
  UNSAT_WATER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_WATER");
  init_pvar(UNSAT_WATER,NULL,'f',0.0);
  UNSAT_WT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_WT_POT");
  init_pvar(UNSAT_WT_POT,NULL,'f',0.0);

  MAC_HEIGHT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_HEIGHT");
  init_pvar(MAC_HEIGHT,NULL,'f',0.0);
  MAC_LAI = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_LAI");
  init_pvar(MAC_LAI,NULL,'f',0.0);
  LAI_eff = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"LAI_eff");
  init_pvar(LAI_eff,NULL,'f',0.0);
  MAC_LIGHT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_LIGHT_CF");
  init_pvar(MAC_LIGHT_CF,NULL,'f',0.0);
  MAC_MAX_BIO = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_MAX_BIO");
  init_pvar(MAC_MAX_BIO,NULL,'f',0.0);
  MAC_NOPH_BIOMAS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_NOPH_BIOMAS");
  init_pvar(MAC_NOPH_BIOMAS,NULL,'f',0.0);
  MAC_NUT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_NUT_CF");
  init_pvar(MAC_NUT_CF,NULL,'f',0.0);
  MAC_PH_BIOMAS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_PH_BIOMAS");
  init_pvar(MAC_PH_BIOMAS,NULL,'f',0.0);
  MAC_PROD_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_PROD_CF");
  init_pvar(MAC_PROD_CF,NULL,'f',0.0);
  MAC_REL_BIOM = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_REL_BIOM");
  init_pvar(MAC_REL_BIOM,NULL,'f',0.0);
  MAC_SALT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_SALT_CF");
  init_pvar(MAC_SALT_CF,NULL,'f',0.0);
  MAC_TEMP_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_TEMP_CF");
  init_pvar(MAC_TEMP_CF,NULL,'f',0.0);
  MAC_TOT_BIOM = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_TOT_BIOM");
  init_pvar(MAC_TOT_BIOM,NULL,'f',0.0);
  MAC_WATER_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_WATER_CF");
  init_pvar(MAC_WATER_CF,NULL,'f',0.0);
  MAC_WATER_AVAIL_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_WATER_AVAIL_CF");
  init_pvar(MAC_WATER_AVAIL_CF,NULL,'f',0.0);
  NPHBIO_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_AVAIL");
  init_pvar(NPHBIO_AVAIL,NULL,'f',0.0);
  NPHBIO_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_MORT");
  init_pvar(NPHBIO_MORT,NULL,'f',0.0);
  NPHBIO_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_MORT_POT");
  init_pvar(NPHBIO_MORT_POT,NULL,'f',0.0);
  NPHBIO_REFUGE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_REFUGE");
  init_pvar(NPHBIO_REFUGE,NULL,'f',0.0);
  NPHBIO_SAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_SAT");
  init_pvar(NPHBIO_SAT,NULL,'f',0.0);
  NPHBIO_TRANSLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_TRANSLOC");
  init_pvar(NPHBIO_TRANSLOC,NULL,'f',0.0);
  NPHBIO_TRANSLOC_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_TRANSLOC_POT");
  init_pvar(NPHBIO_TRANSLOC_POT,NULL,'f',0.0);
  PHBIO_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_AVAIL");
  init_pvar(PHBIO_AVAIL,NULL,'f',0.0);
  PHBIO_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_MORT");
  init_pvar(PHBIO_MORT,NULL,'f',0.0);
  PHBIO_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_MORT_POT");
  init_pvar(PHBIO_MORT_POT,NULL,'f',0.0);
  PHBIO_NPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_NPP");
  init_pvar(PHBIO_NPP,NULL,'f',0.0);
  PHBIO_REFUGE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_REFUGE");
  init_pvar(PHBIO_REFUGE,NULL,'f',0.0);
  PHBIO_SAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_SAT");
  init_pvar(PHBIO_SAT,NULL,'f',0.0); 
  PHBIO_TRANSLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_TRANSLOC");
  init_pvar(PHBIO_TRANSLOC,NULL,'f',0.0);
  PHBIO_TRANSLOC_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_TRANSLOC_POT");
  init_pvar(PHBIO_TRANSLOC_POT,NULL,'f',0.0);

  phbio_npp_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_npp_P");
  init_pvar(phbio_npp_P,NULL,'l',0.0);
  phbio_npp_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_npp_OM");
  init_pvar(phbio_npp_OM,NULL,'l',0.0);
  phbio_mort_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_mort_P");
  init_pvar(phbio_mort_P,NULL,'l',0.0);
  phbio_mort_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_mort_OM");
  init_pvar(phbio_mort_OM,NULL,'l',0.0);
  phbio_transl_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_transl_P");
  init_pvar(phbio_transl_P,NULL,'l',0.0);
  phbio_transl_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_transl_OM");
  init_pvar(phbio_transl_OM,NULL,'l',0.0);
  nphbio_transl_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_transl_P");
  init_pvar(nphbio_transl_P,NULL,'l',0.0);
  nphbio_transl_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_transl_OM");
  init_pvar(nphbio_transl_OM,NULL,'l',0.0);
  nphbio_mort_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_mort_P");
  init_pvar(nphbio_mort_P,NULL,'l',0.0);
  nphbio_mort_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_mort_OM");
  init_pvar(nphbio_mort_OM,NULL,'l',0.0);
  mac_nph_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_P");
  init_pvar(mac_nph_P,NULL,'l',0.0);
  mac_nph_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_PC");
  init_pvar(mac_nph_PC,NULL,'l',0.0);
  mac_nph_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_OM");
  init_pvar(mac_nph_OM,NULL,'l',0.0);
  mac_nph_CtoOM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_CtoOM");
  init_pvar(mac_nph_CtoOM,NULL,'l',0.0);
  mac_ph_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_P");
  init_pvar(mac_ph_P,NULL,'l',0.0);
  mac_ph_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_PC");
  init_pvar(mac_ph_PC,NULL,'l',0.0);
  mac_ph_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_OM");
  init_pvar(mac_ph_OM,NULL,'l',0.0);
  mac_ph_CtoOM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_CtoOM");
  init_pvar(mac_ph_CtoOM,NULL,'l',0.0);

  mac_nph_PC_rep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"mac_nph_PC_rep");
  init_pvar(mac_nph_PC_rep,NULL,'f',0.0);
  mac_ph_PC_rep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"mac_ph_PC_rep");
  init_pvar(mac_ph_PC_rep,NULL,'f',0.0);

  TP_SED_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SED_WT");
  init_pvar(TP_SED_WT,NULL,'l',0.0);
  TP_SED_CONC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SED_CONC");
  init_pvar(TP_SED_CONC,NULL,'l',0.0);
  TP_SEDWT_CONCACT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SEDWT_CONCACT");
  init_pvar(TP_SEDWT_CONCACT,NULL,'l',0.0);
  TP_SF_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SF_WT");
  init_pvar(TP_SF_WT,NULL,'l',0.0);
  TP_SFWT_CONC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SFWT_CONC");
  init_pvar(TP_SFWT_CONC,NULL,'l',0.0);
  TP_SORB = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SORB");
  init_pvar(TP_SORB,NULL,'l',0.0);
  TP_SORBCONC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SORBCONC");
  init_pvar(TP_SORBCONC,NULL,'l',0.0);
  TP_SED_WT_AZ = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SED_WT_AZ");
  init_pvar(TP_SED_WT_AZ,NULL,'l',0.0);
  TP_Act_to_Tot = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_Act_to_Tot");
  init_pvar(TP_Act_to_Tot,NULL,'l',0.0);

  TP_Act_to_TotRep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_Act_to_TotRep");
  init_pvar(TP_Act_to_TotRep,NULL,'f',0.0);
  TP_SORBCONC_rep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SORBCONC_rep");
  init_pvar(TP_SORBCONC_rep,NULL,'f',0.0);

  TP_DNFLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_DNFLOW");
  init_pvar(TP_DNFLOW,NULL,'f',0.0);
  TP_DNFLOW_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_DNFLOW_POT");
  init_pvar(TP_DNFLOW_POT,NULL,'f',0.0);
  TP_Atm_Depos = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_Atm_Depos");
  init_pvar(TP_Atm_Depos,NULL,'f',0.0);
/* v2.6 added map of atmospheric TP deposition rate */
  TP_AtmosDepos = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_AtmosDepos");
  init_pvar(TP_AtmosDepos,NULL,'f',0.0);
  TP_K = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_K");
  init_pvar(TP_K,NULL,'f',0.0);
  TP_SED_MINER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SED_MINER");
  init_pvar(TP_SED_MINER,NULL,'f',0.0);
  TP_SEDWT_CONCACTMG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SEDWT_CONCACTMG");
  init_pvar(TP_SEDWT_CONCACTMG,NULL,'f',0.0);
  TP_SEDWT_UPTAKE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SEDWT_UPTAKE");
  init_pvar(TP_SEDWT_UPTAKE,NULL,'f',0.0);
  TP_SFWT_CONC_MG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SFWT_CONC_MG");
  init_pvar(TP_SFWT_CONC_MG,NULL,'f',0.0);
  TP_SFWT_MINER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SFWT_MINER");
  init_pvar(TP_SFWT_MINER,NULL,'f',0.0);
  TP_SFWT_UPTAK = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SFWT_UPTAK");
  init_pvar(TP_SFWT_UPTAK,NULL,'f',0.0);
  TP_settl = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_settl");
  init_pvar(TP_settl,NULL,'f',0.0);
  TP_SORB_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SORB_POT");
  init_pvar(TP_SORB_POT,NULL,'f',0.0);
  TP_SORBTION = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SORBTION");
  init_pvar(TP_SORBTION,NULL,'f',0.0);
  TP_UPFLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_UPFLOW");
  init_pvar(TP_UPFLOW,NULL,'f',0.0);
  TP_UPFLOW_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_UPFLOW_POT");
  init_pvar(TP_UPFLOW_POT,NULL,'f',0.0);
  P_SUM_CELL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"P_SUM_CELL");
  init_pvar(P_SUM_CELL,NULL,'f',0.0);

  DINdummy = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DINdummy");
  init_pvar(DINdummy,NULL,'l',0.0);

  WQMsettlVel = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"WQMsettlVel");
  init_pvar(WQMsettlVel,NULL,'f',0.0);
  TP_settlDays = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_settlDays");
  init_pvar(TP_settlDays,NULL,'f',0.0);

  SAL_SED_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SAL_SED_WT");
  init_pvar(SAL_SED_WT,NULL,'l',0.0);
  SAL_SF_WT_mb = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SAL_SF_WT_mb");
  init_pvar(SAL_SF_WT_mb,NULL,'l',0.0);
  SALT_SED_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SALT_SED_WT");
  init_pvar(SALT_SED_WT,NULL,'l',0.0);
  SALT_SURF_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SALT_SURF_WT");
  init_pvar(SALT_SURF_WT,NULL,'l',0.0);
        
  SAL_SF_WT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAL_SF_WT");
  init_pvar(SAL_SF_WT,NULL,'f',0.0);
  SALT_Atm_Depos = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_Atm_Depos");
  init_pvar(SALT_Atm_Depos,NULL,'f',0.0);
  SALT_SED_TO_SF_FLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_SED_TO_SF_FLOW");
  init_pvar(SALT_SED_TO_SF_FLOW,NULL,'f',0.0);
  SALT_SFWAT_DOWNFL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_SFWAT_DOWNFL");
  init_pvar(SALT_SFWAT_DOWNFL,NULL,'f',0.0);
  SALT_SFWAT_DOWNFL_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_SFWAT_DOWNFL_POT");
  init_pvar(SALT_SFWAT_DOWNFL_POT,NULL,'f',0.0);

  FLOC_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DECOMP");
  init_pvar(FLOC_DECOMP,NULL,'f',0.0);
  FLOC_DECOMP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DECOMP_POT");
  init_pvar(FLOC_DECOMP_POT,NULL,'f',0.0);
  FLOC_DECOMP_QUAL_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DECOMP_QUAL_CF");
  init_pvar(FLOC_DECOMP_QUAL_CF,NULL,'f',0.0);
  FLOC_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_Z");
  init_pvar(FLOC_Z,NULL,'f',0.0);
  FLOC_DEPO = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DEPO");
  init_pvar(FLOC_DEPO,NULL,'f',0.0);
  FLOC_DEPO_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DEPO_POT");
  init_pvar(FLOC_DEPO_POT,NULL,'f',0.0);
  FLOC_FR_ALGAE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_FR_ALGAE");
  init_pvar(FLOC_FR_ALGAE,NULL,'f',0.0);
  FLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC");
  init_pvar(FLOC,NULL,'f',0.0);
  FlocP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP");
  init_pvar(FlocP,NULL,'f',0.0);
  FlocP_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_DECOMP");
  init_pvar(FlocP_DECOMP,NULL,'f',0.0);
  FlocP_FR_ALGAE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_FR_ALGAE");
  init_pvar(FlocP_FR_ALGAE,NULL,'f',0.0);
  FlocP_PhBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_PhBio");
  init_pvar(FlocP_PhBio,NULL,'f',0.0);
  FlocP_DEPO = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_DEPO");
  init_pvar(FlocP_DEPO,NULL,'f',0.0);
  FlocP_OMrep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_OMrep");
  init_pvar(FlocP_OMrep,NULL,'f',0.0);
  FlocP_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"FlocP_OM");
  init_pvar(FlocP_OM,NULL,'l',0.0);

  SOLRADGRD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SOLRADGRD");
  init_pvar(SOLRADGRD,NULL,'f',0.0);
  SOLRAD274 = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SOLRAD274");
  init_pvar(SOLRAD274,NULL,'f',0.0);
  AIR_TEMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"AIR_TEMP");
  init_pvar(AIR_TEMP,NULL,'f',0.0);
  H2O_TEMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"H2O_TEMP");
  init_pvar(H2O_TEMP,NULL,'f',0.0);

/* habitat-specific parameter arrays */
  HP_ALG_MAX = NULL;

  HP_DOM_MAXDEPTH = NULL;
  HP_DOM_AEROBTHIN = NULL;

  HP_TP_CONC_GRAD = NULL;

  HP_SALT_ICSEDWAT = NULL;
  HP_SALT_ICSFWAT = NULL;

  HP_PHBIO_MAX = NULL;
  HP_NPHBIO_MAX = NULL;
  HP_MAC_MAXHT = NULL;
  HP_NPHBIO_ROOTDEPTH = NULL;
  HP_MAC_MAXROUGH = NULL;
  HP_MAC_MINROUGH = NULL;
  HP_MAC_MAXLAI = NULL;
  HP_MAC_MAXCANOPCOND = NULL;
  HP_MAC_CANOPDECOUP = NULL;
  HP_MAC_TEMPOPT = NULL;
  HP_MAC_LIGHTSAT = NULL;
  HP_MAC_KSP = NULL;
  HP_PHBIO_RCNPP = NULL;
  HP_PHBIO_RCMORT = NULL;
  HP_MAC_WAT_TOLER = NULL;
  HP_MAC_SALIN_THRESH = NULL;
  HP_PHBIO_IC_CTOOM = NULL;
  HP_NPHBIO_IC_CTOOM = NULL;
  HP_PHBIO_IC_PC = NULL;
  HP_NPHBIO_IC_PC = NULL;
  HP_MAC_TRANSLOC_RC = NULL;

  HP_HYD_RCINFILT = NULL;
  HP_HYD_SPEC_YIELD = NULL;
  HP_HYD_POROSITY = NULL;

  HP_FLOC_IC = NULL;
  HP_FLOC_IC_CTOOM = NULL;
  HP_FLOC_IC_PC = NULL;

  HP_SfDepthLo = NULL;
  HP_SfDepthHi = NULL;
  HP_SfDepthInt = NULL;
  HP_PhosLo = NULL;
  HP_PhosHi = NULL;
  HP_PhosInt = NULL;
  HP_FireInt = NULL;
  HP_SalinLo = NULL;
  HP_SalinHi = NULL;
  HP_SalinInt = NULL;
  
  usrErr("Done."); /* console message */

}


/*******
the remaining functions are time-series graph interpolaters
******/

float graph7(unsigned char y, float x)
{
  float s;
  int ig=0, Np=10;
  while(1) {
  if (x <= g7[ig][0]) { if(ig>0) ig=ig-1; else return(g7[0][1+y]);}
  else if (x > g7[ig][0] && x <= g7[ig+1][0]) {
         s =   (g7[ig+1][1+y]-g7[ig][1+y])/
            (g7[ig+1][0]-g7[ig][0]);
         return(s * (x-g7[ig][0]) + g7[ig][1+y]); }
    else if (x > g7[ig+1][0]) { if(ig<Np) ig=ig+1; else return(g7[Np][1+y]);}
  }
}

float graph8(unsigned char y, float x)
{
  float s;
  int ig=0, Np=10;
  while(1) {
  if (x <= g8[ig][0]) { if(ig>0) ig=ig-1; else return(g8[0][1+y]);}
  else if (x > g8[ig][0] && x <= g8[ig+1][0]) {
         s =   (g8[ig+1][1+y]-g8[ig][1+y])/
            (g8[ig+1][0]-g8[ig][0]);
         return(s * (x-g8[ig][0]) + g8[ig][1+y]); }
    else if (x > g8[ig+1][0]) { if(ig<Np) ig=ig+1; else return(g8[Np][1+y]);}
  }
}


/***** end UnitMod.c source ***/