FAO Penman-Monteith公式（彭曼公式）计算参考蒸散量ET0的Python代码

        1998年，联合国粮农组织（FAO）正式提出用Penman-Monteith公式作为计算ET0的唯一标准方法。作者最近做模型需要计算这个参数，用Python写了计算的相关代码。

        根据FAO推荐的Penman-Monteith方法，ET0计算公式如下：

        以日为步长，则式中，ET0为参考蒸散量，mm/day；Rn为作物表面上的净辐射，MJ/（m2 day）；G为土壤热通量，MJ/（m2 day）；T为2米高处日平均气温，℃；u2为2米高处的风速，m/s；es为饱合水汽压，kPa；ea为实际水汽压，kPa；es-ea为饱和水汽压差，kPa；∆为饱和水汽压曲线的斜率；r为湿度计常数，kPa/℃。

def PM_ET0(Tmax,Tmin,Tmean,ea,P,u2,Rn): 
    #Tmax、Tmin、Tmean为最高、最低、平均气温；Tdew为露点温度；P为本站气压;u2为2m风速;Rn为净辐射
    Tmean = 0.5*(Tmax+Tmin)
    #土壤热通量
    G = 0
    #饱和水汽压
    es = 0.5*( e0(Tmax)+e0(Tmin) )
    #实际水汽压
    # ea = ea
    #饱和水汽压曲线斜率
    delta = 4098*e0(Tmean)/( (Tmean+237.3)*(Tmean+237.3) )
    #湿度计常数
    r = 0.665*0.001*P
    # ET0 = ( 0.408*delta*(Rn-G) + r*900*u2*(es-ea)/(Tmean+273) )/( delta + r*(1 + 0.34*u2) )
    
    ET0_1 = 0.408 * delta * (Rn-G)
    ET0_2 = r * 900* u2 * (es-ea)/(Tmean+273)
    ET0_3 = delta + r * (1 + 0.34*u2)
    ET0 = (ET0_1 + ET0_2)/ET0_3
    
    return ET0

（1）土壤热通量（G）

        土壤热通量可用复杂的模型描述。因为土壤热通量相对于Rn来说很小，尤其是在表面被植物覆盖和计算时间段是24小时或更长的时候。土壤温度随大气温度变化，对于长时段步长，计算G可用简单公式：

        式中：G为土壤热通量，MJ/（m2 day）；Cs为土壤热容量，MJ/（m3 ℃）；Ti为第i时刻大气温度，℃；Ti-1为第i-1时刻大气温度，℃；∆t为时间步长，day；∆Z为有效土壤深度，m。因为参考作物表面之下以1天或10天为时段的土壤热通量值相对小，因此它可以忽略，所以有Gday≈0。

（2）2米高处日平均气温（T）

        在温度变化对气象参数值影响很小的情况下，如果没有日平均气温观测值，日平均日气温(Tmean)按FAO Penman-Monteith公式计算饱和水汽压关系曲线(∆)的斜率和平均空气密度的影响(Pa)。为了标准化，以24小时为时段的T被定义为日最高温度(Tmax)和日最低温度(Tmin)的均值,而不是以小时观测温度的平均值，即：

（3）饱和水汽压（es）

        用平均气温代替日最高和最低气温，会使饱和水汽压估计值偏低，相应的水汽压差(表示大气蒸发能力的参数)也变小，也会出现低估参照腾发的情况。因此，平均饱和水汽压必须用与日最高和最低气温对应的饱和水汽压之间的平均值来计算。由于饱和水汽压与空气温度有关,它可用空气温度计算。关系式为：

def e0(T):
    e = 0.6108*(math.e**( (17.27*T/(T+237.3))))
    return e

        式中，e0（T）为空气温度T时的水汽压，kPa。因为上式是非线性函数，必须以日平均最高和最低气温的平均值来计算：        

 （4）饱和水汽压曲线的斜率（∆）

（5）露点温度计算实际水汽压（ea）

        露点温度是空气冷却使空气内的水汽含量达到饱和时的温度。实际水汽压(ea)在露点温度时是饱和水汽压。

（6）湿度计常数（r）

        式中，r为湿度计常数，kPa/℃；p为大气压，kPa；λ为汽化潜热，2.45MJ/kg；Cp为常压下的比热，1.013×10-3MJ/（kg ℃）；ε为水蒸汽分子量与干燥空气分子量的比，其值为0.662。 

