<commit_before>/*
Notes from <NAME>, Feb 2015

baseline on does not result in N loss.... This might be a useful way to deal
with nonvascular N input, alternately it might be indicate a problem in the
nitrogen cycle.

A fix to this problem is proposed, ie, offfpen ncycle..., so right now, the
minimum change will be best....

TEM6 does this, but dealing with nfixation, emissions, etc, is upcoming

Continue to deal with this as a closed N cycle, give nonvascular plants priority
in N uptake over soil.
*/


/*
 * Vegetation_Bgc.cpp
 *
 * Purpose: Calculating (IN)GPP, (IN)NPP, RM, RG, N UPTAKE
 *            and resulted C & N changes
 *
 * History:
 *   June 28, 2011, by <NAME>:
 *     (1) Recoding based on DOS-TEM's code;
 *     (2) Multiple vegetation C&N pools added for DVM
 *
 * Important:
 *   (1) Parameters are read from 'CohortLookup.cpp', and set to
 *       'bgcpar' (struct:: vegpar_bgc)
 *
 *   (2) Calibrated Parameters are also read from 'CohortLookup.cpp'
 *       initially, and set to 'calpar' (strut:: vegpar_cal)
 *
 *   (3) The calculation is for ONE PFT only, so when calling it must
 *       be set the important index:
 *          pfttype
 *
 *   (4) FIVE (5) data pointers must be initialized by calling
 *       corresponding 'set...' methods
 *          chtlu, ed, bd, sb, fd
 *
 *   (5) all calculations are NOT 'fpc' adjusted - so, the
 *       unit is PER projected ground
 *
 */

#include "../TEMLogger.h"
#include "../TEMUtilityFunctions.h"

#include "Vegetation_Bgc.h"

extern src::severity_logger< severity_level > glg;

Vegetation_Bgc::Vegetation_Bgc() {
  bgcpar.dc2n = 0.000519;
  //increase of c:n = 0.176 if the nitrogen of veg is
  //  decreased by 15% per doubling of co2
  //proportional change in c:n per change co2
  //0.15 increase in C:N with double co2
};

Vegetation_Bgc::~Vegetation_Bgc() {
};

//Yuan: the parameterization and initialization is done for one PFT
//set the bgc parameters from inputs
void Vegetation_Bgc::initializeParameter() {
  BOOST_LOG_SEV(glg, info) << "Initializing parameters in Vegetation_Bgc from chtlu (CohortLookup) values.";
  bgcpar.kc      = chtlu->kc[ipft];
  bgcpar.ki      = chtlu->ki[ipft];
  bgcpar.tmin    = chtlu->tmin[ipft];
  bgcpar.toptmin = chtlu->toptmin[ipft];
  bgcpar.toptmax = chtlu->toptmax[ipft];
  bgcpar.tmax    = chtlu->tmax[ipft];
  bgcpar.raq10a0 = chtlu->raq10a0[ipft];
  bgcpar.raq10a1 = chtlu->raq10a1[ipft];
  bgcpar.raq10a2 = chtlu->raq10a2[ipft];
  bgcpar.raq10a3 = chtlu->raq10a3[ipft];
  bgcpar.knuptake = chtlu->knuptake[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    bgcpar.cpart[i]   = chtlu->cpart[i][ipft];
  }

  bgcpar.c2na = chtlu->c2na[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    bgcpar.c2nb[i]   = chtlu->c2nb[i][ipft];
    bgcpar.c2neven[i] = chtlu->initc2neven[i][ipft];
    bgcpar.c2nmin[i] = chtlu->c2nmin[i][ipft];
  }

  bgcpar.labncon = chtlu->labncon[ipft];
  // Calibrated parameters for vegetation BGC
  calpar.cmax = chtlu->cmax[ipft];
  calpar.nmax = chtlu->nmax[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    calpar.cfall[i] = chtlu->cfall[i][ipft];
    calpar.nfall[i] = chtlu->nfall[i][ipft];
  }

  calpar.kra = chtlu->kra[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    calpar.krb[i] = chtlu->krb[i][ipft];
  }

  calpar.frg = chtlu->frg[ipft];
};

