FAOPenman-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）
meta viewport
式中，r为湿度计常数，kPa/℃；p为⼤⽓压，kPa；λ为汽化潜热，2.45MJ/kg；Cp为常压下的⽐热，1.013×10-3MJ/（kg ℃）；ε为⽔蒸汽分⼦量与⼲燥空⽓分⼦量的⽐，其值为0.662。
净辐射缺测时的公式计算
⼤部分情况下没有净辐射的观测资料，根据FAO推荐的Penman-Monteith⽅法，净辐射可由以下公式计算。其中，Rn是进来的净短波辐射Rns与出去的净长波辐射RnL的差值：
净短波辐射Rns是由射⼊进来和反射出去的太阳辐射的平衡关系得出的：
正则表达式php
其中，α为冠层反射系数，以草为假想作物时取值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)javaswing在那个包下
# w = w_F(lat,DEC)
# N = N_F(w)
Q0 = Q0_F(dn,lat,DEC)
return Q0
"""
太阳总辐射 Q计算公式
"""
气象python零基础入门教程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]
linux下tomcat安装及配置教程
if x%4 == 0 and x%100 != 0 :  #判断闰年
runnian = 1
elif x%400 == 0:
runnian = 1
else:
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.s(a) + 0.070257*math.sin(a) - 0.s(2*a)
+ 0.000907*math.sin(2*a) - 0.s(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):
sql导出insert语句TR = 2 * math.pi * (dn - 1) / 365 #⽇⾓
#轨道订正系数
OCF =( 1.00011 + 0.034221 * s(TR) + 0.00128 * math.sin(TR) + 0.s(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) + s(lat) * s(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
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

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。