#include <cstdio>
#include <cmath>
#include <iostream>
#include <string> 


#include "config.h"
#include "controller.h"
#include "input_files.h"
#include "grid_cell.h"
#include "zone.h"
#include "grid_info.h"

using namespace std;
using namespace glm;

/******************************************************************************/
GridCell::GridCell () {
   memset(this, 0, sizeof(GridCell));
   zdis = NZ - 1;
   Ldis = NL - 1;
}


/******************************************************************************/
void GridCell::init(size_t now, size_t y, size_t x, Zone* zone,
                    GridInfo& grid, InputFiles& inputs, Config& config) {
   
   this->y = y;
   this->x = x;
   lat = grid.lat[y];
   lon = grid.lon[x];
   cellArea = grid.cellArea[y][x];

   zoneTotal = zone;

   water = grid.wtr_->get(y, x);
   ice   = grid.ice_->get(y, x);
   if (config.shft_cult == 1) { 
      shft = inputs.shftFile_->get(y, x);
   }
   else {
      shft = grid.shiftCult[y][x];
   }
   
   if (config.BEST_CASE) {
      if (config.BEST_CASE_MIN_FLOWS_T5 && shft) {
         shiftingCultivation = true;
      } else {
         shiftingCultivation = false;
      }
   } else {
      if (config.adjust_smart_flow_option == 5
          && shft) {
         shiftingCultivation = true;
      } else {
         shiftingCultivation = false;
      }
   }
   
   potentialNPP = grid.vrgnNPPAboveGround[y][x];
   potentialBiomass = grid.vrgnBiomassAboveGround[y][x];
   rsForestLoss = grid.rs_forest_loss_2012[y][x];

   if (config.run_management == 1) {  
      croptype_ratio[C3ANN] = grid.C3_annual_input_2000[y][x];
      croptype_ratio[C4ANN] = grid.C4_annual_input_2000[y][x];
      croptype_ratio[C3PER] = grid.C3_perennial_input_2000[y][x];
      croptype_ratio[C4PER] = grid.C4_perennial_input_2000[y][x];
      croptype_ratio[C3NFX] = grid.N_fixing_input_2000[y][x];
      croptype_ratio_total = croptype_ratio[C3ANN] + croptype_ratio[C4ANN]
                           + croptype_ratio[C3PER] + croptype_ratio[C4PER]
                           + croptype_ratio[C3NFX];
                        
   }
      
   if (potentialBiomass > 2.0) {
      potentialForest = true;
   } else {
      potentialForest = false;
   }
   
   loadInputs(now, inputs, config);

   if (config.futureRun | !config.isInitialRun | config.extensionRun ){
     loadRestarts(now,grid,config);
   }
   
   disaggregateInitialLandUse(now, grid, config);

   if (config.run_management == 1) {  
    if (!config.futureRun & !config.extensionRun){
      computeIrrigatedArea (now);
      computeFertilizerAmount (now);
      computeWHBiofuelFractions (now);
      computeFloodedFraction (now);
      computeCroplandBiofuelFractions (now);
      }
      else if (config.extensionRun){
      computeIrrigatedAreaExtension (now);
      computeFertilizerAmountExtension (now);
      }
      computeTillageRate (now);
      computeFutureCroplandBiofuelFractions (now);
      computeBiomassHarvestFraction (now);

   }

}


/******************************************************************************/
const char* glm::getLandUseName(lu_types lu) {
   switch (lu) {
      case glm::VRGN: return "vrgn"; break;
      case glm::SCND: return "scnd"; break;
      case glm::CROP: return "crop"; break;
      case glm::PSTR: return "pstr"; break;
      case glm::URBN: return "urbn"; break;
      default: std::cerr << "unknown landuse type\n"; return "";
   }
}


/******************************************************************************/
string glm::getFlowName(size_t source, size_t target) {
   return string() + getLandUseName((lu_types)source)
                   + "2" + getLandUseName((lu_types)target);
}

/******************************************************************************/
const char* glm::getSubLandUseName(sub_lu_types lu) {
   switch (lu) {
      case glm::PRIMF:     return "primf"; break;
      case glm::PRIMN:     return "primn"; break;
      case glm::SECDF:     return "secdf"; break;
      case glm::SECDN:     return "secdn"; break;
      case glm::URBAN:     return "urban"; break;
      case glm::C3ANN:     return "c3ann"; break;
      case glm::C4ANN:     return "c4ann"; break;
      case glm::C3PER:     return "c3per"; break;
      case glm::C4PER:     return "c4per"; break;
      case glm::C3NFX:     return "c3nfx"; break;
      case glm::PASTR:     return "pastr"; break;
      case glm::RANGE:     return "range"; break;
      default: std::cerr << "unknown landuse type\n"; return "";
   }
}


/******************************************************************************/
const char* glm::getSubLandUseLongName(sub_lu_types lu) {
   switch (lu) {
      case glm::PRIMF:     return "primary vegetation on LUH2 potential forest land"; break;
      case glm::PRIMN:     return "primary vegetation on LUH2 potential non-forest land"; break;
      case glm::SECDF:     return "secondary vegetation on LUH2 potential forest land"; break;
      case glm::SECDN:     return "secondary vegetation on LUH2 potential non-forest land"; break;
      case glm::URBAN:     return "urban land"; break;
      case glm::C3ANN:     return "C3 annual crops"; break;
      case glm::C4ANN:     return "C4 annual crops"; break;
      case glm::C3PER:     return "C3 perennial crops"; break;
      case glm::C4PER:     return "C4 perennial crops"; break;
      case glm::C3NFX:     return "C3 nitrogen-fixing crops"; break;
      case glm::PASTR:     return "managed pasture"; break;
      case glm::RANGE:     return "rangeland"; break;
      default: std::cerr << "unknown landuse type\n"; return "";
   }
}


/******************************************************************************/
string glm::getSubFlowName(size_t source, size_t target) {
   return string() + getSubLandUseName((sub_lu_types)source)
                   + "_to_" + getSubLandUseName((sub_lu_types)target);
}

/******************************************************************************/
void GridCell::Flows::resetNegatives(bool warn, const char* caller) {

   for (size_t from=VRGN; from<N_LU_TYPES; ++from) {
      for (size_t to=VRGN; to<N_LU_TYPES; ++to) {
         if (transitions[from][to] < 0.0) {
            if (warn) {
               std::cerr << caller << " flow bug: "
                         << getFlowName(from, to)
                         << " " << transitions[from][to]
                         << "<0 -- resetting to 0" << std::endl;
            }
            transitions[from][to] = 0.0;
         }
      }
   }

   resetNegativeToZero (vrgnHarvArea1, warn, "flow bug: vrgnHarvArea1", caller);
   resetNegativeToZero (vrgnHarvArea2, warn, "flow bug: vrgnHarvArea2", caller);
   resetNegativeToZero (scndHarvArea1, warn, "flow bug: scndHarvArea1", caller);
   resetNegativeToZero (scndHarvArea2, warn, "flow bug: scndHarvArea1", caller);
   resetNegativeToZero (scndHarvArea3, warn, "flow bug: scndHarvArea1", caller);
}


/******************************************************************************/
void GridCell::Flows::resetNegativeSubs(bool warn, const char* caller) {

   for (size_t from=0; from<N_SUB_LU_TYPES; ++from) {
      for (size_t to=0; to<N_SUB_LU_TYPES; ++to) {
         if (subTransitions[from][to] < 0.0) {
            if (warn) {
               std::cerr << caller << " flow bug: "
                         << getSubFlowName(from, to)
                         << " " << subTransitions[from][to]
                         << "<0 -- resetting to 0" << std::endl;
            }
            transitions[from][to] = 0.0;
         }
      }
   }

   resetNegativeToZero (vrgnHarvArea1, warn, "flow bug: vrgnHarvArea1", caller);
   resetNegativeToZero (vrgnHarvArea2, warn, "flow bug: vrgnHarvArea2", caller);
   resetNegativeToZero (scndHarvArea1, warn, "flow bug: scndHarvArea1", caller);
   resetNegativeToZero (scndHarvArea2, warn, "flow bug: scndHarvArea1", caller);
   resetNegativeToZero (scndHarvArea3, warn, "flow bug: scndHarvArea1", caller);
}


/******************************************************************************/
double GridCell::harvestProbability (double biomass, InputFiles& inputs) {
   
   bool flagger = true;
   double p = 0.0;

   for (int i=0; i<N_PROB-1; i++) {
      if (biomass > (i * 0.5) && biomass <= ((i + 1) * 0.5)) {
         p = inputs.harvProbabilities[i];
         flagger = false;
         break;
      }
      if (i == 0 && biomass <= (i * 1.0)) {
         flagger = false;
         break;
      }
   }
   if (flagger) {
      cerr << "BUG: prob_harv .not. set. biomass: " << biomass;
      p = inputs.harvProbabilities[N_PROB-1];
   }

   return p;
}


/******************************************************************************/
void GridCell::initNextStatesAndFlows (size_t next) {
   memset(&frac[next], 0, sizeof(State));
   memset(&flow, 0, sizeof(Flows));
}


/******************************************************************************/
void GridCell::loadInputs (size_t t, InputFiles& inputs, Config& config) {
   
   frac[t].lu[CROP] = inputs.cropFile_->get(y, x);
   frac[t].lu[PSTR] = inputs.pstrFile_->get(y, x);
   frac[t].lu[VRGN] = inputs.othrFile_->get(y, x);
   frac[t].lu[URBN] = inputs.urbnFile_->get(y, x);
   if ((!config.futureRun)){
   frac[t].subLU[PASTR] = inputs.mpstFile_->get(y, x);
   frac[t].subLU[RANGE] = inputs.rangFile_->get(y, x);
   }
   
   if (config.shft_cult == 1) {
      frac[t].shftCult = inputs.shftFile_->get(y, x);
   }
   if (config.run_management == 1){
      frac[t].irrgData = inputs.irrgFile_->get(y, x);
      frac[t].riceData = inputs.riceFile_->get(y, x);
   }   
   if (config.futureRun == 1 | config.extensionRun==1) {
      frac[t].iamFlood = inputs.iamFloodFile_->get(y, x);
      frac[t].iamC4P = inputs.c4pFile_->get(y, x);
      frac[t].iamC3P = inputs.c3pFile_->get(y, x);
      frac[t].iamC3A = inputs.c3aFile_->get(y, x);
      frac[t].iamC4A = inputs.c4aFile_->get(y, x);
      frac[t].iamC3N = inputs.c3nFile_->get(y, x);
      frac[t].iamFert[C3ANN] = inputs.iamC3AFertValsFile_->get(y, x);
      frac[t].iamFert[C3PER] = inputs.iamC3PFertValsFile_->get(y, x);
      frac[t].iamFert[C4ANN] = inputs.iamC4AFertValsFile_->get(y, x);
      frac[t].iamFert[C4PER] = inputs.iamC4PFertValsFile_->get(y, x);
      frac[t].iamFert[C3NFX] = inputs.iamC3NFertValsFile_->get(y, x);
      frac[t].iamIrrig[C3ANN] = inputs.iamC3AIrrigValsFile_->get(y, x);
      frac[t].iamIrrig[C3PER] = inputs.iamC3PIrrigValsFile_->get(y, x);
      frac[t].iamIrrig[C4ANN] = inputs.iamC4AIrrigValsFile_->get(y, x);
      frac[t].iamIrrig[C4PER] = inputs.iamC4PIrrigValsFile_->get(y, x);
      frac[t].iamIrrig[C3NFX] = inputs.iamC3NIrrigValsFile_->get(y, x);
      
      frac[t].iamCropBiofuel = inputs.iamFutureCropBiofuelFile_->get(y, x);
      frac[t].iamC3PBiofuel = inputs.iamC3PBiofuelFile_->get(y, x);
      frac[t].iamC4PBiofuel = inputs.iamC4PBiofuelFile_->get(y, x);
   }
 
}