// set the initial states from inputs
void Vegetation_Bgc::initializeState() {
  double totvegn = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.c[i]    = chtlu->initvegc[i][ipft];
    bd->m_vegs.strn[i] = chtlu->initvegn[i][ipft]*0.95;
    totvegn += chtlu->initvegn[i][ipft];
  }

  bd->m_vegs.labn  = totvegn*0.05;
  bd->m_vegs.deadc = chtlu->initdeadc[ipft];
  bd->m_vegs.deadn = chtlu->initdeadn[ipft];
  //
  bd->m_vegs.call    = 0.;
  bd->m_vegs.strnall = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.call += bd->m_vegs.c[i];

    if(nfeed) {
      bd->m_vegs.strnall += bd->m_vegs.strn[i];
    }
  }

  if (nfeed) {
    bd->m_vegs.nall = bd->m_vegs.strnall + bd->m_vegs.labn;
  }
};

//set the initial states from restart inputs
void Vegetation_Bgc::set_state_from_restartdata(const RestartData & rdata) {
  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.c[i]    = rdata.vegc[i][ipft];
    bd->m_vegs.strn[i] = rdata.strn[i][ipft];
  }

  bd->m_vegs.labn = rdata.labn[ipft];
  bd->m_vegs.deadc = rdata.deadc[ipft];
  bd->m_vegs.deadn = rdata.deadn[ipft];
  //
  bd->m_vegs.call    = 0.;
  bd->m_vegs.strnall = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.call += bd->m_vegs.c[i];

    if(nfeed) {
      bd->m_vegs.strnall += bd->m_vegs.strn[i];
    }
  }

  if (nfeed) {
    bd->m_vegs.nall = bd->m_vegs.strnall + bd->m_vegs.labn;
  }
};

