def adj_for_low_temp(self, param, Tk, lower_bound=0.0, upper_bound=10.0):
        """ 
        Function allowing Jmax/Vcmax to be forced linearly to zero at low T
        """
        Tc = Tk - const.DEG_TO_KELVIN

        if float_lt(Tc, lower_bound):
            param = 0.0
        elif float_lt(Tc, upper_bound):
            param *= (Tc - lower_bound) / (upper_bound - lower_bound)

        return param

    def decay_in_dry_soils(self, decay_rate, decay_rate_dry):
        """Decay rates (e.g. leaf litterfall) can increase in dry soil, adjust
        decay param

        Parameters:
        -----------
        decay_rate : float
            default model parameter decay rate [tonnes C/ha/day]
        decay_rate_dry : float
            default model parameter dry deacy rate [tonnes C/ha/day]

        Returns:
        --------
        decay_rate : float
            adjusted deacy rate if the soil is dry [tonnes C/ha/day]

        """
        # turn into fraction...
        smc_root = self.state.pawater_root / self.params.wcapac_root
        
        new_decay_rate = (decay_rate_dry - (decay_rate_dry - decay_rate) * 
                         (smc_root - self.params.watdecaydry) /
                         (self.params.watdecaywet - self.params.watdecaydry))

        if float_lt(new_decay_rate, decay_rate):
            new_decay_rate = decay_rate

        if float_gt(new_decay_rate, decay_rate_dry):
            new_decay_rate = decay_rate_dry

        return new_decay_rate

def nc_ratio(carbon_val, nitrogen_val, pool):
    """Calculate nitrogen:carbon ratios

    Parameters:
    ----------
    carbon_val : float
        C value
    nitrogen_val: float
        N value#

    Returns:
    --------
    value : float
        N:C ratio
    """
    if float_lt(carbon_val, 0.0):
        # Note, previously the else branch was set to 1E6...presumably this 
        # was a hack to deal with the scenario for example where 
        # self.state.metabsurf and self.state.metabsurfn both start at zero.  
        # This was fine as this ratio isn't used in the code. Since I have 
        # commented out these diagnostics we shouldn't end up here unless there 
        # really is an error!!
        msg = "Dianostic for %s pool N:C has invalid values C:%s, N:%s" % \
                (pool, carbon_val, nitrogen_val)
        raise ValueError(msg)
    return nitrogen_val / carbon_val

    def soil_temp_factor(self, project_day):
        """Soil-temperature activity factor (A9). Fit to Parton's fig 2a 

        Parameters:
        -----------
        project_day : int
            current simulation day (index)

        Returns:
        --------
        tfac : float
            soil temperature factor [degC]

        """
        tsoil = self.met_data['tsoil'][project_day]

        if float_gt(tsoil, 0.0):
            self.fluxes.tfac_soil_decomp = (0.0326 + 0.00351 * tsoil**1.652 - 
                                            (tsoil / 41.748)**7.19)
            if float_lt(self.fluxes.tfac_soil_decomp, 0.0):
                self.fluxes.tfac_soil_decomp = 0.0
        else:
            # negative number cannot be raised to a fractional power
            # number would need to be complex
            self.fluxes.tfac_soil_decomp = 0.0

        return self.fluxes.tfac_soil_decomp

    def nc_limit(self, cpool, npool, ncmin, ncmax):
        """ Release N to 'Inorgn' pool or fix N from 'Inorgn', in order to keep
        the  N:C ratio of a litter pool within the range 'ncmin' to 'ncmax'.

        Parameters:
        -----------
        cpool : float
            various C pool (state)
        npool : float
            various N pool (state)
        ncmin : float
            maximum N:C ratio
        ncmax : float
            minimum N:C ratio

        Returns:
        --------
        fix/rel : float
            amount of N to be added/released from the inorganic pool

        """
        nmax = cpool * ncmax
        nmin = cpool * ncmin
    
        if float_gt(npool, nmax):  #release
            rel = npool - nmax
            self.fluxes.nlittrelease += rel 
            return -rel
        elif float_lt(npool, nmin):   #fix
            fix = nmin - npool
            self.fluxes.nlittrelease -= fix
            return fix
        else:
            return 0.0 

    def soil_temp_factor(self, project_day):
        """Soil-temperature activity factor (A9).

        Parameters:
        -----------
        project_day : int
            current simulation day (index)

        Returns:
        --------
        tfac : float
            soil temperature factor [degC]

        """
        tsoil = self.met_data["tsoil"][project_day]

        if float_gt(tsoil, 0.0):
            tfac = 0.0326 + 0.00351 * tsoil ** 1.652 - (tsoil / 41.748) ** 7.19
            if float_lt(tfac, 0.0):
                tfac = 0.0
        else:
            # negative number cannot be raised to a fractional power
            # number would need to be complex
            tfac = 0.0

        return tfac

    def assim(self, ci, gamma_star, a1, a2):
        """Morning and afternoon calcultion of photosynthesis with the 
        limitation defined by the variables passed as a1 and a2, i.e. if we 
        are calculating vcmax or jmax limited.
        
        Parameters:
        ----------
        ci : float
            intercellular CO2 concentration.
        gamma_star : float
            CO2 compensation point in the abscence of mitochondrial respiration
        a1 : float
            variable depends on whether the calculation is light or rubisco 
            limited.
        a2 : float
            variable depends on whether the calculation is light or rubisco 
            limited.

        Returns:
        -------
        assimilation_rate : float
            assimilation rate assuming either light or rubisco limitation.
        """
        if float_lt(ci, gamma_star):
            return 0.0
        else:
            return a1 * (ci - gamma_star) / (a2 + ci)

