Python utilities.float_lt函数代码示例

    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