净辐射缺测时的公式计算

        大部分情况下没有净辐射的观测资料，根据FAO推荐的Penman-Monteith方法，净辐射可由以下公式计算。其中，Rn是进来的净短波辐射Rns与出去的净长波辐射RnL的差值：

        净短波辐射Rns是由射入进来和反射出去的太阳辐射的平衡关系得出的： 

        其中，α为冠层反射系数，以草为假想作物时取值0.23，无量纲；Rs为天文辐射，MJ /m2 day。Stefan-Boltzmann定律定量揭示了长波能量发射的速率与表面的绝对温度的4次方成比例。然而，由于物质的吸收和天空的向下辐射，从地球表面上离去的净能量通量小于由Stefan-Boltzmann定律给出的发射能通量。水蒸气，云团，二氧化碳和尘埃是长波辐射的吸收者和发射者，当估计从地球表面的净输出通量时我们应该知道它们的浓度。由于湿度和云团起着很重要的作用，因此在评价长波辐射的净输出通量时，对Stefan-Boltzmann定律中用这两个因素作校正。而假定其它长波辐射的吸收者浓度是常数，则有：

        其中，RnL为净长波辐射，MJ /m2 day；σ为Stefan-Boltzmann常数，取4.903*10-9 MJ/K4 m2 day；Tmax，k为24小时内最高绝对温度值，K；Tmin，k为24小时内最低绝对温度值，K；ea为实际水汽压，kPa；Rs为太阳总辐射，MJ /m2 day；Rso为晴空辐射，MJ /m2 day。在没有修正值情况下，晴空辐射可由以下公式计算：

        其中，Z为测站海拔，m；Ra为天文辐射，MJ /m2 day。 

        下面是计算各种辐射的代码，有部分刚接触python时候写的，方法比较蠢hh，可以在改进下。

"""
天文辐射 Q0计算公式
"""
def cal_Q0(year,month,day,lat):  
    #year,month,day分别为年、月、日；lat为纬度;
    dn = rixu(year,month,day)
    DEC = dec(dn)
    # w = w_F(lat,DEC)
    # N = N_F(w)
    Q0 = Q0_F(dn,lat,DEC)
    return Q0

"""
太阳总辐射 Q计算公式
"""
def cal_Q(Q0,Q_fact,N):
     #Q0为天文辐射；Q_fact为日照时数；N为可照时数
    n = Q_fact
    S1 = s1(n,N)  

    a = 0.248
    b = 0.752
    Q = Q0 * (a + b * S1)
    return Q


"""
定义列表第i个数之前所有元素的和
"""
def sum_list(items,i):  
    sum_numbers = 0
    c = 0
    if i == 1:
        return 0
    else:
        for x in items: 
            if c!= i-1:
                sum_numbers += x
                c += 1
            else:
                break
        return sum_numbers

"""
定义日序计算方法
"""   
def rixu(x,y,z):   
    #x为年,y为月,z为日
    run = [31,29,31,30,31,30,31,31,30,31,30,31]
    flat = [31,28,31,30,31,30,31,31,30,31,30,31]
    if x%4 == 0 and x%100 != 0 :   #判断闰年
        runnian = 1
    elif x%400 == 0:
        runnian = 1
    else:
        runnian = 0
    if runnian == 1:
        d = sum_list(run,y) + z
    else:
        d = sum_list(flat,y) + z
    return d

"""
定义太阳赤纬计算方法
"""
def dec(rixu):
    int(rixu)
    a = 2 * math.pi * (rixu - 1) / 365.2422
    d = ( 0.006918 - 0.399912*math.cos(a) + 0.070257*math.sin(a) - 0.006758*math.cos(2*a)
                     + 0.000907*math.sin(2*a) - 0.002697*math.cos(3*a) + 0.00148*math.sin(3*a) )
    float('%.2f' % d)
    declination = d
    #d = math.degrees(declination)
    return declination

"""
定义时角w的算法
"""
def w_F(lat,dec):
    t = -1 * math.tan((31.56 / 360 * 2 * math.pi)) * math.tan(dec)
    w = math.acos(t)
    return w

"""
定义N的算法
"""
def N_F(w):
    N = 24 * w / math.pi
    return N

"""
定义Q0的算法
"""
def Q0_F(dn,lat,dec):

    TR = 2 * math.pi * (dn - 1) / 365 #日角
    #轨道订正系数
    OCF =( 1.00011 + 0.034221 * math.cos(TR) + 0.00128 * math.sin(TR) + 0.000719*math.cos(2*TR)
            + 0.000077 * math.sin(2*TR) )
    t = -1 * math.tan((31.56 / 360 * 2 * math.pi)) * math.tan(dec)
    w = math.acos(t)
    lat = 31.56  / 360 * 2 * math.pi
    Q = 86400 * 1367 * OCF * (w * math.sin(lat)*math.sin(dec) + math.cos(lat) * math.cos(dec) 
                        * math.sin(w)) / (math.pi)
    return Q

"""
定义S1的算法
"""
def  s1(n,N):
    s = n / N
    return s

"""
净长波辐射计算公式
"""
def R_net_long(Tmax,Tmin,ea,Q,Q0,Z):        # Q=Rs Q0=Ra
    Tmax = 273.6 + Tmax
    Tmin = 273.6 + Tmin
    Rso = (0.75 + 0.00002*Z)*Q0
    Rs = Q
    RnL = 4.903*10e-9*0.25*(Tmax**4-Tmin**4)*(0.34-0.14*ea**0.5)*(1.35*Rs/Rso-0.35)
    return RnL

"""
净辐射计算公式
"""
def R_net(Q,RnL,a=0.23):
    Rs = (1-a)*Q
    Rns = Rs - RnL
    return Rns

来源：SWYTC