def main():

    # sweep the cmd line
    options, args = cmdline_parser()

    from file_parser import initialise_model_data
    # pylint: disable=C0103
    # pylint: disable=C0324
    # pylint: disable=C0103

    fname = "gday"
    fdir = "/Users/mdekauwe/src/python/GDAY_model/params"
    (adj_control, adj_params,
        adj_state, adj_files,
        adj_fluxes, met_data) = initialise_model_data(fname, default_dir=fdir)

    # figure out photosynthesis
    P = PlantProdModel(adj_control, adj_params, adj_state, adj_fluxes, met_data)

    adj_state.lai = (adj_params.slainit * const.M2_AS_HA /
                            const.KG_AS_TONNES / adj_params.cfracts *
                            adj_state.shoot)

    # Specific LAI (m2 onesided/kg DW)
    adj_state.sla = adj_params.slainit



    adj_control.assim_model = 3

    num_days = len(met_data['doy'])
    for i in xrange(num_days):

        if float_lt(adj_state.lai, adj_params.lai_cover):
            gcover = adj_state.lai / adj_params.lai_cover
        else:
            gcover = 1.0

        adj_state.fapar = ((1.0 - exp(-adj_params.kext * adj_state.lai /
                    gcover)) * gcover)

        adj_state.shootnc = adj_state.shootn / adj_state.shoot

        P.run_sim(i)


        print adj_fluxes.gpp / const.HA_AS_M2 * const.TONNES_AS_G

        #print adj_state.lai


        # this is done in derive so do here
        # Specific LAI (m2 onesided/kg DW)
        adj_state.sla = (adj_state.lai / const.M2_AS_HA *
                            const.KG_AS_TONNES *
                            adj_params.cfracts / adj_state.shoot)


    return

    def carbon_production(self, project_day, daylen):
        """ Calculate GPP, NPP and plant respiration

        Parameters:
        -----------
        project_day : integer
            simulation day
        daylen : float
            daytime length (hrs)

        References:
        -----------
        * Jackson, J. E. and Palmer, J. W. (1981) Annals of Botany, 47, 561-565.
        """

        if self.state.lai > 0.0:
            # average leaf nitrogen content (g N m-2 leaf)
            leafn = (self.state.shootnc * self.params.cfracts /
                     self.state.sla * const.KG_AS_G)
            
            # total nitrogen content of the canopy
            self.state.ncontent = leafn * self.state.lai
        else:
            self.state.ncontent = 0.0
         
        # fractional ground cover.
        if float_lt(self.state.lai, self.params.lai_cover):
            frac_gcover = self.state.lai / self.params.lai_cover
        else:
            frac_gcover = 1.0

        # Radiance intercepted by the canopy, accounting for partial closure
        # Jackson and Palmer (1981), derived from beer's law
        if self.state.lai > 0.0:
            self.state.light_interception = ((1.0 - exp(-self.params.kext *
                                             self.state.lai / frac_gcover)) *
                                             frac_gcover)
        else:
            self.state.light_interception = 0.0
        
        if self.control.water_stress:
            # Calculate the soil moisture availability factors [0,1] in the 
            # topsoil and the entire root zone
            (self.state.wtfac_tsoil, 
             self.state.wtfac_root) = self.sm.calculate_soil_water_fac()
        else:
            # really this should only be a debugging option!
            self.state.wtfac_tsoil = 1.0
            self.state.wtfac_root = 1.0
            
        # Estimate photosynthesis 
        if self.control.assim_model == "BEWDY":
            self.bw.calculate_photosynthesis(frac_gcover, project_day, daylen)
        elif self.control.assim_model == "MATE":
            self.mt.calculate_photosynthesis(project_day, daylen)
        else:
            raise AttributeError('Unknown assimilation model')