/******************************************************************************/
void GridCell::loadRestarts (size_t t, GridInfo& grid, Config& config) {

  frac[t].lu[CROP] = grid.restartC3Ann_->get(y,x) + grid.restartC3Per_->get(y,x) + grid.restartC4Ann_->get(y,x) + grid.restartC4Per_->get(y,x) + grid.restartC3NFx_->get(y,x);
  if (frac[t].lu[CROP] > 2.0 | frac[t].lu[CROP] < 0.0){
     frac[t].lu[CROP] = 0;
  }
  frac[t].lu[PSTR] = grid.restartPastr_->get(y,x) + grid.restartRange_->get(y,x);
  if (frac[t].lu[PSTR] > 2.0 | frac[t].lu[PSTR] < 0.0){
     frac[t].lu[PSTR] = 0;
  }
  frac[t].lu[SCND] = grid.restartScndF_->get(y,x) + grid.restartScndN_->get(y,x);
  if (frac[t].lu[SCND] > 2.0 | frac[t].lu[SCND] < 0.0){
     frac[t].lu[SCND] = 0;
  }
  frac[t].lu[VRGN] = grid.restartPrimF_->get(y,x) + grid.restartPrimN_->get(y,x);
  if (frac[t].lu[VRGN] > 2.0 | frac[t].lu[VRGN] < 0.0){
     frac[t].lu[VRGN] = 0;
  }
  frac[t].lu[URBN] = grid.restartUrban_->get(y,x);
  if (frac[t].lu[URBN] > 2.0 | frac[t].lu[URBN] < 0.0){
     frac[t].lu[URBN] = 0;
  }
  frac[t].scndMeanAge = grid.restartSMA_->get(y,x);
  if (frac[t].scndMeanAge > 1e19 | frac[t].scndMeanAge < 0.0){
     frac[t].scndMeanAge = 0;
  }
  frac[t].scndMeanBiomass = grid.restartSMB_->get(y,x);
  if (frac[t].scndMeanBiomass > 1e19 | frac[t].scndMeanBiomass < 0.0){
     frac[t].scndMeanBiomass = 0;
  }
  
  frac[t].fert[C3ANN] = grid.restartC3AnnFert_->get(y,x);
  if (frac[t].fert[C3ANN] > 1e19 | frac[t].fert[C3ANN] < 0.0){
     frac[t].fert[C3ANN] = 0;
  }
  frac[t].fert[C4ANN] = grid.restartC4AnnFert_->get(y,x);
  if (frac[t].fert[C4ANN] > 1e19 | frac[t].fert[C4ANN] < 0.0){
     frac[t].fert[C4ANN] = 0;
  }
  frac[t].fert[C3PER] = grid.restartC3PerFert_->get(y,x);
  if (frac[t].fert[C3PER] > 1e19 | frac[t].fert[C3PER] < 0.0){
     frac[t].fert[C3PER] = 0;
  }
  frac[t].fert[C4PER] = grid.restartC4PerFert_->get(y,x);
  if (frac[t].fert[C4PER] > 1e19 | frac[t].fert[C4PER] < 0.0){
     frac[t].fert[C4PER] = 0;
  }
  frac[t].fert[C3NFX] = grid.restartC3NfxFert_->get(y,x);
  if (frac[t].fert[C3NFX] > 1e19 | frac[t].fert[C3NFX] < 0.0){
     frac[t].fert[C3NFX] = 0;
  }
  
  //irrigation
  frac[t].irrig[C3ANN] = grid.restartC3AnnIrrig_->get(y,x);
  if (frac[t].irrig[C3ANN] > 1e19 | frac[t].irrig[C3ANN] < 0.0){
     frac[t].irrig[C3ANN] = 0;
  }
  frac[t].irrig[C4ANN] = grid.restartC4AnnIrrig_->get(y,x);
  if (frac[t].irrig[C4ANN] > 1e19 | frac[t].irrig[C4ANN] < 0.0){
     frac[t].irrig[C4ANN] = 0;
  }
  frac[t].irrig[C3PER] = grid.restartC3PerIrrig_->get(y,x);
  if (frac[t].irrig[C3PER] > 1e19 | frac[t].irrig[C3PER] < 0.0){
     frac[t].irrig[C3PER] = 0;
  }
  frac[t].irrig[C4PER] = grid.restartC4PerIrrig_->get(y,x);
  if (frac[t].irrig[C4PER] > 1e19 | frac[t].irrig[C4PER] < 0.0){
     frac[t].irrig[C4PER] = 0;
  }
  frac[t].irrig[C3NFX] = grid.restartC3NfxIrrig_->get(y,x);
  if (frac[t].irrig[C3NFX] > 1e19 | frac[t].irrig[C3NFX] < 0.0){
     frac[t].irrig[C3NFX] = 0;
  }
  frac[t].flooded = grid.restartFlood_->get(y,x);
  if (frac[t].flooded > 1e19 | frac[t].flooded < 0.0){
     frac[t].flooded = 0;
  }
  
  //biofuels
  frac[t].cropBiofuels[C3ANN] = grid.restartCropbfC3Ann_->get(y,x);
  if (frac[t].cropBiofuels[C3ANN] > 1e19 | frac[t].cropBiofuels[C3ANN] < 0.0){
     frac[t].cropBiofuels[C3ANN] = 0;
  }
  frac[t].cropBiofuels[C4ANN] = grid.restartCropbfC4Ann_->get(y,x);
  if (frac[t].cropBiofuels[C4ANN] > 1e19 | frac[t].cropBiofuels[C4ANN] < 0.0){
     frac[t].cropBiofuels[C4ANN] = 0;
  }
  frac[t].cropBiofuels[C3PER] = grid.restartCropbfC3Per_->get(y,x);
  if (frac[t].cropBiofuels[C3PER] > 1e19 | frac[t].cropBiofuels[C3PER] < 0.0){
     frac[t].cropBiofuels[C3PER] = 0;
  }
  frac[t].cropBiofuels[C4PER] = grid.restartCropbfC4Per_->get(y,x);
  if (frac[t].cropBiofuels[C4PER] > 1e19 | frac[t].cropBiofuels[C4PER] < 0.0){
     frac[t].cropBiofuels[C4PER] = 0;
  }
  frac[t].cropBiofuels[C3NFX] = grid.restartCropbfC3Nfx_->get(y,x);
  if (frac[t].cropBiofuels[C3NFX] > 1e19 | frac[t].cropBiofuels[C3NFX] < 0.0){
     frac[t].cropBiofuels[C3NFX] = 0;
  }
  
  //harvested fraction
  frac[t].biomassHarvestFraction[C3PER] = grid.restartFharvC3Per_->get(y,x);
  if (frac[t].biomassHarvestFraction[C3PER] > 1e19 | frac[t].biomassHarvestFraction[C3PER] < 0.0){
     frac[t].biomassHarvestFraction[C3PER] = 0;
  }
  frac[t].biomassHarvestFraction[C4PER] = grid.restartFharvC4Per_->get(y,x);
  if (frac[t].biomassHarvestFraction[C4PER] > 1e19 | frac[t].biomassHarvestFraction[C4PER] < 0.0){
     frac[t].biomassHarvestFraction[C4PER] = 0;
  }
  
  //fate of wood harvest
  frac[t].industrialRoundwoodFraction = grid.restartRndwd_->get(y,x);
  if (frac[t].industrialRoundwoodFraction > 1e19 | frac[t].industrialRoundwoodFraction < 0.0){
     frac[t].industrialRoundwoodFraction = 0;
  }
  frac[t].fuelwoodFraction = grid.restartFulwd_->get(y,x);
  if (frac[t].fuelwoodFraction > 1e19 | frac[t].fuelwoodFraction < 0.0){
     frac[t].fuelwoodFraction = 0;
  }
  
  frac[t].commercialBiofuelsFraction = grid.restartCombf_->get(y,x);
  if (frac[t].commercialBiofuelsFraction > 1e19 | frac[t].commercialBiofuelsFraction < 0.0){
     frac[t].commercialBiofuelsFraction = 0;
  }
  
  
  
  
}

/******************************************************************************/
void GridCell::transitions (size_t now, size_t next,
                            InputFiles& inputs, Config& config,
                            clock_t& sft, clock_t& asft, int year, GridInfo& grid) {

   clock_t start = clock();
   initNextStatesAndFlows(next);

   loadInputs(next, inputs, config);

   urbanFlows(now, next);

   alternativeSmartFlow(now, next, config.SMART_FLOW_BUG_PRINT);

   if (shiftingCultivation) {
      adjustSmartFlow(now, config.shft_cult, grid, config);
   }

   sft += clock() - start;
   start = clock();
   asft += clock() - start;
   
   updateForestLoss (next, now, year);

   updateBiomassTallies (now, inputs, year);

}

/******************************************************************************/
void GridCell::cropRotations (size_t now, size_t next) {
   // Compute crop rotations i.e. area transitioning from one crop functional type to another. 
   // Crop rotations do not involve creation of new croplands or abandonment of existing rotations.

   double next_rotations_c3nfx = frac[now].subLU[C3NFX] - flow.getSub(C3NFX,SECDF) - flow.getSub(C3NFX,SECDN) - 
                                 flow.getSub(C3NFX,PASTR) - flow.getSub(C3NFX,RANGE) - flow.getSub(C3NFX,URBAN) -
                                 flow.getSub(C3NFX,C4PER);
   double next_rotations_c4ann = frac[now].subLU[C4ANN] - flow.getSub(C4ANN,SECDF) - flow.getSub(C4ANN,SECDN) - 
                                 flow.getSub(C4ANN,PASTR) - flow.getSub(C4ANN,RANGE) - flow.getSub(C4ANN,URBAN) -
                                 flow.getSub(C4ANN,C4PER);
   double next_rotations_c3ann = frac[now].subLU[C3ANN] - flow.getSub(C3ANN,SECDF) - flow.getSub(C3ANN,SECDN) - 
                                 flow.getSub(C3ANN,PASTR) - flow.getSub(C3ANN,RANGE) - flow.getSub(C3ANN,URBAN) -
                                 flow.getSub(C3ANN,C4PER);

   // Lower area CFT determines total amount that will undergo rotation
   double min_area = std::min(next_rotations_c4ann, next_rotations_c3nfx);
   flow.setSub(C4ANN, C3NFX, zoneTotal->c4ann_to_c3nfx * min_area);
   flow.setSub(C3NFX, C4ANN, zoneTotal->c4ann_to_c3nfx * min_area);
   
   min_area = std::min(next_rotations_c4ann, next_rotations_c3ann);
   flow.setSub(C3ANN, C4ANN, zoneTotal->c4ann_to_c3ann * min_area);
   flow.setSub(C4ANN, C3ANN, zoneTotal->c4ann_to_c3ann * min_area);   

   min_area = std::min(next_rotations_c3ann, next_rotations_c3nfx); 
   flow.setSub(C3ANN, C3NFX, zoneTotal->c3ann_to_c3nfx * min_area);
   flow.setSub(C3NFX, C3ANN, zoneTotal->c3ann_to_c3nfx * min_area);      
}

/******************************************************************************/
void GridCell::expandC4Perennials (size_t now, size_t next) {
   // Compute crop rotations i.e. area transitioning from one crop functional type to another. 
   // Crop rotations do not involve creation of new croplands or abandonment of existing rotations.

   double next_rotations_c3nfx = frac[now].subLU[C3NFX] - flow.getSub(C3NFX,SECDF) - flow.getSub(C3NFX,SECDN) - 
                                 flow.getSub(C3NFX,PASTR) - flow.getSub(C3NFX,RANGE) - flow.getSub(C3NFX,URBAN);
   double next_rotations_c4ann = frac[now].subLU[C4ANN] - flow.getSub(C4ANN,SECDF) - flow.getSub(C4ANN,SECDN) - 
                                 flow.getSub(C4ANN,PASTR) - flow.getSub(C4ANN,RANGE) - flow.getSub(C4ANN,URBAN);
   double next_rotations_c3ann = frac[now].subLU[C3ANN] - flow.getSub(C3ANN,SECDF) - flow.getSub(C3ANN,SECDN) - 
                                 flow.getSub(C3ANN,PASTR) - flow.getSub(C3ANN,RANGE) - flow.getSub(C3ANN,URBAN);
   double next_rotations_c4per = frac[now].subLU[C4PER] - flow.getSub(C4PER,SECDF) - flow.getSub(C4PER,SECDN) - 
                                 flow.getSub(C4PER,PASTR) - flow.getSub(C4PER,RANGE) - flow.getSub(C4PER,URBAN) +
                                 flow.getSub(SECDF,C4PER) + flow.getSub(SECDN,C4PER) + flow.getSub(PRIMF,C4PER) + 
                                 flow.getSub(PRIMN,C4PER) + flow.getSub(RANGE,C4PER) + flow.getSub(PASTR,C4PER) +
                                 flow.getSub(URBAN,C4PER);

   double LUHDeltaC4P = next_rotations_c4per - frac[now].subLU[C4PER];
   double IAMDeltaC4P = frac[next].iamC4P - frac[now].iamC4P;
   
   double c4annLoss = 0.0, c3annLoss = 0.0, c3nfxLoss = 0.0;
   
   if (IAMDeltaC4P>0.000001 & croptype_ratio_total>0){
   if (LUHDeltaC4P>=IAMDeltaC4P){
   // do nothing
   }
   else {
      double targetDelta = IAMDeltaC4P - LUHDeltaC4P;
      if ( (next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx) >= targetDelta ){
      // enough crops to match IAMDeltaC4P
         c3annLoss = targetDelta*next_rotations_c3ann/(next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx);
         c4annLoss = targetDelta*next_rotations_c4ann/(next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx);
         c3nfxLoss = targetDelta*next_rotations_c3nfx/(next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx);
         
         flow.setSub(C3ANN,C4PER,c3annLoss);
         flow.setSub(C4ANN,C4PER,c4annLoss);
         flow.setSub(C3NFX,C4PER,c3nfxLoss);
      }
      else {
      // convert all of other annuals to C4P, but will be less than IAMDeltaC4P
         flow.setSub(C3ANN,C4PER,next_rotations_c3ann);
         flow.setSub(C4ANN,C4PER,next_rotations_c4ann);
         flow.setSub(C3NFX,C4PER,next_rotations_c3nfx);
      }
   }
   }
     
}

/******************************************************************************/
void GridCell::expandC4Annuals (size_t now, size_t next) {
   // Compute crop rotations i.e. area transitioning from one crop functional type to another. 
   // Crop rotations do not involve creation of new croplands or abandonment of existing rotations.

   double next_rotations_c3nfx = frac[now].subLU[C3NFX] - flow.getSub(C3NFX,SECDF) - flow.getSub(C3NFX,SECDN) - 
                                 flow.getSub(C3NFX,PASTR) - flow.getSub(C3NFX,RANGE) - flow.getSub(C3NFX,URBAN);
   double next_rotations_c4ann = frac[now].subLU[C4ANN] - flow.getSub(C4ANN,SECDF) - flow.getSub(C4ANN,SECDN) - 
                                 flow.getSub(C4ANN,PASTR) - flow.getSub(C4ANN,RANGE) - flow.getSub(C4ANN,URBAN)+
                                 flow.getSub(SECDF,C4ANN) + flow.getSub(SECDN,C4ANN) + flow.getSub(PRIMF,C4ANN) + 
                                 flow.getSub(PRIMN,C4ANN) + flow.getSub(RANGE,C4ANN) + flow.getSub(PASTR,C4ANN) +
                                 flow.getSub(URBAN,C4ANN);
   double next_rotations_c3ann = frac[now].subLU[C3ANN] - flow.getSub(C3ANN,SECDF) - flow.getSub(C3ANN,SECDN) - 
                                 flow.getSub(C3ANN,PASTR) - flow.getSub(C3ANN,RANGE) - flow.getSub(C3ANN,URBAN);
  
   double LUHDeltaC4A = next_rotations_c4ann - frac[now].subLU[C4ANN];
   double IAMDeltaC4A = frac[next].iamC4A - frac[now].iamC4A;
   
   double c3annLoss = 0.0, c3nfxLoss = 0.0;
   
   if (IAMDeltaC4A>0.000001 & croptype_ratio_total>0){
   if (LUHDeltaC4A>=IAMDeltaC4A){
   // do nothing
   }
   else {
      double targetDelta = IAMDeltaC4A - LUHDeltaC4A;
      if ( (next_rotations_c3ann + next_rotations_c3nfx) >= targetDelta ){
      // enough crops to match IAMDeltaC4P
         c3annLoss = targetDelta*next_rotations_c3ann/(next_rotations_c3ann + next_rotations_c3nfx);
         c3nfxLoss = targetDelta*next_rotations_c3nfx/(next_rotations_c3ann + next_rotations_c3nfx);
         
         flow.setSub(C3ANN,C4ANN,c3annLoss);
         flow.setSub(C3NFX,C4ANN,c3nfxLoss);
      }
      else {
      // convert all of other annuals to C4P, but will be less than IAMDeltaC4P
         flow.setSub(C3ANN,C4ANN,next_rotations_c3ann);
         flow.setSub(C3NFX,C4ANN,next_rotations_c3nfx);
      }
   }
   }
   
}