//Yuan: the calculation is done for PFTs one by one
void Vegetation_Bgc::prepareIntegration(const bool &nfeedback) {
  //option of N module
  this->set_nfeed(nfeedback);
  //canopy conductance for GPP
  double pet = ed->m_l2a.pet;    //

  if (pet<=0.0) {
    pet = 0.0001;
  }

  double eet = ed->m_l2a.eet;    //
  bd->m_vegd.gv = getGV(eet,  pet);
  // previous 'topt' is the average of previous 10 years
  // which is used as Topt in GPP tfactor AND updated every year in 'phenology'
  double prvtopt = 0.;
  deque<double> toptdeque = cd->toptque[ipft];
  int recnum=toptdeque.size();

  if (recnum>0) {
    for(int i=0; i<recnum; i++) {
      prvtopt += toptdeque[i]/recnum;
    }
  } else {  //no previous year 'topt'
    prvtopt = (bgcpar.tmin+bgcpar.tmax)/2.;
  }

  // temperature factor for GPP
  bd->m_vegd.ftemp = getTempFactor4GPP(ed->m_atms.ta, prvtopt);

  // temperature factor for plant respiration
  bd->m_vegd.raq10 = getRaq10(ed->m_atms.ta);

  // dleafc is for "delta" or "difference"
  // this is the difference between the yearly maximum leaf biomass
  // (potentially affected by drought) and the (current) monthly value of
  // leaf biomass.

  // leaf C requirement for growth, used to dynamically determine C allocation
  // yearly max. leaf C, adjusted by seasonal foliage growth index 'fleaf'
  dleafc = cd->m_vegd.maxleafc[ipft] * cd->m_vegd.fleaf[ipft];

  // C requirement of foliage growth at current timestep
  dleafc -= bd->m_vegs.c[I_leaf];
  dleafc = fmax(0.0, dleafc);
  // litter-falling seasonal adjustment

  //previous 10 year mean of growing season degree-day, using for
  //  normalizing current growing period needed for litterfalling
  double prvttime = 0.0;

  deque<double> ttimedeque = cd->prvgrowingttimeque[ipft];
  recnum=ttimedeque.size();

  if (recnum>0) {
    for(int i=0; i<recnum; i++) {
      prvttime += ttimedeque[i]/recnum;
    }
  } else {  //no previous year 'ttime'
    prvttime = (bgcpar.tmin+bgcpar.tmax)/2.;
  }

  /*TUCKER FEB 2015: some of the logic behind the functional type and temporal 
    adjustment of litterfall appears to be faulty. For now, simplifying model
    so that litterfall == cfall * vegc, for annual plants, cfall for all 
    components = 0.083 (= 1 per year), for deciduous plants, cfall for leaves
    = 0.083. These values are set in cmt_calparbgc.txt. No values should exceed 0.083.
  for (int i=0; i<NUM_PFT_PART; i++) {
    //assuming 'calpar.cfall' is the max. monthly fraction,
    //  we need to modify them for specific PFT types
    if (cd->m_veg.ifdeciwoody[ipft]) {
      //for deciduous woody species, leaf max. falling fraction is 1.0;
      //  stem/root will take the calibrated.
      calpar.cfall[I_leaf] = 1.0;
      calpar.nfall[I_leaf] = 1.0;
    }

    if (!cd->m_veg.ifperenial[ipft]) {
      //for annual species, leaf/stem/root max. falling/dying fraction is 1.0,
      //  no matter what calibrated
      calpar.cfall[i] = 1.0;
      calpar.nfall[i] = 1.0;
    }
  }


  //assuming 'calpar.cfall' is the max. monthly fraction,
  //  and allowing the following seasonal variation
  fltrfall = 1.; // non-growing season, max. litterfall assumed

  if (cd->m_vegd.growingttime[ipft]>0. && prvttime>0.) {
    fltrfall = fmin(1., cd->m_vegd.growingttime[ipft]/prvttime);
  }
  */

  // dead standing C falling
  //if(cd->yrsdist < 10.0) {//cd->mthsdist<108
  if(cd->mthsdist < 108) {
    if(bd->m_vegs.deadc0 > 0) {
      bd->m_v2soi.d2wdebrisc = bd->m_vegs.deadc0/9.0/12.0;

      if(bd->m_v2soi.d2wdebrisc > bd->m_vegs.deadc){
        bd->m_v2soi.d2wdebrisc = bd->m_vegs.deadc; 
        bd->m_vegs.deadc = 0;
      }
      else{
        bd->m_vegs.deadc -= bd->m_v2soi.d2wdebrisc;
      }

      if (nfeed==1) {
        bd->m_v2soi.d2wdebrisn = bd->m_vegs.deadn0/9.0/12.0;

        if(bd->m_v2soi.d2wdebrisn > bd->m_vegs.deadn){
          bd->m_v2soi.d2wdebrisn = bd->m_vegs.deadn;
          bd->m_vegs.deadn = 0;
        }
        else{
          bd->m_vegs.deadn -= bd->m_v2soi.d2wdebrisn;
        }

      }
    }
  } else {
    bd->m_vegs.deadc = 0.0;

    if (nfeed==1) {
      bd->m_vegs.deadn = 0.0;
    }
  }

  //assign states to temporary state variable
  tmp_vegs.deadc            = bd->m_vegs.deadc;

  if (nfeed) {
    tmp_vegs.deadn = bd->m_vegs.deadn;
  }

  tmp_vegs.call    = 0.;

  if (nfeed) {
    tmp_vegs.strnall = 0.;
  }

  for (int i=0; i<NUM_PFT_PART; i++) {
    tmp_vegs.c[i]    = bd->m_vegs.c[i];
    tmp_vegs.call   += bd->m_vegs.c[i];

    if (nfeed) {
      tmp_vegs.strn[i] = bd->m_vegs.strn[i];
      tmp_vegs.strnall+= bd->m_vegs.strn[i];
    }
  }

  if (nfeed) {
    tmp_vegs.labn = bd->m_vegs.labn;
    tmp_vegs.nall = bd->m_vegs.labn+bd->m_vegs.strnall;
  }
};

