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

UnitMod.c

Go to the documentation of this file.
00001 
00038 /*###############
00039 To change among model grid scales of different resolution/extent: 
00040   no code changes are needed here or elsewhere
00041 ###############
00042 
00043  ***
00044  NOTE: The init_eqns module is not fully updated for properly
00045  turning off individual **interacting** *biological/chemical* modules at runtime.  
00046  If one *biological/chemical* module is turned off, 
00047  they all need to be turned off. 
00048 
00049  The following *biological/chemical* modules must be ON or OFF as a group:
00050  (cell_dyn2 + cell_dyn4 + cell_dyn8 + cell_dyn9 + cell_dyn12)
00051  cell_dyn13, the net settling rate module, can be turned on only when above module group is off
00052  ***
00053 
00054 */
00055 
00056   
00057 /* General NOTES on revisions to this source file:
00058         april 2000 VERSION 2.1 - output available for application for water quality performance
00059         july 2002 VERSIOIN 2.1a - first widely-available public release of code -
00060               - misc modifications, such as revising method and types of
00061               parameters that are read in (reduced number of parms in "habitat-specific" dbase);
00062               - added some more code documentation, etc.
00063         July 2003 VERSION 2.2.0 - moved declarations into header file, other similar changes
00064               to help reviewers better understand organization of codes.
00065               - Incorporated new rainfall module (highly revised rain_inp.c) that generalizes the
00066               code to better accomodate implementations of all scale; NO DATA CHANGE - passes identical data
00067               from SFWMM data file to the ELM.
00068               - Included are some specific, minor code changes (that do not alter calculations)  
00069         July 2004 VERSION 2.2.1 - improved code documentation
00070 ******        - AN ESSENTIAL COMPANION TO THIS CODE DOCUMENTATION IS NOW THE "ModelOutlist_creator_version#.xls"
00071 ******            workbook (OpenOffice/Excel).  In this workbook, all of the variables
00072 ******            that are spatial arrays (and other types) are defined/described (including their units).
00073               - In order to avoid the "clutter" of comments interspersing virtually every equation,
00074               this workbook was used to fully define the use of the variables in the model.
00075               - The code here STILL contains comments that describe the variables that may be most 
00076               important or that are better understood by explicitly stating underlying assumptions. 
00077               - Included are several specific, minor code changes (still maintaining v2.1a calcs). 
00078         Aug 2004 VERSION 2.3.0 - added bathymetry map input (dynamic boundary conditions in other src files -> version upgrade)
00079                 - added bathymetry map input
00080                 - added new spatial time series data input for dynamic boundary conditions
00081                 - version increment to 2.3 due to changes in other source files
00082         Oct 2004 VERSION 2.3.1 - internal release 
00083                 - checked that results reasonable, not fully checked
00084         Nov 2004 VERSION 2.3.2 - documentation upgrade
00085                 - removed standing detritus (was operative in v2.3.0 & earlier), consumer, and fire modules completely from code
00086                 - removed extraneous comments, inoperative variables, inoperative parameters
00087                 - added Doxygen tags in comments
00088         Apr 2005 v2.4.4: changed init maps: icMult removed, icMacBio added
00089         May 2006 v2.5.0: added some debug-level output for boundary condition stage
00090 
00091 */
00092 
00093 #include "unitmod.h"
00094 #include "unitmod_vars.h"
00095 #include "unitmod_habparms.h"
00096 #include "unitmod_globparms.h"
00097 #include "budgstats_statvars.h"  
00098 
00119 int call_cell_dyn(int sector,int step)
00120  {
00121   int rv=0;
00122 
00123   switch(sector) {
00124 
00125     case 99: { stats(step); rv=1; } break;
00126     case 0: { horizFlow(step); rv=1; } break;
00127     case 1: { cell_dyn1(step); rv=1; } break;
00128     case 2: { cell_dyn2(step); rv=1; } break;
00129     case 4: { cell_dyn4(step); rv=1; } break;
00130     case 7: { cell_dyn7(step); rv=1; } break;
00131     case 8: { cell_dyn8(step); rv=1; } break;
00132     case 9: { cell_dyn9(step); rv=1; } break;
00133     case 10: { cell_dyn10(step); rv=1; } break;
00134     case 12: { cell_dyn12(step); rv=1; } break;
00135     case 13: { cell_dyn13(step); rv=1; } break;
00136             default:  printf("Warning, undefined sector number:%d\n",sector);
00137   }
00138   return rv;
00139 }
00140 
00141 
00142 /*###############
00143 The following modules (cell_dyn*) are the dynamic calcs for the model.
00144 They are called in the order determined by the call_cell_dyn function.
00145 ################*/
00146 
00147 
00148 /*******/
00157 /* Parameter definitions:
00158       global parameters in GlobalParms.xls OpenOffice/Excel sheet (text export=GlobalParms)
00159       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
00160 ******/
00161 
00162 void cell_dyn1(int step)
00163 {
00164    int fail = -1;
00165    int ix, iy, cellLoc, stat=1; 
00166    
00167    /* v2.5 may be temporary? */
00168    FILE *boundaryFile; 
00169    char filename[40];
00170 
00171      /* SimTime.TIME is a (float) julian day counter since initialization (calc'd in main function, Generic_Driver.c source).  
00172          (SimTime.TIME includes fractions of days if the vertical DT<1.0, but it is
00173          unlikely that the vertical DT will deviate from 1 day). */
00174 /* daily re-mapping of coarse-grid input data */
00175      if (FMOD(SimTime.TIME,1.0) < 0.001) {
00176  
00177       /* remap sfwmm (or other grid) rain data to ELM grid scale */
00178        stat=rain_data_wmm(boundcond_rain); 
00179        if(stat==fail)  
00180        { 
00181          sprintf(msgStr,"Problem with rainfall data, time = %f.\n",SimTime.TIME); 
00182          WriteMsg(msgStr,1); 
00183          exit(fail);    
00184        }
00185 
00186      /* remap sfwmm (or other grid) stage/water depth data to ELM grid scale */
00187       stat=stage_data_wmm(boundcond_depth);  
00188      if(stat==-1)
00189      {
00190        sprintf(msgStr,"Problem with stage data, time = %f.\n",SimTime.TIME);
00191        WriteMsg(msgStr,1);
00192        exit(fail);
00193      } 
00194      
00195     /* ***** */
00196     /* v2.5 should below become permanent? Output the sfwmm boundary stage data */
00197     if ( debug > 4 )    
00198          {
00199                 sprintf( filename, "%s/%s/Output/Debug/wmm_depth_%d.txt", OutputPath, ProjName, (int)SimTime.TIME );
00200    
00201 
00202          if ( ( boundaryFile = fopen( filename, "w" ) ) ==  NULL )
00203          {
00204           printf( "Can't open the %s boundary condition depth debug file!\n", filename );
00205            exit(-1) ;
00206          }
00207          /* fprintf ( boundaryFile, "ROWS=%d\nCOLUMNS=%d\nNotes\"Debug data: SFWMM output, ponded depths (pos/neg). \n", s0, s1 ); */
00208 
00209          for ( ix = 1; ix <= s0 ; ix++ ) 
00210            { for ( iy = 1 ; iy <= s1 ; iy++ )
00211                 fprintf( boundaryFile, "%f\t", boundcond_depth[T(ix,iy)]) ;
00212             fprintf ( boundaryFile, "\n" );
00213           }   
00214          fclose( boundaryFile ) ;
00215          if (SimTime.TIME>60) { 
00216                 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); 
00217                 exit(-1); 
00218          }
00219          }
00220          /* v2.5 should above become permanent? Output the sfwmm boundary stage data */
00221      /* ***** */
00222 
00223 
00224      /* remap sfwmm (or other grid) potential ET data to ELM grid scale */
00225      stat=evap_data_wmm(boundcond_ETp); 
00226      if(stat==-1)
00227      {
00228        sprintf(msgStr,"Problem with ETp data, time = %f.\n",SimTime.TIME);
00229        WriteMsg(msgStr,1);
00230        exit(fail);
00231      } 
00232      
00233      } /* end of mapping multi-grid data */
00234 
00235 
00236             
00237     /* the "julian" day counter within a 365 day "year"  */
00238     DAYJUL = ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 365.0);
00239     /* DAYLENGTH is not used */
00240     /* DAYLENGTH = AMPL*Sin((DAYJUL-79.0)*0.01721)+12.0; */ /* length of daylight (hours) */
00241 
00242                 
00243 /*  Nikolov & Zeller (1992) solar radiation algorithm 
00244    the algorithm and all parameters are published in the (Ecol. Mod., 61:149-168) manuscript,
00245    and the only modifiable parameters (in ELM database system) are local latitude and altitude */
00246     SOLDEC1 = 0.39785*Sin(4.868961+0.017203*DAYJUL   
00247                           +0.033446*Sin(6.224111+0.017202*DAYJUL));
00248     SOLCOSDEC = sqrt(1.0-SOLDEC1*SOLDEC1);
00249     SOLELEV_SINE = Sin(LATRAD)*SOLDEC1+Cos(LATRAD)*SOLCOSDEC;
00250     SOLALTCORR = (1.0-Exp(-0.014*(GP_ALTIT-274.0)/(SOLELEV_SINE*274.0)));
00251     SOLDEC = Arctan(SOLDEC1/sqrt(1.0-SOLDEC1*SOLDEC1));
00252     SOLRISSET_HA1 = -Tan(LATRAD)*Tan(SOLDEC);
00253     SOLRISSET_HA = ( (SOLRISSET_HA1==0.0) ) ? ( PI*0.5 ) : (   ( (SOLRISSET_HA1<0.0) ) ? 
00254                                                                ( PI+Arctan(sqrt(1.0-SOLRISSET_HA1*SOLRISSET_HA1)/SOLRISSET_HA1)  ) : 
00255                                                                (    Arctan(sqrt(1.0-SOLRISSET_HA1*SOLRISSET_HA1)/SOLRISSET_HA1)));
00256     SOLRADATMOS = 458.37*2.0*(1.0+0.033*Cos(360.0/365.0*PI/180.0*DAYJUL))
00257         * ( Cos(LATRAD)*Cos(SOLDEC)*Sin(SOLRISSET_HA) 
00258             + SOLRISSET_HA*180.0/(57.296*PI)*Sin(LATRAD)*Sin(SOLDEC));
00259             
00260   
00261         /* daily habitat switching */
00262     if ( (SimTime.TIME - (int)SimTime.TIME) < DT/2 && HabSwitchOn ) /* HabSwitchOn flag to invoke switching */
00263         for(ix=1; ix<=s0; ix++)  
00264             for(iy=1; iy<=s1; iy++)  
00265                 if(ON_MAP[cellLoc= T(ix,iy)]) { /* TPtoSoil == kg/kg */
00266                     HAB [cellLoc] = HabSwitch (ix, iy, SURFACE_WAT, TPtoSOIL, (int*)FIREdummy, HAB); 
00267                 } 
00268 }
00269 
00270 
00271 
00272 /*******/
00279 /*
00280   NC_ALG[cellLoc] == carbon mass of the NonCalcareous (more appropriately, the eutrophic, or non-native) periphyton community (gC/m^2)
00281   C_ALG[cellLoc] == carbon mass of the Calcareous (more appropriately, the oligotrophic, or native) periphyton community (gC/m^2)
00282 
00283   NC_ALG_P[cellLoc] == phosphorus mass of NonCalcareous (more appropriately, the eutrophic, or non-native) periphyton community (gP/m^2)
00284   C_ALG_P[cellLoc] == phosphorus mass of Calcareous (more appropriately, the oligotrophic, or native) periphyton community (gP/m^2)
00285   
00286 Parameter definitions:
00287       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
00288       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
00289 */
00290 void cell_dyn2(int step)
00291  {
00292  int ix, iy, cellLoc; 
00293  float reduc, min_litTemp,I_ISat, Z_extinct; 
00294  double PO4Pconc, PO4P;
00295  float C_ALG_thresh_CF, mortPot;
00296  
00297   for(ix=1; ix<=s0; ix++) {
00298   for(iy=1; iy<=s1; iy++) {
00299 
00300     if(ON_MAP[cellLoc= T(ix,iy)])  {
00301 
00302 /* these thresholds need updating when a habitat type of a grid cell changes */
00303      ALG_REFUGE[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*GP_ALG_REF_MULT;
00304      ALG_SAT[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*0.9;
00305 
00306      NC_ALG_AVAIL_MORT[cellLoc]  = Max(NC_ALG[cellLoc]-ALG_REFUGE[cellLoc],0);
00307      C_ALG_AVAIL_MORT[cellLoc]  = Max(C_ALG[cellLoc]-ALG_REFUGE[cellLoc],0);
00308 /* bio-avail P (PO4) is calc'd from TP, using pre-processed regression for predicting PO4 from TP */
00309 /* assume that periphyton (microbial) alkaline phosphotase activity keeps PO4 at least 10% of TP conc */
00310      PO4Pconc = Max(TP_SFWT_CONC_MG[cellLoc]*GP_PO4toTP + GP_PO4toTPint,
00311                                 0.10 * TP_SFWT_CONC_MG[cellLoc]); /* mg/L */
00312 
00313 /* light, water, temperature controls apply to both calc and non-calc */
00314      ALG_LIGHT_EXTINCT[cellLoc]  = GP_alg_light_ext_coef; 
00315                      /* algal self-shading implicit in density-dependent constraint function later */
00316      ALG_INCID_LIGHT[cellLoc]  = SOLRADGRD[cellLoc]*Exp(-MAC_LAI[cellLoc]*GP_ALG_SHADE_FACTOR);
00317      Z_extinct = SURFACE_WAT[cellLoc]*ALG_LIGHT_EXTINCT[cellLoc];
00318      I_ISat = ALG_INCID_LIGHT[cellLoc]/GP_ALG_LIGHT_SAT;
00319                      /* averaged over whole water column (based on Steele '65) */
00320      ALG_LIGHT_CF[cellLoc]  = ( Z_extinct > 0.0 ) ? 
00321                      ( 2.718/Z_extinct * (Exp(-I_ISat * Exp(-Z_extinct)) - Exp(-I_ISat)) ) :
00322                      (I_ISat*Exp(1.0-I_ISat));
00323                      /* low-water growth constraint ready for something better based on data */
00324      ALG_WAT_CF[cellLoc]  = ( SURFACE_WAT[cellLoc]>0.0 ) ? ( 1.0 ) : ( 0.0);
00325      ALG_TEMP_CF[cellLoc]  = tempCF(0, 0.20, GP_ALG_TEMP_OPT, 5.0, 40.0, H2O_TEMP[cellLoc]);
00326 
00327       min_litTemp = Min(ALG_LIGHT_CF[cellLoc],ALG_TEMP_CF[cellLoc]);
00328 
00329 /* the 2 communities have same form of growth response to avail phosphorus */
00330      NC_ALG_NUT_CF[cellLoc]  =
00331                      Exp(-GP_alg_uptake_coef * Max(GP_NC_ALG_KS_P-PO4Pconc, 0.0)/GP_NC_ALG_KS_P) ; /* mg/L */
00332      C_ALG_NUT_CF[cellLoc]  =
00333                      Exp(-GP_alg_uptake_coef * Max(GP_C_ALG_KS_P-PO4Pconc, 0.0)/GP_C_ALG_KS_P) ; /* mg/L */
00334 
00335          /* the form of the control function assumes that at very low
00336              P conc, the alkaline phosphotase activity of the microbial assemblage scavenges P, maintaining a minimum nutrient availability to community */
00337      NC_ALG_PROD_CF[cellLoc]  = Min(min_litTemp,ALG_WAT_CF[cellLoc])*Max(NC_ALG_NUT_CF[cellLoc], GP_alg_alkP_min);
00338      C_ALG_PROD_CF[cellLoc]   = Min(min_litTemp,ALG_WAT_CF[cellLoc])*Max(C_ALG_NUT_CF[cellLoc], GP_alg_alkP_min);
00339 
00340      NC_ALG_RESP_POT[cellLoc]  = 
00341          ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
00342          ( 0.0) :
00343          ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc]*NC_ALG_AVAIL_MORT[cellLoc] ); 
00344      C_ALG_RESP_POT[cellLoc]  = 
00345          ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
00346          ( 0.0) :
00347          ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc] *C_ALG_AVAIL_MORT[cellLoc] ); 
00348 
00349      NC_ALG_RESP[cellLoc]  =  
00350          ( NC_ALG_RESP_POT[cellLoc]*DT>NC_ALG_AVAIL_MORT[cellLoc] ) ? 
00351          ( NC_ALG_AVAIL_MORT[cellLoc]/DT ) : 
00352          ( NC_ALG_RESP_POT[cellLoc]);
00353      C_ALG_RESP[cellLoc]  =  
00354          ( C_ALG_RESP_POT[cellLoc]*DT>C_ALG_AVAIL_MORT[cellLoc] ) ? 
00355          ( C_ALG_AVAIL_MORT[cellLoc]/DT ) : 
00356          ( C_ALG_RESP_POT[cellLoc]);
00357                  
00358           /* this is the threshold control function that increases
00359               calcareous/native periph mortality (likely due to loss of calcareous sheath) as P conc. increases */
00360      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);
00361 
00362      NC_ALG_MORT_POT[cellLoc]  = 
00363          ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
00364          ( NC_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT_DRY ) : 
00365          ( NC_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT);
00366      C_ALG_MORT_POT[cellLoc]  = 
00367          ( UNSAT_DEPTH[cellLoc]>GP_algMortDepth ) ? 
00368          ( C_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT_DRY ) : 
00369          ( C_ALG_thresh_CF * C_ALG_AVAIL_MORT[cellLoc]*GP_ALG_RC_MORT);
00370 
00371      NC_ALG_MORT[cellLoc]  = 
00372          ( (NC_ALG_MORT_POT[cellLoc]+NC_ALG_RESP[cellLoc])*DT>NC_ALG_AVAIL_MORT[cellLoc] ) ? 
00373          ( (NC_ALG_AVAIL_MORT[cellLoc]-NC_ALG_RESP[cellLoc]*DT)/DT ) : 
00374          ( NC_ALG_MORT_POT[cellLoc]);
00375      C_ALG_MORT[cellLoc]  = 
00376          ( (C_ALG_MORT_POT[cellLoc]+C_ALG_RESP[cellLoc])*DT>C_ALG_AVAIL_MORT[cellLoc] ) ? 
00377          ( (C_ALG_AVAIL_MORT[cellLoc]-C_ALG_RESP[cellLoc]*DT)/DT ) : 
00378          ( C_ALG_MORT_POT[cellLoc]);
00379                  
00380 /* gross production of the 2 communities (gC/m2/d, NOT kgC/m2/d) */
00381 /* density constraint on both noncalc and calc, competition effect accentuated by calc algae */
00382 
00383      NC_ALG_GPP[cellLoc]  =  NC_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*NC_ALG[cellLoc]       
00384              *Max( (1.0-(GP_AlgComp*C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
00385      C_ALG_GPP[cellLoc]  =  C_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*C_ALG[cellLoc] 
00386        *Max( (1.0-(        C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
00387        
00388                      
00389 /* P uptake is dependent on available P and is relative to a maximum P:C ratio for the tissue
00390    (g C/m^2/d * P:Cmax * dimless * dimless = gP/m2/d (NOT kg) )*/
00391      NC_ALG_GPP_P[cellLoc] = NC_ALG_GPP[cellLoc] *GP_ALG_PC * NC_ALG_NUT_CF[cellLoc]
00392                      * Max(1.0-NC_ALG_PC[cellLoc]/GP_ALG_PC, 0.0); 
00393      C_ALG_GPP_P[cellLoc] = C_ALG_GPP[cellLoc] * GP_ALG_PC * C_ALG_NUT_CF[cellLoc]
00394                      * Max(1.0-C_ALG_PC[cellLoc]/GP_ALG_PC, 0.0); 
00395                  
00396 /* check for available P mass (the nutCF does not) */
00397      PO4P = ramp(Min(PO4Pconc * SFWT_VOL[cellLoc],  conv_kgTOg*TP_SF_WT[cellLoc]) ); /*g P available; v2.5 put in ramp (constrain non-neg) */
00398      reduc = ( (NC_ALG_GPP_P[cellLoc]+C_ALG_GPP_P[cellLoc]) > 0) ? 
00399                      (PO4P / ( (NC_ALG_GPP_P[cellLoc]+C_ALG_GPP_P[cellLoc])*CELL_SIZE*DT) ) :
00400                      (1.0);
00401     /* can have high conc, but low mass of P avail, in presence of high peri biomass and high demand */
00402     /* reduce the production proportionally if excess demand is found */
00403     if (reduc < 1.0) {
00404                      NC_ALG_GPP[cellLoc] *= reduc;   
00405                      NC_ALG_GPP_P[cellLoc] *= reduc;   
00406                      C_ALG_GPP[cellLoc] *= reduc; 
00407                      C_ALG_GPP_P[cellLoc] *= reduc; 
00408                  }
00409 
00410 /* state variables calc'd (gC/m2, NOT kgC/m2) */
00411      NC_ALG[cellLoc] =  NC_ALG[cellLoc] 
00412          + (NC_ALG_GPP[cellLoc]
00413           - NC_ALG_RESP[cellLoc] - NC_ALG_MORT[cellLoc]) * DT;
00414 
00415      C_ALG[cellLoc] =  C_ALG[cellLoc] 
00416          + (C_ALG_GPP[cellLoc] 
00417          - C_ALG_RESP[cellLoc]  - C_ALG_MORT[cellLoc]) * DT;
00418 
00419 /* carbon NPP not currently used elsewhere, only for output */
00420      NC_ALG_NPP[cellLoc]  = NC_ALG_GPP[cellLoc]-NC_ALG_RESP[cellLoc]; 
00421      C_ALG_NPP[cellLoc]  = C_ALG_GPP[cellLoc]-C_ALG_RESP[cellLoc];                 
00422 /* carbon total biomass of both communities, not currently used elsewhere, only for output */
00423      ALG_TOT[cellLoc] = NC_ALG[cellLoc] + C_ALG[cellLoc];
00424                  
00425                  
00426 /*  now calc the P  associated with the C fluxes (GPP_P already calc'd) */
00427      mortPot = (double) NC_ALG_MORT[cellLoc] * NC_ALG_PC[cellLoc];
00428      NC_ALG_MORT_P[cellLoc] = (mortPot*DT>NC_ALG_P[cellLoc]) ?
00429                      (NC_ALG_P[cellLoc]/DT) :
00430                      (mortPot);
00431      mortPot = (double) C_ALG_MORT[cellLoc] * C_ALG_PC[cellLoc];
00432      C_ALG_MORT_P[cellLoc] = (mortPot*DT>C_ALG_P[cellLoc]) ?
00433                      (C_ALG_P[cellLoc]/DT) :
00434                      (mortPot);
00435                  
00436 
00437 /* state variables calc'd (gP/m2, NOT kgP/m2) */
00438      NC_ALG_P[cellLoc] =  NC_ALG_P[cellLoc] 
00439          + (NC_ALG_GPP_P[cellLoc] - NC_ALG_MORT_P[cellLoc]) * DT;
00440 
00441      C_ALG_P[cellLoc] =  C_ALG_P[cellLoc] 
00442          + (C_ALG_GPP_P[cellLoc] - C_ALG_MORT_P[cellLoc]) * DT;
00443 
00444      NC_ALG_PC[cellLoc] = (NC_ALG[cellLoc]>0.0) ?
00445                      (NC_ALG_P[cellLoc]/ NC_ALG[cellLoc]) :
00446                      ( GP_ALG_PC * 0.03); /* default to 3% of  max P:C */
00447      C_ALG_PC[cellLoc] = (C_ALG[cellLoc]>0.0) ?
00448                      (C_ALG_P[cellLoc]/ C_ALG[cellLoc]) :
00449                      ( GP_ALG_PC * 0.03 ); /* default to 3% of max P:C */
00450      NC_ALG_PCrep[cellLoc] = (float)NC_ALG_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */
00451      C_ALG_PCrep[cellLoc] = (float)C_ALG_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */
00452 
00453      if (debug > 0 && NC_ALG[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NC_ALG C biomass (%f m) in cell (%d,%d)!", 
00454                               SimTime.TIME, NC_ALG[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
00455                  if (debug > 0 && C_ALG[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg C_ALG C biomass (%f m) in cell (%d,%d)!", 
00456                               SimTime.TIME, C_ALG[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
00457                  if (debug > 0 && NC_ALG_P[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NC_ALG_P P biomass (%f m) in cell (%d,%d)!", 
00458                               SimTime.TIME, NC_ALG_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
00459                  if (debug > 0 && C_ALG_P[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg C_ALG_P P biomass (%f m) in cell (%d,%d)!", 
00460                               SimTime.TIME, C_ALG_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
00461 
00462                      
00463      TP_SFWT_UPTAK[cellLoc]  = (NC_ALG_GPP_P[cellLoc]+C_ALG_GPP_P[cellLoc])
00464                      *0.001*CELL_SIZE; /* gP/m2 => kg P */
00465 /* recalc P in surface water state variable (kg P) */
00466      TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] - TP_SFWT_UPTAK[cellLoc] * DT;
00467      TP_SFWT_CONC[cellLoc]  = 
00468          ( SFWT_VOL[cellLoc] > 0.0 ) ? 
00469          ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
00470          ( 0.0); /* used in P fluxes for mass balance */
00471      TP_SFWT_CONC_MG[cellLoc]  = 
00472          ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
00473          (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
00474          (0.0); /* (g/m3==mg/L) used for reporting and other modules to evaluate P conc when water is present */
00475 
00476           }
00477   }
00478   }
00479 }
00480 
00481 
00482 
00483 
00484 /*******/
00491 /*     DEPOS_ORG_MAT[cellLoc] == mass of Deposited Organic Matter (DOM/AOM) of soil zone, w/o floc layer (kg OM/m^2)
00492 
00493      DOP[cellLoc] ==  mass of Deposited Organic Phosphorus of soil zone, w/o floc layer, w/o P sorbed (kg P/m^2)
00494   
00495 Parameter definitions:
00496       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
00497       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
00498 */
00499 void cell_dyn4(int step)
00500  {
00501 int ix, iy, cellLoc;
00502  float TPsoil, TP_sedMin, TP_sfMin;
00503   
00504  
00505 
00506   for(ix=1; ix<=s0; ix++) {
00507   for(iy=1; iy<=s1; iy++) {
00508 
00509     if(ON_MAP[cellLoc= T(ix,iy)])  {
00510 
00511 
00512 /* inputs of organic matter (kg OM/m2)*/
00513      DOM_fr_nphBio[cellLoc] = nphbio_mort_OM[cellLoc];
00514      DOM_FR_FLOC[cellLoc]  =  FLOC_DEPO[cellLoc] ; 
00515 
00516 /* losses of organic matter (kg OM/m2) */
00517 
00518      DOM_QUALITY_CF[cellLoc]  =
00519           Min(Exp(-GP_DOM_decomp_coef * Max(GP_DOM_DECOMP_POPT-TP_SEDWT_CONCACTMG[cellLoc], 0.0)
00520           /GP_DOM_DECOMP_POPT),1.0) ; /* mg/L */
00521           DOM_TEMP_CF[cellLoc] = tempCF(0, 0.20, GP_DOM_DECOMP_TOPT, 5.0, 40.0, H2O_TEMP[cellLoc]);
00522      DOM_SED_AEROB_Z[cellLoc]  = Min(Max(UNSAT_DEPTH[cellLoc],HP_DOM_AEROBTHIN[HAB[cellLoc]]),HP_DOM_MAXDEPTH[HAB[cellLoc]]);
00523      DOM_SED_ANAEROB_Z[cellLoc]  = HP_DOM_MAXDEPTH[HAB[cellLoc]]-DOM_SED_AEROB_Z[cellLoc];
00524 
00525                      /* GP_calibDecomp is an adjustable calib parm  */
00526      DOM_DECOMP_POT[cellLoc]  = (double) GP_calibDecomp*GP_DOM_RCDECOMP*DOM_QUALITY_CF[cellLoc]*DOM_TEMP_CF[cellLoc]*DEPOS_ORG_MAT[cellLoc]
00527          *(Min(DOM_SED_AEROB_Z[cellLoc]/HP_DOM_MAXDEPTH[HAB[cellLoc]],1.0)*soil_MOIST_CF[cellLoc]
00528          +GP_DOM_DECOMPRED*Min(DOM_SED_ANAEROB_Z[cellLoc]/HP_DOM_MAXDEPTH[HAB[cellLoc]],1.0) );
00529      DOM_DECOMP[cellLoc]  =  
00530          ( (DOM_DECOMP_POT[cellLoc])*DT>DEPOS_ORG_MAT[cellLoc] ) ? 
00531          ( (DEPOS_ORG_MAT[cellLoc])/DT ) : 
00532          ( DOM_DECOMP_POT[cellLoc]);
00533 /* calc state var (kg OM/m2) */
00534      DEPOS_ORG_MAT[cellLoc] =  DEPOS_ORG_MAT[cellLoc] + 
00535          ( DOM_fr_nphBio[cellLoc] + DOM_FR_FLOC[cellLoc] 
00536          - DOM_DECOMP[cellLoc] ) * DT;
00537 
00538 /* soil elevation */
00539      DOM_Z[cellLoc] = (double) DEPOS_ORG_MAT[cellLoc] / DOM_BD[cellLoc] ; /* (m) organic depth  */
00540      SED_ELEV[cellLoc]  = DOM_Z[cellLoc]+Inorg_Z[cellLoc]+SED_INACT_Z[cellLoc];  /* total land surface elevation, including model GP_DATUM_DISTANCE below NGVD 1929 (m) */
00541 
00542 /* P DOM stoich (kg P /m2) */
00543      DOP_nphBio[cellLoc] = nphbio_mort_P[cellLoc];    
00544      DOP_FLOC[cellLoc] = FlocP_DEPO[cellLoc]; 
00545 
00546      DOP_DECOMP[cellLoc] = (double) DOM_DECOMP[cellLoc] * DOM_P_OM[cellLoc]; 
00547 
00548 /* calc state var of total "organic" P in soil (NOT including dissolved in pore water or sorbed) (kgP/m2) */
00549      DOP[cellLoc] =  DOP[cellLoc]  
00550          + ( DOP_nphBio[cellLoc]  + DOP_FLOC[cellLoc] 
00551          - DOP_DECOMP[cellLoc]) * DT; /* kgP/m2 */
00552 
00553 /* now the P ratio */
00554      DOM_P_OM[cellLoc] = (DEPOS_ORG_MAT[cellLoc]>0.0) ? ( DOP[cellLoc] / DEPOS_ORG_MAT[cellLoc]) : (0.0); /* kgP/kgOM */
00555      TPsoil = DOP[cellLoc]*CELL_SIZE + TP_SORB[cellLoc]; /* kg TP in soil */
00556      TPtoSOIL[cellLoc] = ((DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc])>0.0) ?
00557          (  TPsoil / (DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc]) ) : (0.0); /* kgP/kgsoil */
00558      TPtoVOL[cellLoc] =  (CELL_SIZE * DOM_Z[cellLoc]>0.0) ?
00559          (TPsoil / (CELL_SIZE * DOM_Z[cellLoc]) ) :
00560          (0.0); /* kgP/m3 soil */
00561          
00562      TPtoSOIL_rep[cellLoc] = TPtoSOIL[cellLoc] * conv_kgTOmg; /* reporting purposes only (kg/kg->mg/kg) */
00563      TPtoVOL_rep[cellLoc] = TPtoVOL[cellLoc] * conv_kgTOg; /* reporting purposes only (kg/m3->g/m3 == ug/cm3) */
00564 
00565 /* now the P gain in sed water with decomp
00566  a small proportion goes into surface water P (below) */
00567      TP_sedMin  =  (1.0 - HP_DOM_AEROBTHIN[HAB[cellLoc]]  / HP_DOM_MAXDEPTH[HAB[cellLoc]] )
00568           * DOP_DECOMP[cellLoc] * CELL_SIZE;
00569                  
00570 /* calc P in sed water state variables (kg P) */
00571      TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] + TP_sedMin * DT;
00572              /* this is the active zone, where uptake, sorption, and mineralization take place */
00573      TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] + TP_sedMin * DT;
00574 
00575      TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
00576           (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
00577           (0.0);
00578      TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
00579           ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
00580           (TP_SED_CONC[cellLoc]); /* g/L */
00581      TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m3==mg/L */
00582               
00583 /* now store the ratio of the conc in the active zone relative to total, prior to horiz fluxes
00584 ***** very simple constant, code in transition **** */              
00585      TP_Act_to_Tot[cellLoc] = 1.0 / HP_TP_CONC_GRAD[HAB[cellLoc]];
00586      TP_Act_to_TotRep[cellLoc] =  (float) TP_Act_to_Tot[cellLoc];
00587                  
00588 /* now the P gain in surface water with decomp in the very thin upper layer of the soil */
00589                      /* if there is no surface water present, assume that this
00590                         relative contribution will be an additional sorbed component that
00591                         is introduced to surface water column immediately upon hydration
00592                         with surface water */
00593      TP_sfMin  =  HP_DOM_AEROBTHIN[HAB[cellLoc]] / HP_DOM_MAXDEPTH[HAB[cellLoc]]
00594            * DOP_DECOMP[cellLoc] * CELL_SIZE;
00595                  
00596 /* calc P in surface water state variable (kg P) */
00597      TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] + TP_sfMin * DT;
00598      TP_SFWT_CONC[cellLoc]  = 
00599            ( SFWT_VOL[cellLoc] > 0.0 ) ? 
00600            ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
00601            ( 0.0); /* used in P fluxes for mass balance */
00602      TP_SFWT_CONC_MG[cellLoc]  = 
00603            ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
00604            (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
00605            (0.0); /* (g/m3==mg/L) used for reporting and other modules to evaluate P conc when water is present */
00606 
00607 /* for reporting only: calc sum of all P storages in grid cells (budget calcs do same for Basin/Indicator-Regions) (g P /m^2) */
00608      P_SUM_CELL[cellLoc] = ( (C_ALG_P[cellLoc] + NC_ALG_P[cellLoc]) * 0.001 * CELL_SIZE + /* gP/m2 => kgP */
00609            (mac_nph_P[cellLoc] + mac_ph_P[cellLoc] )* CELL_SIZE + /* kgP/m2 => kgP */
00610            TP_SORB[cellLoc] +
00611            ( FlocP[cellLoc] + DOP[cellLoc]  ) * CELL_SIZE + /* kgP/m2 => kgP */
00612            TP_SED_WT[cellLoc] + TP_SF_WT[cellLoc] ) /* kgP */
00613            /CELL_SIZE * conv_kgTOg; /* kg P/m^2 => g P/m^2 */
00614                  
00615     }
00616   }
00617   }
00618 }
00619 
00620 
00621 /*******/
00628 void horizFlow(int step)
00629 {
00630   int it;
00631   float celWid_to_ref;
00632 
00633   celWid_to_ref = 1.0 - GP_dispLenRef/celWid; /* relative diff between the current model cell width and a reference grid cell that has numerical dispersion == actual dispersion */
00634   dispParm_scaled = celWid_to_ref * GP_dispParm;
00635 
00636   for ( it = 0; it < hyd_iter; it++ )
00637   {
00638     sprintf(msgStr,"\b\b\b%3d",it); 
00639     usrErr0(msgStr);
00640        
00641     /* horizontal raster-vector canal fluxes and water management functions
00642        Water Management switch set at runtime in Driver.parm
00643        this routine also integrates surface, unsaturated, and saturated exchanges
00644      */
00645     /* Nitrogen (DINdummy) argument is dummy argument placeholder */
00646         if (WatMgmtOn  ) {
00647       Run_Canal_Network(SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N,SAT_WATER,HP_HYD_POROSITY,
00648                         HYD_RCCONDUCT, DINdummy,TP_SF_WT,SALT_SURF_WT,DINdummy,TP_SED_WT,SALT_SED_WT,
00649                         UNSAT_WATER,HP_HYD_SPEC_YIELD ); 
00650     }
00651     
00652 
00653 
00654     /* Function for horiz fluxing of surface water, no exchange with sat/unsat water */
00655     /* if the order of the solute is changed, be sure to change the mass bal info in FluxSWstuff fcn */
00656     /* Nitrogen (DINdummy) argument is dummy argument placeholder */
00657     Flux_SWater(it,SURFACE_WAT,SED_ELEV,HYD_MANNINGS_N,SALT_SURF_WT,DINdummy,TP_SF_WT);       
00658   
00659     /* Function for horiz fluxing of ground water and its vertical itegration with unsat and surface water 
00660        It is only called every other hyd_iter iteration, and passes the realized number of gwat iterations to the function.  
00661        If the order of the solutes is changed, be sure to change the mass bal info in FluxGWstuff fcn 
00662      */
00663     /* Nitrogen (DINdummy) argument is dummy argument placeholder */
00664     if ( it%2 ) 
00665       Flux_GWater((it-1)/2, SAT_WATER, UNSAT_WATER, SURFACE_WAT, HYD_RCCONDUCT, HP_HYD_POROSITY,
00666                    HP_HYD_SPEC_YIELD, SED_ELEV, SALT_SED_WT, DINdummy, TP_SED_WT, SALT_SURF_WT, DINdummy, TP_SF_WT);  
00667 
00668   }  /* end of hyd_iter */  
00669 
00670 }
00671 
00672 
00673 
00674 /*******/
00684 /*
00685   SURFACE_WAT[cellLoc] == water storage above the land surface elevation (meters)
00686   UNSAT_WATER[cellLoc] == water storage in the pore spaces of the unsaturated zone (meters)
00687   SAT_WATER[cellLoc] == water storage in the pore spaces of the saturated zone (meters)
00688   
00689 
00690 Parameter definitions:
00691       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
00692       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
00693 ******/
00694 void cell_dyn7(int step)
00695  {
00696     int ix, iy, cellLoc; 
00697     float SatWat_Root_CF, field_cap; 
00698     float mann_height, N_density, f_LAI_eff, sfwat_pr1;
00699     float cloudy=0.0;
00700     
00701 /* the horizontal raster and vector fluxes are always called before this cell_dyn */
00702     for(ix=1; ix<=s0; ix++) {
00703         for(iy=1; iy<=s1; iy++) {
00704             if(ON_MAP[cellLoc= T(ix,iy)])  {
00705 
00706               if (debug > 3) { /* these are old, relatively coarse checks - surface-ground water integration occurs in Fluxes.c */
00707                   if (SAT_WT_HEAD[cellLoc] - 0.01 > SED_ELEV[cellLoc] ) {
00708                       sprintf(msgStr,"Day %6.1f. Warning - SAT_WT_HEAD exceeds elev in cell (%d,%d) by %f.", 
00709                               SimTime.TIME, ix,iy,(SAT_WT_HEAD[cellLoc] - SED_ELEV[cellLoc]) ); 
00710                       WriteMsg( msgStr,True ); }
00711                   if (SURFACE_WAT[cellLoc] > 0.2 && UNSAT_DEPTH[cellLoc] > 0.1 ) {
00712                       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).", 
00713                               SimTime.TIME, SURFACE_WAT[cellLoc], UNSAT_DEPTH[cellLoc],UNSAT_MOIST_PRP[cellLoc], ix,iy ); 
00714                       WriteMsg( msgStr,True );  }
00715                   if (SAT_WATER[cellLoc] < -0.01) { /* this seems unnecessary but... */
00716                       sprintf(msgStr,"Day %6.1f: capacityERR - neg SAT_WATER (%f m) in cell (%d,%d) before cell_dyn7!", 
00717                               SimTime.TIME, SAT_WATER[cellLoc], ix,iy ); 
00718                       WriteMsg( msgStr,True );  dynERRORnum++; }
00719                   if (SURFACE_WAT[cellLoc] < -0.01) {
00720                       sprintf(msgStr,"Day %6.1f: capacityERR - neg SURFACE_WAT (%f m) in cell (%d,%d) before cell_dyn7!", 
00721                               SimTime.TIME, SURFACE_WAT[cellLoc], ix,iy ); 
00722                       WriteMsg( msgStr,True );  dynERRORnum++; }
00723               }
00724 
00725 /*  note that rainfall during a time step is added to surface water storage and available */
00726 /* for runoff (horizFlow) before the calc of infiltration & ET associated with that new input */
00727 /* (infiltration/ET etc will be of avail water the next time step after a rainfall event and horiz flows) */
00728                   /* NSM/SFWMM rainfall input data, created in cell_dyn1, convert here from tenths of mm to meters */
00729               SF_WT_FROM_RAIN[cellLoc] = boundcond_rain[cellLoc]*0.0001;  /*  tenths of mm *0.0001 = m */
00730 
00731               /* solar radiation at altitude of 274m in atmosphere cal/cm2/d) */
00732               /* v2.2+, CLOUDY (cloudiness) spatiotemporal data "temporarily" unavailable, is constant in space and time, local var "cloudy" */
00733               SOLRAD274[cellLoc] = SOLRADATMOS*(SOLBETA-GP_SOLOMEGA* ( ( cloudy>0.0 ) ? ( cloudy ) : ( 0.0) ) ) -SOLALPHA;
00734               SOLRADGRD[cellLoc] = SOLRAD274[cellLoc]+((SOLRADATMOS+1.0)-SOLRAD274[cellLoc])*SOLALTCORR;
00735               H2O_TEMP[cellLoc] = AIR_TEMP[cellLoc]; /* v2.2+, temperature data "temporarily" unavailable, is constant in space and time */
00736               
00737 /******** determine new unsat potentials */
00738               SAT_WT_HEAD[cellLoc]  = SAT_WATER[cellLoc]/HP_HYD_POROSITY[HAB[cellLoc]];
00739               UNSAT_DEPTH[cellLoc]  = SED_ELEV[cellLoc]-SAT_WT_HEAD[cellLoc];
00740               UNSAT_CAP[cellLoc]  =  UNSAT_DEPTH[cellLoc]*HP_HYD_POROSITY[HAB[cellLoc]];
00741           
00742               UNSAT_MOIST_PRP[cellLoc]  = 
00743                   ( UNSAT_CAP[cellLoc]>0.0 ) ? 
00744                   ( Min(UNSAT_WATER[cellLoc]/UNSAT_CAP[cellLoc],1.0) ) : 
00745                   ( 1.0); 
00746                   /* determining the pathway of flow of surface water depending on depth
00747                      of an unsat zone relative to the surface water  */ 
00748               SAT_VS_UNSAT[cellLoc]  = 1/Exp(100.0*Max((SURFACE_WAT[cellLoc]-UNSAT_DEPTH[cellLoc]),0.0)); 
00749      /* 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) */
00750               UNSAT_HYD_COND_CF[cellLoc]  = graph7(0x0,UNSAT_MOIST_PRP[cellLoc] ); 
00751                      /* field capacity = porosity - specific yield; spec yield== proportion of total soil vol
00752                             that represents water that can be moved by gravity */
00753               field_cap = (HP_HYD_POROSITY[HAB[cellLoc]]-HP_HYD_SPEC_YIELD[HAB[cellLoc]]);
00754                   /* unsat_avail is proportion of water in pore space available for gravitational flow (above field capacity) */
00755                   /* e.g., when moisture prop in pore space <= field_cap/pore_space, no percolation */
00756                   /* using old moisture proportion (hasn't changed unless unsat zone was replaced by sat water) */
00757               UNSAT_AVAIL[cellLoc]  = Max(UNSAT_MOIST_PRP[cellLoc]
00758                                           -(field_cap)/HP_HYD_POROSITY[HAB[cellLoc]],0.0);
00759               UNSAT_WT_POT[cellLoc]  = Max(UNSAT_CAP[cellLoc]-UNSAT_WATER[cellLoc],0.0);
00760 
00761 /******** now determine the potential total transpiration and evaporation  */
00762 /* Potential ET is input data used in SFWMM v5.4 */
00763 /* GP_calibET is an adjustable calibration parameter (close to 1.0, adjusted in global parameter input file)  */
00764              HYD_EVAP_CALC[cellLoc]  =  boundcond_ETp[cellLoc] * 0.0001*GP_calibET;  /*  tenths of mm *0.0001 = m */
00765 
00766     /*  Leaf Area Index (LAI) of emergent macrophytes: this effective LAI estimates leaf area index that is above ponded surface water  */
00767               LAI_eff[cellLoc] =  
00768                   (MAC_HEIGHT[cellLoc]>0.0) ? 
00769                   (Max(1.0 - SURFACE_WAT[cellLoc]/MAC_HEIGHT[cellLoc], 0.0)*MAC_LAI[cellLoc]) : 
00770                   (0.0)  ;
00771  
00772           /* control function (0-1) of relative magnitude of the effective Leaf Area Index  */
00773               f_LAI_eff = exp(-LAI_eff[cellLoc]); 
00774               
00775               
00776               HYD_TOT_POT_TRANSP[cellLoc]  = HYD_EVAP_CALC[cellLoc] * (1.0-f_LAI_eff); 
00777 
00778      /* 0-1 control function of saturated water available to roots - capillary draw when roots don't quite reach down to water table */
00779               SatWat_Root_CF =  Exp(-10.0* Max(UNSAT_DEPTH[cellLoc]-HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]],0.0) ); 
00780      /*  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 */
00781               HYD_WATER_AVAIL[cellLoc]  = (UNSAT_DEPTH[cellLoc] > HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]] ) ? 
00782                   ( Max(UNSAT_MOIST_PRP[cellLoc], SatWat_Root_CF) ) :
00783                   ( 1.0 ); 
00784     /* 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) */
00785               MAC_WATER_AVAIL_CF[cellLoc]  = graph8(0x0,HYD_WATER_AVAIL[cellLoc]); 
00786 
00787 /******** next calc the potential and actual flows */
00788 /* unsatLoss(1) unsat to sat percolation */
00789   /*unsat to sat flow here only includes percolation (rising water table accomodated in update after horiz fluxes) */ 
00790               UNSAT_PERC[cellLoc]  = Min(HP_HYD_RCINFILT[HAB[cellLoc]]*UNSAT_HYD_COND_CF[cellLoc],UNSAT_AVAIL[cellLoc]*UNSAT_WATER[cellLoc]);
00791               UNSAT_TO_SAT_FL[cellLoc]  =  
00792                   ( (UNSAT_PERC[cellLoc])*DT > UNSAT_WATER[cellLoc] ) ? 
00793                   ( UNSAT_WATER[cellLoc]/DT ) : 
00794                   ( UNSAT_PERC[cellLoc]);
00795 /* unsatLoss(2) unsat to atmosph transpiration */
00796               HYD_UNSAT_POT_TRANS[cellLoc]  = (UNSAT_DEPTH[cellLoc] > HP_NPHBIO_ROOTDEPTH[HAB[cellLoc]] ) ?
00797                   ( HYD_TOT_POT_TRANSP[cellLoc]*MAC_WATER_AVAIL_CF[cellLoc] ) :
00798                   (0.0); /* no unsat transp if roots are in saturated zone */
00799               UNSAT_TRANSP[cellLoc]  = 
00800                   ( (HYD_UNSAT_POT_TRANS[cellLoc]+UNSAT_TO_SAT_FL[cellLoc])*DT>UNSAT_WATER[cellLoc] ) ? 
00801                   ( (UNSAT_WATER[cellLoc]-UNSAT_TO_SAT_FL[cellLoc]*DT)/DT ) : 
00802                   ( HYD_UNSAT_POT_TRANS[cellLoc]);
00803 
00804 /* satLoss(1) sat to deep aquifer recharge **RATE parameter IS ALWAYS SET to ZERO  *****/
00805               SAT_WT_RECHG[cellLoc]  = 
00806                   ( GP_HYD_RCRECHG*HP_HYD_SPEC_YIELD[HAB[cellLoc]]/HP_HYD_POROSITY[HAB[cellLoc]]*DT>SAT_WATER[cellLoc] ) ? 
00807                   ( SAT_WATER[cellLoc]/DT ) : 
00808                   ( GP_HYD_RCRECHG*HP_HYD_SPEC_YIELD[HAB[cellLoc]]/HP_HYD_POROSITY[HAB[cellLoc]]); 
00809                  
00810 /* sat to surf upflow  when gwater exceeds soil capacity due to lateral inflows
00811    accomodated in gwFluxes */
00812 
00813 /* satLoss(2) sat to unsat with lowering water table due to recharge loss ONLY (horiz outflows accomodated in gwFluxes)
00814    (leaves field capacity amount in unsat zone)*/
00815               SAT_TO_UNSAT_FL[cellLoc]  =  
00816                   ( SAT_WT_RECHG[cellLoc]*field_cap/HP_HYD_POROSITY[HAB[cellLoc]]*DT > SAT_WATER[cellLoc] ) ? 
00817                   ( (SAT_WATER[cellLoc])/DT ) : 
00818                   ( SAT_WT_RECHG[cellLoc]*field_cap/HP_HYD_POROSITY[HAB[cellLoc]]) ;
00819 /* satLoss(3) sat to atmosph */
00820               HYD_SAT_POT_TRANS[cellLoc]  = HYD_TOT_POT_TRANSP[cellLoc]*SatWat_Root_CF; 
00821               SAT_WT_TRANSP[cellLoc]  =  
00822                   ( (HYD_SAT_POT_TRANS[cellLoc]+SAT_TO_UNSAT_FL[cellLoc])*DT > SAT_WATER[cellLoc] ) ? 
00823                   ( (SAT_WATER[cellLoc]-SAT_TO_UNSAT_FL[cellLoc]*DT)/DT ) : 
00824                   ( HYD_SAT_POT_TRANS[cellLoc]);
00825 
00826 /* sfwatLoss(1) sf to sat */
00827                      /* downflow to saturated assumed to occur in situations with small
00828                         unsat layer overlain by significant surface water (SAT_VS_UNSAT very small),
00829                         with immediate hydraulic connectivity betweent the two storages */
00830               SF_WT_TO_SAT_DOWNFLOW[cellLoc]  = 
00831                   ( (1.0-SAT_VS_UNSAT[cellLoc])*UNSAT_WT_POT[cellLoc]*DT>SURFACE_WAT[cellLoc] ) ? 
00832                   ( SURFACE_WAT[cellLoc]/DT ) : 
00833                   ( (1.0-SAT_VS_UNSAT[cellLoc])*UNSAT_WT_POT[cellLoc]); 
00834 /* sfwatLoss(2) sf to unsat infiltration (sat_vs_unsat splits these losses to groundwater, but downflow to sat is given priority) */
00835               SF_WT_POT_INF[cellLoc]  = 
00836                   ( (SAT_VS_UNSAT[cellLoc]*HP_HYD_RCINFILT[HAB[cellLoc]]+SF_WT_TO_SAT_DOWNFLOW[cellLoc])*DT>SURFACE_WAT[cellLoc] ) ? 
00837                   ( (SURFACE_WAT[cellLoc]-SF_WT_TO_SAT_DOWNFLOW[cellLoc]*DT)/DT ) : 
00838                   ( SAT_VS_UNSAT[cellLoc]*HP_HYD_RCINFILT[HAB[cellLoc]]);
00839               SF_WT_INFILTRATION[cellLoc]  = 
00840                   ( SF_WT_POT_INF[cellLoc]*DT > (UNSAT_WT_POT[cellLoc]-SF_WT_TO_SAT_DOWNFLOW[cellLoc]*DT) ) ? 
00841                   ((UNSAT_WT_POT[cellLoc]-SF_WT_TO_SAT_DOWNFLOW[cellLoc]*DT)/DT ) : 
00842                   ( SF_WT_POT_INF[cellLoc]);
00843               sfwat_pr1 = SF_WT_INFILTRATION[cellLoc]+SF_WT_TO_SAT_DOWNFLOW[cellLoc];
00844 /* sfwatLoss(3) sf to atmosph */
00845               SF_WT_EVAP[cellLoc]  =  
00846                   ( (f_LAI_eff*HYD_EVAP_CALC[cellLoc]+sfwat_pr1 )*DT>SURFACE_WAT[cellLoc] ) ? 
00847                   ( (SURFACE_WAT[cellLoc]-sfwat_pr1*DT)/DT ) : 
00848                   ( f_LAI_eff*HYD_EVAP_CALC[cellLoc]); 
00849 
00850 
00851 /******** then update the state variable storages */
00852 
00853 /* unsat loss priority:  percolation,  transpiration */
00854 /* calc unsaturated storage state var (m) */
00855               UNSAT_WATER[cellLoc] = UNSAT_WATER[cellLoc] 
00856                   + (SAT_TO_UNSAT_FL[cellLoc] + SF_WT_INFILTRATION[cellLoc] 
00857                      - UNSAT_TO_SAT_FL[cellLoc] - UNSAT_TRANSP[cellLoc]) * DT;
00858 
00859 /* sat loss priority:  recharge to deep aquifer, re-assign to unsat with lowered water table, transpiration */
00860 /* calc saturated storage state var (m) */
00861               SAT_WATER[cellLoc] =  SAT_WATER[cellLoc] 
00862                   + (UNSAT_TO_SAT_FL[cellLoc] + SF_WT_TO_SAT_DOWNFLOW[cellLoc] 
00863                      - SAT_WT_TRANSP[cellLoc] - SAT_TO_UNSAT_FL[cellLoc] - SAT_WT_RECHG[cellLoc]) * DT;
00864 
00865 /* sfwat loss priority: downflow to replace groundwater loss, infiltration to unsat, evaporation */
00866 /* calc surface storage state var (m) */
00867               SURFACE_WAT[cellLoc] = SURFACE_WAT[cellLoc] 
00868                   + (SF_WT_FROM_RAIN[cellLoc]  
00869                      - SF_WT_EVAP[cellLoc] - SF_WT_INFILTRATION[cellLoc] - SF_WT_TO_SAT_DOWNFLOW[cellLoc]) * DT;
00870 
00871 /******** lastly, update of head calcs, unsat capacity, moisture proportion, etc.
00872  in order to calc water in diff storages for solute concentrations */
00873               SAT_WT_HEAD[cellLoc]  = SAT_WATER[cellLoc]/HP_HYD_POROSITY[HAB[cellLoc]];
00874               UNSAT_DEPTH[cellLoc]  = Max(SED_ELEV[cellLoc]-SAT_WT_HEAD[cellLoc],0.0);
00875               UNSAT_CAP[cellLoc]  =  UNSAT_DEPTH[cellLoc]*HP_HYD_POROSITY[HAB[cellLoc]];
00876 
00877               UNSAT_MOIST_PRP[cellLoc]  = 
00878                   ( UNSAT_CAP[cellLoc]>0.0 ) ? 
00879                   ( Min(UNSAT_WATER[cellLoc]/UNSAT_CAP[cellLoc],1.0) ) : 
00880                   ( 1.0); 
00881               HYD_DOM_ACTWAT_VOL[cellLoc]  = ( Min(HP_DOM_MAXDEPTH[HAB[cellLoc]],UNSAT_DEPTH[cellLoc])*UNSAT_MOIST_PRP[cellLoc] +
00882                                                Max(HP_DOM_MAXDEPTH[HAB[cellLoc]]-UNSAT_DEPTH[cellLoc], 0.0)*HP_HYD_POROSITY[HAB[cellLoc]] )
00883                   *CELL_SIZE; 
00884                   /* flag for presence of small amount of water storage in the active zone must be present */ 
00885               HYD_DOM_ACTWAT_PRES[cellLoc]  = 
00886                   ( HYD_DOM_ACTWAT_VOL[cellLoc] > CELL_SIZE*0.01 ) ? 
00887                   ( 1.0 ) : ( 0.0); 
00888               HYD_SED_WAT_VOL[cellLoc]  = (SAT_WATER[cellLoc]+UNSAT_WATER[cellLoc])*CELL_SIZE;
00889               SFWT_VOL[cellLoc]  = SURFACE_WAT[cellLoc]*CELL_SIZE;
00890 
00891               HydTotHd[cellLoc]  = SAT_WT_HEAD[cellLoc]+SURFACE_WAT[cellLoc]; /* only used for reporting purposes */
00892 
00893                   /* at the square of xx% of the macrophyte's height, the manning's n
00894                      calc will indicate the macrophyte *starting* to bend over,
00895                      starting to offer increasingly less resistance */
00896               mann_height = (GP_mann_height_coef*MAC_HEIGHT[cellLoc])*(GP_mann_height_coef*MAC_HEIGHT[cellLoc]); 
00897               N_density = Max(HP_MAC_MAXROUGH[HAB[cellLoc]] * MAC_REL_BIOM[cellLoc], HP_MAC_MINROUGH[HAB[cellLoc]] );
00898                   /* manning's n for overland (horiz) flow */
00899               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 */   
00900               HYD_MANNINGS_N[cellLoc]  = Max(-Abs((N_density-HP_MAC_MINROUGH[HAB[cellLoc]])
00901                                                   *(pow(2.0,(1.0-SURFACE_WAT[cellLoc]/mann_height))-1.0) ) 
00902                                              + N_density,HP_MAC_MINROUGH[HAB[cellLoc]]);
00903 
00904                   /* sum of transpiration for output only */
00905               HYD_TRANSP[cellLoc]  = UNSAT_TRANSP[cellLoc]+SAT_WT_TRANSP[cellLoc];
00906               HYD_ET[cellLoc]  = HYD_TRANSP[cellLoc]+SF_WT_EVAP[cellLoc];
00907 
00908             }
00909   }
00910     }
00911  }
00912 
00913 
00914 /*******/
00920 /*   MAC_NOPH_BIOMAS[cellLoc] == carbon mass of live Non-Photosynthetic tissues of macrophytes (kgC/m^2)
00921    MAC_PH_BIOMAS[cellLoc] == carbon mass of live Photosynthetic tissues of macrophytes (kgC/m^2)
00922 
00923    mac_nph_P[cellLoc] == phosphorus mass of live Non-Photosynthetic tissues of macrophytes (kgP/m^2)
00924    mac_ph_P[cellLoc] == phosphorus mass of live Photosynthetic tissues of macrophytes (kgP/m^2)
00925 
00926    mac_nph_OM[cellLoc] == organic matter mass of live Non-Photosynthetic tissues of macrophytes (kgOM/m^2)
00927    mac_ph_OM[cellLoc] == organic matter mass of live Photosynthetic tissues of macrophytes (kgOM/m^2)
00928 
00929 spatial (cell) location defines habitat ( arrayOf[HAB[cellLoc]] );
00930 habitat type can switch when global dynamic function (in cell_dyn1) calls the succession function
00931 
00932 Parameter definitions:
00933       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
00934       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
00935 */
00936 void cell_dyn8(int step)
00937  {
00938 int ix, iy, cellLoc; 
00939 float reduc, NPP_P, min_litTemp, nphbio_ddep, phbio_ddep, MAC_PHtoNPH, MAC_PHtoNPH_Init;
00940 
00941   for(ix=1; ix<=s0; ix++) {
00942   for(iy=1; iy<=s1; iy++) {
00943 
00944     if(ON_MAP[cellLoc= T(ix,iy)])  { 
00945               
00946 /* these thresholds need updating when a habitat type of a grid cell changes */
00947       MAC_MAX_BIO[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]+HP_PHBIO_MAX[HAB[cellLoc]];
00948       NPHBIO_REFUGE[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
00949       NPHBIO_SAT[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]*0.9;
00950       PHBIO_REFUGE[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
00951       PHBIO_SAT[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]]*0.9;
00952 /* various control functions for production calc */
00953      MAC_LIGHT_CF[cellLoc]  = SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]
00954          *Exp(1.0-SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]);
00955      MAC_TEMP_CF[cellLoc]  = tempCF(0, 0.20, HP_MAC_TEMPOPT[HAB[cellLoc]], 5.0, 40.0, AIR_TEMP[cellLoc]);
00956      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 */
00957      MAC_WATER_CF[cellLoc]  = Min(MAC_WATER_AVAIL_CF[cellLoc], 
00958          Max(1.0-Max( (SURFACE_WAT[cellLoc]-HP_MAC_WAT_TOLER[HAB[cellLoc]])
00959          /HP_MAC_WAT_TOLER[HAB[cellLoc]],0.0),0.0));
00960      MAC_NUT_CF[cellLoc]  = 
00961                      Exp(-GP_mac_uptake_coef * Max(HP_MAC_KSP[HAB[cellLoc]]-TP_SEDWT_CONCACTMG[cellLoc], 0.0)
00962                          /HP_MAC_KSP[HAB[cellLoc]]) ; /* mg/L */
00963      MAC_SALT_CF[cellLoc]  = 1.0; /* made inactive here; can set up relation to HP_MAC_SALIN_THRESH[HAB[cellLoc]] parm */
00964      min_litTemp = Min(MAC_LIGHT_CF[cellLoc], MAC_TEMP_CF[cellLoc]);
00965      MAC_PROD_CF[cellLoc]  = Min(min_litTemp,MAC_WATER_CF[cellLoc])
00966           *MAC_NUT_CF[cellLoc] *MAC_SALT_CF[cellLoc];
00967 /* net primary production, kg C/m2/d */
00968      PHBIO_NPP[cellLoc]  = HP_PHBIO_RCNPP[HAB[cellLoc]]*MAC_PROD_CF[cellLoc]*MAC_PH_BIOMAS[cellLoc]
00969          * Max( (1.0-MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc]), 0.0);
00970 /* 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 ) */
00971      NPP_P = PHBIO_NPP[cellLoc]  * HP_PHBIO_IC_PC[HAB[cellLoc]]  * Max(MAC_NUT_CF[cellLoc]*2.0,1.0)
00972                      * Max(1.0-mac_ph_PC[cellLoc]/HP_PHBIO_IC_PC[HAB[cellLoc]], 0.0);
00973 /* 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 */
00974      reduc = (NPP_P > 0.0) ? 
00975                      (Max(TP_SED_WT_AZ[cellLoc],0.0) / ( NPP_P*CELL_SIZE*DT) ) : /* within-plant variable stoichiometry */
00976                      (1.0);
00977     /* reduce the production proportionally if excess demand is found */
00978                 /* can have high conc, but low mass of P avail, in presence of high plant biomass and high demand */
00979      if (reduc < 1.0) {
00980                      PHBIO_NPP[cellLoc] *= reduc;
00981                      NPP_P *= reduc;
00982                  }
00983                  
00984 /* losses from photobio */
00985      phbio_ddep = Max(1.0-Max( (PHBIO_SAT[cellLoc]-MAC_PH_BIOMAS[cellLoc]) 
00986                                            /(PHBIO_SAT[cellLoc]-PHBIO_REFUGE[cellLoc]),0.0),0.0);
00987      PHBIO_AVAIL[cellLoc]  = MAC_PH_BIOMAS[cellLoc]*phbio_ddep;
00988      MAC_PHtoNPH_Init = HP_PHBIO_MAX[HAB[cellLoc]] / HP_NPHBIO_MAX[HAB[cellLoc]] ; /*if habitat type changes */
00989      MAC_PHtoNPH = (MAC_NOPH_BIOMAS[cellLoc]>0.0) ? ( MAC_PH_BIOMAS[cellLoc] / MAC_NOPH_BIOMAS[cellLoc]) : (MAC_PHtoNPH_Init);
00990                  
00991      NPHBIO_TRANSLOC_POT[cellLoc]  = (MAC_PHtoNPH>MAC_PHtoNPH_Init) ?
00992                      (exp(HP_MAC_TRANSLOC_RC[HAB[cellLoc]] *(MAC_PHtoNPH-MAC_PHtoNPH_Init)) - 1.0) :
00993                      (0.0); 
00994      NPHBIO_TRANSLOC[cellLoc]  =  
00995          ( (NPHBIO_TRANSLOC_POT[cellLoc])*DT >PHBIO_AVAIL[cellLoc] ) ? 
00996          ( (PHBIO_AVAIL[cellLoc])/DT ) : 
00997          ( NPHBIO_TRANSLOC_POT[cellLoc]);
00998 
00999      PHBIO_MORT_POT[cellLoc]  = HP_PHBIO_RCMORT[HAB[cellLoc]] *PHBIO_AVAIL[cellLoc] 
01000          *(1.0 + (1.0-MAC_WATER_AVAIL_CF[cellLoc]) )/2.0;
01001      PHBIO_MORT[cellLoc]  =
01002                                 ( (PHBIO_MORT_POT[cellLoc]+NPHBIO_TRANSLOC[cellLoc])*DT>PHBIO_AVAIL[cellLoc] ) ? 
01003          ( (PHBIO_AVAIL[cellLoc]-NPHBIO_TRANSLOC[cellLoc]*DT)/DT ) : 
01004          ( PHBIO_MORT_POT[cellLoc]);
01005 
01006 
01007 /* losses from non-photobio  */
01008      nphbio_ddep = Max(1.0-Max((NPHBIO_SAT[cellLoc]-MAC_NOPH_BIOMAS[cellLoc])
01009                                           /(NPHBIO_SAT[cellLoc]-NPHBIO_REFUGE[cellLoc]),0.0),0.0);
01010      NPHBIO_AVAIL[cellLoc]  = MAC_NOPH_BIOMAS[cellLoc]*nphbio_ddep; 
01011 
01012                  PHBIO_TRANSLOC_POT[cellLoc]  = (MAC_PHtoNPH<MAC_PHtoNPH_Init) ?
01013                      (exp(HP_MAC_TRANSLOC_RC[HAB[cellLoc]] *(MAC_PHtoNPH_Init-MAC_PHtoNPH)) - 1.0) :
01014                      (0.0); 
01015      PHBIO_TRANSLOC[cellLoc]  =  
01016          ( (PHBIO_TRANSLOC_POT[cellLoc])*DT >NPHBIO_AVAIL[cellLoc] ) ? 
01017          ( (NPHBIO_AVAIL[cellLoc])/DT ) : 
01018          ( PHBIO_TRANSLOC_POT[cellLoc]);
01019      NPHBIO_MORT_POT[cellLoc]  = NPHBIO_AVAIL[cellLoc]*HP_PHBIO_RCMORT[HAB[cellLoc]]
01020                      * (1.0 + Max(1.0-MAC_PH_BIOMAS[cellLoc]/HP_PHBIO_MAX[HAB[cellLoc]],0.0) )/2.0; /* decreased mort w/ increasing photobio */
01021      NPHBIO_MORT[cellLoc]  =
01022                                 ( (PHBIO_TRANSLOC[cellLoc]+NPHBIO_MORT_POT[cellLoc])*DT>NPHBIO_AVAIL[cellLoc] ) ? 
01023          ( (NPHBIO_AVAIL[cellLoc]-PHBIO_TRANSLOC[cellLoc]*DT)/DT ) : 
01024          ( NPHBIO_MORT_POT[cellLoc]);
01025                  
01026 
01027 /* calc nonphotosynthetic biomass state var (kgC/m2) */
01028      MAC_NOPH_BIOMAS[cellLoc] =  MAC_NOPH_BIOMAS[cellLoc] 
01029                      + (NPHBIO_TRANSLOC[cellLoc] - NPHBIO_MORT[cellLoc] 
01030                         - PHBIO_TRANSLOC[cellLoc]  ) * DT;
01031                  
01032 /* calc nonphotosynthetic biomass state var (kgC/m2) */
01033      MAC_PH_BIOMAS[cellLoc] = MAC_PH_BIOMAS[cellLoc] 
01034          + (PHBIO_TRANSLOC[cellLoc] + PHBIO_NPP[cellLoc] 
01035            - PHBIO_MORT[cellLoc] - NPHBIO_TRANSLOC[cellLoc]) * DT;
01036 
01037 /* total biomass */
01038      MAC_TOT_BIOM[cellLoc]  = MAC_PH_BIOMAS[cellLoc]+MAC_NOPH_BIOMAS[cellLoc];
01039 /* book-keeping calcs used in other modules */
01040      MAC_REL_BIOM[cellLoc]  = 
01041          ( MAC_TOT_BIOM[cellLoc] > 0.0 ) ? 
01042          MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc] : 
01043          0.0001;
01044      MAC_HEIGHT[cellLoc]  = pow(MAC_REL_BIOM[cellLoc],0.33)*HP_MAC_MAXHT[HAB[cellLoc]];
01045      MAC_LAI[cellLoc]  = MAC_REL_BIOM[cellLoc]*HP_MAC_MAXLAI[HAB[cellLoc]];
01046 /* note that an "effective" LAI that accounts for water depth is calc'd in hydro module */
01047 
01048 /*  now calc the P and organic matter associated with the C fluxes */
01049                  phbio_npp_P[cellLoc] = NPP_P;     /* within-plant variable stoichiometry */
01050                  phbio_npp_OM[cellLoc] = PHBIO_NPP[cellLoc] / HP_PHBIO_IC_CTOOM[HAB[cellLoc]]; /* habitat-specfic stoichiometry */
01051 
01052                  phbio_mort_P[cellLoc] = PHBIO_MORT[cellLoc] * mac_ph_PC[cellLoc];
01053                  phbio_mort_OM[cellLoc] = PHBIO_MORT[cellLoc] / mac_ph_CtoOM[cellLoc];
01054 
01055                  phbio_transl_P[cellLoc] = PHBIO_TRANSLOC[cellLoc] * mac_nph_PC[cellLoc];
01056                  phbio_transl_OM[cellLoc] = PHBIO_TRANSLOC[cellLoc] / mac_nph_CtoOM[cellLoc];
01057 
01058                  nphbio_transl_P[cellLoc] = NPHBIO_TRANSLOC[cellLoc] * mac_ph_PC[cellLoc];
01059                  nphbio_transl_OM[cellLoc] = NPHBIO_TRANSLOC[cellLoc] / mac_ph_CtoOM[cellLoc];
01060                  
01061                  nphbio_mort_P[cellLoc] = NPHBIO_MORT[cellLoc] * mac_nph_PC[cellLoc];
01062                  nphbio_mort_OM[cellLoc] = NPHBIO_MORT[cellLoc] / mac_nph_CtoOM[cellLoc];
01063 
01064 
01065 /* state vars: now calc the P and OM associated with those C state vars */
01066                  mac_nph_P[cellLoc] = mac_nph_P[cellLoc]
01067                      + (nphbio_transl_P[cellLoc] - nphbio_mort_P[cellLoc]
01068                         - phbio_transl_P[cellLoc]  ) * DT;
01069                  mac_nph_PC[cellLoc] = ( (MAC_NOPH_BIOMAS[cellLoc] > 0.0) && (mac_nph_P[cellLoc] > 0.0)) ?
01070                      (mac_nph_P[cellLoc] / MAC_NOPH_BIOMAS[cellLoc]) : /* these second mass checks not needed */
01071                      0.3 * HP_NPHBIO_IC_PC[HAB[cellLoc]]; /* default to 0.3 of max for habitat */
01072                  mac_nph_PC_rep[cellLoc] = (float)mac_nph_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */
01073                  
01074 
01075                  mac_nph_OM[cellLoc] = mac_nph_OM[cellLoc]
01076                      + (nphbio_transl_OM[cellLoc] - nphbio_mort_OM[cellLoc]
01077                         - phbio_transl_OM[cellLoc] ) * DT;
01078                  mac_nph_CtoOM[cellLoc] = ( (mac_nph_OM[cellLoc] > 0.0) && (MAC_NOPH_BIOMAS[cellLoc]>0.0)) ?
01079                      (MAC_NOPH_BIOMAS[cellLoc] / mac_nph_OM[cellLoc]) :
01080                      HP_NPHBIO_IC_CTOOM[HAB[cellLoc]];
01081 
01082                  mac_ph_P[cellLoc] = mac_ph_P[cellLoc]
01083                      + (phbio_transl_P[cellLoc] + phbio_npp_P[cellLoc] - phbio_mort_P[cellLoc]
01084                         - nphbio_transl_P[cellLoc] ) * DT;
01085                  mac_ph_PC[cellLoc] = ( (MAC_PH_BIOMAS[cellLoc] > 0.0) && (mac_ph_P[cellLoc]>0.0)) ?
01086                      (mac_ph_P[cellLoc] / MAC_PH_BIOMAS[cellLoc]) :
01087                      0.3 * HP_PHBIO_IC_PC[HAB[cellLoc]]; /* default to 0.3 of max for habitat */
01088                  mac_ph_PC_rep[cellLoc] = (float)mac_ph_PC[cellLoc] * conv_kgTOmg; /* variable for output _rep-orting only */
01089 
01090                  mac_ph_OM[cellLoc] = mac_ph_OM[cellLoc]
01091                      + (phbio_transl_OM[cellLoc] + phbio_npp_OM[cellLoc] - phbio_mort_OM[cellLoc]
01092                         - nphbio_transl_OM[cellLoc]  ) * DT;
01093                  mac_ph_CtoOM[cellLoc] = ( (mac_ph_OM[cellLoc] > 0.0) && (MAC_PH_BIOMAS[cellLoc]>0.0)) ?
01094                      (MAC_PH_BIOMAS[cellLoc] / mac_ph_OM[cellLoc]) :
01095                      HP_PHBIO_IC_CTOOM[HAB[cellLoc]];
01096 
01097                  if (debug > 0 && MAC_NOPH_BIOMAS[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NPhoto C biomass (%f m) in cell (%d,%d)!", 
01098                               SimTime.TIME, MAC_NOPH_BIOMAS[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
01099                  if (debug > 0 && MAC_PH_BIOMAS[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg Photo C biomass (%f m) in cell (%d,%d)!", 
01100                               SimTime.TIME, MAC_PH_BIOMAS[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
01101                  if (debug > 0 && mac_nph_P[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg NPhoto P biomass (%f m) in cell (%d,%d)!", 
01102                               SimTime.TIME, mac_nph_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
01103                  if (debug > 0 && mac_ph_P[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg Photo P biomass (%f m) in cell (%d,%d)!", 
01104                               SimTime.TIME, mac_ph_P[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
01105 
01106 
01107 /******** phosphorus removed from water here *************/
01108      TP_SEDWT_UPTAKE[cellLoc]  = phbio_npp_P[cellLoc]*CELL_SIZE; 
01109 /* recalc P in sed water state variable (kg P) */
01110      TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] - (TP_SEDWT_UPTAKE[cellLoc]) * DT;
01111      TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
01112                      (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
01113                   (0.0);
01114                  
01115                       /* this is the active zone, where uptake, sorption, and mineralization take place */
01116                 TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] - (TP_SEDWT_UPTAKE[cellLoc]) * DT;
01117                  TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
01118                      ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
01119                      (TP_SED_CONC[cellLoc]); /* g/L */
01120                  TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m3==mg/L */
01121               
01122           }
01123   }
01124   }
01125 }
01126 
01127 
01128 /*******/
01135 /*
01136 TP_SF_WT[cellLoc] == mass of P in surface water (kg P)
01137 TP_SED_WT[cellLoc] == mass of P in sediment/soil (pore) water (kg P)
01138 TP_SORB[cellLoc] == mass of P sorbed to soil (kg P)
01139 
01140 Parameter definitions:
01141       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
01142       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
01143 */
01144 void cell_dyn9(int step)
01145  {
01146  int ix, iy, cellLoc; 
01147  float TPpartic, TPsettlRat, TP_settl_pot;
01148  double PO4Pconc, nonPO4Pconc;
01149  
01150   for(ix=1; ix<=s0; ix++) {
01151   for(iy=1; iy<=s1; iy++) {
01152 
01153     if(ON_MAP[cellLoc= T(ix,iy)])  {
01154 /* calc concentrations after horiz fluxes */
01155               TP_SFWT_CONC[cellLoc]  = 
01156                   ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01157                   ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01158                   ( 0.0); /* used in P fluxes for mass balance */
01159              /* using regression for predicting PO4 from TP */
01160               PO4Pconc =  Max( TP_SFWT_CONC[cellLoc]*GP_PO4toTP + 0.001*GP_PO4toTPint,0.0);  /* g/l  (note that intercept may be <0)*/
01161       /* after spatial (horiz) fluxes, recalc the active zone P mass based on old active/total ratio */
01162               TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
01163                   (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
01164                   (0.0);
01165             /* this is the active zone, where uptake, sorption, and mineralization take place */
01166              TP_SED_WT_AZ[cellLoc] = TP_SED_CONC[cellLoc] * TP_Act_to_Tot[cellLoc] * HYD_DOM_ACTWAT_VOL[cellLoc];
01167               TP_SEDWT_CONCACT[cellLoc] =(HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
01168                   ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
01169                   ( TP_SED_CONC[cellLoc]);
01170 
01171 /* inputs to surf  P (kg P) (rain conc in mg/L) */
01172               TP_FR_RAIN[cellLoc]  = SF_WT_FROM_RAIN[cellLoc]*CELL_SIZE*GP_TP_IN_RAIN*0.001;
01173 
01174 /* calc various loss and/or vertical exchange potentials (kg P) */
01175               /* diffusive flux */
01176               TP_UPFLOW_POT[cellLoc]  = /*  advective upflow has been handled in surf-ground integration in fluxes.c  */
01177                   Max((TP_SEDWT_CONCACT[cellLoc]-PO4Pconc)
01178                         *GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0); 
01179               TP_UPFLOW[cellLoc]  = 
01180                   ( (TP_UPFLOW_POT[cellLoc])*DT>TP_SED_WT_AZ[cellLoc] ) ? 
01181                   ( (TP_SED_WT_AZ[cellLoc])/DT ) : 
01182                   ( TP_UPFLOW_POT[cellLoc]);
01183                   /* units = mgP/L */
01184               TP_K[cellLoc]  = Max(GP_TP_K_SLOPE*TP_SORBCONC[cellLoc]+GP_TP_K_INTER,0.0);
01185 /*fix rate */
01186               TP_SORB_POT[cellLoc]  = 
01187                   ( HYD_DOM_ACTWAT_PRES[cellLoc]>0.0 ) ? 
01188                   ( (double) 0.001 
01189                     *(TP_K[cellLoc]*(pow(Max(TP_SEDWT_CONCACT[cellLoc],0.0),0.8) )
01190                       *0.001*(DEPOS_ORG_MAT[cellLoc]*CELL_SIZE+DIM[cellLoc])-TP_SORB[cellLoc] ) ) : 
01191                   ( 0.0);
01192               if (TP_SORB_POT[cellLoc]>0.0) {
01193                   TP_SORBTION[cellLoc]  =  
01194                       ( (TP_SORB_POT[cellLoc]+TP_UPFLOW[cellLoc])*DT>TP_SED_WT_AZ[cellLoc] ) ? 
01195                       ( (TP_SED_WT_AZ[cellLoc]-TP_UPFLOW[cellLoc]*DT)/DT ) : 
01196                       ( TP_SORB_POT[cellLoc]);
01197               }
01198               else { /* neg sorption, loss from sorb variable*/
01199                   TP_SORBTION[cellLoc]  =  
01200                       ( (-TP_SORB_POT[cellLoc])*DT>TP_SORB[cellLoc] ) ? 
01201                       ( (-TP_SORB[cellLoc])/DT ) : 
01202                       ( TP_SORB_POT[cellLoc]);
01203               }
01204               
01205               /* diffusive + advective flux */
01206               TP_DNFLOW_POT[cellLoc]  = (SF_WT_INFILTRATION[cellLoc]+SF_WT_TO_SAT_DOWNFLOW[cellLoc])
01207                   *CELL_SIZE*TP_SFWT_CONC[cellLoc]   
01208                   + Max((PO4Pconc-TP_SEDWT_CONCACT[cellLoc])
01209                         *GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0) ;
01210               TP_DNFLOW[cellLoc]  =  
01211                   ( ( TP_DNFLOW_POT[cellLoc])*DT > TP_SF_WT[cellLoc] ) ? 
01212                   ( ( TP_SF_WT[cellLoc])/DT ) : 
01213                   ( TP_DNFLOW_POT[cellLoc]);
01214 /* calc P in sed water state variables (kg P) */
01215               TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] +
01216                   ( TP_DNFLOW[cellLoc] - TP_UPFLOW[cellLoc] - TP_SORBTION[cellLoc]) * DT;
01217              /* this is the active zone, where uptake, sorption, and mineralization take place */
01218               TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] +
01219                   (TP_DNFLOW[cellLoc] - TP_UPFLOW[cellLoc] - TP_SORBTION[cellLoc]) * DT;
01220 
01221 /* calc P in surface water state variable (kg P) */
01222               TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] + 
01223                   (TP_UPFLOW[cellLoc] + TP_FR_RAIN[cellLoc] 
01224                    - TP_DNFLOW[cellLoc]) * DT;
01225 
01226 /* calc P sorbed-to-soil state variable (kg P) */
01227               TP_SORB[cellLoc] = TP_SORB[cellLoc] + (TP_SORBTION[cellLoc]) * DT;
01228 
01229               TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
01230                   (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
01231                   (0.0); /* g/L */
01232               TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
01233                   ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
01234                   (TP_SED_CONC[cellLoc]); /* g/L */
01235               TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m^3==mg/L */
01236 
01237               TP_SORBCONC[cellLoc] = ((DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc])>0.0) ?
01238                   ( (double) TP_SORB[cellLoc]*conv_kgTOg / (DEPOS_ORG_MAT[cellLoc]*CELL_SIZE + DIM[cellLoc]) ) :
01239                   (0.0); /* gP/kgsoil */
01240 
01241               TP_SORBCONC_rep[cellLoc] = TP_SORBCONC[cellLoc] * conv_gTOmg; /* reporting purposes only (g/kg->mg/kg)*/
01242 
01243               TP_SFWT_CONC[cellLoc]  = 
01244                   ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01245                   ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01246                   ( 0.0); /* g/L used in P fluxes for mass balance */
01247               TP_SFWT_CONC_MG[cellLoc]  = 
01248                   ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
01249                   (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
01250                   (0.0); /* g/m^3==mg/L used for reporting and other modules to evaluate P conc when water is present */
01251 
01252 /* the following is an empirical method to calculate settling of particulate P out of the water column
01253    into the sediments.  It uses the fixed ratio of PO4 to TP, but allows for a decreasing proportion of
01254    TP to be in (large) particulate form as TP concentration drops below a chosen threshold - the sum of
01255    the TP is considered to be dissolved plus large particulate plus small particulate (that cannot settle out) */
01256                   /* mass (kg) of particulate fraction of TP, available for settling to sediments */
01257                   /* using regression for predicting PO4 from TP */
01258               PO4Pconc =  Max(TP_SFWT_CONC_MG[cellLoc]*GP_PO4toTP + GP_PO4toTPint,0.0);  /* mg/l (note that intercept may be <0)*/
01259               nonPO4Pconc = Max(TP_SFWT_CONC_MG[cellLoc]-PO4Pconc,0.0); /* non-PO4 conc, mg/l */
01260               TPpartic = nonPO4Pconc * (1.0-exp(-nonPO4Pconc/GP_TPpart_thresh)) *0.001 * SFWT_VOL[cellLoc] ; /* kg particulate P */
01261 
01262 
01263                   /* settling rate, 1/d */
01264               TPsettlRat = ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ?
01265                   (GP_settlVel/SURFACE_WAT[cellLoc]) :
01266                   0.0;
01267               
01268                   /* potential settling of particulate TP (kg/d) */
01269               TP_settl_pot = TPsettlRat * TPpartic;
01270               TP_settl[cellLoc]  =  
01271                   ( ( TP_settl_pot)*DT > TPpartic ) ? 
01272                   ( (TPpartic)/DT ) : 
01273                   ( TP_settl_pot);
01274 /* calc P in surface water state variable (kg P) */
01275               TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] - TP_settl[cellLoc] * DT;
01276 
01277 /* various book-keeping calcs used in other modules */
01278 /* conc surf and sed wat = kgP/m3==gP/L, another var calc for mg/L */
01279               TP_SFWT_CONC[cellLoc]  = 
01280                   ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01281                   ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01282                   ( 0.0); /* used in P fluxes for mass balance */
01283               TP_SFWT_CONC_MG[cellLoc]  = 
01284                   ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
01285                   (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
01286                   (0.0); /* g/m3==mg/L used for reporting and other modules to evaluate P conc when water is present */
01287               
01288               if (debug > 0 && TP_SF_WT[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg TP_SF_WT (%f m) in cell (%d,%d)!", 
01289                                                                     SimTime.TIME, TP_SF_WT[cellLoc], ix,iy ); WriteMsg( msgStr,True ); usrErr( msgStr ); dynERRORnum++; }
01290               if (debug > 0 && TP_SED_WT_AZ[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg TP_SED_WT_AZ (%f m) in cell (%d,%d)!", 
01291                                                                       SimTime.TIME, TP_SED_WT_AZ[cellLoc], ix,iy ); WriteMsg( msgStr,True ); usrErr( msgStr ); dynERRORnum++; }
01292 
01293     }
01294   }
01295   }
01296 }
01297 
01298 /*******/
01308 /*
01309 Parameter definitions:
01310       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
01311       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
01312 */
01313 void cell_dyn13(int step)
01314  {
01315  int ix, iy, cellLoc; 
01316  float TPsettlRat, TP_settl_pot;
01317  
01318   for(ix=1; ix<=s0; ix++) {
01319   for(iy=1; iy<=s1; iy++) { 
01320 
01321     if(ON_MAP[cellLoc= T(ix,iy)])  {
01322 /* concentration of P in surface water used in P fluxes for mass transfer (kgP/m3==gP/L) */
01323               TP_SFWT_CONC[cellLoc]  = 
01324                   ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01325                   ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01326                   ( 0.0); 
01327 /* rainfall inputs to surface water P (m * m^2 * g/m^3 * kg/g == kg P) */
01328               TP_FR_RAIN[cellLoc]  = SF_WT_FROM_RAIN[cellLoc]
01329                   *CELL_SIZE*GP_TP_IN_RAIN*0.001;              
01330 /* TP settling rate calculation (m/d) = 0 if water depth (m) less than depth threshold parm */
01331               if (SURFACE_WAT[cellLoc] > GP_WQMthresh ) {
01332                   TPsettlRat = WQMsettlVel[cellLoc];
01333                   TP_settlDays[cellLoc]++;
01334               }
01335               else
01336                   TPsettlRat = 0.0;
01337               
01338 /* before we had to put in the day accumulator*/
01339 /*               TPsettlRat = ( SURFACE_WAT[cellLoc] > GP_WQMthresh ) ?  */
01340 /*                   (WQMsettlVel[cellLoc]) : 0.0; */
01341 /* potential settling of particulate TP (m/d * m^2 * kg/m^3 = kg/d) */
01342               TP_settl_pot = TPsettlRat * CELL_SIZE * TP_SFWT_CONC[cellLoc];
01343 
01344 /*  like EWQM, shows mass bal error when ->   TP_settl[cellLoc]  =   TP_settl_pot; */
01345               
01346 /* constrain settling to be no more than kg P available in water column */
01347               TP_settl[cellLoc]  =   
01348                   ( ( TP_settl_pot)*DT > TP_SF_WT[cellLoc] ) ?  
01349                   ( (TP_SF_WT[cellLoc])/DT ) :  
01350                   ( TP_settl_pot); 
01351 /* calc P in surface water state variable (kg P) */
01352               TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] +
01353                   ( TP_FR_RAIN[cellLoc] - TP_settl[cellLoc] ) * DT;
01354 
01355 /* this was in EWQM!!! (mass balance error!):  if (TP_SF_WT[cellLoc]<0.0) TP_SF_WT[cellLoc]=0.0; */
01356               
01357 /* concentration of P in surface water used in P fluxes for mass transfer (kgP/m3==gP/L) */
01358               TP_SFWT_CONC[cellLoc]  = 
01359                   ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01360                   ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01361                   ( 0.0); 
01362 /* concentration used for reporting (e.g., in maps) when water is present. 
01363    The hydraulic detention depth parm also used for all
01364    concentration reporting thresholds (currently 1 cm)
01365 */
01366               TP_SFWT_CONC_MG[cellLoc]  = 
01367                   ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ?  
01368                   (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
01369                   (0.0);  /* g/m3==mg/L */
01370     }
01371   }
01372   }
01373 }
01374 
01375 
01376 
01377 /*******/
01385 /*
01386 SALT_SURF_WT[cellLoc] == mass of salt dissolved in surface water (kg salt)
01387 SALT_SED_WT[cellLoc] == mass of salt dissolved in sediment/soil (pore) water (kg salt)
01388 
01389 
01390 Parameter definitions:
01391       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms)
01392       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
01393 */
01394 void cell_dyn10(int step)
01395  {
01396 int ix, iy, cellLoc;
01397  double SALT_SED_TO_SF_FLOW_pot;
01398  
01399 
01400   for(ix=1; ix<=s0; ix++) {
01401   for(iy=1; iy<=s1; iy++) {
01402 
01403     if(ON_MAP[cellLoc= T(ix,iy)])  {
01404      SAL_SF_WT_mb[cellLoc]  = 
01405                      ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01406                      ( SALT_SURF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01407                      ( 0.0); /* used in salt fluxes for mass balance */
01408      SAL_SED_WT[cellLoc]  = 
01409                      ( HYD_SED_WAT_VOL[cellLoc]>0.0 ) ? 
01410                      ( SALT_SED_WT[cellLoc]/HYD_SED_WAT_VOL[cellLoc] ) : 
01411                      ( 0.0);
01412 
01413               /* diffusive + advective flux */ 
01414                  SALT_SFWAT_DOWNFL_POT[cellLoc]  = (SF_WT_INFILTRATION[cellLoc] + SF_WT_TO_SAT_DOWNFLOW[cellLoc])
01415                      * CELL_SIZE*SAL_SF_WT_mb[cellLoc]
01416                      + Max((SAL_SF_WT_mb[cellLoc]-SAL_SED_WT[cellLoc])
01417                            * GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0)  ; /* TODO: get cl diffusion parm; diffusion parms same as P */
01418      SALT_SFWAT_DOWNFL[cellLoc]  =  
01419                      ( SALT_SFWAT_DOWNFL_POT[cellLoc]*DT>SALT_SURF_WT[cellLoc] ) ? 
01420                      ( SALT_SURF_WT[cellLoc]/DT ) : 
01421                      ( SALT_SFWAT_DOWNFL_POT[cellLoc]);
01422 
01423               /* diffusive flux */  
01424                  SALT_SED_TO_SF_FLOW_pot  =  
01425                     /*  advective upflow has been handled in surf-ground integration in fluxes.c  */
01426                      Max((SAL_SED_WT[cellLoc]-SAL_SF_WT_mb[cellLoc])
01427                             *GP_TP_DIFFCOEF*8.64/GP_TP_DIFFDEPTH*CELL_SIZE,0.0)  ; /* TODO: get cl diffusion parm; diffusion parms same as P */
01428                  SALT_SED_TO_SF_FLOW[cellLoc]  =  
01429                      ( SALT_SED_TO_SF_FLOW_pot*DT>SALT_SED_WT[cellLoc] ) ? 
01430                      ( SALT_SED_WT[cellLoc]/DT ) : 
01431                      ( SALT_SED_TO_SF_FLOW_pot );
01432 /* calc state vars (kg salt) */
01433      SALT_SED_WT[cellLoc] =  SALT_SED_WT[cellLoc]  
01434                      + (SALT_SFWAT_DOWNFL[cellLoc] - SALT_SED_TO_SF_FLOW[cellLoc]) * DT;
01435 
01436      SALT_SURF_WT[cellLoc] = SALT_SURF_WT[cellLoc] 
01437                      + (SALT_SED_TO_SF_FLOW[cellLoc] - SALT_SFWAT_DOWNFL[cellLoc] ) * DT;
01438 
01439 /* book-keeping concentration calcs, (kg/m3==g/L==ppt) used in other modules */
01440      SAL_SF_WT_mb[cellLoc]  = 
01441                      ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01442                      ( SALT_SURF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01443                      ( 0.0); /* used in salt fluxes for mass balance */
01444      SAL_SF_WT[cellLoc]  = 
01445                      ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
01446                      ( SAL_SF_WT_mb[cellLoc]  ) : 
01447                      ( 0.0); /* used for reporting and other modules to evaluate salinity when water is present */
01448      SAL_SED_WT[cellLoc]  = 
01449                      ( HYD_SED_WAT_VOL[cellLoc]>0.0 ) ? 
01450                      ( SALT_SED_WT[cellLoc]/HYD_SED_WAT_VOL[cellLoc] ) : 
01451                      ( 0.0);
01452 
01453     }
01454   }
01455   }
01456                   
01457 }
01458 
01459 
01460 
01467 /*
01468    FLOC[cellLoc] == mass of organic flocculent material at the interface between soil and surface water (kg OM/m^2) \n
01469    FlocP[cellLoc] == mass of phosphorus in the flocculent material at the interface between soil and surface water (kg P/m^2) \n
01470 
01471 This routine originally was Suspended Organic Matter; was modified to instead
01472 represent the organic matter in the flocculent sediment layer. \n
01473 
01474 Parameter definitions: \n
01475       global parameters in GlobalParms.xls Excel sheet (text export=GlobalParms) \n
01476       habitat-specific parameters in HabParms.fmp FileMakerPro database (text export=HabParms)
01477 */
01478 void cell_dyn12(int step)
01479  {
01480  int ix, iy, cellLoc; 
01481  float FlocP_DECOMP_pot, FlocP_DEPO_pot, FlocP_settl, Floc_settl;
01482  
01483   for(ix=1; ix<=s0; ix++) {
01484   for(iy=1; iy<=s1; iy++) {
01485 
01486     if(ON_MAP[cellLoc= T(ix,iy)])  {
01487 /* inputs (kg OM / m2)  */
01488                   /* all periphyton mortality goes to floc */
01489               FLOC_FR_ALGAE[cellLoc]  = (double) (C_ALG_MORT[cellLoc]+NC_ALG_MORT[cellLoc])
01490                   /GP_ALG_C_TO_OM*0.001 ; 
01491                   /* all photobiomass mortality goes to floc */
01492               Floc_fr_phBio[cellLoc]  = phbio_mort_OM[cellLoc];
01493 
01494              /* all settling from water column goes to floc */
01495               FlocP_settl = TP_settl[cellLoc] / CELL_SIZE; /* kg P/m2; */
01496                   /* the particulate P settling is assumed a fixed Redfield P:OM ratio */
01497               Floc_settl =   FlocP_settl / GP_TP_P_OM; 
01498                  
01499           
01500 /* outputs (kg OM / m2) */
01501               FLOC_DECOMP_QUAL_CF[cellLoc]  = /* use the avg conc of sed and sf water here */
01502                   Exp(-GP_DOM_decomp_coef * Max(GP_DOM_DECOMP_POPT-
01503                                  (TP_SFWT_CONC_MG[cellLoc]+TP_SEDWT_CONCACTMG[cellLoc])/2.0, 0.0)
01504                       /GP_DOM_DECOMP_POPT) ; /* mg/L */
01505 /* decomp in surface water has higher rate than in soil,
01506  assuming this stock is of much higer substrate quality than the total soil layer */
01507 /* GP_calibDecomp is an adjustable calib parm */
01508               soil_MOIST_CF[cellLoc]  =  (UNSAT_DEPTH[cellLoc]>HP_DOM_AEROBTHIN[HAB[cellLoc]]) ?
01509                      ( Max(UNSAT_MOIST_PRP[cellLoc],0.0) ) :
01510                      ( 1.0 );
01511               FLOC_DECOMP_POT[cellLoc]  = GP_calibDecomp * 10.0*GP_DOM_RCDECOMP*FLOC[cellLoc]
01512                    *DOM_TEMP_CF[cellLoc] *FLOC_DECOMP_QUAL_CF[cellLoc] * soil_MOIST_CF[cellLoc];
01513               FLOC_DECOMP[cellLoc]  = 
01514                   ( (FLOC_DECOMP_POT[cellLoc])*DT>FLOC[cellLoc] ) ? 
01515                   ( (FLOC[cellLoc])/DT ) : 
01516                   ( FLOC_DECOMP_POT[cellLoc]);
01517 
01518 /* the incorporation of the floc layer into the "true" soil layer occurs
01519    at a rate dependent on the floc depth under flooded conditions, then constant rate under dry conditions */
01520               FLOC_DEPO_POT[cellLoc]  = 
01521                   ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
01522                   ( FLOC_Z[cellLoc]/GP_FlocMax * FLOC[cellLoc]*GP_Floc_rcSoil ) : 
01523                   ( FLOC[cellLoc]*GP_Floc_rcSoil);
01524               FLOC_DEPO[cellLoc]  = 
01525                   ( (FLOC_DEPO_POT[cellLoc]+FLOC_DECOMP[cellLoc])*DT>FLOC[cellLoc] ) ? 
01526                   ( (FLOC[cellLoc]-FLOC_DECOMP[cellLoc]*DT)/DT ) : 
01527                   ( FLOC_DEPO_POT[cellLoc]); 
01528 /* calc the state var (kg OM / m2) */
01529               FLOC[cellLoc] =  FLOC[cellLoc] 
01530                   + ( Floc_settl + Floc_fr_phBio[cellLoc] + FLOC_FR_ALGAE[cellLoc]
01531                       - FLOC_DECOMP[cellLoc] - FLOC_DEPO[cellLoc] ) * DT;
01532 
01533 /* the depth of floc is dependent on a fixed floc bulk density */
01534               FLOC_Z[cellLoc] = (double) FLOC[cellLoc] / GP_Floc_BD ;
01535                  
01536 
01537 /* Floc phosphorus (kg P / m2)  */
01538               FlocP_FR_ALGAE[cellLoc]  = (double) (NC_ALG_MORT_P[cellLoc]
01539                                           + C_ALG_MORT_P[cellLoc])*0.001; /* kg P/m2 */
01540               FlocP_PhBio[cellLoc] = phbio_mort_P[cellLoc] ;    
01541 
01542               FlocP_DECOMP_pot =  FLOC_DECOMP[cellLoc] * FlocP_OM[cellLoc];
01543               FlocP_DECOMP[cellLoc]  = 
01544                   ( (FlocP_DECOMP_pot)*DT>FlocP[cellLoc] ) ? 
01545                   ( (FlocP[cellLoc])/DT ) : 
01546                   ( FlocP_DECOMP_pot); 
01547               FlocP_DEPO_pot =  FLOC_DEPO[cellLoc] * FlocP_OM[cellLoc];
01548               FlocP_DEPO[cellLoc]  = 
01549                   ( (FlocP_DEPO_pot+FlocP_DECOMP[cellLoc])*DT>FlocP[cellLoc] ) ? 
01550                   ( (FlocP[cellLoc]-FlocP_DECOMP[cellLoc]*DT)/DT ) : 
01551                   ( FlocP_DEPO_pot); 
01552               
01553 /* calc the state var (kg P / m2) */
01554               FlocP[cellLoc] =  FlocP[cellLoc] 
01555                   + ( FlocP_settl + FlocP_PhBio[cellLoc] + FlocP_FR_ALGAE[cellLoc]
01556                       - FlocP_DECOMP[cellLoc] - FlocP_DEPO[cellLoc] ) * DT;
01557 
01558               FlocP_OM[cellLoc] = ( FLOC[cellLoc]>0.0) ? (FlocP[cellLoc]/FLOC[cellLoc]) : (0.0); /* kgP/kgOM */
01559               FlocP_OMrep[cellLoc] = (float) FlocP_OM[cellLoc] * conv_kgTOmg; /* mgP/kgOM, variable for output _rep-orting only */
01560               
01561               if (debug > 0 && FLOC[cellLoc] < -GP_MinCheck) { sprintf(msgStr,"Day %6.1f: ERROR - neg FLOC OM biomass (%f m) in cell (%d,%d)!", 
01562                                                                     SimTime.TIME, FLOC[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
01563               if (debug > 0 && FlocP[cellLoc] < -GP_MinCheck)  { sprintf(msgStr,"Day %6.1f: ERROR - neg FLOC P biomass (%f m) in cell (%d,%d)!", 
01564                                                                       SimTime.TIME, FlocP[cellLoc], ix,iy ); WriteMsg( msgStr,True ); dynERRORnum++;}
01565 
01566 /* now the P gain in sediment pore water with decomp - 90% goes to porewater, 10% to sfwat */
01567               TP_SED_MINER[cellLoc]  =  0.90 * FlocP_DECOMP[cellLoc] * CELL_SIZE ; 
01568 /* calc P in sed water state variables (kg P) */
01569               TP_SED_WT[cellLoc] =  TP_SED_WT[cellLoc] + 
01570                   (TP_SED_MINER[cellLoc]) * DT;
01571             /* this is the active zone, where uptake, sorption, and mineralization take place */
01572               TP_SED_WT_AZ[cellLoc] =  TP_SED_WT_AZ[cellLoc] + 
01573                   (TP_SED_MINER[cellLoc]) * DT;
01574 
01575               TP_SED_CONC[cellLoc] = (HYD_SED_WAT_VOL[cellLoc]>0.0) ?
01576                   (TP_SED_WT[cellLoc] / HYD_SED_WAT_VOL[cellLoc]) :
01577                   (0.0);
01578                TP_SEDWT_CONCACT[cellLoc] = ( HYD_DOM_ACTWAT_PRES[cellLoc] > 0.0) ?
01579                   ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
01580                   (TP_SED_CONC[cellLoc]);
01581               TP_SEDWT_CONCACTMG[cellLoc]  = TP_SEDWT_CONCACT[cellLoc]*conv_kgTOg; /* g/m3==mg/L */
01582               
01583               
01584 /* now the P gain in surface water with decomp - 90% goes to porewater, 10% to sfwat */
01585               TP_SFWT_MINER[cellLoc]  = 0.10*FlocP_DECOMP[cellLoc] * CELL_SIZE ;  
01586 /* calc P in surface water state variable (kg P) */
01587               TP_SF_WT[cellLoc] = TP_SF_WT[cellLoc] + 
01588                   (TP_SFWT_MINER[cellLoc]) * DT;
01589               TP_SFWT_CONC[cellLoc]  = 
01590                   ( SFWT_VOL[cellLoc] > 0.0 ) ? 
01591                   ( TP_SF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : 
01592                   ( 0.0); /* used in P fluxes for mass balance */
01593               TP_SFWT_CONC_MG[cellLoc]  = 
01594                   ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? 
01595                   (TP_SFWT_CONC[cellLoc]*conv_kgTOg) : 
01596                   (0.0); /* g/m3==mg/L used for reporting and other modules to evaluate P conc when water is present */
01597     }
01598   }
01599   }
01600 }
01601 
01602 
01613 float tempCF(int form, float c1, float topt, float tmin, float tmax, float tempC) 
01614 {
01615   if (form == 1) {
01616     /* Lassiter et al. 1975, where c1 = curvature parm */
01617     return (Exp(c1*(tempC-topt)) * pow(((tmax-tempC)/(tmax-topt)), (c1*(tmax-topt))) );
01618   }
01619   else {
01620     /* Jorgensen 1976 */
01621     return ( Exp(-2.3 * Abs((tempC-topt)/(topt-tmin))) );
01622   }
01623 }
01624 
01626 void init_static_data(void)
01627 {
01628   usrErr0("Acquiring static spatial data..."); /* console message */
01629   
01630   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 */
01631   read_map_file("BoundCond",BCondFlow,'d',1.0,0.0);       /* boundary condition flow restrictions, dimless integer attributes */
01632             /* ONLY when running the EWQM settling rate version of ELM */
01633   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) */
01634   read_map_file("basins",basn,'d',1.0,0.0);               /* Basins/Indicator-Region map, dimless integer attributes  */
01635 
01636   alloc_mem_stats (); /* allocate memory for budget & stats arrays (need to have read the Basin/Indicator-Region map first) */
01637   BIRinit(); /* Set up the Basin/Indicator-Region (BIR) linkages/inheritances */
01638   BIRoutfiles(); /* Basin/Indicator-Region output files */
01639   
01640   usrErr("Done.");
01641 } 
01642 
01643 
01645 void init_dynam_data(void)
01646 {
01647   usrErr0("Acquiring dynamic spatial data..."); /* console message */
01648   
01649   read_map_file("Elevation",ELEVATION,'f',0.01,0.0);      /*  positive elevation relative to NGVD 1929 datum (positive above NGVD29; all data in cm NGVD29, converted to m NGVD29) */
01650   read_map_file("HAB",HAB,'c',1.0,0.0);                   /* habitat (veg classifications, dimless unsigned char attributes) */
01651   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) */
01652   read_map_file("icSfWt",SURFACE_WAT,'f',0.01,0.0);       /* initial ponded surface water (data in cm, converted to m) */
01653   read_map_file("icUnsat",UNSAT_DEPTH,'f',0.01,0.0);      /* initial depth of unsaturated zone (data in cm, converted to m) */
01654   read_map_file("soilTP",TPtoSOIL,'f',0.000001,0.0);  /* soil TP map (data in mgP/kgSoil, converted to kgP/kgSoil) */
01655   read_map_file("soilBD",BulkD,'f',1.0,0.0);         /* soil bulk dens map (data in mgSoil/cm3 == kgSoil/m3)  */
01656   read_map_file("soil_orgBD",DOM_BD,'f',1.0,0.0);    /* organic soil bulk dens map (data in mgOM/cm3 == kgOM/m3)  */
01657 
01658   /* 2 static mapps need to be here for re-initialing w/ a multiplier (sensitivity analysis) */ 
01659   read_map_file("HydrCond",HYD_RCCONDUCT,'f',1.0,0.0);   /* hydraulic conductivity (no conversion, data in m/d)  */
01660   read_map_file("Bathymetry",Bathymetry,'f',0.01,0.0);      /* v2.3+, positive estuarine bathymetry (depth) relative to NGVD 1929 datum (positive depth below NGVD29; all positive depth data in cm, converted to m) */
01661 
01662   usrErr("Done."); /* console message */
01663 } 
01664 
01665 
01671 void init_eqns(void)
01672  {
01673 int ix,iy, cellLoc;
01674 float tmp; /* used in checking nutrient availability for MichMent kinetics */
01675 float min_litTemp; /* used to determine min of temper, light cf's for alg and macs */
01676 float I_ISat, Z_extinct, PO4Pconc, PO4P;
01677 float MAC_PHtoNPH, MAC_PHtoNPH_Init;
01678 float soil_propOrg; 
01679 #define satDensity 0.9 /* assign the relative proportion (0 - 1) of a population maximum to be the saturation density for a population */
01680 
01681   usrErr0("Initializing unit model equations..."); /* console message */
01682   
01683   SimTime.TIME = 0;
01684   DAYJUL = ( FMOD(SimTime.TIME,365.0) >0.0 ) ? ( FMOD(SimTime.TIME,365.0) ) : ( 365.0);
01685   LATRAD = ((int)(GP_LATDEG)+(GP_LATDEG-(int)(GP_LATDEG))*5.0/3.0)*PI/180.0;
01686   /* AMPL = Exp(7.42+0.045*LATRAD*180.0/PI)/3600.0; */
01687   /* DAYLENGTH = AMPL*Sin((DAYJUL-79.0)*0.01721)+12.0; */
01688   SOLALPHA = 32.9835-64.884*(1.0-1.3614*Cos(LATRAD));
01689   SOLDEC1 = 0.39785*Sin(4.868961+0.017203*DAYJUL   +0.033446*Sin(6.224111+0.017202*DAYJUL));
01690   SOLCOSDEC = sqrt(1.0-pow(SOLDEC1,2.0));
01691   SOLELEV_SINE = Sin(LATRAD)*SOLDEC1+Cos(LATRAD)*SOLCOSDEC;
01692   SOLALTCORR = (1.0-Exp(-0.014*(GP_ALTIT-274.0)/(SOLELEV_SINE*274.0)));
01693   SOLBETA = 0.715-0.3183*(1.0-1.3614*Cos(LATRAD));
01694   SOLDEC = Arctan(SOLDEC1/sqrt(1.0-pow(SOLDEC1,2.0)));
01695   SOLRISSET_HA1 = -Tan(LATRAD)*Tan(SOLDEC);
01696   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)));
01697   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));
01698         
01699         for(ix=0; ix<=s0+1; ix++) {
01700             for(iy=0; iy<=s1+1; iy++) { 
01701                   cellLoc = T(ix,iy);
01702                   
01703                   AIR_TEMP[cellLoc] = 25.0;  /* spatial time series unavailable after 1995; globally constant in v2.2+ */
01704 /* rainfall read from sfwmm data, remapped to ELM resolution */
01705                   boundcond_rain[cellLoc] =  SF_WT_FROM_RAIN[cellLoc] = boundcond_ETp[cellLoc] = boundcond_depth[cellLoc] = 0.0;
01706        
01707                                 /* used to have cloudiness influence on GP_SOLOMEGA, now 0 */
01708                  SOLRAD274[cellLoc] = SOLRADATMOS*(SOLBETA-GP_SOLOMEGA*0.0 )-SOLALPHA;
01709                  SOLRADGRD[cellLoc] = SOLRAD274[cellLoc]+((SOLRADATMOS+1.0)-SOLRAD274[cellLoc])*SOLALTCORR;
01710                  H2O_TEMP[cellLoc] = AIR_TEMP[cellLoc];
01711          
01712                  ALG_REFUGE[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*GP_ALG_REF_MULT;
01713                  ALG_SAT[cellLoc] = HP_ALG_MAX[HAB[cellLoc]]*0.9;
01714 
01715               /* v2.3: with southern everglades topo, put bathymetry back into the model */
01716                  ELEVATION[cellLoc] = ELEVATION[cellLoc] * GP_IC_ELEV_MULT;
01717                  Bathymetry[cellLoc] = Bathymetry[cellLoc] * GP_IC_BATHY_MULT;
01718                  SED_ELEV[cellLoc] =  ELEVATION[cellLoc] - Bathymetry[cellLoc] + GP_DATUM_DISTANCE; 
01719                  SED_INACT_Z[cellLoc] = SED_ELEV[cellLoc]-HP_DOM_MAXDEPTH[HAB[cellLoc]]; 
01720 
01721                  MAC_MAX_BIO[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]+HP_PHBIO_MAX[HAB[cellLoc]];
01722                  NPHBIO_REFUGE[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
01723                  NPHBIO_SAT[cellLoc] = HP_NPHBIO_MAX[HAB[cellLoc]] * satDensity;
01724                  PHBIO_REFUGE[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]]*GP_MAC_REFUG_MULT;
01725                  PHBIO_SAT[cellLoc] = HP_PHBIO_MAX[HAB[cellLoc]] * satDensity;
01726             }
01727         }
01728         for(ix=1; ix<s0+1; ix++) {
01729             for(iy=1; iy<s1+1; iy++) { cellLoc = T(ix,iy);
01730 /*initialization of major state variables */
01731 
01732               /* hydro */
01733                  HYD_RCCONDUCT[cellLoc] = HYD_RCCONDUCT[cellLoc] * GP_calibGWat;
01734 
01735 /* map */        UNSAT_DEPTH[cellLoc] = Max(UNSAT_DEPTH[cellLoc] + GP_HYD_IC_UNSAT_ADD,0.0); /* m */
01736                  SAT_WT_HEAD[cellLoc] = SED_ELEV[cellLoc]- UNSAT_DEPTH[cellLoc]; /* m */
01737                  SAT_WATER[cellLoc] = SAT_WT_HEAD[cellLoc] * HP_HYD_POROSITY[HAB[cellLoc]]; /* m */
01738 
01739                  UNSAT_WATER[cellLoc] = HP_HYD_POROSITY[HAB[cellLoc]]*UNSAT_DEPTH[cellLoc] *GP_HYD_ICUNSATMOIST; /* m */
01740                  UNSAT_CAP[cellLoc] = UNSAT_DEPTH[cellLoc]*HP_HYD_POROSITY[HAB[cellLoc]]; /* m */
01741                  UNSAT_MOIST_PRP[cellLoc]  =  /* dimless proportion, 0-1 */
01742                      ( UNSAT_CAP[cellLoc]>0.0 ) ? 
01743                      ( Min(UNSAT_WATER[cellLoc]/UNSAT_CAP[cellLoc],1.0) ) : 
01744                      ( 1.0); 
01745 
01746 /* map */        SURFACE_WAT[cellLoc] =  Max(SURFACE_WAT[cellLoc] + GP_HYD_IC_SFWAT_ADD, 0.0); /* m */
01747                  SFWT_VOL[cellLoc] = SURFACE_WAT[cellLoc]*CELL_SIZE;
01748 
01749                  HYD_DOM_ACTWAT_VOL[cellLoc] = Max(HP_DOM_MAXDEPTH[HAB[cellLoc]]-UNSAT_DEPTH[cellLoc]*
01750                                                    (1.0-UNSAT_MOIST_PRP[cellLoc]),0.0)*HP_HYD_POROSITY[HAB[cellLoc]]*CELL_SIZE;
01751                  HYD_SED_WAT_VOL[cellLoc] = (SAT_WATER[cellLoc]+UNSAT_WATER[cellLoc])*CELL_SIZE;
01752                  
01753               /* soil */
01754 /* map */        DOM_BD[cellLoc] = DOM_BD[cellLoc] * GP_IC_DOM_BD_MULT;
01755 /* map */        BulkD[cellLoc] = BulkD[cellLoc] * GP_IC_BulkD_MULT;
01756                  soil_propOrg = DOM_BD[cellLoc] / BulkD[cellLoc];
01757                  DIM[cellLoc] = (1.0 - soil_propOrg) * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] * CELL_SIZE; /* kg inorganic matter */
01758                  Inorg_Z[cellLoc] = (1.0 - soil_propOrg) * HP_DOM_MAXDEPTH[HAB[cellLoc]]; /*  fixed inorganic depth (m) */
01759                  DOM_Z[cellLoc] = HP_DOM_MAXDEPTH[HAB[cellLoc]] - Inorg_Z[cellLoc]; /* m */
01760 
01761                  DEPOS_ORG_MAT[cellLoc] = DOM_BD[cellLoc]*DOM_Z[cellLoc]; /* (mgOM/cm3 ==> kgOM/m3) * m = kgOM/m2 */
01762 
01763                  DOM_SED_AEROB_Z[cellLoc] = Min(UNSAT_DEPTH[cellLoc]+HP_DOM_AEROBTHIN[HAB[cellLoc]],HP_DOM_MAXDEPTH[HAB[cellLoc]]); /* m */
01764                  DOM_SED_ANAEROB_Z[cellLoc] = HP_DOM_MAXDEPTH[HAB[cellLoc]]-DOM_SED_AEROB_Z[cellLoc]; /* m */
01765 
01766 /* map */        TPtoSOIL[cellLoc] = TPtoSOIL[cellLoc] * GP_IC_TPtoSOIL_MULT; /* kgP/kgSoil */
01767                  DOP[cellLoc] =  (1.0-GP_sorbToTP)*TPtoSOIL[cellLoc] * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] ; /* kgP/kg_soil * kg_soil/m3 * m == kgP/m2 */
01768 
01769               /* floc layer of soil */
01770                  FLOC[cellLoc] = HP_FLOC_IC[HAB[cellLoc]]; /* kg OM/m2  */
01771                  FlocP[cellLoc] = FLOC[cellLoc]*HP_FLOC_IC_PC[HAB[cellLoc]]*HP_FLOC_IC_CTOOM[HAB[cellLoc]]; /* kg P /m2 */
01772                  FLOC_Z[cellLoc] = (double) FLOC[cellLoc] / GP_Floc_BD ; /* m */
01773 
01774               /* phosphorus */
01775 /* v2.4.4 */       TP_SFWT_CONC[cellLoc]  = GP_TP_ICSFWAT * conv_mgTOg; /* mg/L * g/mg  => g/L */
01776                  TP_SF_WT[cellLoc] =TP_SFWT_CONC[cellLoc] * SFWT_VOL[cellLoc]; /* kg P */
01777                  TP_SFWT_CONC_MG[cellLoc] = TP_SFWT_CONC[cellLoc] * conv_gTOmg; /* mg/L */
01778                       /* using regression for predicting PO4 from TP */
01779                  PO4Pconc = Max(TP_SFWT_CONC_MG[cellLoc]*GP_PO4toTP + GP_PO4toTPint, 0.10 * TP_SFWT_CONC_MG[cellLoc]); /* mg/L */
01780                  PO4P = PO4Pconc * SFWT_VOL[cellLoc] /conv_kgTOg; /*kg P available (from conc. in mg/l) */
01781 
01782 /* v2.4.4 */       TP_SED_CONC[cellLoc] = GP_TP_ICSEDWAT * conv_mgTOg; /* mg/L * g/mg => g/L */
01783                  TP_SED_WT[cellLoc] = TP_SED_CONC[cellLoc] * HYD_SED_WAT_VOL[cellLoc]; /* kg P */
01784                      /* this is the active zone, where uptake, sorption, and mineralization take place */
01785                  TP_Act_to_Tot[cellLoc] = 1.0 / HP_TP_CONC_GRAD[HAB[cellLoc]]; /* ratio of TP conc in active zone relative to total */
01786                  TP_SED_WT_AZ[cellLoc] = TP_SED_CONC[cellLoc] * TP_Act_to_Tot[cellLoc] * HYD_DOM_ACTWAT_VOL[cellLoc]; /* kg P */
01787                  TP_SEDWT_CONCACT[cellLoc] =(HYD_DOM_ACTWAT_VOL[cellLoc] > 0.0) ?
01788                      ( TP_SED_WT_AZ[cellLoc]/HYD_DOM_ACTWAT_VOL[cellLoc] ) :
01789                      ( TP_SED_CONC[cellLoc]); /* g/L */
01790                  TP_SEDWT_CONCACTMG[cellLoc] = TP_SEDWT_CONCACT[cellLoc]*conv_gTOmg; /* mg/L */
01791 
01792                  TP_SORB[cellLoc] = GP_sorbToTP*TPtoSOIL[cellLoc]
01793                      * BulkD[cellLoc] * HP_DOM_MAXDEPTH[HAB[cellLoc]] * CELL_SIZE; /* dimless * kgP/kg_soil * kg_soil/m3 * m * m^2 == kgP */
01794 
01795 
01796               /* salt */
01797                  SALT_SED_WT[cellLoc] = HYD_SED_WAT_VOL[cellLoc]*HP_SALT_ICSEDWAT[HAB[cellLoc]];
01798                  SAL_SED_WT[cellLoc] = ( HYD_SED_WAT_VOL[cellLoc]>0.0 ) ? ( SALT_SED_WT[cellLoc]/HYD_SED_WAT_VOL[cellLoc] ) : ( 0.0);
01799                  SALT_SURF_WT[cellLoc] = SFWT_VOL[cellLoc]*HP_SALT_ICSFWAT[HAB[cellLoc]];
01800                  SAL_SF_WT[cellLoc] = ( SURFACE_WAT[cellLoc] > GP_DetentZ ) ? ( SALT_SURF_WT[cellLoc]/SFWT_VOL[cellLoc] ) : ( 0.0);
01801 
01802               /* periphyton */
01803 /* 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 */
01804 /* 2.4.4 */       C_ALG[cellLoc]  = HP_ALG_MAX[HAB[cellLoc]] * GP_ALG_IC_MULT * (1.0 - GP_ALG_REF_MULT); /* g C/m2 */
01805                  /* ic PC of periph in oligotrophic area is 3% of max P:C, varies across space via (0.1->1.0) map */
01806 /* 2.4.4 */       NC_ALG_PC[cellLoc] = GP_ALG_PC; /* gP/ gC */
01807 /* 2.4.4 */       C_ALG_PC[cellLoc]  = GP_ALG_PC; /* gP/ gC */
01808 
01809                  NC_ALG_P[cellLoc] = NC_ALG[cellLoc]*NC_ALG_PC[cellLoc];   /* g P/m2 */
01810                  C_ALG_P[cellLoc] = C_ALG[cellLoc]*C_ALG_PC[cellLoc];  /* g P/m2 */  
01811 
01812               /* macrophytes */
01813 /* 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 */
01814                    /*  now calc the P and OM associated with that C */
01815 /* 2.4.4 */       mac_ph_PC[cellLoc] = HP_PHBIO_IC_PC[HAB[cellLoc]]; 
01816                  mac_ph_P[cellLoc] = MAC_PH_BIOMAS[cellLoc] * mac_ph_PC[cellLoc]; /* kg P/m2 */
01817                  mac_ph_OM[cellLoc] = MAC_PH_BIOMAS[cellLoc] / HP_PHBIO_IC_CTOOM[HAB[cellLoc]];
01818                  mac_ph_CtoOM[cellLoc] = HP_PHBIO_IC_CTOOM[HAB[cellLoc]];
01819                  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);
01820 
01821 /* 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 */
01822                    /*  now calc the P and OM associated with that C */
01823 /* 2.4.4 */       mac_nph_PC[cellLoc] = HP_NPHBIO_IC_PC[HAB[cellLoc]]; 
01824                  mac_nph_P[cellLoc] = MAC_NOPH_BIOMAS[cellLoc] * mac_nph_PC[cellLoc];  /* kg P/m2 */
01825                  mac_nph_OM[cellLoc] = MAC_NOPH_BIOMAS[cellLoc] / HP_NPHBIO_IC_CTOOM[HAB[cellLoc]];
01826                  mac_nph_CtoOM[cellLoc] = HP_NPHBIO_IC_CTOOM[HAB[cellLoc]];
01827                  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);
01828 
01829                  MAC_REL_BIOM[cellLoc] = ( MAC_TOT_BIOM[cellLoc] > 0 ) ? MAC_TOT_BIOM[cellLoc]/MAC_MAX_BIO[cellLoc] : 0.0;
01830                  MAC_LAI[cellLoc] = MAC_REL_BIOM[cellLoc]*HP_MAC_MAXLAI[HAB[cellLoc]];
01831                  MAC_HEIGHT[cellLoc] = MAC_REL_BIOM[cellLoc]*HP_MAC_MAXHT[HAB[cellLoc]];
01832                  LAI_eff[cellLoc] =  (MAC_HEIGHT[cellLoc]>0.0) ? (Max(1.0 - SURFACE_WAT[cellLoc]/MAC_HEIGHT[cellLoc], 0.0)*MAC_LAI[cellLoc]) : (0.0)  ;                 
01833 
01834 /* end of initialization of major state variables */
01835                  
01836 /* *************************** */
01837 /* These are the multiple calculations used if particular modules are turned off. \n
01838  NOTE: THIS section (of init_eqns() ) is not fully updated for properly
01839  turning off individual **interacting** *biological/chemical* modules.  
01840  If one *biological/chemical* module is turned off, 
01841  they all need to be turned off. (Note that cell_dyn's 3,5,6 should always be off). \n
01842 
01843  *** \n
01844  The following *biological/chemical* modules must be ON or OFF as a group:
01845  (cell_dyn2 + cell_dyn4 + cell_dyn8 + cell_dyn9  + cell_dyn12)
01846  cell_dyn13, the net settling rate module, can be turned on only when above module group is off
01847  *** \n
01848  
01849  In particular, the budget will
01850  not properly reflect actual dynamics if those 
01851  modules are not treated as a group.
01852 */
01853      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);
01854                      /* calcareous periphyton */
01855      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);
01856                  ALG_TEMP_CF[cellLoc]  = tempCF(0, 0.20, GP_ALG_TEMP_OPT, 5.0, 40.0, H2O_TEMP[cellLoc]);
01857      NC_ALG_RESP_POT[cellLoc]  = 
01858          ( UNSAT_DEPTH[cellLoc]<0.05 ) ? 
01859          ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc]*NC_ALG[cellLoc] ) : 
01860          ( 0.0);
01861      NC_ALG_RESP[cellLoc]  =  
01862          ( NC_ALG_RESP_POT[cellLoc]*DT>NC_ALG[cellLoc] ) ? 
01863          ( NC_ALG[cellLoc]/DT ) : 
01864          ( NC_ALG_RESP_POT[cellLoc]);
01865                      /* calcareous periphyton */
01866      C_ALG_RESP_POT[cellLoc]  = 
01867          ( UNSAT_DEPTH[cellLoc]<0.05 ) ? 
01868          ( GP_ALG_RC_RESP*ALG_TEMP_CF[cellLoc]*C_ALG[cellLoc] ) : 
01869          ( 0.0);
01870      C_ALG_RESP[cellLoc]  =  
01871          ( C_ALG_RESP_POT[cellLoc]*DT>C_ALG[cellLoc] ) ? 
01872          ( C_ALG[cellLoc]/DT ) : 
01873          ( C_ALG_RESP_POT[cellLoc]);
01874 
01875      NC_ALG_AVAIL_MORT[cellLoc] = NC_ALG[cellLoc]-ALG_REFUGE[cellLoc];
01876      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]);
01877                      /* calcareous periphyton */
01878      C_ALG_AVAIL_MORT[cellLoc]  = C_ALG[cellLoc]-ALG_REFUGE[cellLoc];
01879      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]);
01880 
01881 /* light, water, temperature controls apply to calc and non-calc periphyton */
01882      ALG_LIGHT_EXTINCT[cellLoc]  =  0.01; /* light extinction coef */
01883                      /* algal self-shading implicit in density-dependent constraint function later */
01884      ALG_INCID_LIGHT[cellLoc]  = SOLRADGRD[cellLoc]*Exp(-MAC_LAI[cellLoc]*GP_ALG_SHADE_FACTOR);
01885                  Z_extinct = SURFACE_WAT[cellLoc]*ALG_LIGHT_EXTINCT[cellLoc];
01886      I_ISat = ALG_INCID_LIGHT[cellLoc]/GP_ALG_LIGHT_SAT;
01887                      /*  averaged over whole water column (based on Steele '65) */
01888      ALG_LIGHT_CF[cellLoc]  = ( Z_extinct > 0.0 ) ? 
01889          ( 2.718/Z_extinct * (Exp(-I_ISat * Exp(-Z_extinct)) - Exp(-I_ISat)) ) :
01890                 (I_ISat*Exp(1.0-I_ISat));
01891                      /*  low-water growth constraint ready for something better based on data */
01892      ALG_WAT_CF[cellLoc]  = ( SURFACE_WAT[cellLoc]>0.0 ) ? ( 1.0 ) : ( 0.0);
01893 /* the 2 communities have same growth response to avail phosphorus - avail P is roughly calc'd from TP */
01894      NC_ALG_NUT_CF[cellLoc]  =  PO4Pconc/(PO4Pconc+GP_NC_ALG_KS_P) ;
01895      C_ALG_NUT_CF[cellLoc]  = PO4Pconc/(PO4Pconc+GP_C_ALG_KS_P); 
01896       min_litTemp = Min(ALG_LIGHT_CF[cellLoc],ALG_TEMP_CF[cellLoc]);
01897       NC_ALG_PROD_CF[cellLoc]  = Min(min_litTemp,ALG_WAT_CF[cellLoc])*NC_ALG_NUT_CF[cellLoc];
01898       C_ALG_PROD_CF[cellLoc]   = Min(min_litTemp,ALG_WAT_CF[cellLoc])*C_ALG_NUT_CF[cellLoc];
01899 /* gross production of the 2 communities (gC/m2, NOT kgC/m2) */
01900                      /* density constraint contains both noncalc and calc, competition effect accentuated by calc algae */
01901                      /* used to increase calc growth by factor of 10 */
01902       NC_ALG_GPP[cellLoc]  =  NC_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*NC_ALG[cellLoc]       
01903              *Max( (1.0-(GP_AlgComp*C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
01904      C_ALG_GPP[cellLoc]  =  C_ALG_PROD_CF[cellLoc]*GP_ALG_RC_PROD*C_ALG[cellLoc] 
01905        *Max( (1.0-(    C_ALG[cellLoc]+NC_ALG[cellLoc])/HP_ALG_MAX[HAB[cellLoc]]),0.0);
01906 /* check for available P mass (the MichMent kinetics nutCF do not) */
01907      tmp = ( (NC_ALG_GPP[cellLoc]+C_ALG_GPP[cellLoc]) > 0) ? 
01908          PO4P / ( (NC_ALG_GPP[cellLoc]+C_ALG_GPP[cellLoc]) 
01909          * 0.001*GP_ALG_PC*CELL_SIZE*DT) :
01910          1.0;
01911      if (tmp < 1.0) {
01912          NC_ALG_GPP[cellLoc] *= tmp;   
01913          C_ALG_GPP[cellLoc] *= tmp; 
01914     /* can have high conc, but low mass of P avail, in presence of high peri biomass and high demand */
01915     /* reduce the production proportionally if excess demand is found */
01916        }
01917 /* total of calc and noncalc algae available and their total NPP */
01918      NC_ALG_NPP[cellLoc]  = NC_ALG_GPP[cellLoc]-NC_ALG_RESP[cellLoc]; 
01919      C_ALG_NPP[cellLoc]  = C_ALG_GPP[cellLoc]-C_ALG_RESP[cellLoc]; 
01920 
01921 
01922      DOM_QUALITY_CF[cellLoc]  = (Max(TP_SFWT_CONC_MG[cellLoc],TP_SEDWT_CONCACTMG[cellLoc]))
01923          /GP_DOM_DECOMP_POPT; /* mg/L */
01924      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)));
01925      soil_MOIST_CF[cellLoc] = pow(Max(UNSAT_MOIST_PRP[cellLoc],0.0),0.75);
01926      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));
01927      DOM_DECOMP[cellLoc] =  ( (DOM_DECOMP_POT[cellLoc])*DT>DEPOS_ORG_MAT[cellLoc] ) ? ( (DEPOS_ORG_MAT[cellLoc])/DT ) : ( DOM_DECOMP_POT[cellLoc]);
01928 /* added for P DOM stoich */
01929      DOP_DECOMP[cellLoc] = DOM_DECOMP[cellLoc] * DOM_P_OM[cellLoc]; 
01930 
01931      SAT_VS_UNSAT[cellLoc] = 1/Exp(100.0*Max((SURFACE_WAT[cellLoc]-UNSAT_DEPTH[cellLoc]),0.0));
01932      UNSAT_WT_POT[cellLoc] = Max(UNSAT_CAP[cellLoc]-UNSAT_WATER[cellLoc],0.0);
01933          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]); 
01934      SAT_WT_RECHG[cellLoc] = ( GP_HYD_RCRECHG*DT>SAT_WATER[cellLoc] ) ? ( SAT_WATER[cellLoc]/DT ) : ( GP_HYD_RCRECHG);
01935      SF_WT_POT_INF[cellLoc] = ( (SURFACE_WAT[cellLoc]<HP_HYD_RCINFILT[HAB[cellLoc]]*DT) ) ? ( SURFACE_WAT[cellLoc]/DT ) : ( HP_HYD_RCINFILT[HAB[cellLoc]]);
01936      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]);
01937      HYD_DOM_ACTWAT_PRES[cellLoc] = ( HYD_DOM_ACTWAT_VOL[cellLoc] > 0.01 ) ? ( 1.0 ) : ( 0.0);
01938      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)));
01939 
01940      MAC_LIGHT_CF[cellLoc] = SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]*Exp(1.0-SOLRADGRD[cellLoc]/HP_MAC_LIGHTSAT[HAB[cellLoc]]);
01941      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]])));
01942      MAC_WATER_AVAIL_CF[cellLoc] = graph8(0x0,HYD_WATER_AVAIL[cellLoc]);
01943      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));
01944      MAC_NUT_CF[cellLoc] =  TP_SEDWT_CONCACT[cellLoc]/(TP_SEDWT_CONCACT[cellLoc]+HP_MAC_KSP[HAB[cellLoc]]*0.001) ;
01945 
01946      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));
01947      min_litTemp = Min(MAC_LIGHT_CF[cellLoc], MAC_TEMP_CF[cellLoc]);
01948      MAC_PROD_CF[cellLoc]  = Min(min_litTemp,MAC_WATER_CF[cellLoc])
01949          *MAC_NUT_CF[cellLoc]*MAC_SALT_CF[cellLoc];
01950      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);
01951 /* check for available P mass that will be taken up (MichMent kinetics in nutCF does not) */
01952      tmp = (PHBIO_NPP[cellLoc] > 0) ? 
01953          TP_SED_WT[cellLoc] / ( PHBIO_NPP[cellLoc] * HP_NPHBIO_IC_PC[HAB[cellLoc]]*CELL_SIZE*DT) :
01954          1.0;
01955      if (tmp < 1.0) PHBIO_NPP[cellLoc] *= tmp;   
01956     /* reduce the production proportionally if excess demand is found */
01957     /* can have high conc, but low mass of P avail, in presence of high plant biomass and high demand */
01958 /* now add the P and OM associated with that C */
01959      phbio_npp_P[cellLoc] = PHBIO_NPP[cellLoc] * HP_PHBIO_IC_PC[HAB[cellLoc]];     /* habitat-specfic stoichiometry */
01960      phbio_npp_OM[cellLoc] = PHBIO_NPP[cellLoc] / HP_PHBIO_IC_CTOOM[HAB[cellLoc]]; /* habitat-specfic stoichiometry */
01961 
01962 /* init ("target") ph/nph ratio and new transloc algorithm */
01963      MAC_PHtoNPH_Init = HP_PHBIO_MAX[HAB[cellLoc]] / HP_NPHBIO_MAX[HAB[cellLoc]] ;
01964      MAC_PHtoNPH = MAC_PH_BIOMAS[cellLoc] / MAC_NOPH_BIOMAS[cellLoc];
01965 
01966      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) */ 
01967 
01968      PHBIO_TRANSLOC[cellLoc] =  ( (PHBIO_TRANSLOC_POT[cellLoc])*DT >NPHBIO_AVAIL[cellLoc] ) ? ( (NPHBIO_AVAIL[cellLoc])/DT ) : ( PHBIO_TRANSLOC_POT[cellLoc]);
01969 /*  now remove the P and OM associated with that C */
01970      phbio_transl_P[cellLoc] = PHBIO_TRANSLOC[cellLoc] * mac_nph_PC[cellLoc];
01971      phbio_transl_OM[cellLoc] = PHBIO_TRANSLOC[cellLoc] / mac_nph_CtoOM[cellLoc];
01972      NPHBIO_MORT_POT[cellLoc] = NPHBIO_AVAIL[cellLoc]*HP_PHBIO_RCMORT[HAB[cellLoc]]*(1.0-MAC_PH_BIOMAS[cellLoc]/HP_PHBIO_MAX[HAB[cellLoc]]);
01973      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]);
01974      PHBIO_MORT_POT[cellLoc] = HP_PHBIO_RCMORT[HAB[cellLoc]] *PHBIO_AVAIL[cellLoc] *(1.0-MAC_WATER_AVAIL_CF[cellLoc]);
01975 /* now remove the P and OM associated with that C */
01976      nphbio_mort_P[cellLoc] = NPHBIO_MORT[cellLoc] * mac_nph_PC[cellLoc];
01977      nphbio_mort_OM[cellLoc] = NPHBIO_MORT[cellLoc] / mac_nph_CtoOM[cellLoc];
01978 
01979      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) */ 
01980      NPHBIO_TRANSLOC[cellLoc] = ( (NPHBIO_TRANSLOC_POT[cellLoc])*DT >MAC_PH_BIOMAS[cellLoc] ) ? ( (MAC_PH_BIOMAS[cellLoc])/DT ) : ( NPHBIO_TRANSLOC_POT[cellLoc]);
01981 /*  now remove the P and OM associated with that C */
01982      nphbio_transl_P[cellLoc] = NPHBIO_TRANSLOC[cellLoc] * mac_ph_PC[cellLoc];
01983      nphbio_transl_OM[cellLoc] = NPHBIO_TRANSLOC[cellLoc] / mac_ph_CtoOM[cellLoc];
01984      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]);
01985 /*  now remove the P associated with that C */
01986      phbio_mort_P[cellLoc] = PHBIO_MORT[cellLoc] * mac_ph_PC[cellLoc];
01987      phbio_mort_OM[cellLoc] = PHBIO_MORT[cellLoc] / mac_ph_CtoOM[cellLoc];
01988 
01989 
01990      FLOC_DECOMP_QUAL_CF[cellLoc] = Min(TP_SFWT_CONC_MG[cellLoc]/GP_DOM_DECOMP_POPT,1.0) ;
01991      FLOC_DECOMP_POT[cellLoc] = GP_DOM_RCDECOMP*FLOC[cellLoc]*DOM_TEMP_CF[cellLoc] *FLOC_DECOMP_QUAL_CF[cellLoc];
01992      FLOC_DECOMP[cellLoc] = ( (FLOC_DECOMP_POT[cellLoc])*DT>FLOC[cellLoc] ) ? ( (FLOC[cellLoc])/DT ) : ( FLOC_DECOMP_POT[cellLoc]);
01993      FLOC_DEPO_POT[cellLoc] = (  SURFACE_WAT[cellLoc] > GP_DetentZ ) ? ( FLOC[cellLoc]*0.05 ) : ( FLOC[cellLoc]/DT);
01994      FLOC_DEPO[cellLoc] = ( (FLOC_DEPO_POT[cellLoc]+FLOC_DECOMP[cellLoc])*DT>FLOC[cellLoc] ) ? ( (FLOC[cellLoc]-FLOC_DECOMP[cellLoc]*DT)/DT ) : ( FLOC_DEPO_POT[cellLoc]);
01995  
01996      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]]);
01997   } /* spatial loop end */
01998   } /* spatial loop end */
01999   usrErr("Done.");
02000 
02001 } /* end of init_eqns */
02002 
02003 
02007 void init_canals(int irun) 
02008 {
02009    if (irun == 1) {
02010       usrErr("Initializing Water Management...");
02011       Canal_Network_Init(GP_DATUM_DISTANCE,SED_ELEV); /* WatMgmt.c - initialize the canal network topology and data */
02012       usrErr("Done Water Management.");
02013    }
02014    else {
02015       reinitCanals();
02016    }
02017 
02018 }
02019 
02021 void init_succession(void) 
02022 {
02023    HabSwitch_Init( ); 
02024 
02025 }
02026 
02028 void reinitBIR(void) 
02029 {
02030    usrErr0("Re-initializing Basin/Indicator-Region info...");
02031    BIRstats_reset(); 
02032    BIRbudg_reset(); 
02033    Cell_reset_avg(); 
02034    Cell_reset_hydper(); 
02035    usrErr("Done.");
02036 }
02037 
02039 void reinitCanals(void) 
02040 {
02041    usrErr0("Re-initializing Canal depths & constituent masses...");
02042    CanalReInit();  
02043    usrErr("Done.");
02044 }
02045 
02046 
02059 void gen_output(int step, ViewParm *view)
02060 {
02061     #define numOutputs 50000
02062     static int iw[numOutputs];
02063     int oIndex;
02064     ViewParm   *vp;
02065 
02072 /* TODO: the ModelOutlistCreator does not provide the correct argument for non-floats in tgen[] code generator - hand corrected doubles and unsigned chars here */
02073 
02074     static outVar_struct tgen[] = {
02075       { (float**)&TP_settlDays, "TP_settlDays", 'f' },
02076       { (float**)&FLOC, "FLOC", 'f' },
02077       { (float**)&FLOC_DECOMP, "FLOC_DECOMP", 'f' },
02078       { (float**)&FLOC_DECOMP_POT, "FLOC_DECOMP_POT", 'f' },
02079       { (float**)&FLOC_DECOMP_QUAL_CF, "FLOC_DECOMP_QUAL_CF", 'f' },
02080       { (float**)&FLOC_DEPO, "FLOC_DEPO", 'f' },
02081       { (float**)&FLOC_DEPO_POT, "FLOC_DEPO_POT", 'f' },
02082       { (float**)&FLOC_FR_ALGAE, "FLOC_FR_ALGAE", 'f' },
02083       { (float**)&FLOC_Z, "FLOC_Z", 'f' },
02084       { (float**)&FlocP_OMrep, "FlocP_OMrep", 'f' },
02085       { (float**)&soil_MOIST_CF, "soil_MOIST_CF", 'f' },
02086       { (float**)&TP_SED_MINER, "TP_SED_MINER", 'f' },
02087       { (float**)&TP_SFWT_MINER, "TP_SFWT_MINER", 'f' },
02088       { (float**)&AIR_TEMP, "AIR_TEMP", 'f' },
02089       { (unsigned char**)&HAB, "HAB", 'c' },
02090       { (float**)&SOLRAD274, "SOLRAD274", 'f' },
02091       { (float**)&SOLRADGRD, "SOLRADGRD", 'f' },
02092       { (float**)&H2O_TEMP, "H2O_TEMP", 'f' },
02093       { (float**)&HYD_DOM_ACTWAT_PRES, "HYD_DOM_ACTWAT_PRES", 'f' },
02094       { (float**)&HYD_DOM_ACTWAT_VOL, "HYD_DOM_ACTWAT_VOL", 'f' },
02095       { (float**)&HYD_ET, "HYD_ET", 'f' },
02096       { (float**)&HYD_EVAP_CALC, "HYD_EVAP_CALC", 'f' },
02097       { (float**)&HYD_MANNINGS_N, "HYD_MANNINGS_N", 'f' },
02098       { (float**)&HYD_SAT_POT_TRANS, "HYD_SAT_POT_TRANS", 'f' },
02099       { (float**)&HYD_SED_WAT_VOL, "HYD_SED_WAT_VOL", 'f' },
02100       { (float**)&HYD_TOT_POT_TRANSP, "HYD_TOT_POT_TRANSP", 'f' },
02101       { (float**)&HYD_TRANSP, "HYD_TRANSP", 'f' },
02102       { (float**)&HYD_UNSAT_POT_TRANS, "HYD_UNSAT_POT_TRANS", 'f' },
02103       { (float**)&HYD_WATER_AVAIL, "HYD_WATER_AVAIL", 'f' },
02104       { (float**)&HydTotHd, "HydTotHd", 'f' },
02105       { (float**)&LAI_eff, "LAI_eff", 'f' },
02106       { (float**)&MAC_WATER_AVAIL_CF, "MAC_WATER_AVAIL_CF", 'f' },
02107       { (float**)&SAT_TO_UNSAT_FL, "SAT_TO_UNSAT_FL", 'f' },
02108       { (float**)&SAT_VS_UNSAT, "SAT_VS_UNSAT", 'f' },
02109       { (float**)&SAT_WATER, "SAT_WATER", 'f' },
02110       { (float**)&SAT_WT_HEAD, "SAT_WT_HEAD", 'f' },
02111       { (float**)&SAT_WT_RECHG, "SAT_WT_RECHG", 'f' },
02112       { (float**)&SAT_WT_TRANSP, "SAT_WT_TRANSP", 'f' },
02113       { (float**)&SF_WT_EVAP, "SF_WT_EVAP", 'f' },
02114       { (float**)&SF_WT_FROM_RAIN, "SF_WT_FROM_RAIN", 'f' },
02115       { (float**)&SF_WT_INFILTRATION, "SF_WT_INFILTRATION", 'f' },
02116       { (float**)&SF_WT_POT_INF, "SF_WT_POT_INF", 'f' },
02117       { (float**)&SF_WT_TO_SAT_DOWNFLOW, "SF_WT_TO_SAT_DOWNFLOW", 'f' },
02118       { (float**)&SFWT_VOL, "SFWT_VOL", 'f' },
02119       { (float**)&SURFACE_WAT, "SURFACE_WAT", 'f' },
02120       { (float**)&UNSAT_AVAIL, "UNSAT_AVAIL", 'f' },
02121       { (float**)&UNSAT_CAP, "UNSAT_CAP", 'f' },
02122       { (float**)&UNSAT_DEPTH, "UNSAT_DEPTH", 'f' },
02123       { (float**)&UNSAT_HYD_COND_CF, "UNSAT_HYD_COND_CF", 'f' },
02124       { (float**)&UNSAT_MOIST_PRP, "UNSAT_MOIST_PRP", 'f' },
02125       { (float**)&UNSAT_PERC, "UNSAT_PERC", 'f' },
02126       { (float**)&UNSAT_TO_SAT_FL, "UNSAT_TO_SAT_FL", 'f' },
02127       { (float**)&UNSAT_TRANSP, "UNSAT_TRANSP", 'f' },
02128       { (float**)&UNSAT_WATER, "UNSAT_WATER", 'f' },
02129       { (float**)&UNSAT_WT_POT, "UNSAT_WT_POT", 'f' },
02130       { (float**)&ELEVATION, "ELEVATION", 'f' },
02131       { (float**)&HYD_RCCONDUCT, "HYD_RCCONDUCT", 'f' },
02132       { (unsigned char**)&ON_MAP, "ON_MAP", 'c' },
02133       { (float**)&SED_INACT_Z, "SED_INACT_Z", 'f' },
02134       { (float**)&MAC_HEIGHT, "MAC_HEIGHT", 'f' },
02135       { (float**)&MAC_LAI, "MAC_LAI", 'f' },
02136       { (float**)&MAC_LIGHT_CF, "MAC_LIGHT_CF", 'f' },
02137       { (float**)&MAC_MAX_BIO, "MAC_MAX_BIO", 'f' },
02138       { (float**)&MAC_NOPH_BIOMAS, "MAC_NOPH_BIOMAS", 'f' },
02139       { (float**)&mac_nph_PC_rep, "mac_nph_PC_rep", 'f' },
02140       { (float**)&MAC_NUT_CF, "MAC_NUT_CF", 'f' },
02141       { (float**)&MAC_PH_BIOMAS, "MAC_PH_BIOMAS", 'f' },
02142       { (float**)&mac_ph_PC_rep, "mac_ph_PC_rep", 'f' },
02143       { (float**)&MAC_PROD_CF, "MAC_PROD_CF", 'f' },
02144       { (float**)&MAC_REL_BIOM, "MAC_REL_BIOM", 'f' },
02145       { (float**)&MAC_SALT_CF, "MAC_SALT_CF", 'f' },
02146       { (float**)&MAC_TEMP_CF, "MAC_TEMP_CF", 'f' },
02147       { (float**)&MAC_TOT_BIOM, "MAC_TOT_BIOM", 'f' },
02148       { (float**)&MAC_WATER_CF, "MAC_WATER_CF", 'f' },
02149       { (float**)&NPHBIO_AVAIL, "NPHBIO_AVAIL", 'f' },
02150       { (float**)&NPHBIO_MORT, "NPHBIO_MORT", 'f' },
02151       { (float**)&NPHBIO_MORT_POT, "NPHBIO_MORT_POT", 'f' },
02152       { (float**)&NPHBIO_REFUGE, "NPHBIO_REFUGE", 'f' },
02153       { (float**)&NPHBIO_SAT, "NPHBIO_SAT", 'f' },
02154       { (float**)&NPHBIO_TRANSLOC, "NPHBIO_TRANSLOC", 'f' },
02155       { (float**)&NPHBIO_TRANSLOC_POT, "NPHBIO_TRANSLOC_POT", 'f' },
02156       { (float**)&PHBIO_AVAIL, "PHBIO_AVAIL", 'f' },
02157       { (float**)&PHBIO_MORT, "PHBIO_MORT", 'f' },
02158       { (float**)&PHBIO_MORT_POT, "PHBIO_MORT_POT", 'f' },
02159       { (float**)&PHBIO_NPP, "PHBIO_NPP", 'f' },
02160       { (float**)&PHBIO_REFUGE, "PHBIO_REFUGE", 'f' },
02161       { (float**)&PHBIO_SAT, "PHBIO_SAT", 'f' },
02162       { (float**)&PHBIO_TRANSLOC, "PHBIO_TRANSLOC", 'f' },
02163       { (float**)&PHBIO_TRANSLOC_POT, "PHBIO_TRANSLOC_POT", 'f' },
02164       { (float**)&TP_SEDWT_UPTAKE, "TP_SEDWT_UPTAKE", 'f' },
02165       { (float**)&ALG_INCID_LIGHT, "ALG_INCID_LIGHT", 'f' },
02166       { (float**)&ALG_LIGHT_CF, "ALG_LIGHT_CF", 'f' },
02167       { (float**)&ALG_LIGHT_EXTINCT, "ALG_LIGHT_EXTINCT", 'f' },
02168       { (float**)&ALG_REFUGE, "ALG_REFUGE", 'f' },
02169       { (float**)&ALG_SAT, "ALG_SAT", 'f' },
02170       { (float**)&ALG_TEMP_CF, "ALG_TEMP_CF", 'f' },
02171       { (float**)&ALG_TOT, "ALG_TOT", 'f' },
02172       { (float**)&ALG_WAT_CF, "ALG_WAT_CF", 'f' },
02173       { (float**)&C_ALG, "C_ALG", 'f' },
02174       { (float**)&C_ALG_AVAIL_MORT, "C_ALG_AVAIL_MORT", 'f' },
02175       { (float**)&C_ALG_GPP, "C_ALG_GPP", 'f' },
02176       { (float**)&C_ALG_MORT, "C_ALG_MORT", 'f' },
02177       { (float**)&C_ALG_MORT_POT, "C_ALG_MORT_POT", 'f' },
02178       { (float**)&C_ALG_NPP, "C_ALG_NPP", 'f' },
02179       { (float**)&C_ALG_NUT_CF, "C_ALG_NUT_CF", 'f' },
02180       { (float**)&C_ALG_PROD_CF, "C_ALG_PROD_CF", 'f' },
02181       { (float**)&C_ALG_RESP, "C_ALG_RESP", 'f' },
02182       { (float**)&C_ALG_RESP_POT, "C_ALG_RESP_POT", 'f' },
02183       { (float**)&NC_ALG, "NC_ALG", 'f' },
02184       { (float**)&NC_ALG_AVAIL_MORT, "NC_ALG_AVAIL_MORT", 'f' },
02185       { (float**)&NC_ALG_GPP, "NC_ALG_GPP", 'f' },
02186       { (float**)&NC_ALG_MORT, "NC_ALG_MORT", 'f' },
02187       { (float**)&NC_ALG_MORT_POT, "NC_ALG_MORT_POT", 'f' },
02188       { (float**)&NC_ALG_NPP, "NC_ALG_NPP", 'f' },
02189       { (float**)&NC_ALG_NUT_CF, "NC_ALG_NUT_CF", 'f' },
02190       { (float**)&NC_ALG_PROD_CF, "NC_ALG_PROD_CF", 'f' },
02191       { (float**)&NC_ALG_RESP, "NC_ALG_RESP", 'f' },
02192       { (float**)&NC_ALG_RESP_POT, "NC_ALG_RESP_POT", 'f' },
02193       { (float**)&TP_SFWT_UPTAK, "TP_SFWT_UPTAK", 'f' },
02194       { (float**)&TP_Act_to_TotRep, "TP_Act_to_TotRep", 'f' },
02195       { (float**)&TP_DNFLOW, "TP_DNFLOW", 'f' },
02196       { (float**)&TP_DNFLOW_POT, "TP_DNFLOW_POT", 'f' },
02197       { (float**)&TP_FR_RAIN, "TP_FR_RAIN", 'f' },
02198       { (float**)&TP_K, "TP_K", 'f' },
02199       { (double**)&TP_SED_CONC, "TP_SED_CONC", 'l' },
02200       { (double**)&TP_SED_WT, "TP_SED_WT", 'l' },
02201       { (double**)&TP_SED_WT_AZ, "TP_SED_WT_AZ", 'l' },
02202       { (double**)&TP_SEDWT_CONCACT, "TP_SEDWT_CONCACT", 'l' },
02203       { (float**)&TP_SEDWT_CONCACTMG, "TP_SEDWT_CONCACTMG", 'f' },
02204       { (float**)&TP_settl, "TP_settl", 'f' },
02205       { (double**)&TP_SF_WT, "TP_SF_WT", 'l' },
02206       { (double**)&TP_SFWT_CONC, "TP_SFWT_CONC", 'l' },
02207       { (float**)&TP_SFWT_CONC_MG, "TP_SFWT_CONC_MG", 'f' },
02208       { (double**)&TP_SORB, "TP_SORB", 'l' },
02209       { (float**)&TP_SORB_POT, "TP_SORB_POT", 'f' },
02210       { (double**)&TP_SORBCONC, "TP_SORBCONC", 'l' },
02211       { (float**)&TP_SORBCONC_rep, "TP_SORBCONC_rep", 'f' },
02212       { (float**)&TP_SORBTION, "TP_SORBTION", 'f' },
02213       { (float**)&TP_UPFLOW, "TP_UPFLOW", 'f' },
02214       { (float**)&TP_UPFLOW_POT, "TP_UPFLOW_POT", 'f' },
02215       { (double**)&SAL_SED_WT, "SAL_SED_WT", 'l' },
02216       { (float**)&SAL_SF_WT, "SAL_SF_WT", 'f' },
02217       { (float**)&SALT_SED_TO_SF_FLOW, "SALT_SED_TO_SF_FLOW", 'f' },
02218       { (double**)&SALT_SED_WT, "SALT_SED_WT", 'l' },
02219       { (float**)&SALT_SFWAT_DOWNFL, "SALT_SFWAT_DOWNFL", 'f' },
02220       { (float**)&SALT_SFWAT_DOWNFL_POT, "SALT_SFWAT_DOWNFL_POT", 'f' },
02221       { (double**)&SALT_SURF_WT, "SALT_SURF_WT", 'l' },
02222       { (double**)&DEPOS_ORG_MAT, "DEPOS_ORG_MAT", 'l' },
02223       { (float**)&DOM_DECOMP, "DOM_DECOMP", 'f' },
02224       { (float**)&DOM_DECOMP_POT, "DOM_DECOMP_POT", 'f' },
02225       { (float**)&DOM_FR_FLOC, "DOM_FR_FLOC", 'f' },
02226       { (double**)&DOM_P_OM, "DOM_P_OM", 'l' },
02227       { (float**)&DOM_QUALITY_CF, "DOM_QUALITY_CF", 'f' },
02228       { (float**)&DOM_SED_AEROB_Z, "DOM_SED_AEROB_Z", 'f' },
02229       { (float**)&DOM_SED_ANAEROB_Z, "DOM_SED_ANAEROB_Z", 'f' },
02230       { (float**)&DOM_TEMP_CF, "DOM_TEMP_CF", 'f' },
02231       { (float**)&DOM_Z, "DOM_Z", 'f' },
02232       { (double**)&DOP, "DOP", 'l' },
02233       { (float**)&DOP_DECOMP, "DOP_DECOMP", 'f' },
02234       { (float**)&P_SUM_CELL, "P_SUM_CELL", 'f' },
02235       { (float**)&SED_ELEV, "SED_ELEV", 'f' },
02236       { (float**)&TPtoSOIL_rep, "TPtoSOIL_rep", 'f' },
02237       { (float**)&Floc_fr_phBioAvg, "Floc_fr_phBioAvg", 'f' },
02238       { (float**)&TPSfMinAvg, "TPSfMinAvg", 'f' },
02239       { (float**)&ETAvg, "ETAvg", 'f' },
02240       { (float**)&EvapAvg, "EvapAvg", 'f' },
02241       { (float**)&HydPerAnn, "HydPerAnn", 'f' },
02242       { (float**)&LAI_effAvg, "LAI_effAvg", 'f' },
02243       { (float**)&Manning_nAvg, "Manning_nAvg", 'f' },
02244       { (float**)&RainAvg, "RainAvg", 'f' },
02245       { (float**)&SfWatAvg, "SfWatAvg", 'f' },
02246       { (float**)&TotHeadAvg, "TotHeadAvg", 'f' },
02247       { (float**)&TranspAvg, "TranspAvg", 'f' },
02248       { (float**)&UnsatMoistAvg, "UnsatMoistAvg", 'f' },
02249       { (float**)&UnsatZavg, "UnsatZavg", 'f' },
02250       { (float**)&mac_nph_PCAvg, "mac_nph_PCAvg", 'f' },
02251       { (float**)&Mac_nphBioAvg, "Mac_nphBioAvg", 'f' },
02252       { (float**)&Mac_nphMortAvg, "Mac_nphMortAvg", 'f' },
02253       { (float**)&Mac_nppAvg, "Mac_nppAvg", 'f' },
02254       { (float**)&mac_ph_PCAvg, "mac_ph_PCAvg", 'f' },
02255       { (float**)&Mac_phBioAvg, "Mac_phBioAvg", 'f' },
02256       { (float**)&Mac_phMortAvg, "Mac_phMortAvg", 'f' },
02257       { (float**)&Mac_totBioAvg, "Mac_totBioAvg", 'f' },
02258       { (float**)&MacNutCfAvg, "MacNutCfAvg", 'f' },
02259       { (float**)&MacWatCfAvg, "MacWatCfAvg", 'f' },
02260       { (float**)&TPSedUptAvg, "TPSedUptAvg", 'f' },
02261       { (float**)&C_Peri_mortAvg, "C_Peri_mortAvg", 'f' },
02262       { (float**)&C_Peri_nppAvg, "C_Peri_nppAvg", 'f' },
02263       { (float**)&C_Peri_PCAvg, "C_Peri_PCAvg", 'f' },
02264       { (float**)&C_PeriAvg, "C_PeriAvg", 'f' },
02265       { (float**)&C_PeriNutCFAvg, "C_PeriNutCFAvg", 'f' },
02266       { (float**)&C_PeriRespAvg, "C_PeriRespAvg", 'f' },
02267       { (float**)&NC_Peri_mortAvg, "NC_Peri_mortAvg", 'f' },
02268       { (float**)&NC_Peri_nppAvg, "NC_Peri_nppAvg", 'f' },
02269       { (float**)&NC_Peri_PCAvg, "NC_Peri_PCAvg", 'f' },
02270       { (float**)&NC_PeriAvg, "NC_PeriAvg", 'f' },
02271       { (float**)&NC_PeriNutCFAvg, "NC_PeriNutCFAvg", 'f' },
02272       { (float**)&NC_PeriRespAvg, "NC_PeriRespAvg", 'f' },
02273       { (float**)&PeriAvg, "PeriAvg", 'f' },
02274       { (float**)&PeriLiteCFAvg, "PeriLiteCFAvg", 'f' },
02275       { (float**)&TPSfUptAvg, "TPSfUptAvg", 'f' },
02276       { (float**)&TP_settlAvg, "TP_settlAvg", 'f' },
02277       { (float**)&TPSedWatAvg, "TPSedWatAvg", 'f' },
02278       { (float**)&TPSfWatAvg, "TPSfWatAvg", 'f' },
02279       { (float**)&SaltSedAvg, "SaltSedAvg", 'f' },
02280       { (float**)&SaltSfAvg, "SaltSfAvg", 'f' },
02281       { (float**)&SedElevAvg, "SedElevAvg", 'f' },
02282       { (float**)&TPSedMinAvg, "TPSedMinAvg", 'f' },
02283       { (float**)&TPSorbAvg, "TPSorbAvg", 'f' },
02284       { (float**)&TPtoSOILAvg, "TPtoSOILAvg", 'f' },
02285       { (float**)&TPtoVOLAvg, "TPtoVOLAvg", 'f' },
02286 
02287 
02288       { NULL, NULL, '\0' },
02289   };
02290 
02291     outVar_struct *ptable;
02292     int  i;
02293 
02294     for (i = 0, ptable = tgen; ptable->pfvar != NULL; ptable++, i++)
02295     {
02296   vp = view + i;
02297 
02298 /* TODO: develop flexible output of calendar-based and julian-based outsteps (jan 11, 2005) */
02299 if (vp->step > 0)
02300                 
02301             if (strcmp(ptable->pname,"HydPerAnn")!=0) { /* i.e., not the HydPerAnn variable */
02302                 
02303                 if  (step % vp->step == 0 && (vp->step != CalMonOut) ) {   /* standard julian-day outstep interval variables (note: the != CalMonOut needed for step=0) */ 
02304                     oIndex = iw[i]++;
02305                     write_output(oIndex, vp, *(ptable->pfvar), ptable->pname, ptable->ctype, step);
02306                 }
02307                 else if ( (avgPrint) && (vp->step == CalMonOut) ) { /* variables output at the special 1-calendar-month outstep interval */  
02308                     oIndex = iw[i]++;
02309                     write_output(oIndex, vp, *(ptable->pfvar), ptable->pname, ptable->ctype, step);
02310                 }
02311             }
02312         
02313             else
02314                 if (FMOD(DAYJUL, 273.0) ==0) { /* hydroperiod is printed at a special-case time (approximately Oct 1 every year) */
02315                     oIndex = iw[i]++;
02316                     write_output(oIndex, vp, *(ptable->pfvar), ptable->pname, ptable->ctype, step);
02317                 }
02318         
02319      
02320     }
02321 
02322         /* after printing, zero the arrays holding averages or hydroperiods (v2.2.1note using avgPrint var)*/
02323     if (avgPrint) {
02324         Cell_reset_avg();
02325     }
02326 
02327     if (FMOD(DAYJUL, 273.0) ==0) {
02328         Cell_reset_hydper();
02329     }
02330 
02331 
02332 } /* end of gen_output routine */
02333 
02334 
02348 void ReadHabParms ( char* s_parm_name, int s_parm_relval )
02349 {
02350 
02351  /* Periphyton/Algae */
02352   HP_ALG_MAX = get_hab_parm(s_parm_name, s_parm_relval, "HP_ALG_MAX"); 
02353 
02354 /* Soil */
02355   HP_DOM_MAXDEPTH = get_hab_parm(s_parm_name, s_parm_relval,"HP_DOM_MAXDEPTH");  
02356   HP_DOM_AEROBTHIN = get_hab_parm(s_parm_name, s_parm_relval,"HP_DOM_AEROBTHIN"); 
02357 
02358 /* Phosphorus */
02359   HP_TP_CONC_GRAD = get_hab_parm(s_parm_name, s_parm_relval,"HP_TP_CONC_GRAD"); 
02360 
02361 /* Salt/tracer */
02362   HP_SALT_ICSEDWAT = get_hab_parm(s_parm_name, s_parm_relval,"HP_SALT_ICSEDWAT"); 
02363   HP_SALT_ICSFWAT = get_hab_parm(s_parm_name, s_parm_relval,"HP_SALT_ICSFWAT"); 
02364         
02365 /* Macrophytes */
02366   HP_PHBIO_MAX = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_MAX"); 
02367   HP_NPHBIO_MAX = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_MAX"); 
02368   HP_MAC_MAXHT = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXHT"); 
02369   HP_NPHBIO_ROOTDEPTH = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_ROOTDEPTH"); 
02370   HP_MAC_MAXROUGH = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXROUGH"); 
02371   HP_MAC_MINROUGH = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MINROUGH"); 
02372   HP_MAC_MAXLAI = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXLAI"); 
02373   HP_MAC_MAXCANOPCOND = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_MAXCANOPCOND"); /* unused in ELMv2.2, 2.3 */ 
02374   HP_MAC_CANOPDECOUP = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_CANOPDECOUP"); /* unused in ELMv2.2, 2.3 */ 
02375   HP_MAC_TEMPOPT = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_TEMPOPT"); 
02376   HP_MAC_LIGHTSAT = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_LIGHTSAT"); 
02377   HP_MAC_KSP = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_KSP"); 
02378   HP_PHBIO_RCNPP = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_RCNPP"); 
02379   HP_PHBIO_RCMORT = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_RCMORT"); 
02380   HP_MAC_WAT_TOLER = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_WAT_TOLER"); 
02381   HP_MAC_SALIN_THRESH = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_SALIN_THRESH"); /* unused in ELMv2.2, 2.3 */ 
02382   HP_PHBIO_IC_CTOOM = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_IC_CTOOM"); 
02383   HP_NPHBIO_IC_CTOOM = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_IC_CTOOM"); 
02384   HP_PHBIO_IC_PC = get_hab_parm(s_parm_name, s_parm_relval,"HP_PHBIO_IC_PC"); 
02385   HP_NPHBIO_IC_PC = get_hab_parm(s_parm_name, s_parm_relval,"HP_NPHBIO_IC_PC"); 
02386   HP_MAC_TRANSLOC_RC = get_hab_parm(s_parm_name, s_parm_relval,"HP_MAC_TRANSLOC_RC"); 
02387 
02388 /* Hydrology */
02389   HP_HYD_RCINFILT = get_hab_parm(s_parm_name, s_parm_relval,"HP_HYD_RCINFILT"); 
02390   HP_HYD_SPEC_YIELD = get_hab_parm(s_parm_name, s_parm_relval,"HP_HYD_SPEC_YIELD");  
02391   HP_HYD_POROSITY = get_hab_parm(s_parm_name, s_parm_relval,"HP_HYD_POROSITY"); 
02392 
02393 /* Floc */
02394   HP_FLOC_IC = get_hab_parm(s_parm_name, s_parm_relval,"HP_FLOC_IC"); 
02395   HP_FLOC_IC_CTOOM = get_hab_parm(s_parm_name, s_parm_relval,"HP_FLOC_IC_CTOOM"); 
02396   HP_FLOC_IC_PC = get_hab_parm(s_parm_name, s_parm_relval,"HP_FLOC_IC_PC"); 
02397 
02398 /* Habitat succession */
02399   HP_SfDepthLo = get_hab_parm(s_parm_name, s_parm_relval,"HP_SfDepthLo"); 
02400   HP_SfDepthHi = get_hab_parm(s_parm_name, s_parm_relval,"HP_SfDepthHi"); 
02401   HP_SfDepthInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_SfDepthInt"); 
02402   HP_PhosLo = get_hab_parm(s_parm_name, s_parm_relval,"HP_PhosLo"); 
02403   HP_PhosHi = get_hab_parm(s_parm_name, s_parm_relval,"HP_PhosHi"); 
02404   HP_PhosInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_PhosInt"); 
02405   HP_FireInt = get_hab_parm(s_parm_name, s_parm_relval,"HP_FireInt"); 
02406   
02407 }
02408 
02409 
02414 void ReadGlobalParms( char* s_parm_name, int s_parm_relval )
02415  {
02416      
02417 /* Geography, hydrology */
02418     GP_SOLOMEGA = get_global_parm(s_parm_name, s_parm_relval,"GP_SOLOMEGA"); 
02419     GP_ALTIT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALTIT"); 
02420     GP_LATDEG = get_global_parm(s_parm_name, s_parm_relval,"GP_LATDEG"); 
02421     GP_mannDepthPow = get_global_parm(s_parm_name, s_parm_relval,"GP_mannDepthPow"); 
02422     GP_mannHeadPow = get_global_parm(s_parm_name, s_parm_relval,"GP_mannHeadPow"); 
02423     GP_calibGWat = get_global_parm(s_parm_name, s_parm_relval,"GP_calibGWat"); 
02424     GP_IDW_pow = get_global_parm(s_parm_name, s_parm_relval,"GP_IDW_pow"); 
02425     GP_calibET = get_global_parm(s_parm_name, s_parm_relval,"GP_calibET"); 
02426     GP_DATUM_DISTANCE = get_global_parm(s_parm_name, s_parm_relval,"GP_DATUM_DISTANCE"); 
02427     GP_HYD_IC_SFWAT_ADD = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_IC_SFWAT_ADD"); 
02428     GP_HYD_IC_UNSAT_ADD = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_IC_UNSAT_ADD"); 
02429     GP_HYD_RCRECHG = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_RCRECHG"); 
02430     GP_HYD_ICUNSATMOIST = get_global_parm(s_parm_name, s_parm_relval,"GP_HYD_ICUNSATMOIST"); 
02431     GP_DetentZ = get_global_parm(s_parm_name, s_parm_relval,"GP_DetentZ"); 
02432     GP_MinCheck = get_global_parm(s_parm_name, s_parm_relval,"GP_MinCheck"); 
02433     GP_SLRise = get_global_parm(s_parm_name, s_parm_relval,"GP_SLRise"); 
02434     GP_dispLenRef = get_global_parm(s_parm_name, s_parm_relval,"GP_dispLenRef"); 
02435     GP_dispParm = get_global_parm(s_parm_name, s_parm_relval,"GP_dispParm"); 
02436  
02437  /* Periphyton/Algae */
02438     GP_ALG_IC_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_IC_MULT"); 
02439     GP_alg_uptake_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_uptake_coef"); 
02440     GP_ALG_SHADE_FACTOR = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_SHADE_FACTOR"); 
02441     GP_algMortDepth = get_global_parm(s_parm_name, s_parm_relval,"GP_algMortDepth"); 
02442     GP_ALG_RC_MORT_DRY = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_MORT_DRY"); 
02443     GP_ALG_RC_MORT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_MORT"); 
02444     GP_ALG_RC_PROD = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_PROD"); 
02445     GP_ALG_RC_RESP = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_RC_RESP"); 
02446     GP_alg_R_accel = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_R_accel"); 
02447     GP_AlgComp = get_global_parm(s_parm_name, s_parm_relval,"GP_AlgComp"); 
02448     GP_ALG_REF_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_REF_MULT"); 
02449     GP_NC_ALG_KS_P = get_global_parm(s_parm_name, s_parm_relval,"GP_NC_ALG_KS_P"); 
02450     GP_alg_alkP_min = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_alkP_min"); 
02451     GP_C_ALG_KS_P = get_global_parm(s_parm_name, s_parm_relval,"GP_C_ALG_KS_P"); 
02452     GP_ALG_TEMP_OPT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_TEMP_OPT"); 
02453     GP_C_ALG_threshTP = get_global_parm(s_parm_name, s_parm_relval,"GP_C_ALG_threshTP"); 
02454     GP_ALG_C_TO_OM = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_C_TO_OM"); 
02455     GP_alg_light_ext_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_alg_light_ext_coef"); 
02456     GP_ALG_LIGHT_SAT = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_LIGHT_SAT"); 
02457     GP_ALG_PC = get_global_parm(s_parm_name, s_parm_relval,"GP_ALG_PC"); 
02458 
02459 /* Soil */
02460     GP_DOM_RCDECOMP = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_RCDECOMP"); 
02461     GP_DOM_DECOMPRED = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_DECOMPRED"); 
02462     GP_calibDecomp = get_global_parm(s_parm_name, s_parm_relval,"GP_calibDecomp"); 
02463     GP_DOM_decomp_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_decomp_coef"); 
02464     GP_DOM_DECOMP_POPT = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_DECOMP_POPT"); 
02465     GP_DOM_DECOMP_TOPT = get_global_parm(s_parm_name, s_parm_relval,"GP_DOM_DECOMP_TOPT"); 
02466     GP_sorbToTP = get_global_parm(s_parm_name, s_parm_relval,"GP_sorbToTP"); 
02467     GP_IC_ELEV_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_ELEV_MULT"); 
02468     GP_IC_BATHY_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_BATHY_MULT"); 
02469     GP_IC_DOM_BD_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_DOM_BD_MULT"); 
02470     GP_IC_BulkD_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_BulkD_MULT"); 
02471     GP_IC_TPtoSOIL_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_IC_TPtoSOIL_MULT"); 
02472  
02473 /* Macrophytes */
02474     GP_MAC_IC_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_MAC_IC_MULT"); 
02475     GP_MAC_REFUG_MULT = get_global_parm(s_parm_name, s_parm_relval,"GP_MAC_REFUG_MULT"); 
02476     GP_mac_uptake_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_mac_uptake_coef"); 
02477     GP_mann_height_coef = get_global_parm(s_parm_name, s_parm_relval,"GP_mann_height_coef"); 
02478 
02479 /* Floc */
02480     GP_Floc_BD = get_global_parm(s_parm_name, s_parm_relval,"GP_Floc_BD"); 
02481     GP_FlocMax = get_global_parm(s_parm_name, s_parm_relval,"GP_FlocMax"); 
02482     GP_TP_P_OM = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_P_OM"); 
02483     GP_Floc_rcSoil = get_global_parm(s_parm_name, s_parm_relval,"GP_Floc_rcSoil"); 
02484                 
02485 /* Phosphorus */
02486     GP_TP_DIFFCOEF = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_DIFFCOEF"); 
02487     GP_TP_K_INTER = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_K_INTER"); 
02488     GP_TP_K_SLOPE = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_K_SLOPE"); 
02489     GP_WQMthresh = get_global_parm(s_parm_name, s_parm_relval,"GP_WQMthresh"); 
02490     GP_PO4toTP = get_global_parm(s_parm_name, s_parm_relval,"GP_PO4toTP"); 
02491     GP_TP_IN_RAIN = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_IN_RAIN"); 
02492     GP_PO4toTPint = get_global_parm(s_parm_name, s_parm_relval,"GP_PO4toTPint"); 
02493     GP_TP_ICSFWAT = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_ICSFWAT"); 
02494     GP_TP_ICSEDWAT = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_ICSEDWAT"); 
02495     GP_TPpart_thresh = get_global_parm(s_parm_name, s_parm_relval,"GP_TPpart_thresh"); 
02496     GP_TP_DIFFDEPTH = get_global_parm(s_parm_name, s_parm_relval,"GP_TP_DIFFDEPTH"); 
02497     GP_settlVel = get_global_parm(s_parm_name, s_parm_relval,"GP_settlVel"); 
02498     
02499 }
02500 
02504 void get_map_dims() 
02505 {
02506   read_map_dims("Elevation");
02507 }
02508 
02514 void alloc_memory() 
02515 {
02516   usrErr0("Allocating Memory...");  /* console message */
02517   
02518   ON_MAP = (unsigned char *) nalloc(sizeof(unsigned char)*(s0+2)*(s1+2),"ON_MAP");
02519   init_pvar(ON_MAP,NULL,'c',0.0);
02520 
02521   BCondFlow = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"BCondFlow");
02522   init_pvar(BCondFlow,NULL,'d',0.0);
02523   HAB = (unsigned char *) nalloc(sizeof(unsigned char)*(s0+2)*(s1+2),"HAB");
02524   init_pvar(HAB,NULL,'c',0.0);
02525   basn = (int *) nalloc(sizeof(int)*(s0+2)*(s1+2),"basn"); /* Basin/Indicator-Region map */
02526   init_pvar(basn,NULL,'d',0.0);
02527 
02528   ALG_INCID_LIGHT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_INCID_LIGHT");
02529   init_pvar(ALG_INCID_LIGHT,NULL,'f',0.0);
02530   ALG_LIGHT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_LIGHT_CF");
02531   init_pvar(ALG_LIGHT_CF,NULL,'f',0.0);
02532   ALG_LIGHT_EXTINCT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_LIGHT_EXTINCT");
02533   init_pvar(ALG_LIGHT_EXTINCT,NULL,'f',0.0);
02534   ALG_WAT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_WAT_CF");
02535   init_pvar(ALG_WAT_CF,NULL,'f',0.0);
02536   ALG_TEMP_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_TEMP_CF");
02537   init_pvar(ALG_TEMP_CF,NULL,'f',0.0);
02538   ALG_REFUGE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_REFUGE");
02539   init_pvar(ALG_REFUGE,NULL,'f',0.0);
02540   ALG_SAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_SAT");
02541   init_pvar(ALG_SAT,NULL,'f',0.0);
02542   ALG_TOT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ALG_TOT");
02543   init_pvar(ALG_TOT,NULL,'f',0.0);
02544 
02545   NC_ALG_AVAIL_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_AVAIL_MORT");
02546   init_pvar(NC_ALG_AVAIL_MORT,NULL,'f',0.0);
02547   NC_ALG_GPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_GPP");
02548   init_pvar(NC_ALG_GPP,NULL,'f',0.0);
02549   NC_ALG_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_MORT");
02550   init_pvar(NC_ALG_MORT,NULL,'f',0.0);
02551   NC_ALG_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_MORT_POT");
02552   init_pvar(NC_ALG_MORT_POT,NULL,'f',0.0);
02553   NC_ALG_NPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_NPP");
02554   init_pvar(NC_ALG_NPP,NULL,'f',0.0);
02555   NC_ALG_NUT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_NUT_CF");
02556   init_pvar(NC_ALG_NUT_CF,NULL,'f',0.0);
02557   NC_ALG_PROD_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_PROD_CF");
02558   init_pvar(NC_ALG_PROD_CF,NULL,'f',0.0);
02559   NC_ALG_RESP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_RESP");
02560   init_pvar(NC_ALG_RESP,NULL,'f',0.0);
02561   NC_ALG_RESP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_RESP_POT");
02562   init_pvar(NC_ALG_RESP_POT,NULL,'f',0.0);
02563   NC_ALG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG");
02564   init_pvar(NC_ALG,NULL,'f',0.0);
02565   C_ALG_AVAIL_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_AVAIL_MORT");
02566   init_pvar(C_ALG_AVAIL_MORT,NULL,'f',0.0);
02567   C_ALG_GPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_GPP");
02568   init_pvar(C_ALG_GPP,NULL,'f',0.0);
02569   C_ALG_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_MORT");
02570   init_pvar(C_ALG_MORT,NULL,'f',0.0);
02571   C_ALG_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_MORT_POT");
02572   init_pvar(C_ALG_MORT_POT,NULL,'f',0.0);
02573   C_ALG_NPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_NPP");
02574   init_pvar(C_ALG_NPP,NULL,'f',0.0);
02575   C_ALG_NUT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_NUT_CF");
02576   init_pvar(C_ALG_NUT_CF,NULL,'f',0.0);
02577   C_ALG_PROD_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_PROD_CF");
02578   init_pvar(C_ALG_PROD_CF,NULL,'f',0.0);
02579   C_ALG_RESP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_RESP");
02580   init_pvar(C_ALG_RESP,NULL,'f',0.0);
02581   C_ALG_RESP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_RESP_POT");
02582   init_pvar(C_ALG_RESP_POT,NULL,'f',0.0);
02583   C_ALG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG");
02584   init_pvar(C_ALG,NULL,'f',0.0);
02585   NC_ALG_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_P");
02586   init_pvar(NC_ALG_P,NULL,'f',0.0);
02587   C_ALG_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_P");
02588   init_pvar(C_ALG_P,NULL,'f',0.0);
02589   NC_ALG_GPP_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_GPP_P");
02590   init_pvar(NC_ALG_GPP_P,NULL,'f',0.0);
02591   C_ALG_GPP_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_GPP_P");
02592   init_pvar(C_ALG_GPP_P,NULL,'f',0.0);
02593   NC_ALG_MORT_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_MORT_P");
02594   init_pvar(NC_ALG_MORT_P,NULL,'f',0.0);
02595   C_ALG_MORT_P = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"C_ALG_MORT_P");
02596   init_pvar(C_ALG_MORT_P,NULL,'f',0.0);
02597   NC_ALG_PCrep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NC_ALG_PCrep");
02598   init_pvar(NC_ALG_PCrep,NULL,'f',0.0);
02599   C_ALG_PCrep = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"C_ALG_PCrep");
02600   init_pvar(C_ALG_PCrep,NULL,'f',0.0);
02601   NC_ALG_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"NC_ALG_PC");
02602   init_pvar(NC_ALG_PC,NULL,'l',0.0);
02603   C_ALG_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"C_ALG_PC");
02604   init_pvar(C_ALG_PC,NULL,'l',0.0);
02605 
02606   DEPOS_ORG_MAT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DEPOS_ORG_MAT");
02607   init_pvar(DEPOS_ORG_MAT,NULL,'l',0.0);
02608         
02609   DOM_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_Z");
02610   init_pvar(DOM_Z,NULL,'f',0.0);
02611   DOM_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_DECOMP");
02612   init_pvar(DOM_DECOMP,NULL,'f',0.0);
02613   DOM_DECOMP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_DECOMP_POT");
02614   init_pvar(DOM_DECOMP_POT,NULL,'f',0.0);
02615   DOM_fr_nphBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_fr_nphBio");
02616   init_pvar(DOM_fr_nphBio,NULL,'f',0.0);
02617   
02618   Floc_fr_phBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"Floc_fr_phBio");
02619   init_pvar(Floc_fr_phBio,NULL,'f',0.0);
02620   DOM_FR_FLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_FR_FLOC");
02621   init_pvar(DOM_FR_FLOC,NULL,'f',0.0);
02622   soil_MOIST_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"soil_MOIST_CF");
02623   init_pvar(soil_MOIST_CF,NULL,'f',0.0);
02624   DOM_QUALITY_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_QUALITY_CF");
02625   init_pvar(DOM_QUALITY_CF,NULL,'f',0.0);
02626   DOM_SED_AEROB_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_SED_AEROB_Z");
02627   init_pvar(DOM_SED_AEROB_Z,NULL,'f',0.0);
02628   DOM_SED_ANAEROB_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_SED_ANAEROB_Z");
02629   init_pvar(DOM_SED_ANAEROB_Z,NULL,'f',0.0);
02630   DOM_TEMP_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_TEMP_CF");
02631   init_pvar(DOM_TEMP_CF,NULL,'f',0.0);
02632   ELEVATION = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"ELEVATION");
02633   init_pvar(ELEVATION,NULL,'f',0.0);
02634   Bathymetry = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"Bathymetry");
02635   init_pvar(Bathymetry,NULL,'f',0.0);
02636   SED_ELEV = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SED_ELEV");
02637   init_pvar(SED_ELEV,NULL,'f',0.0);
02638   SED_INACT_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SED_INACT_Z");
02639   init_pvar(SED_INACT_Z,NULL,'f',0.0);
02640   DOP_FLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOP_FLOC");
02641   init_pvar(DOP_FLOC,NULL,'f',0.0);
02642   DOP_nphBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOP_nphBio");
02643   init_pvar(DOP_nphBio,NULL,'f',0.0);
02644 
02645   DOM_P_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DOM_P_OM");
02646   init_pvar(DOM_P_OM,NULL,'l',0.0);
02647 
02648   TPtoSOIL = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TPtoSOIL");
02649   init_pvar(TPtoSOIL,NULL,'l',0.0);
02650   TPtoSOIL_rep = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TPtoSOIL_rep");
02651   init_pvar(TPtoSOIL_rep,NULL,'l',0.0);
02652   TPtoVOL = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TPtoVOL");
02653   init_pvar(TPtoVOL,NULL,'l',0.0);
02654   TPtoVOL_rep = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TPtoVOL_rep");
02655   init_pvar(TPtoVOL_rep,NULL,'l',0.0);
02656   BulkD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"BulkD");
02657   init_pvar(BulkD,NULL,'f',0.0);
02658   DOM_BD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOM_BD");
02659   init_pvar(DOM_BD,NULL,'f',0.0);
02660   Inorg_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"Inorg_Z");
02661   init_pvar(Inorg_Z,NULL,'f',0.0);
02662   DIM = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DIM");
02663   init_pvar(DIM,NULL,'f',0.0);
02664   DOP_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"DOP_DECOMP");
02665   init_pvar(DOP_DECOMP,NULL,'f',0.0);
02666 
02667   DOP = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DOP");
02668   init_pvar(DOP,NULL,'l',0.0);
02669 
02670 /* placeholder for fire */
02671   FIREdummy = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FIREdummy");
02672   init_pvar(FIREdummy,NULL,'f',0.0);
02673         
02674   HYD_DOM_ACTWAT_PRES = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_DOM_ACTWAT_PRES");
02675   init_pvar(HYD_DOM_ACTWAT_PRES,NULL,'f',0.0);
02676   HYD_DOM_ACTWAT_VOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_DOM_ACTWAT_VOL");
02677   init_pvar(HYD_DOM_ACTWAT_VOL,NULL,'f',0.0);
02678   HYD_EVAP_CALC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_EVAP_CALC");
02679   init_pvar(HYD_EVAP_CALC,NULL,'f',0.0);
02680   HYD_MANNINGS_N = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_MANNINGS_N");
02681   init_pvar(HYD_MANNINGS_N,NULL,'f',0.0);
02682   HYD_RCCONDUCT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_RCCONDUCT");
02683   init_pvar(HYD_RCCONDUCT,NULL,'f',0.0);
02684   HYD_SAT_POT_TRANS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_SAT_POT_TRANS");
02685   init_pvar(HYD_SAT_POT_TRANS,NULL,'f',0.0);
02686   HYD_SED_WAT_VOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_SED_WAT_VOL");
02687   init_pvar(HYD_SED_WAT_VOL,NULL,'f',0.0);
02688   HYD_TOT_POT_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_TOT_POT_TRANSP");
02689   init_pvar(HYD_TOT_POT_TRANSP,NULL,'f',0.0);
02690   HydTotHd = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HydTotHd");
02691   init_pvar(HydTotHd,NULL,'f',0.0);
02692   HYD_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_TRANSP");
02693   init_pvar(HYD_TRANSP,NULL,'f',0.0);
02694   HYD_ET = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_ET");
02695   init_pvar(HYD_ET,NULL,'f',0.0);
02696   HYD_UNSAT_POT_TRANS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_UNSAT_POT_TRANS");
02697   init_pvar(HYD_UNSAT_POT_TRANS,NULL,'f',0.0);
02698   HYD_WATER_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"HYD_WATER_AVAIL");
02699   init_pvar(HYD_WATER_AVAIL,NULL,'f',0.0);
02700 
02701 /* sfwmm rainfall, stage, and pET mapped to elm grid */
02702   boundcond_rain = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"boundcond_rain");
02703   init_pvar(boundcond_rain,NULL,'f',0.0);
02704   boundcond_depth = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"boundcond_depth");
02705   init_pvar(boundcond_depth,NULL,'f',0.0);
02706   boundcond_ETp = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"boundcond_ETp");
02707   init_pvar(boundcond_ETp,NULL,'f',0.0);
02708 
02709   SAT_TO_UNSAT_FL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_TO_UNSAT_FL");
02710   init_pvar(SAT_TO_UNSAT_FL,NULL,'f',0.0);
02711   SAT_VS_UNSAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_VS_UNSAT");
02712   init_pvar(SAT_VS_UNSAT,NULL,'f',0.0);
02713   SAT_WATER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WATER");
02714   init_pvar(SAT_WATER,NULL,'f',0.0);
02715   SAT_WT_HEAD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WT_HEAD");
02716   init_pvar(SAT_WT_HEAD,NULL,'f',0.0);
02717   SAT_WT_RECHG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WT_RECHG");
02718   init_pvar(SAT_WT_RECHG,NULL,'f',0.0);
02719   SAT_WT_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAT_WT_TRANSP");
02720   init_pvar(SAT_WT_TRANSP,NULL,'f',0.0);
02721   SF_WT_EVAP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_EVAP");
02722   init_pvar(SF_WT_EVAP,NULL,'f',0.0);
02723   SF_WT_FROM_RAIN = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_FROM_RAIN");
02724   init_pvar(SF_WT_FROM_RAIN,NULL,'f',0.0);
02725   SF_WT_INFILTRATION = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_INFILTRATION");
02726   init_pvar(SF_WT_INFILTRATION,NULL,'f',0.0);
02727   SF_WT_POT_INF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_POT_INF");
02728   init_pvar(SF_WT_POT_INF,NULL,'f',0.0);
02729   SF_WT_TO_SAT_DOWNFLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SF_WT_TO_SAT_DOWNFLOW");
02730   init_pvar(SF_WT_TO_SAT_DOWNFLOW,NULL,'f',0.0);
02731   SFWT_VOL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SFWT_VOL");
02732   init_pvar(SFWT_VOL,NULL,'f',0.0);
02733   SURFACE_WAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SURFACE_WAT");
02734   init_pvar(SURFACE_WAT,NULL,'f',0.0);
02735   UNSAT_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_AVAIL");
02736   init_pvar(UNSAT_AVAIL,NULL,'f',0.0);
02737   UNSAT_CAP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_CAP");
02738   init_pvar(UNSAT_CAP,NULL,'f',0.0);
02739   UNSAT_DEPTH = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_DEPTH");
02740   init_pvar(UNSAT_DEPTH,NULL,'f',0.0);
02741   UNSAT_HYD_COND_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_HYD_COND_CF");
02742   init_pvar(UNSAT_HYD_COND_CF,NULL,'f',0.0);
02743   UNSAT_MOIST_PRP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_MOIST_PRP");
02744   init_pvar(UNSAT_MOIST_PRP,NULL,'f',0.0);
02745   UNSAT_PERC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_PERC");
02746   init_pvar(UNSAT_PERC,NULL,'f',0.0);
02747   UNSAT_TO_SAT_FL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_TO_SAT_FL");
02748   init_pvar(UNSAT_TO_SAT_FL,NULL,'f',0.0);
02749   UNSAT_TRANSP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_TRANSP");
02750   init_pvar(UNSAT_TRANSP,NULL,'f',0.0);
02751   UNSAT_WATER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_WATER");
02752   init_pvar(UNSAT_WATER,NULL,'f',0.0);
02753   UNSAT_WT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"UNSAT_WT_POT");
02754   init_pvar(UNSAT_WT_POT,NULL,'f',0.0);
02755 
02756   MAC_HEIGHT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_HEIGHT");
02757   init_pvar(MAC_HEIGHT,NULL,'f',0.0);
02758   MAC_LAI = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_LAI");
02759   init_pvar(MAC_LAI,NULL,'f',0.0);
02760   LAI_eff = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"LAI_eff");
02761   init_pvar(LAI_eff,NULL,'f',0.0);
02762   MAC_LIGHT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_LIGHT_CF");
02763   init_pvar(MAC_LIGHT_CF,NULL,'f',0.0);
02764   MAC_MAX_BIO = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_MAX_BIO");
02765   init_pvar(MAC_MAX_BIO,NULL,'f',0.0);
02766   MAC_NOPH_BIOMAS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_NOPH_BIOMAS");
02767   init_pvar(MAC_NOPH_BIOMAS,NULL,'f',0.0);
02768   MAC_NUT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_NUT_CF");
02769   init_pvar(MAC_NUT_CF,NULL,'f',0.0);
02770   MAC_PH_BIOMAS = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_PH_BIOMAS");
02771   init_pvar(MAC_PH_BIOMAS,NULL,'f',0.0);
02772   MAC_PROD_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_PROD_CF");
02773   init_pvar(MAC_PROD_CF,NULL,'f',0.0);
02774   MAC_REL_BIOM = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_REL_BIOM");
02775   init_pvar(MAC_REL_BIOM,NULL,'f',0.0);
02776   MAC_SALT_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_SALT_CF");
02777   init_pvar(MAC_SALT_CF,NULL,'f',0.0);
02778   MAC_TEMP_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_TEMP_CF");
02779   init_pvar(MAC_TEMP_CF,NULL,'f',0.0);
02780   MAC_TOT_BIOM = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_TOT_BIOM");
02781   init_pvar(MAC_TOT_BIOM,NULL,'f',0.0);
02782   MAC_WATER_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_WATER_CF");
02783   init_pvar(MAC_WATER_CF,NULL,'f',0.0);
02784   MAC_WATER_AVAIL_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"MAC_WATER_AVAIL_CF");
02785   init_pvar(MAC_WATER_AVAIL_CF,NULL,'f',0.0);
02786   NPHBIO_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_AVAIL");
02787   init_pvar(NPHBIO_AVAIL,NULL,'f',0.0);
02788   NPHBIO_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_MORT");
02789   init_pvar(NPHBIO_MORT,NULL,'f',0.0);
02790   NPHBIO_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_MORT_POT");
02791   init_pvar(NPHBIO_MORT_POT,NULL,'f',0.0);
02792   NPHBIO_REFUGE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_REFUGE");
02793   init_pvar(NPHBIO_REFUGE,NULL,'f',0.0);
02794   NPHBIO_SAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_SAT");
02795   init_pvar(NPHBIO_SAT,NULL,'f',0.0);
02796   NPHBIO_TRANSLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_TRANSLOC");
02797   init_pvar(NPHBIO_TRANSLOC,NULL,'f',0.0);
02798   NPHBIO_TRANSLOC_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"NPHBIO_TRANSLOC_POT");
02799   init_pvar(NPHBIO_TRANSLOC_POT,NULL,'f',0.0);
02800   PHBIO_AVAIL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_AVAIL");
02801   init_pvar(PHBIO_AVAIL,NULL,'f',0.0);
02802   PHBIO_MORT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_MORT");
02803   init_pvar(PHBIO_MORT,NULL,'f',0.0);
02804   PHBIO_MORT_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_MORT_POT");
02805   init_pvar(PHBIO_MORT_POT,NULL,'f',0.0);
02806   PHBIO_NPP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_NPP");
02807   init_pvar(PHBIO_NPP,NULL,'f',0.0);
02808   PHBIO_REFUGE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_REFUGE");
02809   init_pvar(PHBIO_REFUGE,NULL,'f',0.0);
02810   PHBIO_SAT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_SAT");
02811   init_pvar(PHBIO_SAT,NULL,'f',0.0); 
02812   PHBIO_TRANSLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_TRANSLOC");
02813   init_pvar(PHBIO_TRANSLOC,NULL,'f',0.0);
02814   PHBIO_TRANSLOC_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"PHBIO_TRANSLOC_POT");
02815   init_pvar(PHBIO_TRANSLOC_POT,NULL,'f',0.0);
02816 
02817   phbio_npp_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_npp_P");
02818   init_pvar(phbio_npp_P,NULL,'l',0.0);
02819   phbio_npp_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_npp_OM");
02820   init_pvar(phbio_npp_OM,NULL,'l',0.0);
02821   phbio_mort_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_mort_P");
02822   init_pvar(phbio_mort_P,NULL,'l',0.0);
02823   phbio_mort_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_mort_OM");
02824   init_pvar(phbio_mort_OM,NULL,'l',0.0);
02825   phbio_transl_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_transl_P");
02826   init_pvar(phbio_transl_P,NULL,'l',0.0);
02827   phbio_transl_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"phbio_transl_OM");
02828   init_pvar(phbio_transl_OM,NULL,'l',0.0);
02829   nphbio_transl_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_transl_P");
02830   init_pvar(nphbio_transl_P,NULL,'l',0.0);
02831   nphbio_transl_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_transl_OM");
02832   init_pvar(nphbio_transl_OM,NULL,'l',0.0);
02833   nphbio_mort_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_mort_P");
02834   init_pvar(nphbio_mort_P,NULL,'l',0.0);
02835   nphbio_mort_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"nphbio_mort_OM");
02836   init_pvar(nphbio_mort_OM,NULL,'l',0.0);
02837   mac_nph_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_P");
02838   init_pvar(mac_nph_P,NULL,'l',0.0);
02839   mac_nph_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_PC");
02840   init_pvar(mac_nph_PC,NULL,'l',0.0);
02841   mac_nph_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_OM");
02842   init_pvar(mac_nph_OM,NULL,'l',0.0);
02843   mac_nph_CtoOM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_CtoOM");
02844   init_pvar(mac_nph_CtoOM,NULL,'l',0.0);
02845   mac_ph_P = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_P");
02846   init_pvar(mac_ph_P,NULL,'l',0.0);
02847   mac_ph_PC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_PC");
02848   init_pvar(mac_ph_PC,NULL,'l',0.0);
02849   mac_ph_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_OM");
02850   init_pvar(mac_ph_OM,NULL,'l',0.0);
02851   mac_ph_CtoOM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_CtoOM");
02852   init_pvar(mac_ph_CtoOM,NULL,'l',0.0);
02853 
02854   mac_nph_PC_rep = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_nph_PC_rep");
02855   init_pvar(mac_nph_PC_rep,NULL,'l',0.0);
02856   mac_ph_PC_rep = (float *) nalloc(sizeof(double)*(s0+2)*(s1+2),"mac_ph_PC_rep");
02857   init_pvar(mac_ph_PC_rep,NULL,'l',0.0);
02858 
02859   TP_SED_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SED_WT");
02860   init_pvar(TP_SED_WT,NULL,'l',0.0);
02861   TP_SED_CONC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SED_CONC");
02862   init_pvar(TP_SED_CONC,NULL,'l',0.0);
02863   TP_SEDWT_CONCACT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SEDWT_CONCACT");
02864   init_pvar(TP_SEDWT_CONCACT,NULL,'l',0.0);
02865   TP_SF_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SF_WT");
02866   init_pvar(TP_SF_WT,NULL,'l',0.0);
02867   TP_SFWT_CONC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SFWT_CONC");
02868   init_pvar(TP_SFWT_CONC,NULL,'l',0.0);
02869   TP_SORB = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SORB");
02870   init_pvar(TP_SORB,NULL,'l',0.0);
02871   TP_SORBCONC = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SORBCONC");
02872   init_pvar(TP_SORBCONC,NULL,'l',0.0);
02873   TP_SED_WT_AZ = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_SED_WT_AZ");
02874   init_pvar(TP_SED_WT_AZ,NULL,'l',0.0);
02875   TP_Act_to_Tot = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"TP_Act_to_Tot");
02876   init_pvar(TP_Act_to_Tot,NULL,'l',0.0);
02877 
02878   TP_Act_to_TotRep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_Act_to_TotRep");
02879   init_pvar(TP_Act_to_TotRep,NULL,'f',0.0);
02880   TP_SORBCONC_rep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SORBCONC_rep");
02881   init_pvar(TP_SORBCONC_rep,NULL,'f',0.0);
02882 
02883   TP_DNFLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_DNFLOW");
02884   init_pvar(TP_DNFLOW,NULL,'f',0.0);
02885   TP_DNFLOW_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_DNFLOW_POT");
02886   init_pvar(TP_DNFLOW_POT,NULL,'f',0.0);
02887   TP_FR_RAIN = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_FR_RAIN");
02888   init_pvar(TP_FR_RAIN,NULL,'f',0.0);
02889   TP_K = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_K");
02890   init_pvar(TP_K,NULL,'f',0.0);
02891   TP_SED_MINER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SED_MINER");
02892   init_pvar(TP_SED_MINER,NULL,'f',0.0);
02893   TP_SEDWT_CONCACTMG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SEDWT_CONCACTMG");
02894   init_pvar(TP_SEDWT_CONCACTMG,NULL,'f',0.0);
02895   TP_SEDWT_UPTAKE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SEDWT_UPTAKE");
02896   init_pvar(TP_SEDWT_UPTAKE,NULL,'f',0.0);
02897   TP_SFWT_CONC_MG = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SFWT_CONC_MG");
02898   init_pvar(TP_SFWT_CONC_MG,NULL,'f',0.0);
02899   TP_SFWT_MINER = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SFWT_MINER");
02900   init_pvar(TP_SFWT_MINER,NULL,'f',0.0);
02901   TP_SFWT_UPTAK = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SFWT_UPTAK");
02902   init_pvar(TP_SFWT_UPTAK,NULL,'f',0.0);
02903   TP_settl = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_settl");
02904   init_pvar(TP_settl,NULL,'f',0.0);
02905   TP_SORB_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SORB_POT");
02906   init_pvar(TP_SORB_POT,NULL,'f',0.0);
02907   TP_SORBTION = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_SORBTION");
02908   init_pvar(TP_SORBTION,NULL,'f',0.0);
02909   TP_UPFLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_UPFLOW");
02910   init_pvar(TP_UPFLOW,NULL,'f',0.0);
02911   TP_UPFLOW_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_UPFLOW_POT");
02912   init_pvar(TP_UPFLOW_POT,NULL,'f',0.0);
02913   P_SUM_CELL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"P_SUM_CELL");
02914   init_pvar(P_SUM_CELL,NULL,'f',0.0);
02915 
02916   DINdummy = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"DINdummy");
02917   init_pvar(DINdummy,NULL,'l',0.0);
02918 
02919   WQMsettlVel = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"WQMsettlVel");
02920   init_pvar(WQMsettlVel,NULL,'f',0.0);
02921   TP_settlDays = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"TP_settlDays");
02922   init_pvar(TP_settlDays,NULL,'f',0.0);
02923 
02924   SAL_SED_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SAL_SED_WT");
02925   init_pvar(SAL_SED_WT,NULL,'l',0.0);
02926   SAL_SF_WT_mb = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SAL_SF_WT_mb");
02927   init_pvar(SAL_SF_WT_mb,NULL,'l',0.0);
02928   SALT_SED_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SALT_SED_WT");
02929   init_pvar(SALT_SED_WT,NULL,'l',0.0);
02930   SALT_SURF_WT = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"SALT_SURF_WT");
02931   init_pvar(SALT_SURF_WT,NULL,'l',0.0);
02932         
02933   SAL_SF_WT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SAL_SF_WT");
02934   init_pvar(SAL_SF_WT,NULL,'f',0.0);
02935   SALT_SED_TO_SF_FLOW = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_SED_TO_SF_FLOW");
02936   init_pvar(SALT_SED_TO_SF_FLOW,NULL,'f',0.0);
02937   SALT_SFWAT_DOWNFL = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_SFWAT_DOWNFL");
02938   init_pvar(SALT_SFWAT_DOWNFL,NULL,'f',0.0);
02939   SALT_SFWAT_DOWNFL_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SALT_SFWAT_DOWNFL_POT");
02940   init_pvar(SALT_SFWAT_DOWNFL_POT,NULL,'f',0.0);
02941 
02942   FLOC_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DECOMP");
02943   init_pvar(FLOC_DECOMP,NULL,'f',0.0);
02944   FLOC_DECOMP_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DECOMP_POT");
02945   init_pvar(FLOC_DECOMP_POT,NULL,'f',0.0);
02946   FLOC_DECOMP_QUAL_CF = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DECOMP_QUAL_CF");
02947   init_pvar(FLOC_DECOMP_QUAL_CF,NULL,'f',0.0);
02948   FLOC_Z = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_Z");
02949   init_pvar(FLOC_Z,NULL,'f',0.0);
02950   FLOC_DEPO = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DEPO");
02951   init_pvar(FLOC_DEPO,NULL,'f',0.0);
02952   FLOC_DEPO_POT = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_DEPO_POT");
02953   init_pvar(FLOC_DEPO_POT,NULL,'f',0.0);
02954   FLOC_FR_ALGAE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC_FR_ALGAE");
02955   init_pvar(FLOC_FR_ALGAE,NULL,'f',0.0);
02956   FLOC = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FLOC");
02957   init_pvar(FLOC,NULL,'f',0.0);
02958   FlocP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP");
02959   init_pvar(FlocP,NULL,'f',0.0);
02960   FlocP_DECOMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_DECOMP");
02961   init_pvar(FlocP_DECOMP,NULL,'f',0.0);
02962   FlocP_FR_ALGAE = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_FR_ALGAE");
02963   init_pvar(FlocP_FR_ALGAE,NULL,'f',0.0);
02964   FlocP_PhBio = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_PhBio");
02965   init_pvar(FlocP_PhBio,NULL,'f',0.0);
02966   FlocP_DEPO = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_DEPO");
02967   init_pvar(FlocP_DEPO,NULL,'f',0.0);
02968   FlocP_OMrep = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"FlocP_OMrep");
02969   init_pvar(FlocP_OMrep,NULL,'f',0.0);
02970   FlocP_OM = (double *) nalloc(sizeof(double)*(s0+2)*(s1+2),"FlocP_OM");
02971   init_pvar(FlocP_OM,NULL,'l',0.0);
02972 
02973   SOLRADGRD = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SOLRADGRD");
02974   init_pvar(SOLRADGRD,NULL,'f',0.0);
02975   SOLRAD274 = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"SOLRAD274");
02976   init_pvar(SOLRAD274,NULL,'f',0.0);
02977   AIR_TEMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"AIR_TEMP");
02978   init_pvar(AIR_TEMP,NULL,'f',0.0);
02979   H2O_TEMP = (float *) nalloc(sizeof(float)*(s0+2)*(s1+2),"H2O_TEMP");
02980   init_pvar(H2O_TEMP,NULL,'f',0.0);
02981 
02982 /* habitat-specific parameter arrays */
02983   HP_ALG_MAX = NULL;
02984 
02985   HP_DOM_MAXDEPTH = NULL;
02986   HP_DOM_AEROBTHIN = NULL;
02987 
02988   HP_TP_CONC_GRAD = NULL;
02989 
02990   HP_SALT_ICSEDWAT = NULL;
02991   HP_SALT_ICSFWAT = NULL;
02992 
02993   HP_PHBIO_MAX = NULL;
02994   HP_NPHBIO_MAX = NULL;
02995   HP_MAC_MAXHT = NULL;
02996   HP_NPHBIO_ROOTDEPTH = NULL;
02997   HP_MAC_MAXROUGH = NULL;
02998   HP_MAC_MINROUGH = NULL;
02999   HP_MAC_MAXLAI = NULL;
03000   HP_MAC_MAXCANOPCOND = NULL;
03001   HP_MAC_CANOPDECOUP = NULL;
03002   HP_MAC_TEMPOPT = NULL;
03003   HP_MAC_LIGHTSAT = NULL;
03004   HP_MAC_KSP = NULL;
03005   HP_PHBIO_RCNPP = NULL;
03006   HP_PHBIO_RCMORT = NULL;
03007   HP_MAC_WAT_TOLER = NULL;
03008   HP_MAC_SALIN_THRESH = NULL;
03009   HP_PHBIO_IC_CTOOM = NULL;
03010   HP_NPHBIO_IC_CTOOM = NULL;
03011   HP_PHBIO_IC_PC = NULL;
03012   HP_NPHBIO_IC_PC = NULL;
03013   HP_MAC_TRANSLOC_RC = NULL;
03014 
03015   HP_HYD_RCINFILT = NULL;
03016   HP_HYD_SPEC_YIELD = NULL;
03017   HP_HYD_POROSITY = NULL;
03018 
03019   HP_FLOC_IC = NULL;
03020   HP_FLOC_IC_CTOOM = NULL;
03021   HP_FLOC_IC_PC = NULL;
03022 
03023   HP_SfDepthLo = NULL;
03024   HP_SfDepthHi = NULL;
03025   HP_SfDepthInt = NULL;
03026   HP_PhosLo = NULL;
03027   HP_PhosHi = NULL;
03028   HP_PhosInt = NULL;
03029   HP_FireInt = NULL;
03030   
03031   usrErr("Done."); /* console message */
03032 
03033 }
03034 
03035 
03036 /*******
03037 the remaining functions are time-series graph interpolaters
03038 ******/
03039 
03047 float graph7(unsigned char y, float x)
03048 {
03049   float s;
03050   int ig=0, Np=10;
03057   while(1) {
03058   if (x <= g7[ig][0]) { if(ig>0) ig=ig-1; else return(g7[0][1+y]);}
03059   else if (x > g7[ig][0] && x <= g7[ig+1][0]) {
03060          s =   (g7[ig+1][1+y]-g7[ig][1+y])/
03061             (g7[ig+1][0]-g7[ig][0]);
03062          return(s * (x-g7[ig][0]) + g7[ig][1+y]); }
03063     else if (x > g7[ig+1][0]) { if(ig<Np) ig=ig+1; else return(g7[Np][1+y]);}
03064   }
03065 }
03066 
03074 float graph8(unsigned char y, float x)
03075 {
03076   float s;
03077   int ig=0, Np=10;
03084   while(1) {
03085   if (x <= g8[ig][0]) { if(ig>0) ig=ig-1; else return(g8[0][1+y]);}
03086   else if (x > g8[ig][0] && x <= g8[ig+1][0]) {
03087          s =   (g8[ig+1][1+y]-g8[ig][1+y])/
03088             (g8[ig+1][0]-g8[ig][0]);
03089          return(s * (x-g8[ig][0]) + g8[ig][1+y]); }
03090     else if (x > g8[ig+1][0]) { if(ig<Np) ig=ig+1; else return(g8[Np][1+y]);}
03091   }
03092 }
03093 
03094 
03095 /***** end UnitMod.c source ***/

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