/******************************************************************************/
void GridCell::expandC3NFixing (size_t now, size_t next) {
   // Compute crop rotations i.e. area transitioning from one crop functional type to another. 
   // Crop rotations do not involve creation of new croplands or abandonment of existing rotations.

   double next_rotations_c3nfx = frac[now].subLU[C3NFX] - flow.getSub(C3NFX,SECDF) - flow.getSub(C3NFX,SECDN) - 
                                 flow.getSub(C3NFX,PASTR) - flow.getSub(C3NFX,RANGE) - flow.getSub(C3NFX,URBAN) +
                                 flow.getSub(SECDF,C3NFX) + flow.getSub(SECDN,C3NFX) + flow.getSub(PRIMF,C3NFX) + 
                                 flow.getSub(PRIMN,C3NFX) + flow.getSub(RANGE,C3NFX) + flow.getSub(PASTR,C3NFX) +
                                 flow.getSub(URBAN,C3NFX) - flow.getSub(C3NFX,C4ANN);
   double next_rotations_c4ann = frac[now].subLU[C4ANN] - flow.getSub(C4ANN,SECDF) - flow.getSub(C4ANN,SECDN) - 
                                 flow.getSub(C4ANN,PASTR) - flow.getSub(C4ANN,RANGE) - flow.getSub(C4ANN,URBAN);
                                 
   double next_rotations_c3ann = frac[now].subLU[C3ANN] - flow.getSub(C3ANN,SECDF) - flow.getSub(C3ANN,SECDN) - 
                                 flow.getSub(C3ANN,PASTR) - flow.getSub(C3ANN,RANGE) - flow.getSub(C3ANN,URBAN) - 
                                 flow.getSub(C3ANN,C4ANN);
  
   double LUHDeltaC3N = next_rotations_c3nfx - frac[now].subLU[C3NFX];
   double IAMDeltaC3N = frac[next].iamC3N - frac[now].iamC3N;
   
   double c3annLoss = 0.0, c4annLoss = 0.0;
   
   if (IAMDeltaC3N>0.000001 & croptype_ratio_total>0){
   if (LUHDeltaC3N>=IAMDeltaC3N){
   // do nothing
   }
   else {
      double targetDelta = IAMDeltaC3N - LUHDeltaC3N;
      if ( (next_rotations_c3ann + next_rotations_c4ann) >= targetDelta ){
      // enough crops to match IAMDeltaC4P
         c3annLoss = targetDelta*next_rotations_c3ann/(next_rotations_c3ann + next_rotations_c3nfx);
         c4annLoss = targetDelta*next_rotations_c4ann/(next_rotations_c3ann + next_rotations_c4ann);
         
         flow.setSub(C3ANN,C3NFX,c3annLoss);
         flow.setSub(C4ANN,C3NFX,c4annLoss);
      }
      else {
      // convert all of other annuals to C4P, but will be less than IAMDeltaC4P
         flow.setSub(C3ANN,C3NFX,next_rotations_c3ann);
         flow.setSub(C4ANN,C3NFX,next_rotations_c4ann);
      }
   }
   }
    
}

/******************************************************************************/
void GridCell::expandC3Perennials (size_t now, size_t next) {
   // Compute crop rotations i.e. area transitioning from one crop functional type to another. 
   // Crop rotations do not involve creation of new croplands or abandonment of existing rotations.

   double next_rotations_c3nfx = frac[now].subLU[C3NFX] - flow.getSub(C3NFX,SECDF) - flow.getSub(C3NFX,SECDN) - 
                                 flow.getSub(C3NFX,PASTR) - flow.getSub(C3NFX,RANGE) - flow.getSub(C3NFX,URBAN) - 
                                 flow.getSub(C3NFX,C4PER);
   double next_rotations_c4ann = frac[now].subLU[C4ANN] - flow.getSub(C4ANN,SECDF) - flow.getSub(C4ANN,SECDN) - 
                                 flow.getSub(C4ANN,PASTR) - flow.getSub(C4ANN,RANGE) - flow.getSub(C4ANN,URBAN) - 
                                 flow.getSub(C4ANN,C4PER);
   double next_rotations_c3ann = frac[now].subLU[C3ANN] - flow.getSub(C3ANN,SECDF) - flow.getSub(C3ANN,SECDN) - 
                                 flow.getSub(C3ANN,PASTR) - flow.getSub(C3ANN,RANGE) - flow.getSub(C3ANN,URBAN) -
                                 flow.getSub(C3ANN,C4PER);
   double next_rotations_c3per = frac[now].subLU[C3PER] - flow.getSub(C3PER,SECDF) - flow.getSub(C3PER,SECDN) - 
                                 flow.getSub(C3PER,PASTR) - flow.getSub(C3PER,RANGE) - flow.getSub(C3PER,URBAN) +
                                 flow.getSub(SECDF,C3PER) + flow.getSub(SECDN,C3PER) + flow.getSub(PRIMF,C3PER) + 
                                 flow.getSub(PRIMN,C3PER) + flow.getSub(RANGE,C3PER) + flow.getSub(PASTR,C3PER) +
                                 flow.getSub(URBAN,C3PER);

   double LUHDeltaC3P = next_rotations_c3per - frac[now].subLU[C3PER];
   double IAMDeltaC3P = frac[next].iamC3P - frac[now].iamC3P;
   
   double c4annLoss = 0.0, c3annLoss = 0.0, c3nfxLoss = 0.0;
   
   if (IAMDeltaC3P>0.000001 & croptype_ratio_total>0){
   if (LUHDeltaC3P>=IAMDeltaC3P){
   // do nothing
   }
   else {
      double targetDelta = IAMDeltaC3P - LUHDeltaC3P;
      if ( (next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx) >= targetDelta ){
      // enough crops to match IAMDeltaC4P
         c3annLoss = targetDelta*next_rotations_c3ann/(next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx);
         c4annLoss = targetDelta*next_rotations_c4ann/(next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx);
         c3nfxLoss = targetDelta*next_rotations_c3nfx/(next_rotations_c4ann + next_rotations_c3ann + next_rotations_c3nfx);
         
         flow.setSub(C3ANN,C3PER,c3annLoss);
         flow.setSub(C4ANN,C3PER,c4annLoss);
         flow.setSub(C3NFX,C3PER,c3nfxLoss);
      }
      else {
      // convert all of other annuals to C4P, but will be less than IAMDeltaC4P
         flow.setSub(C3ANN,C3PER,next_rotations_c3ann);
         flow.setSub(C4ANN,C3PER,next_rotations_c4ann);
         flow.setSub(C3NFX,C3PER,next_rotations_c3nfx);
      }
   }
   }
       
}

/******************************************************************************/
void GridCell::urbanFlows (size_t now, size_t next) {
   // Compute flows between land use types depending on the change in urban area
   // between current year and next
   double deltaU = frac[next].lu[URBN] - frac[now].lu[URBN];
   
   if (fabs(deltaU) < ZEROVALUE_CHECK) {
     // Do nothing if urban area does not change
   }
   else if (deltaU < 0.0) {
   
     // If urban area decreases, determine the land use change from urban to
     // crop, pasture and secondary. The land use change is proportional to the
     // area of the land use class.
     flow.set(URBN, CROP, frac[next].lu[CROP] / (frac[next].lu[CROP] + frac[next].lu[PSTR] + frac[next].lu[SCND] + frac[next].lu[VRGN])* (-1 * deltaU));
     flow.set(URBN, PSTR, frac[next].lu[PSTR] / (frac[next].lu[CROP] + frac[next].lu[PSTR] + frac[next].lu[SCND] + frac[next].lu[VRGN])* (-1 * deltaU));
     flow.set(URBN, SCND, (frac[next].lu[SCND] + frac[next].lu[VRGN]) / (frac[next].lu[CROP] + frac[next].lu[PSTR] + frac[next].lu[SCND] + frac[next].lu[VRGN])* (-1 * deltaU));
   }
   else if (deltaU > 0.0){
     // If urban area increases, check whether existing cropland, pasture and secondary
     // lands can satisfy the increase. If they can, land use change is proportional
     // to their relative areas, if not, take the additional land from virgin primary lands.
     if ((frac[now].lu[CROP] + frac[now].lu[PSTR] + frac[now].lu[SCND]) >= deltaU) {
       flow.set(CROP, URBN, frac[now].lu[CROP] / (frac[now].lu[CROP] + frac[now].lu[PSTR] + frac[now].lu[SCND]) * deltaU);       
       flow.set(PSTR, URBN, frac[now].lu[PSTR] / (frac[now].lu[CROP] + frac[now].lu[PSTR] + frac[now].lu[SCND]) * deltaU);       
       flow.set(SCND, URBN, frac[now].lu[SCND] / (frac[now].lu[CROP] + frac[now].lu[PSTR] + frac[now].lu[SCND]) * deltaU);       
     }
     else{
       flow.set(CROP, URBN, frac[now].lu[CROP]);
       flow.set(PSTR, URBN, frac[now].lu[PSTR]);
       flow.set(SCND, URBN, frac[now].lu[SCND]);
       flow.set(VRGN, URBN, deltaU - frac[now].lu[CROP] - frac[now].lu[PSTR] - frac[now].lu[SCND] );
     }
   }

   frac[now].u_p2u = flow.get(PSTR,URBN);
   frac[now].u_s2u = flow.get(SCND,URBN);
   frac[now].u_v2u = flow.get(VRGN,URBN);
   frac[now].u_u2s = flow.get(URBN,SCND);
   frac[now].u_u2p = flow.get(URBN,PSTR);
   
}


/******************************************************************************/
void GridCell::alternativeSmartFlow (size_t now, size_t next, bool warn) {
   // This function attributes land use change happening from land use types other
   // than urban.
   // Step 1: Compute land use change for cropland, pasture and (virgin + secondary),
   // that is apart from land use change into and out of urban lands.
   double deltaC = frac[next].lu[CROP] - frac[now].lu[CROP] - flow.get(URBN,CROP) + flow.get(CROP,URBN);
   double deltaP = frac[next].lu[PSTR] - frac[now].lu[PSTR] - flow.get(URBN,PSTR) + flow.get(PSTR,URBN);
   double deltaO = frac[next].lu[VRGN] + frac[next].lu[SCND]
                 - frac[now].lu[VRGN] - frac[now].lu[SCND]
                 - flow.get(URBN,SCND) + flow.get(SCND,URBN) + flow.get(VRGN,URBN);
                 
                 double denom = 0.0;

   // Step 2: Compute land use changes for the foll. cases:
   // 1. All zero i.e. do nothing!
   // 2. Two negative or decreasing in fraction of grid cell area: cropland, pasture
   // 3. Two negative or decreasing in fraction of grid cell area: cropland, other
   // 4. Two negative or decreasing in fraction of grid cell area: pasture, other
   // 5. Two positive or increasing in fraction of grid cell area: cropland, pasture
   // 6. Two positive or increasing in fraction of grid cell area: cropland, other
   // 7. Two positive or increasing in fraction of grid cell area: pasture, other   

   // all zero 
   if ((fabs(deltaC) < ZEROVALUE_CHECK) 
       && (fabs(deltaP) < ZEROVALUE_CHECK) 
       && (fabs(deltaO) < ZEROVALUE_CHECK)) {
      // Do nothing
   } 
   // two negative: cropland, pasture
   else if((deltaC <= 0.0) && (deltaP <= 0.0) && (deltaO > 0.0)) {
      // If cropland and pasture change is negative, but other is positive, this
      // means that land use is going from cropland, pasture to the other category.
      flow.set(CROP, SCND, -1 * deltaC);
      flow.set(PSTR, SCND, -1 * deltaP);
   }
   // two negative: cropland, other
   else if((deltaC <= 0.0) && (deltaP > 0.0) && (deltaO <= 0.0)) {
      flow.set(CROP, PSTR, -1 * deltaC);
      denom = frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN);
      if (deltaO< -ZEROVALUE_CHECK & fabs(denom)>ZEROVALUE_CHECK) {
      
     double vrgnRatio = (frac[now].lu[VRGN]-flow.get(VRGN,URBN))/(frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN));
  double scndRatio = (frac[now].lu[SCND]-flow.get(SCND,URBN))/(frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN));

      flow.set(VRGN,PSTR, -1 * vrgnRatio * deltaO);
      flow.set(SCND,PSTR, -1 * scndRatio * deltaO);
      }
      else {
      flow.set(VRGN,PSTR, 0.0);
      flow.set(SCND,PSTR, 0.0);
      }
      
   }
   // two negative: pasture, other
   else if((deltaC > 0.0) && (deltaP <= 0.0) && (deltaO <= 0.0)) {
         flow.set(PSTR, CROP, -1 * deltaP);
         denom = frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN);
	  if (deltaO< -ZEROVALUE_CHECK & fabs(denom)>ZEROVALUE_CHECK) {
      
	 double vrgnRatio = (frac[now].lu[VRGN]-flow.get(VRGN,URBN))/(frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN));
     double scndRatio = (frac[now].lu[SCND]-flow.get(SCND,URBN))/(frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN));

      flow.set(VRGN,CROP, -1 * vrgnRatio * deltaO);
      flow.set(SCND,CROP, -1 * scndRatio * deltaO);
	  }
	  else{
      flow.set(VRGN,CROP, 0.0);
      flow.set(SCND,CROP, 0.0);
	  }
      
   } 
   // two positive: pasture, other    
   else if ((deltaC <= 0.0) && (deltaP > 0.0) && (deltaO > 0.0)) {
      flow.set(CROP, SCND, deltaO);
      flow.set(CROP, PSTR, deltaP);
   }
   // two positive: cropland, pasture
   else if ((deltaC > 0.0) && (deltaP > 0.0) && (deltaO <= 0.0)) {
      denom = frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN);
      if (deltaO< -ZEROVALUE_CHECK & fabs(denom)>ZEROVALUE_CHECK) {
      
     double vrgnRatio = (frac[now].lu[VRGN]-flow.get(VRGN,URBN))/(frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN));
     double scndRatio = (frac[now].lu[SCND]-flow.get(SCND,URBN))/(frac[now].lu[VRGN]+frac[now].lu[SCND]-flow.get(VRGN,URBN)-flow.get(SCND,URBN));

      flow.set(VRGN,CROP, vrgnRatio * deltaC);
      flow.set(SCND,CROP, scndRatio * deltaC);

      flow.set(VRGN,PSTR, vrgnRatio * deltaP);
      flow.set(SCND,PSTR, scndRatio * deltaP);
      }
      else{
      flow.set(VRGN,CROP, 0.0);
      flow.set(SCND,CROP, 0.0);
      flow.set(VRGN,PSTR, 0.0);
      flow.set(SCND,PSTR, 0.0);
      }
     
   }
   // two positive: cropland, other
   else if ((deltaC > 0.0) && (deltaP <= 0.0) && (deltaO > 0.0)) {
      flow.set(PSTR, SCND, deltaO);
      flow.set(PSTR, CROP, deltaC);
   }  
	
	
	
    
   if (warn) {
      flow.resetNegatives(warn, "alternativeSmartFlow");
   }
}