// C and N fluxes without N limitation
void Vegetation_Bgc::delta() {
  /// some environment variables
  double co2 = ed->m_atms.co2;
  double par = ed->m_a2l.par;
  //leaf phenology
  // 1) current EET and previous max. EET controlled
  double fleaf = cd->m_vegd.fleaf[ipft];
  // 2) plant size (biomass C) or age controlled
  double ffoliage = cd->m_vegd.ffoliage[ipft];
  //GPP without N limitation
  double ingppall = getGPP(co2, par, fleaf, ffoliage,
                           bd->m_vegd.ftemp, bd->m_vegd.gv);
  // maintainence respiration first estimation
  double rm = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    if (tmp_vegs.c[i]>0.) {
      bd->m_vegd.kr[i] = getKr(tmp_vegs.c[i], i);
      del_v2a.rm[i] = getRm(tmp_vegs.c[i], bd->m_vegd.raq10, bd->m_vegd.kr[i]);
    } else {
      del_v2a.rm[i]    = 0.;
    }

    rm += del_v2a.rm[i];
  }

  double rmadj = 1.0;  //used below for summarizing ingpp[]

  if (rm>ingppall && rm>0.0) {
    rmadj=ingppall/rm;
  }

  // total available C for allocation
  // if C assimilation available, first goes to maintenance respiration
  double innppall = fmax(0., ingppall-rm)/(1.0+calpar.frg);
  // NPP allocation to leaf estimated first
  double nppl = dleafc;
  // Leaves have the second priority for C assimilation
  del_a2v.innpp[I_leaf] = fmin(nppl, innppall);
  del_v2a.rg[I_leaf] = calpar.frg * del_a2v.innpp[I_leaf];
  double npprgl = del_a2v.innpp[I_leaf] + del_v2a.rg[I_leaf];

  // the rest goes to stem/root, assuming equal priority
  double innpprest = fmax(0., ingppall-rm-npprgl)/(1.0+calpar.frg);
  double cpartrest = 0.;

  for (int i=I_leaf+1; i<NUM_PFT_PART; i++) {
    cpartrest +=bgcpar.cpart[i];
  }

  // Handle stems and roots
  for (int i=I_leaf+1; i<NUM_PFT_PART; i++) {
    del_a2v.innpp[i] = 0.0;
    del_v2a.rg[i]    = 0.0;

    if (cpartrest > 0.0 && innppall > 0.0) {
      del_a2v.innpp[i] = innppall * bgcpar.cpart[i] / cpartrest;
      del_v2a.rg[i]    = calpar.frg * del_a2v.innpp[i];
    }
  }

  // summary of INGPP
  for (int i=0; i<NUM_PFT_PART; i++) {
    del_a2v.ingpp[i] = del_a2v.innpp[i]+del_v2a.rm[i]*rmadj+del_v2a.rg[i];
    //may be over the 'ingppall' due to 'rm', which are biomass-C based,
    //  if not adjusted by 'rmadj'
  }<commit_msg>bug fix: Change math expressions to avoid error with percentage.<commit_after>/*
Notes from <NAME>, Feb 2015

baseline on does not result in N loss.... This might be a useful way to deal
with nonvascular N input, alternately it might be indicate a problem in the
nitrogen cycle.

A fix to this problem is proposed, ie, offfpen ncycle..., so right now, the
minimum change will be best....

TEM6 does this, but dealing with nfixation, emissions, etc, is upcoming

Continue to deal with this as a closed N cycle, give nonvascular plants priority
in N uptake over soil.
*/


/*
 * Vegetation_Bgc.cpp
 *
 * Purpose: Calculating (IN)GPP, (IN)NPP, RM, RG, N UPTAKE
 *            and resulted C & N changes
 *
 * History:
 *   June 28, 2011, by <NAME>:
 *     (1) Recoding based on DOS-TEM's code;
 *     (2) Multiple vegetation C&N pools added for DVM
 *
 * Important:
 *   (1) Parameters are read from 'CohortLookup.cpp', and set to
 *       'bgcpar' (struct:: vegpar_bgc)
 *
 *   (2) Calibrated Parameters are also read from 'CohortLookup.cpp'
 *       initially, and set to 'calpar' (strut:: vegpar_cal)
 *
 *   (3) The calculation is for ONE PFT only, so when calling it must
 *       be set the important index:
 *          pfttype
 *
 *   (4) FIVE (5) data pointers must be initialized by calling
 *       corresponding 'set...' methods
 *          chtlu, ed, bd, sb, fd
 *
 *   (5) all calculations are NOT 'fpc' adjusted - so, the
 *       unit is PER projected ground
 *
 */

#include "../TEMLogger.h"
#include "../TEMUtilityFunctions.h"

#include "Vegetation_Bgc.h"

extern src::severity_logger< severity_level > glg;

Vegetation_Bgc::Vegetation_Bgc() {
  bgcpar.dc2n = 0.000519;
  //increase of c:n = 0.176 if the nitrogen of veg is
  //  decreased by 15% per doubling of co2
  //proportional change in c:n per change co2
  //0.15 increase in C:N with double co2
};

Vegetation_Bgc::~Vegetation_Bgc() {
};