def clip(value, min=None, max=None):
    """ clip value btw defined range """
    if float_lt(value, min):
        value = min
    elif float_gt(value, max):
        value = max
    return value
    

 
    
    

def day_length(date, latitude):
    """ Figure out number of sunlight hours, (hours day-1)
    
    Routine from sdgvm. date is a python object, see datetime library for 
    more info 
    
    Parameters:
    -----------
    date : date format string
        date object, yr/month/day
    latitude : float    
        latitude [degrees]
        
    Returns:
    --------
    dayl : float 
        daylength [hrs]
    
    """
    conv = math.pi / 180.0
    
    # day of year 1-365/366
    doy = int(date.strftime('%j'))
    
    # Total number of days in year
    if calendar.isleap(date.year):
        yr_days = 366.
    else:
        yr_days = 365.
    
    solar_declin = -23.4 * math.cos(conv * yr_days * (doy + 10.0) / yr_days)
    temx = -math.tan(latitude * conv) * math.tan(solar_declin * conv) 
    
    if float_lt(math.fabs(temx), 1.0):
        has = math.acos(temx) / conv
        dayl = 2.0 * has / 15.0
    elif float_gt(temx, 0.0):
        dayl = 0.0
    else:
        dayl = 24.0
    
    return dayl

    def grazer_inputs(self):
        """ Grazer inputs from faeces and urine, flux detd by faeces c:n """
        if self.control.grazing:
            self.params.faecesn = self.fluxes.faecesc / self.params.faecescn
        else:
            self.params.faecesn = 0.0

        # make sure faecesn <= total n input to soil from grazing
        arg = self.fluxes.neaten * self.params.fractosoil
        if float_gt(self.params.faecesn, arg):
            self.params.faecesn = self.fluxes.neaten * self.params.fractosoil

        # urine=total-faeces
        if self.control.grazing:
            self.fluxes.nurine = self.fluxes.neaten * self.params.fractosoil - self.params.faecesn
        else:
            self.fluxes.nurine = 0.0

        if float_lt(self.fluxes.nurine, 0.0):
            self.fluxes.nurine = 0.0

def calc_npp(M, control, params, state, fluxes, met_data):
    
    state.lai = (params.slainit * const.M2_AS_HA /
                            const.KG_AS_TONNES / params.cfracts *
                            state.shoot)
    
    # Specific LAI (m2 onesided/kg DW)
    state.sla = params.slainit
    
    
    year = str(control.startyear)
    month = str(control.startmonth)
    day = str(control.startday)
    datex = datetime.datetime.strptime((year + month + day), "%Y%m%d")
        
    npp = np.zeros(0)
    for project_day in xrange(365):
    
        state.shootnc = state.shootn / state.shoot
        state.ncontent = (state.shootnc * params.cfracts /
                                state.sla * const.KG_AS_G)
        daylen = day_length(datex, params.latitude)
        state.wtfac_root = 1.0
        #state.lai = laidata[project_day]
    
    
        if float_lt(state.lai, params.lai_cover):
            frac_gcover = state.lai / params.lai_cover
        else:
            frac_gcover = 1.0
    
        state.light_interception = ((1.0 - math.exp(-params.kext *
                                            state.lai / frac_gcover)) *
                                            frac_gcover)

        M.calculate_photosynthesis(project_day, daylen)
        
        npp = np.append(npp, fluxes.npp_gCm2)
        datex += datetime.timedelta(days=1)
    return npp.sum()