/******************************************************************************/
void GridCell::adjustSmartFlow (size_t now, int shft_cult, GridInfo& grid, Config& config) {

   double c = frac[now].lu[CROP] - flow.getFlowsAway(CROP);
   double p = frac[now].lu[PSTR] - flow.getFlowsAway(PSTR);
   double v = frac[now].lu[VRGN] - flow.getFlowsAway(VRGN);
   double s = frac[now].lu[SCND] - flow.getFlowsAway(SCND);

   double vsTotal = v + s;

   double c2s = 0.0;
   double p2s = 0.0;
   double v2c = 0.0;
   double v2p = 0.0;
   double s2c = 0.0;
   double s2p = 0.0;

   // Shifting cultivation only happens if cropland fraction of grid
   // cell area is non-negative
   if (c > 0.0) {
   // abandonment of crop
      // best case method, desired abandonment = 1/15 per year
      // determine time_k
      double timeK = 0.0;
      if (shft_cult)  {
         // Read from netCDF file the shifting cultivation rate
	  timeK = frac[now].shftCult;     
      } 
      else {
         timeK = 1.0 / 15.0;
      }
      // Compute the fraction of grid cell area that will be used in shifting cultivation
      if (config.futureRun){
        if (grid.crop2015[y][x] < c){;
           c2s = grid.crop2015[y][x] * timeK;
        }
        else {
           c2s = c * timeK;
        }
           
      }
      else if (config.extensionRun){
       if (grid.crop2015[y][x] < c){;
           c2s = grid.crop2015[y][x] * timeK;
        }
        else {
           c2s = c * timeK;
        }
      }else{
        c2s = c * timeK;
      }
      
      if ((c2s) <= vsTotal) {
	
	if ((zoneTotal->priority == VRGN_PRIORITY) & (s < 10*timeK*c)) {
	          if ((c2s) <= v) {
               v2c = c2s;
            } else { // not enough v, but enough v+s
	       v2c = v;
               s2c = c2s - v2c;
            }
	} else	{
	        if ((c2s) <= s) {
               s2c = c2s;
            } else { // not enough s, but enough v+s
	        s2c = s;
               v2c = c2s - s2c;
            }
         }
      } else { // not enough v+s, take all v+s (everything) that is possible
      
	  v2c =  v;
	   s2c =  s;

         // re-define desired abandonment (flowcs_)
         // to be what is available
         c2s = v2c + s2c;
      }
      
      flow.set(VRGN, CROP, flow.get(VRGN, CROP) + v2c);
      flow.set(VRGN, PSTR, flow.get(VRGN, PSTR) + v2p);
      flow.set(SCND, CROP, flow.get(SCND, CROP) + s2c);
      flow.set(SCND, PSTR, flow.get(SCND, PSTR) + s2p);
      flow.set(CROP, SCND, flow.get(CROP, SCND) + c2s);
      flow.set(PSTR, SCND, flow.get(PSTR, SCND) + p2s);

      frac[now].sc_v2c = v2c;
      frac[now].sc_v2p = v2p;
      frac[now].sc_s2c = s2c;
      frac[now].sc_s2p = s2p;
      frac[now].sc_c2s = c2s;
      frac[now].sc_p2s = p2s;
      frac[now].sc_c = c;
      frac[now].sc_p = p;
      frac[now].sc_v = v;
      frac[now].sc_s = s;
    
      zoneTotal->flowvc_prime += v2c * cellArea / 1000. / 1000.;
      zoneTotal->flowvp_prime += v2p * cellArea / 1000. / 1000.;
      zoneTotal->flowsc_prime += s2c * cellArea / 1000. / 1000.;
      zoneTotal->flowsp_prime += s2p * cellArea / 1000. / 1000.;
      zoneTotal->flowcs_prime += c2s * cellArea / 1000. / 1000.;
      zoneTotal->flowps_prime += p2s * cellArea / 1000. / 1000.;
   }
}


/******************************************************************************/
void GridCell::updateForestLoss (size_t next, size_t now, int year) {
  
  if (potentialForest) {

    if (year >= 2001 & year <=2015) {
     frac[now].RSForestLossRemaining = frac[now].landsat/12;
      frac[now].RSForestLossRemaining = frac[now].RSForestLossRemaining - flow.get(VRGN,CROP) - flow.get(VRGN,PSTR);
      if (frac[now].scndMeanAge > 20) {
        frac[now].RSForestLossRemaining = frac[now].RSForestLossRemaining - flow.get(SCND,CROP) - flow.get(SCND,PSTR);
      }
       if (frac[now].RSForestLossRemaining < 0.0) {
         frac[now].RSForestLossRemaining = 0.0;
       }

     }
     else {
       frac[next].RSForestLossRemaining = frac[now].RSForestLossRemaining;
     }
    
  }


}

/******************************************************************************/
void GridCell::updateBiomassTallies (size_t now, InputFiles& inputs, int year) {

   double fromVrgn = flow.getFlowsAway(VRGN);
   double fromScnd = flow.getFlowsAway(SCND);
   double nextVrgn = frac[now].lu[VRGN] - fromVrgn;
   double nextScnd = frac[now].lu[SCND] - fromScnd;

double vrgnProtect = 0.0;
double scndProtect = 0.0;

if (year >1950 & year <2001 & frac[now].landsat >0.0 & (nextVrgn+nextScnd)>0.0 ) {
  vrgnProtect = nextVrgn/(nextVrgn+nextScnd)*frac[now].landsat;
  scndProtect = nextScnd/(nextVrgn+nextScnd)*frac[now].landsat;
  if (vrgnProtect > nextVrgn){
    vrgnProtect = nextVrgn;
  }
  if (scndProtect > nextScnd){
    scndProtect = nextScnd;
  }
}
   
   if (potentialForest) {
      zoneTotal->converted_forest_land += ((fromVrgn * potentialBiomass)
                                           + (fromScnd * frac[now].scndMeanBiomass))
                                          * cellArea / 1000.0;

 if ((frac[now].RSForestLossRemaining > 0) | (year < 2001) | (year > 2015) ){
      zoneTotal->vrgnForestBiomass[zdis] += (nextVrgn - vrgnProtect) * cellArea
                                          * potentialBiomass / 1000.0;
                                   
      zoneTotal->scndZdisForestBiomass[zdis] += (nextScnd - scndProtect) * cellArea
                                    * frac[now].scndMeanBiomass / 1000.0;
                             
      scndMatureForestBiomass = frac[now].scndMeanBiomass
                              * harvestProbability(frac[now].scndMeanBiomass,
                                                   inputs)
                              * (nextScnd - scndProtect) * cellArea;

      zoneTotal->scndMatureZdisForestBiomass[zdis] += scndMatureForestBiomass / 1000.0;

      scndForestBiomass = (nextScnd - scndProtect) * cellArea * frac[now].scndMeanBiomass; 

 } // RS forest Loss >0 
  else if (frac[now].scndMeanAge < 20){
    scndYoungForestBiomass = (nextScnd - scndProtect) * cellArea * frac[now].scndMeanBiomass
                              - scndMatureForestBiomass;
 
       zoneTotal->scndYoungForestBiomass += scndYoungForestBiomass / 1000.0;
       scndMatureForestBiomass = 0;
       }
       else {
       scndMatureForestBiomass = 0;
       scndYoungForestBiomass = 0;
       }

   } else { // Non forest 

      zoneTotal->vrgnNonForestBiomass += nextVrgn * cellArea
                                       * potentialBiomass / 1000.0;


      zoneTotal->scndNonForestBiomass += nextScnd * cellArea
                                       * frac[now].scndMeanBiomass / 1000.0;
   }

}

/******************************************************************************/
void GridCell::updateCropTypeTallies (size_t now) {
//country total crop type area in ha
        zoneTotal->c3annTotal += frac[now].subLU[C3ANN] * cellArea * 1e-4;
        zoneTotal->c4annTotal += frac[now].subLU[C4ANN] * cellArea * 1e-4;
        zoneTotal->c3perTotal += frac[now].subLU[C3PER] * cellArea * 1e-4;
        zoneTotal->c4perTotal += frac[now].subLU[C4PER] * cellArea * 1e-4;
        zoneTotal->c3nfxTotal += frac[now].subLU[C3NFX] * cellArea * 1e-4;
        zoneTotal->cropTotal += frac[now].lu[CROP] * cellArea * 1e-6; // in km^2
    
    }
      
/******************************************************************************/
void GridCell::updateFertTallies (size_t now) {  
     //   country total fertilizer usage per crop type in kg N 
     if (frac[now].subLU[C3ANN] > 1e-6){
        zoneTotal->c3annFertCurrent += frac[now].fert[C3ANN] * frac[now].subLU[C3ANN] * cellArea *1e-4;
        }
        if (frac[now].subLU[C3ANN] > 1e-6){
        zoneTotal->c4annFertCurrent += frac[now].fert[C4ANN] * frac[now].subLU[C4ANN] * cellArea *1e-4;
        }
        if (frac[now].subLU[C3ANN] > 1e-6){
        zoneTotal->c3perFertCurrent += frac[now].fert[C3PER] * frac[now].subLU[C3PER] * cellArea *1e-4;
        }
        if (frac[now].subLU[C3ANN] > 1e-6){
        zoneTotal->c4perFertCurrent += frac[now].fert[C4PER] * frac[now].subLU[C4PER] * cellArea *1e-4;
        }
        if (frac[now].subLU[C3ANN] > 1e-6){
        zoneTotal->c3nfxFertCurrent += frac[now].fert[C3NFX] * frac[now].subLU[C3NFX] * cellArea *1e-4;  
        }
                                        
                                
}

/******************************************************************************/
void GridCell::updateIrrigTallies (size_t now, size_t next, Config& config) {  
     //   country total fertilizer usage per crop type in km^2 
        zoneTotal->irrigCurrent += frac[now].irrig[C3ANN] * frac[now].subLU[C3ANN] * cellArea *1e-6;
        zoneTotal->irrigCurrent += frac[now].irrig[C4ANN] * frac[now].subLU[C4ANN] * cellArea *1e-6;
        zoneTotal->irrigCurrent += frac[now].irrig[C3PER] * frac[now].subLU[C3PER] * cellArea *1e-6;
        zoneTotal->irrigCurrent += frac[now].irrig[C4PER] * frac[now].subLU[C4PER] * cellArea *1e-6;
        zoneTotal->irrigCurrent += frac[now].irrig[C3NFX] * frac[now].subLU[C3NFX] * cellArea *1e-6;
        
        zoneTotal->irrigCurrentNew += frac[now].irrig[C3ANN] * frac[next].subLU[C3ANN] * cellArea *1e-6;
        zoneTotal->irrigCurrentNew += frac[now].irrig[C4ANN] * frac[next].subLU[C4ANN] * cellArea *1e-6;
        if (config.irrig_flag==1){
           if ((frac[next].subLU[C3PER]-frac[next].subLU[C3PER])<=0){
              zoneTotal->irrigCurrentNew += frac[now].irrig[C3PER] * frac[next].subLU[C3PER] * cellArea *1e-6;
           }
           else if ((frac[next].iamC3PBiofuel-frac[now].iamC3PBiofuel) <= 0){
              zoneTotal->irrigCurrentNew += frac[now].irrig[C3PER] * frac[next].subLU[C3PER] * cellArea *1e-6;
           }
           else {
              zoneTotal->irrigCurrentNew += frac[now].irrig[C3PER] * frac[now].subLU[C3PER] * cellArea *1e-6;
           }
           
           if ((frac[next].subLU[C4PER]-frac[next].subLU[C4PER])<=0){
              zoneTotal->irrigCurrentNew += frac[now].irrig[C4PER] * frac[next].subLU[C4PER] * cellArea *1e-6;
           }
           else if ((frac[next].iamC4PBiofuel-frac[now].iamC4PBiofuel) <= 0){
              zoneTotal->irrigCurrentNew += frac[now].irrig[C4PER] * frac[next].subLU[C4PER] * cellArea *1e-6;
           }
           else {
              zoneTotal->irrigCurrentNew += frac[now].irrig[C4PER] * frac[now].subLU[C4PER] * cellArea *1e-6;
           }
           
        }
        else{
           zoneTotal->irrigCurrentNew += frac[now].irrig[C3PER] * frac[next].subLU[C3PER] * cellArea *1e-6;
           zoneTotal->irrigCurrentNew += frac[now].irrig[C4PER] * frac[next].subLU[C4PER] * cellArea *1e-6;
        }
        zoneTotal->irrigCurrentNew += frac[now].irrig[C3NFX] * frac[next].subLU[C3NFX] * cellArea *1e-6;
        
        
        if (frac[now].irrig[C3ANN] > 0.000001) {
           zoneTotal->irrigAvail += (1-frac[now].irrig[C3ANN]) * frac[next].subLU[C3ANN] * cellArea *1e-6;
        }
        else {
           if (frac[next].subLU[C3ANN] > 0.000001) {
              zoneTotal->irrigNew += frac[next].subLU[C3ANN] * cellArea *1e-6;
           }
        }
        if (frac[now].irrig[C4ANN] > 0.000001) {
           zoneTotal->irrigAvail += (1-frac[now].irrig[C4ANN]) * frac[next].subLU[C4ANN] * cellArea *1e-6;
        }
        else {
           if (frac[next].subLU[C4ANN] > 0.000001) {
              zoneTotal->irrigNew += frac[next].subLU[C4ANN] * cellArea *1e-6;
           }
        }
        if (frac[now].irrig[C3PER] > 0.000001) {
           if (config.irrig_flag==1){
              zoneTotal->irrigAvail += (1-frac[now].irrig[C3PER]) * frac[next].subLU[C3PER] * (1-frac[next].iamC3PBiofuel) * cellArea *1e-6;
           }
           else{
              zoneTotal->irrigAvail += (1-frac[now].irrig[C3PER]) * frac[next].subLU[C3PER] * cellArea *1e-6;
           }
        }
        else {
           if (frac[next].subLU[C3PER] > 0.000001) {
              zoneTotal->irrigNew += frac[next].subLU[C3PER] * cellArea *1e-6;
           }
        }
        if (frac[now].irrig[C4PER] > 0.000001) {
           zoneTotal->irrigAvail += (1-frac[now].irrig[C4PER]) * frac[next].subLU[C4PER] * cellArea *1e-6;
        }
        else {
           if (frac[next].subLU[C4PER] > 0.000001) {
              zoneTotal->irrigNew += frac[next].subLU[C4PER] * cellArea *1e-6;
           }
        }
        if (frac[now].irrig[C3NFX] > 0.000001) {
           zoneTotal->irrigAvail += (1-frac[now].irrig[C3NFX]) * frac[next].subLU[C3NFX] * cellArea *1e-6;   
        }
        else {
           if (frac[next].subLU[C3NFX] > 0.000001) {
              zoneTotal->irrigNew += frac[next].subLU[C3NFX] * cellArea *1e-6;
           }
        }   
        
                                   
                                
}