//Yuan: the parameterization and initialization is done for one PFT
//set the bgc parameters from inputs
void Vegetation_Bgc::initializeParameter() {
  BOOST_LOG_SEV(glg, info) << "Initializing parameters in Vegetation_Bgc from chtlu (CohortLookup) values.";
  bgcpar.kc      = chtlu->kc[ipft];
  bgcpar.ki      = chtlu->ki[ipft];
  bgcpar.tmin    = chtlu->tmin[ipft];
  bgcpar.toptmin = chtlu->toptmin[ipft];
  bgcpar.toptmax = chtlu->toptmax[ipft];
  bgcpar.tmax    = chtlu->tmax[ipft];
  bgcpar.raq10a0 = chtlu->raq10a0[ipft];
  bgcpar.raq10a1 = chtlu->raq10a1[ipft];
  bgcpar.raq10a2 = chtlu->raq10a2[ipft];
  bgcpar.raq10a3 = chtlu->raq10a3[ipft];
  bgcpar.knuptake = chtlu->knuptake[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    bgcpar.cpart[i]   = chtlu->cpart[i][ipft];
  }

  bgcpar.c2na = chtlu->c2na[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    bgcpar.c2nb[i]   = chtlu->c2nb[i][ipft];
    bgcpar.c2neven[i] = chtlu->initc2neven[i][ipft];
    bgcpar.c2nmin[i] = chtlu->c2nmin[i][ipft];
  }

  bgcpar.labncon = chtlu->labncon[ipft];
  // Calibrated parameters for vegetation BGC
  calpar.cmax = chtlu->cmax[ipft];
  calpar.nmax = chtlu->nmax[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    calpar.cfall[i] = chtlu->cfall[i][ipft];
    calpar.nfall[i] = chtlu->nfall[i][ipft];
  }

  calpar.kra = chtlu->kra[ipft];

  for (int i=0; i<NUM_PFT_PART; i++) {
    calpar.krb[i] = chtlu->krb[i][ipft];
  }

  calpar.frg = chtlu->frg[ipft];
};

// set the initial states from inputs
void Vegetation_Bgc::initializeState() {
  double totvegn = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.c[i]    = chtlu->initvegc[i][ipft];
    bd->m_vegs.strn[i] = chtlu->initvegn[i][ipft]*0.95;
    totvegn += chtlu->initvegn[i][ipft];
  }

  bd->m_vegs.labn  = totvegn*0.05;
  bd->m_vegs.deadc = chtlu->initdeadc[ipft];
  bd->m_vegs.deadn = chtlu->initdeadn[ipft];
  //
  bd->m_vegs.call    = 0.;
  bd->m_vegs.strnall = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.call += bd->m_vegs.c[i];

    if(nfeed) {
      bd->m_vegs.strnall += bd->m_vegs.strn[i];
    }
  }

  if (nfeed) {
    bd->m_vegs.nall = bd->m_vegs.strnall + bd->m_vegs.labn;
  }
};

//set the initial states from restart inputs
void Vegetation_Bgc::set_state_from_restartdata(const RestartData & rdata) {
  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.c[i]    = rdata.vegc[i][ipft];
    bd->m_vegs.strn[i] = rdata.strn[i][ipft];
  }

  bd->m_vegs.labn = rdata.labn[ipft];
  bd->m_vegs.deadc = rdata.deadc[ipft];
  bd->m_vegs.deadn = rdata.deadn[ipft];
  //
  bd->m_vegs.call    = 0.;
  bd->m_vegs.strnall = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    bd->m_vegs.call += bd->m_vegs.c[i];

    if(nfeed) {
      bd->m_vegs.strnall += bd->m_vegs.strn[i];
    }
  }

  if (nfeed) {
    bd->m_vegs.nall = bd->m_vegs.strnall + bd->m_vegs.labn;
  }
};

