高纬度地区潜在蒸散发计算（FAO Penman-Monteith方法）

1.Penman-Monteith公式介绍

       FAO Penman-Monteith 方法是世界粮农组织（FAO）推荐计算潜在蒸散量的方法。这里， 定义潜在蒸散量为一种假想参照作物冠层的蒸散速率。假设作物植株高度为0.12m，固定的 作物表面阻力为70m/s，反射率为0.23，非常类似于表面开阔、高度一致、生长旺盛、完全遮盖地面而水分充分适宜的绿色草地的蒸散量。其计算公式如（1）所示：

2.估计潜在蒸散发所用数据

        GSOD数据集中的日最高温、日最低温、日露点温度、10m高处风速（该数据已经过预处理，将单位转化为Penman-Monteith 公式所需单位）。

3.MATLAB计算潜在蒸散发

  本处自定义潜在蒸散发的计算函数，计算是直接输入参数调用本函数即可。

function ET0=PM_eva(Tmax,Tmin,Tdew,uz,z,phi,ele,J)
% Tmax为日最高温，单位℃
% Tmin为日最低温，单位℃
% Tdew为日露点温度，单位℃
% uz为z米处风速，单位m/s
% phi为站点纬度，单位rad(弧度制)
% ele为站点海拔，单位m
% J-日序数，取值范围为1到365或366，1月1日取日序为1 

% 大气参数
P=101.3* ( (293-0.0065*ele) /293 )^5.26 ;  % 大气压
gamma=0.665*10^-3*P;     % 湿度计常数

% 气温-平均气温
Tmean=( Tmax+ Tmin )/2 ;

% 平均饱和水汽压
e0_Tmax= 0.6108*exp( 17.27*Tmax / (Tmax+237.3) );
e0_Tmin= 0.6108*exp( 17.27*Tmin / (Tmin+237.3) );
es=( e0_Tmax + e0_Tmin )/2 ;

% 饱和水汽压曲线斜率(delta)
delta = 4098*(0.6108*exp(17.27*Tmean/(Tmean+237.3) ) )/(Tmean+237.3)^2;

% 实际水汽压
ea=  0.6108*exp( 17.27*Tdew / (Tdew+237.3) );

% 天顶辐射-Ra
% J-日序数，取值范围为1到365或366，1月1日取日序为1 
Gsc = 0.0820; % 太阳常数，取值为0.0820MJ/(m^2×min)
dr = 1 + 0.033*cos(2*pi*J/365); % 日地间相对距离的倒数
Sdelta = 0.408*sin(2*pi*J/365-1.39);  % 太阳磁偏角，单位为弧度rad
ws = acos( -tan(phi)*tan(Sdelta) ); % 日落时角，单位为弧度rad

Ra =  24*60/pi*Gsc*dr*( ws*sin(phi)*sin(Sdelta)+cos(phi)*cos(Sdelta)*sin(ws));

% 太阳辐射-Rs
Krs = 0.162; % 调整系数，陆地取0.162
Rs =Krs*Ra*(Tmax-Tmin)^0.5;

% 晴空太阳辐射-Rso
Rso = (0.75+2*10^(-5)*ele )*Ra;

% 净太阳辐射-Rns
alpha = 0.23;
Rns = (1-alpha)*Rs;

% 净长波辐射
sigma = 4.903*10^(-9);      
TmaxK = Tmax+273.16;  % 一天中的最高绝对温度
TminK = Tmin+273.16;  % 一天中的最低绝对温度
air_correction= (0.34-0.14*ea^(1/2)); % 空气湿度的订正项
cloud_correction= (1.35*Rs/Rso-0.35);

Rnl=sigma*(TmaxK^4+TminK^4)/2 *air_correction*cloud_correction;

% 净辐射-Rn
Rn = Rns -Rnl;

% 土壤热通量-G
G=0;

% 2m高处风速，由z米高处风速转化而来
u2 = uz*4.87/log(67.82*z-5.42);

ET0 = ( 0.408*delta*(Rn-G)+gamma*900/(Tmean+273)*u2*(es-ea) ) / ( delta + gamma*(1+0.34*u2) );
ET0=real(ET0 );
end