/******************************************************************************/
void GridCell::harvest(size_t now, InputFiles& inputs, Config& config, int year) {

   // Priority is taken account when pre-computing ratios.
   // These should be able to be called in any order.

   // Note: getFlowsAway includes VRGN to SCND, but this value should always be
   //       0 when starting this routine. It should also never be set twice in
   //       this routine (one is forest, one is nonforest). Careful if either
   //       of those criteria changes.
   
double vrgnProtect = 0.0;
double scndProtect = 0.0;

double nextVrgn = (frac[now].lu[VRGN] - flow.getFlowsAway(VRGN));
double nextScnd = (frac[now].lu[SCND] - flow.getFlowsAway(SCND));
if (year >1950 & year <2001 & frac[now].landsat >0.0 & (nextVrgn+nextScnd)>0.0 ) {
  vrgnProtect = nextVrgn/(nextVrgn+nextScnd)*frac[now].landsat;
  scndProtect = nextScnd/(nextVrgn+nextScnd)*frac[now].landsat;
  if (vrgnProtect > nextVrgn){
    vrgnProtect = nextVrgn;
  }
  if (scndProtect > nextScnd){
    scndProtect = nextScnd;
  }
}
 
 
   if (config.VIRGIN_HARVEST_SWITCH) {
      if (potentialForest) {

if (frac[now].RSForestLossRemaining>0.0 | year < 2001 | year > 2015){
        flow.set(VRGN, SCND, (frac[now].lu[VRGN] - flow.getFlowsAway(VRGN) - vrgnProtect)
                              * zoneTotal->vrgnForestHarvestRatio[zdis]);

         flow.vrgnHarvBiomass1 = flow.get(VRGN, SCND) * cellArea * potentialBiomass;

         zoneTotal->harvestRemaining -= flow.vrgnHarvBiomass1 / 1000.;
	 	 zoneTotal->totalHarvested += flow.vrgnHarvBiomass1 / 1000.;
	 	 zoneTotal->harvestedVrgnForestBiomass +=flow.vrgnHarvBiomass1 / 1000.;
}

      }
   }

    if (config.SECONDARY_HARVEST_SWITCH && zoneTotal->harvestRemaining) {
      if (potentialForest) {
      if (frac[now].RSForestLossRemaining>0.0 | year < 2001 | year > 2015){

         flow.scndHarvBiomass1 = (frac[now].lu[SCND] - flow.getFlowsAway(SCND) - scndProtect)
                              * zoneTotal->scndZdisForestHarvestRatio[zdis] 
                              * harvestProbability(frac[now].scndMeanBiomass, inputs)* cellArea * frac[now].scndMeanBiomass ;

         zoneTotal->harvestRemaining -= flow.scndHarvBiomass1 / 1000.;
	 	 zoneTotal->totalHarvested += flow.scndHarvBiomass1 / 1000.;
	 	 zoneTotal->harvestedScndMatureForestBiomass +=flow.scndHarvBiomass1 / 1000.;

       }
      }
   }
   
   if (config.FORCE_HARVEST_SWITCH && zoneTotal->harvestRemaining) {
      if (potentialForest) {
      if (frac[now].RSForestLossRemaining>0.0 | year < 2001 | year >2015) {
         flow.scndHarvBiomass2 = (scndForestBiomass - flow.scndHarvBiomass1)
                               * zoneTotal->scndYoungForestHarvestRatio;
         zoneTotal->harvestRemaining -= flow.scndHarvBiomass2 / 1000.;
	  zoneTotal->totalHarvested += flow.scndHarvBiomass2/1000.;
	  }
	   else if (frac[now].scndMeanAge < 20) {
 	  flow.scndHarvBiomass2 = (scndYoungForestBiomass)
                                * zoneTotal->scndYoungForestHarvestRatio;
          zoneTotal->harvestRemaining -= flow.scndHarvBiomass2 / 1000.;
 	  zoneTotal->totalHarvested += flow.scndHarvBiomass2/1000.;
 	  }
	  
      } else {
        
         flow.set(VRGN, SCND, (frac[now].lu[VRGN] - flow.getFlowsAway(VRGN))
                              * zoneTotal->vrgnNonForestHarvestRatio);

         flow.vrgnHarvBiomass2 = flow.get(VRGN, SCND) * cellArea * potentialBiomass;


         zoneTotal->harvestRemaining -= flow.vrgnHarvBiomass2 / 1000.;
	  zoneTotal->totalHarvested += flow.vrgnHarvBiomass2/1000.;

         flow.scndHarvBiomass3 = (frac[now].lu[SCND] - flow.getFlowsAway(SCND))
                               * cellArea * frac[now].scndMeanBiomass
                               * zoneTotal->scndNonForestHarvestRatio;
         zoneTotal->harvestRemaining -= flow.scndHarvBiomass3 / 1000.;
	  zoneTotal->totalHarvested += flow.scndHarvBiomass3/1000.;
      }
   }
   if (potentialBiomass > 0.0) {
      flow.vrgnHarvArea1 = flow.vrgnHarvBiomass1 / potentialBiomass / cellArea;
      flow.vrgnHarvArea2 = flow.vrgnHarvBiomass2 / potentialBiomass / cellArea;
   } else {
      flow.vrgnHarvArea1  = 0.0;
      flow.vrgnHarvArea2  = 0.0;
   }
   if (frac[now].scndMeanBiomass > 0.0) {
      flow.scndHarvArea1  =  flow.scndHarvBiomass1 / frac[now].scndMeanBiomass / cellArea;
      flow.scndHarvArea2  =  flow.scndHarvBiomass2 / frac[now].scndMeanBiomass / cellArea;
      flow.scndHarvArea3  =  flow.scndHarvBiomass3 / frac[now].scndMeanBiomass / cellArea;
   } else {
      flow.scndHarvArea1  = 0.0;
      flow.scndHarvArea2  = 0.0;
      flow.scndHarvArea3  = 0.0;
   }
   
}

   
/******************************************************************************/
void GridCell::doSubTransitions (size_t now, size_t next, Config& config) {

   flow.resetNegatives(config.FLOW_BUG_PRINT, "doSubTransitions");

   // Pasture and rangeland transition areas
   double temp_pastr;
   double temp_range;
   double temp_graz;

   // CFT fractions
   double c3AnnualFraction;
   double c4AnnualFraction;
   double c3PerennialFraction;
   double c4PerennialFraction;
   double NfixingFraction;
   
   sub_lu_types primType;
   sub_lu_types secdType;
            
   c3AnnualFraction = zoneTotal->c3AnnualFraction;
   c4AnnualFraction = zoneTotal->c4AnnualFraction;
   c3PerennialFraction = zoneTotal->c3PerennialFraction;
   c4PerennialFraction = zoneTotal->c4PerennialFraction;
   NfixingFraction = zoneTotal->NfixingFraction;
 
   if (potentialForest) {
      primType = PRIMF;
      secdType = SECDF;
   } else {
      primType = PRIMN;
      secdType = SECDN;
   }
   
   // Copy the 'other' transitions that don't involve crop or pasture
   flow.setSub(primType, secdType, flow.get(VRGN, SCND));
   flow.setSub(primType, URBAN, flow.get(VRGN, URBN));
   flow.setSub(secdType, URBAN, flow.get(SCND, URBN));
   flow.setSub(URBAN, secdType, flow.get(URBN, SCND));

   double IAMDeltaC4P = frac[next].iamC4P - frac[now].iamC4P;
   double IAMDeltaC3P = frac[next].iamC3P - frac[now].iamC3P;
   double IAMDeltaC4A = frac[next].iamC4A - frac[now].iamC4A;
   
   // Create new croplands
   double v2cFrac = flow.get(VRGN, CROP);
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (croptype_ratio_total > 0.0) {
       if (IAMDeltaC4P>=v2cFrac & v2cFrac>0.000001 & IAMDeltaC4P> 0.000001 & config.futureRun==1 & config.c4p_flag==1){
        switch (ct) {
               case glm::C3ANN: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               case glm::C4ANN: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               case glm::C3PER: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               case glm::C4PER: flow.setSub(primType, (sub_lu_types)ct, v2cFrac); break;
               case glm::C3NFX: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
            
         }
        else if (IAMDeltaC4A>=v2cFrac & v2cFrac>0.000001 & IAMDeltaC4A> 0.000001 & config.futureRun==1 & config.c4a_flag==1){

        switch (ct) {
               case glm::C3ANN: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               case glm::C4ANN: flow.setSub(primType, (sub_lu_types)ct, v2cFrac); break;
               case glm::C3PER: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               case glm::C4PER: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               case glm::C3NFX: flow.setSub(primType, (sub_lu_types)ct, 0); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
            }
        else{
         flow.setSub(primType, (sub_lu_types)ct, v2cFrac * croptype_ratio[ct]);
         }
      } else { 
    
         // When GLM thinks a grid cell has agriculture and MONFREDA does not, assign
         // crop functional type fractions based on a pre-computed layer
         if (config.cft_option == 1) {     
            switch (ct) {
               case glm::C3ANN: flow.setSub(primType, (sub_lu_types)ct, v2cFrac * c3AnnualFraction); break;
               case glm::C4ANN: flow.setSub(primType, (sub_lu_types)ct, v2cFrac * c4AnnualFraction); break;
               case glm::C3PER: flow.setSub(primType, (sub_lu_types)ct, v2cFrac * c3PerennialFraction); break;
               case glm::C4PER: flow.setSub(primType, (sub_lu_types)ct, v2cFrac * c4PerennialFraction); break;
               case glm::C3NFX: flow.setSub(primType, (sub_lu_types)ct, v2cFrac * NfixingFraction); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
         } else {
            flow.setSub(primType, (sub_lu_types)ct, v2cFrac / 5.0);
         }
      }
   }
   if (IAMDeltaC4P>=v2cFrac & v2cFrac>0.000001 & IAMDeltaC4P> 0.000001 & config.futureRun==1 & config.c4p_flag==1){
   IAMDeltaC4P = IAMDeltaC4P - v2cFrac;
   }
   if (IAMDeltaC4A>=v2cFrac & v2cFrac>0.000001 & IAMDeltaC4A> 0.000001 & config.futureRun==1 & config.c4a_flag==1){
   IAMDeltaC4A = IAMDeltaC4A - v2cFrac;
   }
      
   double s2cFrac = flow.get(SCND, CROP);
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (croptype_ratio_total > 0.0) {
      if (IAMDeltaC4P>=s2cFrac & s2cFrac>0.000001 & IAMDeltaC4P> 0.000001 & config.futureRun==1 & config.c4p_flag==1){
      switch (ct) {
               case glm::C3ANN: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               case glm::C4ANN: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               case glm::C3PER: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               case glm::C4PER: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac); break;
               case glm::C3NFX: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
        
}
else if (IAMDeltaC4A>=s2cFrac & s2cFrac>0.000001 & IAMDeltaC4A> 0.000001 & config.futureRun==1 & config.c4a_flag==1){
        switch (ct) {
               case glm::C3ANN: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               case glm::C4ANN: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac); break;
               case glm::C3PER: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               case glm::C4PER: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               case glm::C3NFX: flow.setSub(secdType, (sub_lu_types)ct, 0); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
            }

         else{
          flow.setSub(secdType, (sub_lu_types)ct, s2cFrac * croptype_ratio[ct]);
          }
         
      } else {
         if (config.cft_option == 1) {                
            switch (ct) {
               case glm::C3ANN: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac * c3AnnualFraction); break;
               case glm::C4ANN: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac * c4AnnualFraction); break;
               case glm::C3PER: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac * c3PerennialFraction); break;
               case glm::C4PER: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac * c4PerennialFraction); break;
               case glm::C3NFX: flow.setSub(secdType, (sub_lu_types)ct, s2cFrac * NfixingFraction); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
         } else {
            flow.setSub(secdType, (sub_lu_types)ct, s2cFrac / 5.0);
         }
      }
   }
  if (IAMDeltaC4P>=s2cFrac & s2cFrac>0.000001 & IAMDeltaC4P> 0.000001 & config.futureRun==1 & config.c4p_flag==1){
     IAMDeltaC4P = IAMDeltaC4P - s2cFrac;
   }
   if (IAMDeltaC4A>=s2cFrac & s2cFrac>0.000001 & IAMDeltaC4A> 0.000001 & config.futureRun==1 & config.c4a_flag==1){
     IAMDeltaC4A = IAMDeltaC4A - s2cFrac;
   }
   
   double u2cFrac = flow.get(URBN, CROP);
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (croptype_ratio_total > 0.0) {
         flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac * croptype_ratio[ct]);
      } else {
         if (config.cft_option == 1) {                
            switch (ct) {
               case glm::C3ANN: flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac * c3AnnualFraction); break;
               case glm::C4ANN: flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac * c4AnnualFraction); break;
               case glm::C3PER: flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac * c3PerennialFraction); break;
               case glm::C4PER: flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac * c4PerennialFraction); break;
               case glm::C3NFX: flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac * NfixingFraction); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
         } else {
            flow.setSub(URBAN, (sub_lu_types)ct, u2cFrac / 5.0);
         }
      }
   }

   // Abandon/transition-away croplands
   double c2sFrac = flow.get(CROP, SCND);
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (croptype_ratio_total > 0.0) {
         flow.setSub((sub_lu_types)ct, secdType, c2sFrac * croptype_ratio[ct]);
      } else {
         if (config.cft_option == 1) {                
            switch (ct) {
               case glm::C3ANN: flow.setSub((sub_lu_types)ct, secdType, c2sFrac * c3AnnualFraction); break;
               case glm::C4ANN: flow.setSub((sub_lu_types)ct, secdType, c2sFrac * c4AnnualFraction); break;
               case glm::C3PER: flow.setSub((sub_lu_types)ct, secdType, c2sFrac * c3PerennialFraction); break;
               case glm::C4PER: flow.setSub((sub_lu_types)ct, secdType, c2sFrac * c4PerennialFraction); break;
               case glm::C3NFX: flow.setSub((sub_lu_types)ct, secdType, c2sFrac * NfixingFraction); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
         } else {
            flow.setSub((sub_lu_types)ct, secdType, c2sFrac / 5.0);
         }
      }
   }
   double c2uFrac = flow.get(CROP, URBN);
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (croptype_ratio_total > 0.0) {
         flow.setSub((sub_lu_types)ct, URBAN, c2uFrac * croptype_ratio[ct]);
      } else {
         if (config.cft_option == 1) {                
            switch (ct) {
               case glm::C3ANN: flow.setSub((sub_lu_types)ct, URBAN, c2uFrac * c3AnnualFraction); break;
               case glm::C4ANN: flow.setSub((sub_lu_types)ct, URBAN, c2uFrac * c4AnnualFraction); break;
               case glm::C3PER: flow.setSub((sub_lu_types)ct, URBAN, c2uFrac * c3PerennialFraction); break;
               case glm::C4PER: flow.setSub((sub_lu_types)ct, URBAN, c2uFrac * c4PerennialFraction); break;
               case glm::C3NFX: flow.setSub((sub_lu_types)ct, URBAN, c2uFrac * NfixingFraction); break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
         } else {
            flow.setSub((sub_lu_types)ct, URBAN, c2uFrac / 5.0);
         }         
      }
   }

   if (frac[now].u_p2u > 0){
     
     flow.setSub(PASTR, URBAN, frac[now].u_p2u*frac[now].subLU[PASTR]/frac[now].lu[PSTR]);
     flow.setSub(RANGE, URBAN, frac[now].u_p2u*frac[now].subLU[RANGE]/frac[now].lu[PSTR]);
   }
   
   if (frac[now].u_u2p > 0){
     
     flow.setSub(URBAN, PASTR, frac[now].u_u2p*frac[next].subLU[PASTR]/frac[next].lu[PSTR]);
     flow.setSub(URBAN, RANGE, frac[now].u_u2p*frac[next].subLU[RANGE]/frac[next].lu[PSTR]);
   }
   
     
   if (!config.futureRun){
      // temp_pastr, temp_range denote how much pasture and rangeland are
      // going to increase/decrease in the coming time-step
      temp_pastr = frac[next].subLU[PASTR] - frac[now].subLU[PASTR];
      temp_range = frac[next].subLU[RANGE] - frac[now].subLU[RANGE];

      // First take care of urban transitioning into or out of pasture and rangeland
      if (frac[now].u_u2p >0){
         temp_pastr -= frac[now].u_u2p*frac[next].subLU[PASTR]/frac[next].lu[PSTR];
         temp_range -= frac[now].u_u2p*frac[next].subLU[RANGE]/frac[next].lu[PSTR];
      }
      if (frac[now].u_p2u >0){
         temp_pastr += frac[now].u_p2u*frac[now].subLU[PASTR]/frac[now].lu[PSTR];
         temp_range += frac[now].u_p2u*frac[now].subLU[RANGE]/frac[now].lu[PSTR];
      }
    
      temp_graz = frac[next].lu[PSTR] - frac[now].lu[PSTR] - frac[now].u_u2p + frac[now].u_p2u;

      if (temp_graz>0 & temp_pastr<0) {
	flow.setSub(PASTR,RANGE,-1*temp_pastr);
	temp_pastr = 0.0;
	temp_range = temp_graz;
      }
      else if (temp_graz>0 & temp_range<0){
	flow.setSub(RANGE,PASTR,-1*temp_range);
	temp_range = 0.0;
	temp_pastr = temp_graz;
      }
      else if (temp_graz<0 & temp_pastr>0){
	flow.setSub(RANGE,PASTR,temp_pastr);
	temp_pastr = 0.0;
	temp_range = temp_graz;
      }
      else if (temp_graz<0 & temp_range>0){
	flow.setSub(PASTR,RANGE,temp_range);
	temp_range = 0.0;
	temp_pastr = temp_graz;
      }
      else if(temp_graz==0 & temp_pastr<0){
        flow.setSub(PASTR,RANGE,-1*temp_pastr);
	temp_pastr = 0.0;
	temp_range = 0.0;;
      }
      else if (temp_graz==0 & temp_range<0){
	flow.setSub(RANGE,PASTR,-1*temp_range);
	temp_range = 0.0;
	temp_pastr = 0.0;
      }
      else if (fabs(temp_graz)<0.0000009 & fabs(temp_range)<0.0000009 & fabs(temp_pastr)<0.0000009){
	flow.setSub(RANGE,PASTR,-1*temp_range);
	temp_range = 0.0;
	temp_pastr = 0.0;
      }
      else {
      cout << "grazing error \n";
           }
   
      double c2pFrac = flow.get(CROP, PSTR);
      
      for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
        if (fabs(temp_graz)>0.0000000000001){
	  
          if (croptype_ratio_total > 0.0) {
            flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * croptype_ratio[ct] * temp_pastr/(temp_graz));
            flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * croptype_ratio[ct] * temp_range/(temp_graz));

          } else {
               if (config.cft_option == 1) {                
                  switch (ct) {
                     case glm::C3ANN: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * c3AnnualFraction * temp_pastr/(temp_graz));
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * c3AnnualFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C4ANN: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * c4AnnualFraction * temp_pastr/(temp_graz));
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * c4AnnualFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C3PER: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * c3PerennialFraction * temp_pastr/(temp_graz));
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * c3PerennialFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C4PER: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * c4PerennialFraction * temp_pastr/(temp_graz));
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * c4PerennialFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C3NFX: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * NfixingFraction * temp_pastr/(temp_graz));
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * NfixingFraction * temp_range/(temp_graz));
                                      break;
                     default: std::cerr << "Incorrect crop functional type\n";
                  }
               } else {
                  flow.setSub((sub_lu_types)ct, PASTR, c2pFrac / 5.0 * temp_pastr/(temp_graz));
                  flow.setSub((sub_lu_types)ct, RANGE, c2pFrac / 5.0 * temp_range/(temp_graz));
               }            
           }
        }
        else{
          flow.setSub((sub_lu_types)ct, PASTR, 0);
          flow.setSub((sub_lu_types)ct, RANGE, 0);
	}
      }

      // Create new pasture
      
      if (frac[now].sc_v2p <= 0) {
        if (fabs(temp_graz)> 0.0000000000001){
          flow.setSub(primType, PASTR, flow.get(VRGN, PSTR)*temp_pastr/(temp_graz));
          flow.setSub(primType, RANGE, flow.get(VRGN, PSTR)*temp_range/(temp_graz));
        }
        else{
          flow.setSub(primType, PASTR, 0);
          flow.setSub(primType, RANGE, 0);
        }
      }
      else {
        if (fabs(temp_graz)> 0.0000000000001){
          flow.setSub(primType, PASTR, (flow.get(VRGN, PSTR)-frac[now].sc_v2p)*temp_pastr/(temp_graz));
          flow.setSub(primType, RANGE, (flow.get(VRGN, PSTR)-frac[now].sc_v2p)*temp_range/(temp_graz));
        }
        flow.setSub(primType, PASTR, flow.getSub(primType, PASTR) + frac[now].sc_v2p*frac[now].subLU[PASTR]/frac[now].lu[PSTR]);
        flow.setSub(primType, RANGE, flow.getSub(primType, RANGE) + frac[now].sc_v2p*frac[now].subLU[RANGE]/frac[now].lu[PSTR]);
      }

           
      if (frac[now].sc_s2p <= 0) {
        if (fabs(temp_graz) > 0.0000000000001){
          flow.setSub(secdType, PASTR, flow.get(SCND, PSTR)*temp_pastr/(temp_graz));
          flow.setSub(secdType, RANGE, flow.get(SCND, PSTR)*temp_range/(temp_graz));
        }
        else{
          flow.setSub(secdType, PASTR, 0);
          flow.setSub(secdType, RANGE, 0);
        }
      }
      else {
        if (fabs(temp_graz) > 0.0000000000001){
          flow.setSub(secdType, PASTR, (flow.get(SCND, PSTR)-frac[now].sc_s2p)*temp_pastr/(temp_graz));
          flow.setSub(secdType, RANGE, (flow.get(SCND, PSTR)-frac[now].sc_s2p)*temp_range/(temp_graz));
        }
        flow.setSub(secdType, PASTR, flow.getSub(secdType, PASTR) + frac[now].sc_s2p*frac[now].subLU[PASTR]/frac[now].lu[PSTR]);
        flow.setSub(secdType, RANGE, flow.getSub(secdType, RANGE) + frac[now].sc_s2p*frac[now].subLU[RANGE]/frac[now].lu[PSTR]);
      }

      // Abandon/transition-away pasture
      if (fabs(temp_graz) > 0.0000000000001){
         flow.setSub(PASTR, secdType, flow.get(PSTR, SCND)*temp_pastr/(temp_graz));
         flow.setSub(RANGE, secdType, flow.get(PSTR, SCND)*temp_range/(temp_graz));
      }
      else{
        flow.setSub(PASTR, secdType, 0);
        flow.setSub(RANGE, secdType, 0);
      }
      
      double p2cFrac = flow.get(PSTR, CROP);
      for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
        if (fabs(temp_graz) >0.0000000000001){
          if (croptype_ratio_total > 0.0) {
            flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac * croptype_ratio[ct] * temp_pastr/(temp_graz));
            flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac * croptype_ratio[ct] * temp_range/(temp_graz));
          }
          else {
            if (config.cft_option == 1) {                
                  switch (ct) {
                     case glm::C3ANN: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c3AnnualFraction * temp_pastr/(temp_graz));
                                      flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c3AnnualFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C4ANN: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c4AnnualFraction * temp_pastr/(temp_graz));
                                      flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c4AnnualFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C3PER: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c3PerennialFraction * temp_pastr/(temp_graz));
                                      flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c3PerennialFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C4PER: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c4PerennialFraction * temp_pastr/(temp_graz));
                                      flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c4PerennialFraction * temp_range/(temp_graz));
                                      break;
                     case glm::C3NFX: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * NfixingFraction * temp_pastr/(temp_graz));
                                      flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * NfixingFraction * temp_range/(temp_graz));
                                      break;
                     default: std::cerr << "Incorrect crop functional type\n";
                  }
               } else {
                  flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac / 5.0 * temp_pastr/(temp_graz));
                  flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac / 5.0 * temp_range/(temp_graz));
                 }
            }
	}
         else {
            flow.setSub(PASTR, (sub_lu_types)ct,0);
            flow.setSub(RANGE, (sub_lu_types)ct,0);
	 }
	
      }

}
   else { // future run
      flow.setSub(PASTR, URBAN, flow.get(PSTR,URBN)*croptype_ratio[PASTR]);
      flow.setSub(RANGE, URBAN, flow.get(PSTR,URBN)*croptype_ratio[RANGE]); 
  
      flow.setSub(URBAN, PASTR, flow.get(URBN,PSTR)*croptype_ratio[PASTR]);
      flow.setSub(URBAN, RANGE, flow.get(URBN,PSTR)*croptype_ratio[RANGE]);

   
      double c2pFrac = flow.get(CROP, PSTR);
      for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
         if (croptype_ratio_total > 0.0) {
            flow.setSub((sub_lu_types)ct, PASTR, c2pFrac * croptype_ratio[ct] * croptype_ratio[PASTR]);
            flow.setSub((sub_lu_types)ct, RANGE, c2pFrac * croptype_ratio[ct] * croptype_ratio[RANGE]);
         } else {
               if (config.cft_option == 1) {                
                  switch (ct) {
                     case glm::C3ANN: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac  * c3AnnualFraction * croptype_ratio[PASTR]);
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac  * c3AnnualFraction * croptype_ratio[RANGE]);
                                      break;
                     case glm::C4ANN: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac  * c4AnnualFraction * croptype_ratio[PASTR]);
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac  * c4AnnualFraction * croptype_ratio[RANGE]);
                                      break;
                     case glm::C3PER: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac  * c3PerennialFraction * croptype_ratio[PASTR]);
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac  * c3PerennialFraction * croptype_ratio[RANGE]);
                                      break;
                     case glm::C4PER: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac  * c4PerennialFraction * croptype_ratio[PASTR]);
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac  * c4PerennialFraction * croptype_ratio[RANGE]);
                                      break;
                     case glm::C3NFX: flow.setSub((sub_lu_types)ct, PASTR, c2pFrac  * NfixingFraction * croptype_ratio[PASTR]);
                                      flow.setSub((sub_lu_types)ct, RANGE, c2pFrac  * NfixingFraction * croptype_ratio[RANGE]);
                                      break;
                     default: std::cerr << "Incorrect crop functional type\n";
                  }
               } else {
                  flow.setSub((sub_lu_types)ct, PASTR, c2pFrac / 5.0 * croptype_ratio[PASTR]);
                  flow.setSub((sub_lu_types)ct, RANGE, c2pFrac / 5.0 * croptype_ratio[RANGE]);
               }
         }
      }

    
      flow.setSub(primType, PASTR, flow.get(VRGN, PSTR)*croptype_ratio[PASTR]);
      flow.setSub(primType, RANGE, flow.get(VRGN, PSTR)*croptype_ratio[RANGE]);


      flow.setSub(secdType, PASTR, flow.get(SCND, PSTR)*croptype_ratio[PASTR]);
      flow.setSub(secdType, RANGE, flow.get(SCND, PSTR)*croptype_ratio[RANGE]);
     
	 
      // Abandon/transition-away pasture
      flow.setSub(PASTR, secdType, flow.get(PSTR, SCND)*croptype_ratio[PASTR]);
      flow.setSub(RANGE, secdType, flow.get(PSTR, SCND)*croptype_ratio[RANGE]);
  
      double p2cFrac = flow.get(PSTR, CROP);
      for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
         if (croptype_ratio_total > 0.0) {
         if (IAMDeltaC4P>=p2cFrac & p2cFrac>0.000001 & IAMDeltaC4P> 0.000001 & config.futureRun==1 & config.c4p_flag==1){
      switch (ct) {
               case glm::C3ANN: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               case glm::C4ANN: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               case glm::C3PER: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               case glm::C4PER: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac* croptype_ratio[PASTR]);
               flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac* croptype_ratio[RANGE]); break;
               case glm::C3NFX: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
            }
      else if (IAMDeltaC4A>=p2cFrac & p2cFrac>0.000001 & IAMDeltaC4A> 0.000001 & config.futureRun==1 & config.c4a_flag==1){
        switch (ct) {
               case glm::C3ANN: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               case glm::C4ANN: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac* croptype_ratio[PASTR]); 
               flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac* croptype_ratio[RANGE]);break;
               case glm::C3PER: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               case glm::C4PER: flow.setSub(PASTR, (sub_lu_types)ct, 0);
               flow.setSub(RANGE, (sub_lu_types)ct, 0); break;
               case glm::C3NFX: flow.setSub(PASTR, (sub_lu_types)ct, 0); 
               flow.setSub(RANGE, (sub_lu_types)ct, 0);break;
               default: std::cerr << "Incorrect crop functional type\n";
            }
            }
         
         else{
          flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac * croptype_ratio[ct] * croptype_ratio[PASTR]);
          flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac * croptype_ratio[ct] * croptype_ratio[RANGE]);
          }
         
         }
         else {
            if (config.cft_option == 1) {                
               switch (ct) {
                  case glm::C3ANN: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c3AnnualFraction * croptype_ratio[PASTR]);
                                   flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c3AnnualFraction * croptype_ratio[RANGE]);
                                   break;
                  case glm::C4ANN: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c4AnnualFraction * croptype_ratio[PASTR]);
                                   flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c4AnnualFraction * croptype_ratio[RANGE]);
                                   break;
                  case glm::C3PER: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c3PerennialFraction * croptype_ratio[PASTR]);
                                   flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c3PerennialFraction * croptype_ratio[RANGE]);
                                   break;
                  case glm::C4PER: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * c4PerennialFraction * croptype_ratio[PASTR]);
                                   flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * c4PerennialFraction * croptype_ratio[RANGE]);
                                   break;
                  case glm::C3NFX: flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac  * NfixingFraction * croptype_ratio[PASTR]);
                                   flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac  * NfixingFraction * croptype_ratio[RANGE]);
                                   break;
                  default: std::cerr << "Incorrect crop functional type\n";
               }
            } else {
                  flow.setSub(PASTR, (sub_lu_types)ct, p2cFrac / 5.0 * croptype_ratio[PASTR]);
                  flow.setSub(RANGE, (sub_lu_types)ct, p2cFrac / 5.0 * croptype_ratio[RANGE]);
            }
         }
      }  
      
      if (IAMDeltaC4P>=p2cFrac & p2cFrac>0.000001 & IAMDeltaC4P> 0.000001 & config.futureRun==1 & config.c4p_flag==1){
     IAMDeltaC4P = IAMDeltaC4P - p2cFrac;
   }
   if (IAMDeltaC4A>=p2cFrac & p2cFrac>0.000001 & IAMDeltaC4A> 0.000001 & config.futureRun==1 & config.c4a_flag==1){
     IAMDeltaC4A = IAMDeltaC4A - p2cFrac;
   }
     
   }

}
/******************************************************************************/
void GridCell::updateCropTypesRatios (size_t nowxt, Config& config) {

  if (frac[nowxt].lu[CROP]>0.0000001 & croptype_ratio_total>0){
      croptype_ratio[C3ANN] = frac[nowxt].subLU[C3ANN]/frac[nowxt].lu[CROP];
      croptype_ratio[C4ANN] = frac[nowxt].subLU[C4ANN]/frac[nowxt].lu[CROP];
      croptype_ratio[C3PER] = frac[nowxt].subLU[C3PER]/frac[nowxt].lu[CROP];
      croptype_ratio[C4PER] = frac[nowxt].subLU[C4PER]/frac[nowxt].lu[CROP];
      croptype_ratio[C3NFX] = frac[nowxt].subLU[C3NFX]/frac[nowxt].lu[CROP];
      croptype_ratio_total = croptype_ratio[C3ANN] + croptype_ratio[C4ANN]
                           + croptype_ratio[C3PER] + croptype_ratio[C4PER]
                           + croptype_ratio[C3NFX];
                           
      croptype_ratio[C3ANN] = croptype_ratio[C3ANN]/croptype_ratio_total;
      croptype_ratio[C4ANN] = croptype_ratio[C4ANN]/croptype_ratio_total;
      croptype_ratio[C3PER] = croptype_ratio[C3PER]/croptype_ratio_total;
      croptype_ratio[C4PER] = croptype_ratio[C4PER]/croptype_ratio_total;
      croptype_ratio[C3NFX] = croptype_ratio[C3NFX]/croptype_ratio_total;
      croptype_ratio_total = croptype_ratio[C3ANN] + croptype_ratio[C4ANN]
                           + croptype_ratio[C3PER] + croptype_ratio[C4PER]
                           + croptype_ratio[C3NFX];
                                                
                           
    }                       
  }