//Yuan: the calculation is done for PFTs one by one
void Vegetation_Bgc::prepareIntegration(const bool &nfeedback) {
  //option of N module
  this->set_nfeed(nfeedback);
  //canopy conductance for GPP
  double pet = ed->m_l2a.pet;    //

  if (pet<=0.0) {
    pet = 0.0001;
  }

  double eet = ed->m_l2a.eet;    //
  bd->m_vegd.gv = getGV(eet,  pet);
  // previous 'topt' is the average of previous 10 years
  // which is used as Topt in GPP tfactor AND updated every year in 'phenology'
  double prvtopt = 0.;
  deque<double> toptdeque = cd->toptque[ipft];
  int recnum=toptdeque.size();

  if (recnum>0) {
    for(int i=0; i<recnum; i++) {
      prvtopt += toptdeque[i]/recnum;
    }
  } else {  //no previous year 'topt'
    prvtopt = (bgcpar.tmin+bgcpar.tmax)/2.;
  }

  // temperature factor for GPP
  bd->m_vegd.ftemp = getTempFactor4GPP(ed->m_atms.ta, prvtopt);

  // temperature factor for plant respiration
  bd->m_vegd.raq10 = getRaq10(ed->m_atms.ta);

  // dleafc is for "delta" or "difference"
  // this is the difference between the yearly maximum leaf biomass
  // (potentially affected by drought) and the (current) monthly value of
  // leaf biomass.

  // leaf C requirement for growth, used to dynamically determine C allocation
  // yearly max. leaf C, adjusted by seasonal foliage growth index 'fleaf'
  dleafc = cd->m_vegd.maxleafc[ipft] * cd->m_vegd.fleaf[ipft];

  // C requirement of foliage growth at current timestep
  dleafc -= bd->m_vegs.c[I_leaf];
  dleafc = fmax(0.0, dleafc);
  // litter-falling seasonal adjustment

  //previous 10 year mean of growing season degree-day, using for
  //  normalizing current growing period needed for litterfalling
  double prvttime = 0.0;

  deque<double> ttimedeque = cd->prvgrowingttimeque[ipft];
  recnum=ttimedeque.size();

  if (recnum>0) {
    for(int i=0; i<recnum; i++) {
      prvttime += ttimedeque[i]/recnum;
    }
  } else {  //no previous year 'ttime'
    prvttime = (bgcpar.tmin+bgcpar.tmax)/2.;
  }

  /*TUCKER FEB 2015: some of the logic behind the functional type and temporal 
    adjustment of litterfall appears to be faulty. For now, simplifying model
    so that litterfall == cfall * vegc, for annual plants, cfall for all 
    components = 0.083 (= 1 per year), for deciduous plants, cfall for leaves
    = 0.083. These values are set in cmt_calparbgc.txt. No values should exceed 0.083.
  for (int i=0; i<NUM_PFT_PART; i++) {
    //assuming 'calpar.cfall' is the max. monthly fraction,
    //  we need to modify them for specific PFT types
    if (cd->m_veg.ifdeciwoody[ipft]) {
      //for deciduous woody species, leaf max. falling fraction is 1.0;
      //  stem/root will take the calibrated.
      calpar.cfall[I_leaf] = 1.0;
      calpar.nfall[I_leaf] = 1.0;
    }

    if (!cd->m_veg.ifperenial[ipft]) {
      //for annual species, leaf/stem/root max. falling/dying fraction is 1.0,
      //  no matter what calibrated
      calpar.cfall[i] = 1.0;
      calpar.nfall[i] = 1.0;
    }
  }


  //assuming 'calpar.cfall' is the max. monthly fraction,
  //  and allowing the following seasonal variation
  fltrfall = 1.; // non-growing season, max. litterfall assumed

  if (cd->m_vegd.growingttime[ipft]>0. && prvttime>0.) {
    fltrfall = fmin(1., cd->m_vegd.growingttime[ipft]/prvttime);
  }
  */

  // dead standing C falling
  //if(cd->yrsdist < 10.0) {//cd->mthsdist<108
  if(cd->mthsdist < 108) {
    if(bd->m_vegs.deadc0 > 0) {
      bd->m_v2soi.d2wdebrisc = bd->m_vegs.deadc0/9.0/12.0;

      if(bd->m_v2soi.d2wdebrisc > bd->m_vegs.deadc){
        bd->m_v2soi.d2wdebrisc = bd->m_vegs.deadc; 
        bd->m_vegs.deadc = 0;
      }
      else{
        bd->m_vegs.deadc -= bd->m_v2soi.d2wdebrisc;
      }

      if (nfeed==1) {
        bd->m_v2soi.d2wdebrisn = bd->m_vegs.deadn0/9.0/12.0;

        if(bd->m_v2soi.d2wdebrisn > bd->m_vegs.deadn){
          bd->m_v2soi.d2wdebrisn = bd->m_vegs.deadn;
          bd->m_vegs.deadn = 0;
        }
        else{
          bd->m_vegs.deadn -= bd->m_v2soi.d2wdebrisn;
        }

      }
    }
  } else {
    bd->m_vegs.deadc = 0.0;

    if (nfeed==1) {
      bd->m_vegs.deadn = 0.0;
    }
  }

  //assign states to temporary state variable
  tmp_vegs.deadc            = bd->m_vegs.deadc;

  if (nfeed) {
    tmp_vegs.deadn = bd->m_vegs.deadn;
  }

  tmp_vegs.call    = 0.;

  if (nfeed) {
    tmp_vegs.strnall = 0.;
  }

  for (int i=0; i<NUM_PFT_PART; i++) {
    tmp_vegs.c[i]    = bd->m_vegs.c[i];
    tmp_vegs.call   += bd->m_vegs.c[i];

    if (nfeed) {
      tmp_vegs.strn[i] = bd->m_vegs.strn[i];
      tmp_vegs.strnall+= bd->m_vegs.strn[i];
    }
  }

  if (nfeed) {
    tmp_vegs.labn = bd->m_vegs.labn;
    tmp_vegs.nall = bd->m_vegs.labn+bd->m_vegs.strnall;
  }
};