def mate_day_length(date, latitude):
    """ Calculate number of sunlight hours (units = h d-1)
    
    Routine comes from MATE, though not sure how right this is, are the 
    hemispheres inverted? Check
    
    Parameters:
    -----------
    date : date format string
        date object, yr/month/day
    latitude : float    
        latitude [degrees]
        
    Returns:
    --------
    dayl : float 
        daylength [hrs]
    
    """
    conv = math.pi / 180.0
    
    # day of year 1-365/366
    doy = int(date.strftime('%j'))
    
    # Total number of days in year
    if calendar.isleap(date.year):
        yr_days = 366.
    else:
        yr_days = 365.
    
    solar_dec = (23.4 * math.pi / 180.0 * math.cos(2.0 * math.pi / 
                    yr_days * (doy + 10.0)))
    
    if float_lt(latitude, 0.0):
        solar_dec *= -1.0
	 
    dayl = (math.acos(-math.tan(latitude * conv) * math.tan(solar_dec)) * 
	            24.0 / math.pi)

    return dayl

    def decay_in_dry_soils(self, decay_rate, decay_rate_dry):
        """Decay rates (e.g. leaf litterfall) can increase in dry soil, adjust
        decay param. This is based on field measurements by F. J. Hingston 
        (unpublished) cited in Corbeels.

        Parameters:
        -----------
        decay_rate : float
            default model parameter decay rate [tonnes C/ha/day]
        decay_rate_dry : float
            default model parameter dry deacy rate [tonnes C/ha/day]

        Returns:
        --------
        decay_rate : float
            adjusted deacy rate if the soil is dry [tonnes C/ha/day]
        
        Reference:
        ----------
        Corbeels et al. (2005) Ecological Modelling, 187, 449-474.
        
        """
        # turn into fraction...
        smc_root = self.state.pawater_root / self.params.wcapac_root

        new_decay_rate = decay_rate_dry - (decay_rate_dry - decay_rate) * (smc_root - self.params.watdecaydry) / (
            self.params.watdecaywet - self.params.watdecaydry
        )

        if float_lt(new_decay_rate, decay_rate):
            new_decay_rate = decay_rate

        if float_gt(new_decay_rate, decay_rate_dry):
            new_decay_rate = decay_rate_dry

        return new_decay_rate

    def carbon_production(self, date, day, daylen):
        """ Calculate GPP, NPP and plant respiration

        Parameters:
        -----------
        day : intefer
            simulation day
        date : date string object
            date object string (yr/mth/day)
        daylen : float
            daytime length (hrs)

        References:
        -----------
        * Jackson, J. E. and Palmer, J. W. (1981) Annals of Botany, 47, 561-565.
        """

        # leaf nitrogen content
        self.state.ncontent = (self.state.shootnc * self.params.cfracts /
                                self.state.sla * const.KG_AS_G)

        # fractional ground cover.
        if float_lt(self.state.lai, self.params.lai_cover):
            frac_gcover = self.state.lai / self.params.lai_cover
        else:
            frac_gcover = 1.0

        # Radiance intercepted by the canopy, accounting for partial closure
        # Jackson and Palmer (1981), derived from beer's law
        self.state.light_interception = ((1.0 - math.exp(-self.params.kext *
                                            self.state.lai / frac_gcover)) *
                                            frac_gcover)
        
        
        
        sys.path.append("/Users/mdekauwe/src/fortran/maestra")
        import physiol
        import utils
        import maestcom as m
        
        RDFIPT = 7.29966089E-02
        TUIPT = 0.71114331
        TDIPT = 0.16205148
        RNET = 17.933075
        WIND = 8.20037797E-02
        PAR = 17.746254
        TAIR = -0.55100000
        
        RH = 0.73264003
        VPD = 157.64043
        VMFD = 1.5811477
        PRESS = 99700.000
        SOILMOIST = 0.0000000
        
        IECO = 1
        
        
        
        TVJUP = -100.00000
        TVJDN = -100.00000
        THETA = 0.94999999
        AJQ = 0.30000001
        RD0 = 0.50000006
        Q10F = 0.10260000
        RTEMP = 28.000000
        DAYRESP = 0.60000002
        TBELOW = -100.00000
        MODELGS = 4
        GSREF = 0.0000000
        I0 = 0.0000000
        D0 = 2.66481364E-21
        VK1 = 0.0000000
        VK2 = 0.0000000
        VPD1 = 0.0000000
        VPD2 = 0.0000000
        VMFD0 = 0.0000000
        GSJA = 0.0000000
        GSJB = 0.0000000
        T0 = 0.0000000
        TREF = 0.0000000
        TMAX = 0.0000000
        SMD1 = 0.0000000
        SMD2 = 0.0000000
        SOILDATA = 0
        SWPEXP = 0.0000000
        G0 = 0.0000000
        D0L = 2.66246708E-44
        GAMMA = 0.0000000
        G1 = 3.0573249
        GK = 0.0000000
        
        
        
        GSC = 3.10860965E-02
        RD = 2.08401568E-02
        
        
        print 
        
        