/******************************************************************************/
void GridCell::updateStates (size_t now, size_t next, Config& config) {

  double othr = frac[next].lu[VRGN] + frac[next].lu[SCND];
  double pstr = frac[next].lu[PSTR];
  double crop = frac[next].lu[CROP];
  double urbn = frac[next].lu[URBN];
  
  frac[next].priorWH = frac[now].priorWH;
  if (flow.vrgnHarvArea1 >0.0){
  frac[next].priorWH = 1;
  }
  else if (flow.vrgnHarvArea2 > 0.0) {
  frac[next].priorWH = 1;
  }
  else if (flow.scndHarvArea1 >0.0){
  frac[next].priorWH = 1;
  }
  else if (flow.scndHarvArea2 > 0.0) {
  frac[next].priorWH = 1;
  }
  else if (flow.scndHarvArea3 > 0.0) {
  frac[next].priorWH = 1;
  }
  
  frac[next].landsat = frac[now].landsat;
  
   for (size_t target=VRGN; target<N_LU_TYPES; ++target) {
      frac[next].lu[target] = frac[now].lu[target];
      for (size_t source=VRGN; source<N_LU_TYPES; ++source) {
         if (target != source) {
            if (target != VRGN) {
            
               frac[next].lu[target] += flow.get(source, target);
	       
            }
            frac[next].lu[target] -= flow.get(target, source);
         }
      }
   }

   if ((frac[next].lu[VRGN] + frac[next].lu[SCND] + frac[next].lu[PSTR] + frac[next].lu[CROP] + frac[next].lu[URBN]) > 1.00001 ){
     cout << "fractions exceed 1 \n";
   }
    
   
   for (size_t target=PRIMF; target<N_SUB_LU_TYPES; ++target) {
      frac[next].subLU[target] = frac[now].subLU[target];
      for (size_t source=PRIMF; source<N_SUB_LU_TYPES; ++source) {
         if (target != source) {
            if (target != PRIMF && target != PRIMN) {
            
               frac[next].subLU[target] += flow.getSub(source, target);
	       
            }
            frac[next].subLU[target] -= flow.getSub(target, source);
         }
      }
   }
   if ((frac[next].subLU[PRIMF] + frac[next].subLU[PRIMN]+frac[next].subLU[SECDF] + frac[next].subLU[SECDN] + frac[next].subLU[PASTR] +frac[next].subLU[RANGE] + frac[next].subLU[C3ANN]+ frac[next].subLU[C4ANN]+ frac[next].subLU[C3PER]+ frac[next].subLU[C4PER]+ frac[next].subLU[C3NFX] + frac[next].subLU[URBAN]) > 1.00001 ){
     cout << "sub fractions exceed 1 \n";
     int i=1;
      for (size_t target=PRIMF; target<N_SUB_LU_TYPES; ++target) {
	cout << "sublu " << i << " next: " << frac[next].subLU[target] << "\n";
	cout << "sublu " << i << " now: " << frac[now].subLU[target] << "\n";
	i++;
       }
      i=1;
       for (size_t target=VRGN; target<N_LU_TYPES; ++target) {
	cout << "lu " << i << " next: " << frac[next].lu[target] << "\n";
	cout << "lu " << i << " now: " << frac[now].lu[target] << "\n";
	i++;
       }
   }

        if (fabs(othr - frac[next].lu[VRGN] - frac[next].lu[SCND])>0.00001){
	  cout << "other hyde mismatch \n";
	    cout << "othr: " << othr << " computed othr: " << frac[next].lu[VRGN] + frac[next].lu[SCND] << "\n";
	}
        if (fabs(crop - frac[next].lu[CROP])>0.00001){
cout << "cropland hyde mismatch \n";
	}
	if (fabs(pstr - frac[next].lu[PSTR])>0.00001){
cout << "pasture hyde mismatch \n";
	}
	if (fabs(urbn - frac[next].lu[URBN])>0.00001){
cout << "urban hyde mismatch \n";
	}
	
        if (fabs(frac[next].subLU[SECDF] + frac[next].subLU[SECDN] - frac[next].lu[SCND])>0.0001){
	  cout << "secondary mismatch \n";
	}
	if (fabs(frac[next].subLU[PRIMF] + frac[next].subLU[PRIMN] - frac[next].lu[VRGN])>0.0001){
cout << "primary mismatch \n";
	}
        if (fabs(frac[next].subLU[C3ANN] + frac[next].subLU[C4ANN] + frac[next].subLU[C3NFX] + frac[next].subLU[C4PER] + frac[next].subLU[C3PER] - frac[next].lu[CROP])>0.0001){
cout << "cropland mismatch \n\n" << fabs(frac[next].subLU[C3ANN] + frac[next].subLU[C4ANN] + frac[next].subLU[C3NFX] + frac[next].subLU[C4PER] + frac[next].subLU[C3PER] - frac[next].lu[CROP]) << " \n";
 cout << "cropland now: " << frac[now].lu[CROP] << "\n";
 cout << "cropland next: " << frac[next].lu[CROP] << "\n";

	}
	
	if (fabs(frac[next].subLU[PASTR] + frac[next].subLU[RANGE] - frac[next].lu[PSTR])>0.0001){
cout << "pasture mismatch \n";
 
	}
	if (fabs(frac[next].subLU[URBAN] - frac[next].lu[URBN])>0.0001){
cout << "urban mismatch \n";
	}

	
     if (fabs(frac[next].subLU[SECDF] + frac[next].subLU[SECDN] - frac[next].lu[SCND])>0.00001){
     cout << "secondary mismatch. scnd: " <<  frac[next].lu[SCND] << " secdf: " << frac[next].subLU[SECDF] << " secdn: " << frac[next].subLU[SECDN] << "\n";
    }
   
   bool warn = config.STATE_BUG_PRINT;
   for (size_t lu=VRGN; lu<N_LU_TYPES; ++lu) {
      if (resetNegativeToZero(frac[next].lu[lu], warn,
                              (string("updateStates bug: ")
                               + getLandUseName((lu_types)lu)).c_str()) ) {
         printState(lu, now, next);
      }
   }

   for (size_t lu=PRIMF; lu<N_SUB_LU_TYPES; ++lu) {
      if (resetNegativeToZero(frac[next].subLU[lu], warn,
                              (string("updateStates bug: ")
                               + getSubLandUseName((sub_lu_types)lu)).c_str()) ) {
         printSubState(lu, now, next);
      }
   }

   if (config.STATE_PRINT) {
      printAllStates(now, next, "updateStates");
   }

}


