def jmax_and_vcmax_func(self, temp):
        """ Maximum rate of electron transport (jmax) and of rubisco activity

        Parameters:
        -----------
        temp : float
            air temperature

        Returns:
        --------
        jmax : float
            maximum rate of electron transport
        vcmax : float
            maximum rate of Rubisco activity
        """
        if float_gt(temp, 10.0):
            jmax = (self.params.jmaxna * (1.0 + (temp - 25.0) * (0.05 +
                    (temp - 25.0) * (-1.81 * 1E-3 + (temp - 25.0) *
                    (-1.37 * 1E-4)))))
            vcmax = (self.params.vcmaxna * (1.0 + (temp - 25.0) *
                    (0.0485 + (temp - 25.0) * (-6.93 * 1E-4 + (temp - 25.0) *
                    (-3.9 * 1E-5)))))
        elif float_gt(temp, 0.0):
            jmax = self.params.jmaxna * 0.0305 * temp
            vcmax = self.params.vcmaxna * 0.0238 * temp
        else:
            jmax = 0.0
            vcmax = 0.0

        return jmax, vcmax

    def partition_plant_litter_n(self, nsurf, nsoil):
        """ Partition litter N from the plant (surface) and roots into metabolic
        and structural pools  
        
        Parameters:
        -----------
        nsurf : float
            N input from surface pool
        nsoil : float
            N input from soil pool
        """
        # constant structural input n:c as per century
        if not self.control.strfloat:
        
            # dead plant litter -> structural pool

            # n flux -> surface structural pool
            self.fluxes.n_surf_struct_litter = (self.fluxes.surf_struct_litter / 
                                                self.params.structcn)
            # n flux -> soil structural pool
            self.fluxes.n_soil_struct_litter = (self.fluxes.soil_struct_litter / 
                                                self.params.structcn)                         
                                     
            # if not enough N for structural, all available N goes to structural
            if float_gt( self.fluxes.n_surf_struct_litter, nsurf):
                 self.fluxes.n_surf_struct_litter = nsurf
            if float_gt(self.fluxes.n_soil_struct_litter, nsoil):
                self.fluxes.n_soil_struct_litter = nsoil
        
        # structural input n:c is a fraction of metabolic
        else:
            c_surf_struct_litter = (self.fluxes.surf_struct_litter * 
                                    self.params.structrat +
                                    self.fluxes.surf_metab_litter)
            
            if float_eq(c_surf_struct_litter, 0.0):
                 self.fluxes.n_surf_struct_litter = 0.0
            else:
                 self.fluxes.n_surf_struct_litter = (nsurf * 
                                                     self.fluxes.surf_struct_litter *
                                                     self.params.structrat / 
                                                     c_surf_struct_litter)
            
            c_soil_struct_litter = (self.fluxes.soil_struct_litter * 
                                    self.params.structrat +
                                    self.fluxes.soil_metab_litter)
           
            if float_eq(c_soil_struct_litter, 0.0):
                self.fluxes.n_soil_struct_litter = 0.
            else:
                self.fluxes.n_soil_struct_litter = (nsurf * 
                                                    self.fluxes.soil_struct_litter *
                                                    self.params.structrat / 
                                                    c_soil_struct_litter)
        
        # remaining N goes to metabolic pools
        self.fluxes.n_surf_metab_litter = (nsurf - 
                                           self.fluxes.n_surf_struct_litter)
        self.fluxes.n_soil_metab_litter = (nsoil - 
                                           self.fluxes.n_soil_struct_litter)

    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 epsilon(self, asat, par, daylen, alpha):
        """ Canopy scale LUE using method from Sands 1995, 1996. 
        
        Sands derived daily canopy LUE from Asat by modelling the light response
        of photosysnthesis as a non-rectangular hyperbola with a curvature 
        (theta) and a quantum efficiency (alpha). 
        
        Assumptions of the approach are:
         - horizontally uniform canopy
         - PAR varies sinusoidally during daylight hours
         - extinction coefficient is constant all day
         - Asat and incident radiation decline through the canopy following 
           Beer's Law.
         - leaf transmission is assumed to be zero.
           
        * Numerical integration of "g" is simplified to 6 intervals. 

        Parameters:
        ----------
        asat : float
            Light-saturated photosynthetic rate at the top of the canopy
        par : float
            photosyntetically active radiation (umol m-2 d-1)
        daylen : float
            length of day (hrs).
        theta : float
            curvature of photosynthetic light response curve 
        alpha : float
            quantum yield of photosynthesis (mol mol-1)
            
        Returns:
        -------
        lue : float
            integrated light use efficiency over the canopy (umol C umol-1 PAR)

        References:
        -----------
        See assumptions above...
        * Sands, P. J. (1995) Australian Journal of Plant Physiology, 
          22, 601-14.

        """
        delta = 0.16666666667  # subintervals scaler, i.e. 6 intervals
        h = daylen * const.SECS_IN_HOUR  # number of seconds of daylight

        if float_gt(asat, 0.0):
            # normalised daily irradiance
            q = pi * self.params.kext * alpha * par / (2.0 * h * asat)
            integral_g = 0.0
            for i in xrange(1, 13, 2):
                sinx = sin(pi * i / 24.0)
                arg1 = sinx
                arg2 = 1.0 + q * sinx
                arg3 = sqrt((1.0 + q * sinx) ** 2.0 - 4.0 * self.params.theta * q * sinx)
                integral_g += arg1 / (arg2 + arg3) * delta
            lue = alpha * integral_g * pi
        else:
            lue = 0.0

        return lue

    def calc_infiltration(self, rain):
        """ Estimate "effective" rain, or infiltration I guess.

        Simple assumption that infiltration relates to leaf area
        and therefore canopy storage capacity (wetloss). Interception is
        likely to be ("more") erroneous if a canopy is subject to frequent daily
        rainfall I would suggest.

        Parameters:
        -------
        rain : float
            rainfall [mm d-1]

        """

        if float_gt(rain, 0.0):
            self.fluxes.interception = self.state.lai * self.params.wetloss
            self.fluxes.interception = clip(self.fluxes.interception,
                                            min=self.fluxes.interception,
                                            max=rain)

            self.fluxes.erain = (rain * self.params.rfmult -
                                    self.fluxes.interception)
        else:
            self.fluxes.interception = 0.0
            self.fluxes.erain = 0.0

    def calc_transpiration(self):
        """ units mm/day """

        if float_gt(self.fluxes.wue, 0.0):
            self.fluxes.transpiration = self.fluxes.gpp_gCm2 / self.fluxes.wue
        else:
            self.fluxes.transpiration = 0.0

    def calc_soil_evaporation(self, tavg, net_rad, press, daylen, sw_rad):
        """ Use Penman eqn to calculate top soil evaporation flux at the
        potential rate.

        Soil evaporation is dependent upon soil wetness and plant cover. The net
        radiation term is scaled for the canopy cover passed to this func and
        the impact of soil wetness is accounted for in the wtfac term. As the
        soil dries the evaporation component reduces significantly.

        Key assumptions from Ritchie...

        * When plant provides shade for the soil surface, evaporation will not
        be the same as bare soil evaporation. Wind speed, net radiation and VPD
        will all belowered in proportion to the canopy density. Following
        Ritchie role ofwind, VPD are assumed to be negligible and are therefore
        ignored.

        These assumptions are based on work with crops and whether this holds
        for tree shading where the height from the soil to the base of the
        crown is larger is questionable.

        units = (mm/day)

        References:
        -----------
        * Ritchie, 1972, Water Resources Research, 8, 1204-1213.

        Parameters:
        -----------
        tavg : float
            average daytime temp [degC]
        net_rad : float
            net radiation [mj m-2 day-1]
        press : float
            average daytime pressure [kPa]

        Returns:
        --------
        soil_evap : float
            soil evaporation [mm d-1]

        """
        P = Penman()
        soil_evap = P.calc_evaporation(net_rad, tavg, press)
        
        # Surface radiation is reduced by overstory LAI cover. This empirical
        # fit comes from Ritchie (1972) and is formed by a fit between the LAI
        # of 5 crops types and the fraction of observed net radiation at the
        # surface. Whilst the LAI does cover a large range, nominal 0–6, there
        # are only 12 measurements and only three from LAI > 3. So this might
        # not hold as well for a forest canopy?
        # Ritchie 1972, Water Resources Research, 8, 1204-1213.
        if float_gt(self.state.lai, 0.0):
            soil_evap *= exp(-0.398 * self.state.lai)
        
        # reduce soil evaporation if top soil is dry
        soil_evap *= self.state.wtfac_topsoil
        tconv = 60.0 * 60.0 * daylen # seconds to day
        
        return soil_evap * tconv

    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 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 inputs_from_structrual_pool(self, nsurf, nsoil):
        """structural pool input fluxes

        Parameters:
        -----------
        nsurf : float
            N input from surface pool
        nsoil : float
            N input from soil pool
        """
        # constant structural input n:c as per century
        
        if not self.control.strfloat:
            # dead plant -> structural

            # surface
            self.fluxes.nresid[0] = self.fluxes.cresid[0] / self.params.structcn

            # soil
            self.fluxes.nresid[1] = self.fluxes.cresid[1] / self.params.structcn

            # if not enough N for structural, all available N goes to structural
            if float_gt(self.fluxes.nresid[0], nsurf):
                self.fluxes.nresid[0] = nsurf
            if float_gt(self.fluxes.nresid[1], nsoil):
                self.fluxes.nresid[1] = nsoil
        else:
            # structural input n:c is a fraction of metabolic
            cwgtsu = (self.fluxes.cresid[0] * self.params.structrat +
                      self.fluxes.cresid[2])
            if float_eq(cwgtsu, 0.0):
                self.fluxes.nresid[0] = 0.0
            else:
                self.fluxes.nresid[0] = (nsurf * self.fluxes.cresid[0] *
                                         self.params.structrat / cwgtsu)

            cwgtsl = (self.fluxes.cresid[1] * self.params.structrat +
                      self.fluxes.cresid[3])
            if float_eq(cwgtsl, 0.0):
                self.fluxes.nresid[1] = 0.
            else:
                self.fluxes.nresid[1] = (nsurf * self.fluxes.cresid[1] *
                                         self.params.structrat / cwgtsl)

    def calculate_top_of_canopy_n(self):
        """ Calculate the canopy N at the top of the canopy (g N m-2), N0.
        See notes and Chen et al 93, Oecologia, 93,63-69. 
        """

        if float_gt(self.state.lai, 0.0):
            # calculation for canopy N content at the top of the canopy
            N0 = self.state.ncontent * self.params.kext / (1.0 - exp(-self.params.kext * self.state.lai))
        else:
            N0 = 0.0
        return N0

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 adjust_cproduction(self, option):
        """ select model?

        It seems hybrid is the same as biomass, so check this and remove. I
        have set it up so that a call to hybrid just calls biomass

        Parameters:
        -----------
        option : integer
            model option

        """
        if option == 1 or option == 3:
            self.monteith_rescap_model()
        elif option == 2 and float_gt(self.params.fwpmax, self.params.fwpmin):
            self.biomass_model()
        elif option == 3 and float_gt(self.params.fwpmax, self.params.fwpmin):
            self.hybrid_model() # same as calling biomass model!
        else:
            err_msg = "Unknown water bal model (try 1-3): %s\n" % option
            raise RuntimeError, err_msg

    def calc_evaporation(self, vpd, wind, gs, net_rad, tavg, press, canht=None, 
                         ga=None):

        """
        Parameters:
        -----------
        vpd : float
            vapour pressure def [kPa]
        wind : float
            average daytime wind speed [m s-1]
        gs : float
            stomatal conductance [m s-1]
        net_rad : float
            net radiation [mj m-2 s-1] 
        tavg : float
            daytime average temperature [degC]
        press : float
            average daytime pressure [kPa]

        Returns:
        --------
        et : float
            evapotranspiration [mm d-1]

        """
        # if not read from met file calculate atmospheric pressure from sea lev
        if press == None:
            press = self.calc_atmos_pressure()
        
        lambdax = self.calc_latent_heat_of_vapourisation(tavg)
        gamma = self.calc_pyschrometric_constant(lambdax, press)
        slope = self.calc_slope_of_saturation_vapour_pressure_curve(tavg)
        rho = self.calc_density_of_air(tavg)
        if ga is None:
            ga = self.canopy_boundary_layer_conductance(wind, canht)
       
        if float_gt(gs, 0.0):
            # decoupling coefficent, Jarvis and McNaughton, 1986
            # when omega is close to zero, it is said to be well coupled and
            # gs is the dominant controller of water loss (gs<ga).
            e = slope / gamma # chg of latent heat relative to sensible heat of air
            omega = (e + 1.0) / (e + 1.0 + (ga / gs))
            
            arg1 = ((slope * net_rad ) + (rho * self.cp * vpd * ga))
            arg2 = slope + gamma * (1.0 + (ga / gs))
            et = (arg1 / arg2) / lambdax
        else:
            et = 0.0
            omega = 0.0
        
        return et, omega

    def calc_evaporation(self, vpd, wind, gs, net_rad, tavg, press):

        """
        Parameters:
        -----------
        vpd : float
            vapour pressure def [kPa]
        wind : float
            average daytime wind speed [m s-1]
        gs : float
            stomatal conductance [m s-1]
        net_rad : float
            net radiation [mj m-2 s-1] 
        tavg : float
            daytime average temperature [degC]
        press : float
            average daytime pressure [kPa]

        Returns:
        --------
        et : float
            evapotranspiration [mm d-1]

        """
        # if not read from met file calculate atmospheric pressure from sea lev
        if press == None:
            press = self.calc_atmos_pressure()
        
        lambdax = self.calc_latent_heat_of_vapourisation(tavg)
        gamma = self.calc_pyschrometric_constant(lambdax, press)
        slope = self.calc_slope_of_saturation_vapour_pressure_curve(tavg)
        rho = self.calc_density_of_air(tavg)
        
        ga = self.calc_atmos_boundary_layer_conductance(wind)
       
        # our model is a big leaf, so canopy conductance, gc = gs
        gc = gs
        
        if float_gt(gc, 0.0):
            # decoupling coefficent, Jarvis and McNaughton, 1986
            e = slope / gamma # chg of latent heat relative to sensible heat of air
            omega = (e + 1.0) / (e + 1.0 + (ga / gc))
            
            arg1 = ((slope * net_rad ) + (rho * self.cp * vpd * ga))
            arg2 = slope + gamma * (1.0 + ga / gc)
            et = (arg1 / arg2) / lambdax
        else:
            et = 0.0
            omega = 0.0
        
        return et

    def epsilon(self, amax, par, daylen, alpha):
        """ Canopy scale LUE using method from Sands 1995, 1996. 
        
        Numerical integration of g is simplified to 6 intervals. Leaf 
        transmission is assumed to be zero.

        Parameters:
        ----------
        amax : float
            photosynthetic rate at the top of the canopy
        par : float
            incident photosyntetically active radiation
        daylen : float
            length of day in hours.
        theta : float
            curvature of photosynthetic light response curve 
        alpha : float
            quantum yield of photosynthesis (mol mol-1)
            
        Returns:
        -------
        lue : float
            integrated light use efficiency over the canopy

        References:
        -----------
        See assumptions above...
        * Sands, P. J. (1995) Australian Journal of Plant Physiology, 
          22, 601-14.
        * LUE stuff comes from Sands 1996

        """
        delta = 0.16666666667 # subintervals scaler, i.e. 6 intervals
        h = daylen * const.HRS_TO_SECS 
        theta = self.params.theta # local var
        
        if float_gt(amax, 0.0):
            q = pi * self.params.kext * alpha * par / (2.0 * h * amax)
            integral = 0.0
            for i in xrange(1, 13, 2):
                sinx = sin(pi * i / 24.)
                arg1 = sinx
                arg2 = 1.0 + q * sinx 
                arg3 = sqrt((1.0 + q * sinx)**2.0 - 4.0 * theta * q * sinx)
                integral += arg1 / (arg2 + arg3) * delta
            lue = alpha * integral * pi
        else:
            lue = 0.0
        
        return lue

    def calc_wue(self, vpd, ca, amb_co2):
        """water use efficiency

        Not sure of units conversions here, have to ask BM

        Parameters:
        -----------
        vpd : float
            average daily vpd [kPa]
        ca : float
            atmospheric co2, depending on flag set in param file this will be
            ambient or elevated. [umol mol-1]

        """
        if self.control.wue_model == 0:
            # Gday original implementation
            # (gC / kg H20)
            if float_gt(vpd, 0.0):
                # WUE Power law dependence on co2, Pepper et al 2005.
                
                co2_ratio = (ca / amb_co2)
                co2_adjustment = co2_ratio**self.params.co2_effect_on_wue

                # wue inversely proportional to daily mean vpd
                self.fluxes.wue = self.params.wue0 * co2_adjustment / vpd
            else:
                self.fluxes.wue = 0.0
        elif self.control.wue_model == 1 and self.control.assim_model == 7:
            conv = const.MOL_C_TO_GRAMS_C / const.MOL_WATER_TO_GRAMS_WATER
            self.fluxes.wue = (conv * 1000.0 * (ca * const.UMOL_TO_MOL *
                                (1.0 - self.fluxes.cica_avg) /
                                (1.6 * vpd / 101.0)))
                        #if self.fluxes.wue > 20.0: self.fluxes.wue = 20.0    # FIX THIS!!!
            # what is this? ask BM

        elif self.control.wue_model == 2:
            self.fluxes.wue = (self.params.wue0 * 0.27273 / vpd *
                                ca / amb_co2)
        elif self.control.wue_model == 3 :
            if float_eq(self.fluxes.transpiration, 0.0):
                self.fluxes.wue = 0.0
            else:
                self.fluxes.wue = (self.fluxes.gpp_gCm2 /
                                    self.fluxes.transpiration)
        else:
            raise AttributeError('Unknown WUE calculation option')

    def epsilon(self, amax, par, daylen, alpha):
        """ Canopy scale LUE using method from Sands 1995, 1996.

        Parameters:
        ----------
        amax : float
            photosynthetic rate at the top of the canopy
        par : float
            incident photosyntetically active radiation
        daylen : float
            length of day in hours.
        theta : float
            curvature of photosynthetic light response curve 
        alpha : float
            quantum yield of photosynthesis (mol mol-1)
            
        Returns:
        -------
        lue : float
            integrated light use efficiency over the canopy

        References:
        -----------
        See assumptions above...
        * Sands, P. J. (1995) Australian Journal of Plant Physiology, 22, 601-14.

        """
        
        if float_gt(amax, 0.0):
            q = (math.pi * self.params.kext * alpha * par /
                    (2.0 * daylen * const.HRS_TO_SECS * amax))

            # check sands but shouldn't it be 2 * q * sin x on the top?
            f = (lambda x: x / (1.0 + q * x + math.sqrt((1.0 + q * x)**2.0 -
                            4.0 * self.params.theta * q * x)))
            g = [f(math.sin(math.pi * i / 24.)) for i in xrange(1, 13, 2)]
            
            #Trapezoidal rule - seems more accurate
            gg = 0.16666666667 * sum(g)

            lue = alpha * gg * math.pi
        else:
            lue = 0.0

        return lue

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