// C and N fluxes without N limitation
void Vegetation_Bgc::delta() {
  /// some environment variables
  double co2 = ed->m_atms.co2;
  double par = ed->m_a2l.par;
  //leaf phenology
  // 1) current EET and previous max. EET controlled
  double fleaf = cd->m_vegd.fleaf[ipft];
  // 2) plant size (biomass C) or age controlled
  double ffoliage = cd->m_vegd.ffoliage[ipft];
  //GPP without N limitation
  double ingppall = getGPP(co2, par, fleaf, ffoliage,
                           bd->m_vegd.ftemp, bd->m_vegd.gv);
  // maintainence respiration first estimation
  double rm = 0.;

  for (int i=0; i<NUM_PFT_PART; i++) {
    if (tmp_vegs.c[i]>0.) {
      bd->m_vegd.kr[i] = getKr(tmp_vegs.c[i], i);
      del_v2a.rm[i] = getRm(tmp_vegs.c[i], bd->m_vegd.raq10, bd->m_vegd.kr[i]);
    } else {
      del_v2a.rm[i]    = 0.;
    }

    rm += del_v2a.rm[i];
  }

  double rmadj = 1.0;  //used below for summarizing ingpp[]

  if (rm>ingppall && rm>0.0) {
    rmadj=ingppall/rm;
  }

  // total available C for allocation
  // if C assimilation available, first goes to maintenance respiration

  double innppall = (ingppall - rm) * (1.0 - calpar.frg);
  if (innppall < 0.0) { innppall = 0.0; }

  // NPP allocation to leaf estimated first
  double nppl = dleafc;
  // Leaves have the second priority for C assimilation
  del_a2v.innpp[I_leaf] = fmin(nppl, innppall);
  del_v2a.rg[I_leaf] = calpar.frg * del_a2v.innpp[I_leaf];
  double npprgl = del_a2v.innpp[I_leaf] + del_v2a.rg[I_leaf];

  // Assign remainder of C (after respiration, and leaf growth) to
  // stems and roots.
  double innpprest = fmax(0.0, ingppall-rm-npprgl) * (1.0 - calpar.frg);
  double cpartrest = 0.;

  for (int i=I_leaf+1; i<NUM_PFT_PART; i++) {
    cpartrest +=bgcpar.cpart[i];
  }

  // Handle stems and roots
  for (int i=I_leaf+1; i<NUM_PFT_PART; i++) {
    del_a2v.innpp[i] = 0.0;
    del_v2a.rg[i]    = 0.0;

    if (cpartrest > 0.0 && innppall > 0.0) {
      del_a2v.innpp[i] = innppall * bgcpar.cpart[i] / cpartrest;
      del_v2a.rg[i]    = calpar.frg * del_a2v.innpp[i];
    }
  }

  // summary of INGPP
  for (int i=0; i<NUM_PFT_PART; i++) {
    del_a2v.ingpp[i] = del_a2v.innpp[i]+del_v2a.rm[i]*rmadj+del_v2a.rg[i];
    //may be over the 'ingppall' due to 'rm', which are biomass-C based,
    //  if not adjusted by 'rmadj'
  }
{'fraction_non_alphanumeric': 0.1486993603411514, 'fraction_numerical': 0.02183368869936034, 'mean_word_length': 2.8343688685415303, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 269, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6794980', 'n_tokens_mistral': 10433, 'n_tokens_neox': 9849, 'n_words': 2179}