/******************************************************************************/
void GridCell::updateSecondaryLands(size_t now, size_t next) {
   double sma_area_gained = 0.0, sma_area_lost = 0.0, sma_area_notlost = 0.0;
   for (size_t source=VRGN; source!=N_LU_TYPES; ++source) {

      sma_area_gained += flow.get(source, SCND);
      sma_area_lost += flow.get(SCND, source);
   }
   
   
   if (frac[now].scndMeanBiomass > ZEROVALUE_CHECK) {
      sma_area_gained += (flow.scndHarvBiomass1 + flow.scndHarvBiomass2
                          + flow.scndHarvBiomass3) / frac[now].scndMeanBiomass / cellArea;

      sma_area_lost += (flow.scndHarvBiomass1 + flow.scndHarvBiomass2 +
                        flow.scndHarvBiomass3) / frac[now].scndMeanBiomass / cellArea;
   }

   sma_area_notlost = frac[now].lu[SCND] - sma_area_lost;

   if ((frac[now].lu[SCND] + sma_area_gained - sma_area_lost) > ZEROVALUE_CHECK) {
      frac[next].scndMeanAge = (sma_area_notlost * (frac[now].scndMeanAge + 1.0)
                                + sma_area_gained) / (frac[now].lu[SCND] + sma_area_gained
                                                      - sma_area_lost);

   } else {
      frac[next].scndMeanAge = 1.0;
   }

   if (potentialBiomass > 0.0) {
      // 0.38 is wood fraction of NPP
      // 0.75 is the aboveground fraction of NPP
      frac[next].scndMeanBiomass = potentialBiomass * (1.0
                                                       - exp(-(potentialNPP * 0.75 * 0.38
                                                               * frac[next].scndMeanAge)
                                                             / potentialBiomass));
   } else {
      frac[next].scndMeanBiomass = 0.0;
   }
   
   if (resetNegativeToZero (frac[next].scndMeanAge, true, "Bug: sma")) {
      cerr << " s(t)=" << frac[now].lu[SCND]
           << " sma(t)=" << frac[now].scndMeanAge
           << " sma(t+1)=" << frac[next].scndMeanAge
           << " smb(t)=" << frac[now].scndMeanBiomass
           << " sma_area_notlost=" << sma_area_notlost
           << " sma_area_lost=" << sma_area_lost
           << " sma_area_gained=" << sma_area_gained << "\n";
   }

}



/********************************************************************************/
void GridCell::disaggregateInitialLandUse(size_t now, GridInfo& grid, Config& config) {

  double temp_pastr, temp_range;
  
   // Split VRGN into PRIMF and PRIMN
   if (potentialForest) {
      frac[now].subLU[PRIMF] = frac[now].lu[VRGN];
      frac[now].subLU[SECDF] = frac[now].lu[SCND];
   } else {
      frac[now].subLU[PRIMN] = frac[now].lu[VRGN];
      frac[now].subLU[SECDN] = frac[now].lu[SCND];
   }

   if (config.futureRun | config.extensionRun){
     frac[now].subLU[C3ANN] = grid.restartC3Ann_->get(y,x);
     if (frac[now].subLU[C3ANN] > 2.0 | frac[now].subLU[C3ANN] < 0.0){
        frac[now].subLU[C3ANN] = 0;
     }
     frac[now].subLU[C3PER] = grid.restartC3Per_->get(y,x);
     if (frac[now].subLU[C3PER] > 2.0 | frac[now].subLU[C3PER] < 0.0){
        frac[now].subLU[C3PER] = 0;
     }
     frac[now].subLU[C4ANN] = grid.restartC4Ann_->get(y,x);
     if (frac[now].subLU[C4ANN] > 2.0 | frac[now].subLU[C4ANN] < 0.0){
        frac[now].subLU[C4ANN] = 0;
     }
     frac[now].subLU[C4PER] = grid.restartC4Per_->get(y,x);
     if (frac[now].subLU[C4PER] > 2.0 | frac[now].subLU[C4PER] < 0.0){
        frac[now].subLU[C4PER] = 0;
     }
     frac[now].subLU[C3NFX] = grid.restartC3NFx_->get(y,x);
     if (frac[now].subLU[C3NFX] > 2.0 | frac[now].subLU[C3NFX] < 0.0){
        frac[now].subLU[C3NFX] = 0;
     }
     
   }
   else{
     for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
       if (croptype_ratio_total > 0.0) {
	 frac[now].subLU[ct] = frac[now].lu[CROP] * croptype_ratio[ct];
       } else {
            if (config.cft_option == 1) {
               switch (ct) {
                  case glm::C3ANN: frac[now].subLU[ct] = frac[now].lu[CROP] * zoneTotal->c3AnnualFraction; break;
                  case glm::C4ANN: frac[now].subLU[ct] = frac[now].lu[CROP] * zoneTotal->c4AnnualFraction; break;
                  case glm::C3PER: frac[now].subLU[ct] = frac[now].lu[CROP] * zoneTotal->c3PerennialFraction; break;
                  case glm::C4PER: frac[now].subLU[ct] = frac[now].lu[CROP] * zoneTotal->c4PerennialFraction; break;
                  case glm::C3NFX: frac[now].subLU[ct] = frac[now].lu[CROP] * zoneTotal->NfixingFraction; break;
                  default: std::cerr << "Incorrect crop functional type\n";
               }
            } else {
               frac[now].subLU[ct] = frac[now].lu[CROP] * 1.0 / 5.0;
            }
       }
     }
   }

   frac[now].subLU[URBAN] = frac[now].lu[URBN];
   

   if (config.futureRun | config.extensionRun){
     frac[now].subLU[PASTR] = grid.restartPastr_->get(y,x);
     if (frac[now].subLU[PASTR] > 2.0 | frac[now].subLU[PASTR] < 0.0){
        frac[now].subLU[PASTR] = 0;
     }
     frac[now].subLU[RANGE] = grid.restartRange_->get(y,x);
     if (frac[now].subLU[RANGE] > 2.0 | frac[now].subLU[RANGE] < 0.0){
        frac[now].subLU[RANGE] = 0;
     }

         croptype_ratio[PASTR] = grid.future_pasture_fracs[y][x];
         if (croptype_ratio[PASTR]>1.000000000000){
            croptype_ratio[PASTR]=1.000000000000;
         }
         if (croptype_ratio[PASTR]<0.000000000000){
            croptype_ratio[PASTR]=0.000000000000;
         }
         croptype_ratio[RANGE] = 1 - croptype_ratio[PASTR];

   }
   else{
     temp_pastr =  frac[now].subLU[PASTR];
     temp_range =  frac[now].subLU[RANGE];
     if ((temp_pastr + temp_range) !=0){
     frac[now].subLU[PASTR] = temp_pastr/(temp_pastr + temp_range)*frac[now].lu[PSTR];
     frac[now].subLU[RANGE] = temp_range/(temp_pastr + temp_range)*frac[now].lu[PSTR];
     }
     else{
     frac[now].subLU[PASTR] = 0;
     frac[now].subLU[RANGE] = 0;
     }
    
   }

}