#.........这里部分代码省略.........

    def update_plant_state(self, fdecay, rdecay):
        """ Daily change in C content

        Parameters:
        -----------
        fdecay : float
            foliage decay rate
        rdecay : float
            fine root decay rate

        """
        self.state.shoot += (self.fluxes.cpleaf - self.fluxes.deadleaves -
                                self.fluxes.ceaten)
        self.state.root += self.fluxes.cproot - self.fluxes.deadroots
        self.state.branch += self.fluxes.cpbranch - self.fluxes.deadbranch
        self.state.stem += self.fluxes.cpstem - self.fluxes.deadstems
        
        if self.control.deciduous_model:
            self.state.shootn += (self.fluxes.npleaf - 
                                 (self.fluxes.deadleafn - self.fluxes.neaten))
        else:
            self.state.shootn += (self.fluxes.npleaf - 
                                  fdecay * self.state.shootn - 
                                  self.fluxes.neaten)
        self.state.shootn = max(0.0, self.state.shootn)
        
        self.state.rootn += self.fluxes.nproot - rdecay * self.state.rootn
        self.state.branchn += (self.fluxes.npbranch - self.params.bdecay *
                                self.state.branchn)
        self.state.stemnimm += (self.fluxes.npstemimm - self.params.wdecay *
                                self.state.stemnimm)
        self.state.stemnmob += (self.fluxes.npstemmob - self.params.wdecay *
                                self.state.stemnmob -
                                self.params.retransmob * self.state.stemnmob)
        self.state.stemn = self.state.stemnimm + self.state.stemnmob
        
        
        
        
        # maximum leaf n:c ratio is function of stand age
        #  - switch off age effect by setting ncmaxfyoung = ncmaxfold
        age_effect = ((self.state.age - self.params.ageyoung) / 
                        (self.params.ageold - self.params.ageyoung))
        
        ncmaxf = (self.params.ncmaxfyoung - (self.params.ncmaxfyoung -
                    self.params.ncmaxfold) * age_effect)
               
        if float_lt(ncmaxf, self.params.ncmaxfold):
            ncmaxf = self.params.ncmaxfold

        if float_gt(ncmaxf, self.params.ncmaxfyoung):
            ncmaxf = self.params.ncmaxfyoung

        # if foliage or root n:c ratio exceeds its max, then nitrogen uptake is
        # cut back n.b. new ring n/c max is already set because it is related
        # to leaf n:c
        extrar = 0.
        extras = 0.
        if float_gt(self.state.shootn, (self.state.shoot * ncmaxf)):
            extras = self.state.shootn - self.state.shoot * ncmaxf

            #n uptake cannot be reduced below zero.
            if float_gt(extras, self.fluxes.nuptake):
                extras = self.fluxes.nuptake

            self.state.shootn -= extras
            self.fluxes.nuptake -= extras

        ncmaxr = ncmaxf * self.params.ncrfac  # max root n:c

        if float_gt(self.state.rootn, (self.state.root * ncmaxr)):
            extrar = self.state.rootn - self.state.root * ncmaxr

            #n uptake cannot be reduced below zero.
            if float_gt((extras + extrar), self.fluxes.nuptake):
                extrar = self.fluxes.nuptake - extras

            self.state.rootn -= extrar
            self.fluxes.nuptake -= extrar 
        
        if self.control.deciduous_model:
            # update annual fluxes - store for next year
            self.state.clabile_store += self.fluxes.npp
            self.state.aroot_uptake += self.fluxes.nuptake
            
            self.state.aretrans += self.fluxes.retrans
            self.state.anloss += self.fluxes.nloss
            
            # update N:C of plant pools
            if float_eq(self.state.shoot, 0.0):
                self.state.shootnc = 0.0
            else:
                self.state.shootnc = max(self.state.shootn / self.state.shoot, 
                                         self.params.ncfmin)
                
            # N:C of stuctural wood as function of N:C of foliage 
            # (kgN kg-1C)(Medlyn et al 2000)  
            self.state.nc_ws = (self.params.nc_wsa + self.params.nc_wsb * 
                                self.state.shootnc)    
           
