根据gsod数据集计算潜在蒸发量（Matlab）

     gsod数据集来自USAF DATSAV3地表数据和联邦气候综合地表小时(ISH)数据集中包含的天气观测数据。数据集包括全球超过9000多个站点的数据，数据类型包括：日最高温、日最低温、日平均温度、日露点温度、平均风速、最大持续风速、日降水量、积雪深度等。其中温度数据的单位为华氏摄氏度，精确到小数点后一位；风速数据单位为节，精确到小数点后一位；降水量单位为英寸，精确到小数点后两位。

   数据范围从1929年开始，其中1973年后的数据质量最高。数据下载地址：Global Surface Summary of the Day - GSOD (noaa.gov)

  潜在蒸发量的计算采用FAO Penman-Monteith公式，公式如下：

需要用到gsod数据集中的日最高温、日最低温、站点纬度、站点海拔、日露点温度、10m高处风速。具体的计算程序如下：

function ET0=PM_eva(Tmax,Tmin,Tdew,uz,z,phi,ele,J)
% 变量说明
% Tmax为日最高温，单位为℃；
% Tmin为日最低温，单位为℃；
% Tdew为日露点温度，单位为℃；
% uz为z米高处风速，单位为m/s；
% Phi为站点的纬度位置，单位为弧度；
% ele为站点的海拔，单位为m；
% J为一天在一年中的日序，如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高处风速，由zm高处风速转化而来
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