/********************************************************************************/
void GridCell::computeIrrigatedArea (size_t nowxt) {
   // Compute country and crop specific irrigated area
   
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (frac[nowxt].lu[CROP] > ZEROVALUE_CHECK) {
         if (frac[nowxt].lu[CROP]>frac[nowxt].irrgData) {
            frac[nowxt].irrig[ct] = frac[nowxt].irrgData / frac[nowxt].lu[CROP];
         } else {
            frac[nowxt].irrig[ct] = 1;
         }
      } else {
         frac[nowxt].irrig[ct] = 0;
      }
    
   }
}


/********************************************************************************/
void GridCell::computeIrrigatedAreaExtension (size_t nowxt) {
   // Compute country and crop specific irrigated area
   
for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      // Only assign fertilizer usage to CFT if grid cell has any CFT
      if (frac[nowxt].subLU[ct] > 1e-6) {
         switch (ct) {
            case glm::C3ANN: frac[nowxt].irrig[ct] = frac[nowxt].iamIrrig[C3ANN]; break;
            case glm::C4ANN: frac[nowxt].irrig[ct] = frac[nowxt].iamIrrig[C4ANN]; break;
            case glm::C3PER: frac[nowxt].irrig[ct] = frac[nowxt].iamIrrig[C3PER]; break;
            case glm::C4PER: frac[nowxt].irrig[ct] = frac[nowxt].iamIrrig[C4PER]; break;
            case glm::C3NFX: frac[nowxt].irrig[ct] = frac[nowxt].iamIrrig[C3NFX]; break;
         }
      }
      else {
         frac[nowxt].irrig[ct] = 0.0;
      }
   }
   
   
}


/********************************************************************************/
void GridCell::computeIrrigatedAreaFuture (size_t now, size_t next) {
   // Compute country and crop specific irrigated area

for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
  
       if (zoneTotal->irrigNeededChange <=0) {
          frac[next].irrig[ct] = frac[now].irrig[ct] * (1 + zoneTotal->irrigDecFrac);
       }
       else {
          if (frac[now].irrig[ct] > 0.000001){
             frac[next].irrig[ct] = frac[now].irrig[ct] + zoneTotal->irrigIncFrac1*(1-frac[now].irrig[ct]);
             
          }
          else {
             if (frac[next].subLU[ct] > 0.000001){
                frac[next].irrig[ct] = zoneTotal->irrigIncFrac2;
             }
             else {
                frac[next].irrig[ct] = frac[now].irrig[ct];
             }
          }
       }
       
    
   }

}


/********************************************************************************/
void GridCell::computeFertilizerAmount (size_t nowxt) {
   // Compute country and crop functional type-specific fertilizer usage amount
   
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      // Only assign fertilizer usage to CFT if grid cell has any CFT
      if (frac[nowxt].subLU[ct] > 1e-6) {
         switch (ct) {
            case glm::C3ANN: frac[nowxt].fert[ct] = zoneTotal->fertDataC3ANN; break;
            case glm::C4ANN: frac[nowxt].fert[ct] = zoneTotal->fertDataC4ANN; break;
            case glm::C3PER: frac[nowxt].fert[ct] = zoneTotal->fertDataC3PER; break;
            case glm::C4PER: frac[nowxt].fert[ct] = zoneTotal->fertDataC4PER; break;
            case glm::C3NFX: frac[nowxt].fert[ct] = zoneTotal->fertDataC3NFX; break;
         }
      }
      else {
         frac[nowxt].fert[ct] = 0.0;
      }
   }
}


/********************************************************************************/
void GridCell::computeFertilizerAmountExtension (size_t nowxt) {
   // Compute country and crop functional type-specific fertilizer usage amount
   
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      // Only assign fertilizer usage to CFT if grid cell has any CFT
      if (frac[nowxt].subLU[ct] > 1e-6) {
         switch (ct) {
            case glm::C3ANN: frac[nowxt].fert[ct] = frac[nowxt].iamFert[C3ANN]; break;
            case glm::C4ANN: frac[nowxt].fert[ct] = frac[nowxt].iamFert[C4ANN]; break;
            case glm::C3PER: frac[nowxt].fert[ct] = frac[nowxt].iamFert[C3PER]; break;
            case glm::C4PER: frac[nowxt].fert[ct] = frac[nowxt].iamFert[C4PER]; break;
            case glm::C3NFX: frac[nowxt].fert[ct] = frac[nowxt].iamFert[C3NFX]; break;
         }
      }
      else {
         frac[nowxt].fert[ct] = 0.0;
      }
   }
}


/********************************************************************************/
void GridCell::computeFertilizerAmountFuture (size_t now, size_t next, Config& config) {
   // Compute country and crop functional type-specific fertilizer usage amount
   
   double regFert;
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      // Only assign fertilizer usage to CFT if grid cell has any CFT
       switch (ct) {
             case glm::C3ANN: regFert = zoneTotal->fertDataC3ANN; break;
             case glm::C4ANN: regFert = zoneTotal->fertDataC4ANN; break;
             case glm::C3PER: regFert = zoneTotal->fertDataC3PER; break;
             case glm::C4PER: regFert = zoneTotal->fertDataC4PER; break;
             case glm::C3NFX: regFert = zoneTotal->fertDataC3NFX; break;
          }
          
       if ( (frac[next].subLU[ct] > 1e-6) & (frac[now].fert[ct]<1) & (frac[now].iamFert[ct]>0) & (regFert>1e-6) & (config.fert_flag==1) ) {
         switch (ct) {
             case glm::C3ANN: frac[next].fert[ct] = frac[next].iamFert[ct]; break;
             case glm::C4ANN: frac[next].fert[ct] = frac[next].iamFert[ct]; break;
             case glm::C3PER: frac[next].fert[ct] = frac[next].iamFert[ct]; break;
             case glm::C4PER: frac[next].fert[ct] = frac[next].iamFert[ct]; break;
             case glm::C3NFX: frac[next].fert[ct] = frac[next].iamFert[ct]; break;
          }
       }
      
     else if (frac[next].subLU[ct] > 1e-6) {
         switch (ct) {
            case glm::C3ANN: frac[next].fert[ct] = frac[now].fert[ct] * (1 + zoneTotal->fertDataC3ANN/100); break;
            case glm::C4ANN: frac[next].fert[ct] = frac[now].fert[ct] * (1 + zoneTotal->fertDataC4ANN/100); break;
            case glm::C3PER: frac[next].fert[ct] = frac[now].fert[ct] * (1 + zoneTotal->fertDataC3PER/100); break;
            case glm::C4PER: frac[next].fert[ct] = frac[now].fert[ct] * (1 + zoneTotal->fertDataC4PER/100); break;
            case glm::C3NFX: frac[next].fert[ct] = frac[now].fert[ct] * (1 + zoneTotal->fertDataC3NFX/100); break;
         }
         
         if ((frac[next].fert[ct] > frac[next].iamFert[ct]) & (frac[next].fert[ct]>frac[now].fert[ct]) & (config.fert_flag==1) ){
          
frac[next].fert[ct] = frac[now].fert[ct];
         }
        
         if ((frac[next].fert[ct] < 0.8*frac[next].iamFert[ct])&(regFert>1e-6)&(frac[next].subLU[ct] > 1e-6)& (frac[now].iamFert[ct]>=0) & (config.fert_flag==1)) {
            frac[next].fert[ct] = frac[next].fert[ct] + 0.01*(frac[next].iamFert[ct]-frac[next].fert[ct]);
         }
         if (frac[next].fert[ct] < 0.0){
          frac[next].fert[ct] = 0.0;
         }
         if (frac[next].fert[ct]>500 & config.fert_flag==0){
            frac[next].fert[ct] = frac[now].fert[ct];
         }
      }
      else {
         frac[next].fert[ct] = 0.0;
      }
      
      if (frac[next].fert[ct]>500){
      }
   }
}


/********************************************************************************/
void GridCell::computeTillageRate (size_t now) {
   // 1. Check if there is any cropland
   // 2. Only assign tillage fraction to C3 Annuals, C4 Annuals and C3 N-Fixing

   double tilled_crops = frac[now].subLU[C3ANN] + frac[now].subLU[C4ANN] + frac[now].subLU[C3NFX];
   
   if (tilled_crops > ZEROVALUE_CHECK) {
      frac[now].tillage = zoneTotal->tillage;
   } else {
      frac[now].tillage = 0;
   }
}

/********************************************************************************/
void GridCell::precomputeFutureFloodedFraction (size_t now, size_t next) {
   // Management layer for flooded rice (subset of C3ANN crops)
  
  
    frac[next].riceData = frac[now].flooded*frac[now].subLU[C3ANN] + (frac[next].iamFlood - frac[now].iamFlood);
   
   
   
}
/********************************************************************************/
void GridCell::computeFloodedFraction (size_t nowxt) {
   // Management layer for flooded rice (subset of C3ANN crops)
  
      if (frac[nowxt].subLU[C3ANN] > 0.0) {
         frac[nowxt].flooded = frac[nowxt].riceData/frac[nowxt].subLU[C3ANN];
         // Amount of flooded rice cannot exceed C3Annual crops
         if (frac[nowxt].flooded > 1) {
            frac[nowxt].flooded = 1;
         }
         if (frac[nowxt].flooded < 0) {
            frac[nowxt].flooded = 0.0;
         }
      }
      else {
         frac[nowxt].flooded = 0.0;
      }
   
   
}

/********************************************************************************/
void GridCell::computeFloodedFractionExtension (size_t nowxt) {
   // Management layer for flooded rice (subset of C3ANN crops)
  
      if (frac[nowxt].subLU[C3ANN] > 0.0) {
         frac[nowxt].flooded = frac[nowxt].iamFlood;
      }
      else {
         frac[nowxt].flooded = 0.0;
      }
   
   
}
/********************************************************************************/
void GridCell::computeNationalFloodedFraction (size_t now, size_t next) {
   // Management layer for flooded rice (subset of C3ANN crops)
  
      if (frac[next].subLU[C3ANN] > 0.0) {
         frac[next].flooded = frac[now].flooded * frac[now].subLU[C3ANN]* (1+zoneTotal->floodedChange)/frac[next].subLU[C3ANN];
         // Amount of flooded rice cannot exceed C3Annual crops
         if (frac[next].flooded > 1) {
            frac[next].flooded = 1;
         }
         if (frac[next].flooded < 0) {
            frac[next].flooded = 0.0;
         }
      }
      else {
         frac[next].flooded = 0.0;
      }
   
   
}
/********************************************************************************/
void GridCell::computeWHBiofuelFractions (size_t now) {
   // Read country-specific data on the amount of fuelwood, use it to compute industrial roundwood fraction
   
   frac[now].fuelwoodFraction = zoneTotal->fuelwoodFraction;
   frac[now].commercialBiofuelsFraction = 0.0;
   frac[now].industrialRoundwoodFraction = 1.0 - frac[now].fuelwoodFraction - frac[now].commercialBiofuelsFraction;
}


/********************************************************************************/
void GridCell::computeCroplandBiofuelFractions (size_t now) {
   // Read country-specific information on the amount of crop functional type used to produce biofuel 
   
   for (size_t ct=C3ANN; ct<=C3NFX; ++ct) {
      if (frac[now].subLU[ct] > 0) {
         switch (ct) {
            case glm::C3ANN: frac[now].cropBiofuels[ct] = zoneTotal->cropBiofuelsC3ANN; break;
            case glm::C4ANN: frac[now].cropBiofuels[ct] = zoneTotal->cropBiofuelsC4ANN; break;
            case glm::C3PER: frac[now].cropBiofuels[ct] = zoneTotal->cropBiofuelsC3PER; break;
            case glm::C4PER: frac[now].cropBiofuels[ct] = zoneTotal->cropBiofuelsC4PER; break;
            case glm::C3NFX: frac[now].cropBiofuels[ct] = zoneTotal->cropBiofuelsC3NFX; break;
         }
      }
      else {
         frac[now].cropBiofuels[ct] = 0.0;
      }
   }
}

/********************************************************************************/
void GridCell::computeFutureCroplandBiofuelFractions (size_t now) {
   // Read country-specific information on the amount of crop functional type used to produce biofuel 
   
  
      if (frac[now].lu[CROP] > 0) {
        frac[now].totalCropBiofuel = frac[now].iamCropBiofuel;
      }
      else {
         frac[now].totalCropBiofuel = 0.0;
      }
      
      
}


/********************************************************************************/
void GridCell::computeBiomassHarvestFraction (size_t now) {
   // Using Harvest Index values from EPIC model
   // Reference: http://epicapex.tamu.edu/files/2015/10/EPIC.0810-User-Manual-Sept-15.pdf
   // Currently LUH2 1.0 historical release, sugarcane is the only c4 Perennial crop
   // Many more C3 Perennial crops, but they have a wide range of harvest indices, so taking average value
   frac[now].biomassHarvestFraction[C3PER] = 0.1;
   frac[now].biomassHarvestFraction[C4PER] = 0.9;
}