#.........这里部分代码省略.........

    # days in each year
    years = uniq(met_data["year"])
    # years = years[:-1] # dump last year as missing "site" LAI

    days_in_year = [met_data["year"].count(yr) for yr in years]

    for i, yr in enumerate(years):
        daylen = calculate_daylength(days_in_year[i], params.latitude)
        for doy in xrange(days_in_year[i]):
            state.wtfac_root = 1.0
            # state.lai = lai_data[project_day]
            state.lai = 2.0
            if state.lai > 0.0:
                state.shootnc = 0.03  # state.shootn / state.shoot
                state.ncontent = state.shootnc * params.cfracts / params.sla * const.KG_AS_G
            else:
                state.ncontent = 0.0

            if float_lt(state.lai, params.lai_cover):
                frac_gcover = state.lai / params.lai_cover
            else:
                frac_gcover = 1.0

            state.fipar = (1.0 - exp(-params.kext * state.lai / frac_gcover)) * frac_gcover

            M.calculate_photosynthesis(project_day, daylen[doy])

            print fluxes.gpp_gCm2  # , state.lai

            project_day += 1

    def carbon_production(self, project_day, daylen):
        """ Calculate GPP, NPP and plant respiration

        Parameters:
        -----------
        project_day : integer
            simulation day
        daylen : float
            daytime length (hrs)

        References:
        -----------
        * Jackson, J. E. and Palmer, J. W. (1981) Annals of Botany, 47, 561-565.
        """

        # leaf nitrogen content
        if self.state.lai > 0.0:
            # Leaf N content (g m-2)                       
            self.state.ncontent = (self.state.shootnc * self.params.cfracts /
                                   self.state.sla * const.KG_AS_G)
        else:
            self.state.ncontent = 0.0
            
        # fractional ground cover.
        if float_lt(self.state.lai, self.params.lai_cover):
            frac_gcover = self.state.lai / self.params.lai_cover
        else:
            frac_gcover = 1.0

        # Radiance intercepted by the canopy, accounting for partial closure
        # Jackson and Palmer (1981), derived from beer's law
        if self.state.lai > 0.0:
            self.state.light_interception = ((1.0 - exp(-self.params.kext *
                                             self.state.lai / frac_gcover)) *
                                             frac_gcover)
        else:
            self.state.light_interception = 0.0
        
        if self.control.water_stress:
            # Calculate the soil moisture availability factors [0,1] in the 
            # topsoil and the entire root zone
            (self.state.wtfac_tsoil, 
                self.state.wtfac_root) = self.sm.calculate_soil_water_fac()
        else:
            # really this should only be a debugging option!
            self.state.wtfac_tsoil = 1.0
            self.state.wtfac_root = 1.0
            
        # Estimate photosynthesis using an empirical model
        if self.control.assim_model >=0 and self.control.assim_model <= 4:
            self.pp.calculate_photosynthesis(project_day)
        # Estimate photosynthesis using the mechanistic BEWDY model
        elif self.control.assim_model >=5 and self.control.assim_model <= 6:
            # calculate plant C uptake using bewdy
            self.bw.calculate_photosynthesis(frac_gcover, project_day, daylen)
        # Estimate photosynthesis using the mechanistic MATE model. Also need to
        # calculate a water availability scalar to determine Ci:Ca reln.
        elif self.control.assim_model ==7:
            self.mt.calculate_photosynthesis(project_day, daylen)
        else:
            raise AttributeError('Unknown assimilation model')

    def update_plant_state(self, fdecay, rdecay, project_day):
        """ Daily change in C content

        Parameters:
        -----------
        fdecay : float
            foliage decay rate
        rdecay : float
            fine root decay rate

        """
        self.state.shoot += (self.fluxes.cpleaf - self.fluxes.deadleaves -
                                self.fluxes.ceaten)
        self.state.root += self.fluxes.cproot - self.fluxes.deadroots
        self.state.branch += self.fluxes.cpbranch - self.fluxes.deadbranch
        self.state.stem += self.fluxes.cpstem - self.fluxes.deadstems
        
        if self.control.deciduous_model:
            self.state.shootn += (self.fluxes.npleaf - 
                                 (self.fluxes.deadleafn - self.fluxes.neaten))
        else:
            self.state.shootn += (self.fluxes.npleaf - 
                                  fdecay * self.state.shootn - 
                                  self.fluxes.neaten)
        self.state.shootn = max(0.0, self.state.shootn)
        
        self.state.rootn += self.fluxes.nproot - rdecay * self.state.rootn
        self.state.branchn += (self.fluxes.npbranch - self.params.bdecay *
                                self.state.branchn)
        self.state.stemnimm += (self.fluxes.npstemimm - self.params.wdecay *
                                self.state.stemnimm)
        self.state.stemnmob += (self.fluxes.npstemmob - self.params.wdecay *
                                self.state.stemnmob -
                                self.params.retransmob * self.state.stemnmob)
        self.state.stemn = self.state.stemnimm + self.state.stemnmob
        
        self.state.exu_pool += self.fluxes.cprootexudate
        self.fluxes.microbial_resp = self.calc_microbial_resp(project_day)
        self.state.exu_pool -= self.fluxes.microbial_resp
        #print self.fluxes.microbial_resp, self.fluxes.cprootexudate
        
        self.calculate_cn_store()
        
        if not self.control.deciduous_model:
            
            # maximum leaf n:c ratio is function of stand age
            #  - switch off age effect by setting ncmaxfyoung = ncmaxfold
            age_effect = ((self.state.age - self.params.ageyoung) / 
                            (self.params.ageold - self.params.ageyoung))
            
            ncmaxf = (self.params.ncmaxfyoung - (self.params.ncmaxfyoung -
                        self.params.ncmaxfold) * age_effect)
                   
            if float_lt(ncmaxf, self.params.ncmaxfold):
                ncmaxf = self.params.ncmaxfold
    
            if float_gt(ncmaxf, self.params.ncmaxfyoung):
                ncmaxf = self.params.ncmaxfyoung
    
            # if foliage or root n:c ratio exceeds its max, then nitrogen 
            # uptake is cut back n.b. new ring n/c max is already set because 
            # it is related to leaf n:c
            extrar = 0.
            extras = 0.
            if float_gt(self.state.shootn, (self.state.shoot * ncmaxf)):
                extras = self.state.shootn - self.state.shoot * ncmaxf
    
                #n uptake cannot be reduced below zero.
                if float_gt(extras, self.fluxes.nuptake):
                    extras = self.fluxes.nuptake
    
                self.state.shootn -= extras
                self.fluxes.nuptake -= extras
    
            ncmaxr = ncmaxf * self.params.ncrfac  # max root n:c
    
            if float_gt(self.state.rootn, (self.state.root * ncmaxr)):
                extrar = self.state.rootn - self.state.root * ncmaxr
    
                #n uptake cannot be reduced below zero.
                if float_gt((extras + extrar), self.fluxes.nuptake):
                    extrar = self.fluxes.nuptake - extras
    
                self.state.rootn -= extrar
                self.fluxes.nuptake -= extrar