close?all;
clear;
clc

%%??input

MetDy_pt?=?‘G:\outputdata\test‘;

[yr1yr2]?=?deal（19612016）;
outpt?=?‘G:\outputdata\ET‘;
outhd?=?‘SURF_CLI_CHN_MUL_ID_‘;
%%??operate
hd?=?‘SURF_CLI_CHN_MUL_ID_‘;
A=xlsread（‘G:\selectstations.xlsx‘）;
%?SURF_CLI_CHN_MUL_DAY_STATION為氣象臺站表
S=size（A1）;
%?StnIDs?=?A（:1）;??%?提取第一列站點(diǎn)ID
for?r=1:S
????S_ID=A（r1）;
?for?yr?=?yr1:yr2
????fl=dir（[MetDy_ptfilesephdnum2str（S_ID）‘_‘num2str（yr）‘.txt‘]）;
????tmp=dlmread（[MetDy_ptfilesepfl.name]）;
????m=size（tmp1）;
??for?R=1:m
????Y=tmp（R2）;?????%?緯度
????X=tmp（R3）;?????%?經(jīng)度
????Z=tmp（R4）;?????%?海拔
????h=10;?????????????%?hight?of?the?wind?speed?（風(fēng)速?高度）
????
????%?define?variables???（輸入氣象要素?變量）
????year=tmp（R5）;
????month=tmp（R6）;
????Tmean=tmp（R9）;
????Tmax=tmp（R10）;
????Tmin=tmp（R11）;
????RH=tmp（R13）;
????WS=tmp（R14）/10;?????
????SSD=tmp（R15）/10;????????
????u=WS*0.7480;?%?風(fēng)速?
%?????******************大氣壓*************
????Pre=101.3*（（293-0.065*Z）/293）^5.26;
????
%?????******************飽和水汽壓曲線斜率*************

????delta=（4098*0.6108*exp（17.27*Tmean/（Tmean+273.3）））/（（Tmean+273.3）^2）;
%?????pre=Pre;????????%??air?pressure
????
????
????%?********step?5:psychrometric?constant干濕常數(shù)********
????gamma=6.65E-4*Pre;
?????
????
????%********step?6:Delta?Term（DT）（auxiliary?for?calculation?for?Radiation?Term）********
????DT=delta/（delta+gamma*（1+0.34*u））;
????
????
????%?********step?7:Psi?Term（PT）（auxiliary?for?calculation?for?Wind?Term）********
????PT=gamma/（delta+gamma*（1+0.34*u））;
????%?********step?8:Tempreture?Term（TT）（auxiliary?for?calculation?for?Wind
????%?Term）********
????TT=900*（Tmean+273）^-1*u;
??????
????%?********step?9:mean?saturation?vapor?pressure?derived?from?air?tempreture（eS）********
????eTmax=0.6108*exp（17.27*Tmax/（Tmax+237.3））;
????eTmin=0.6108*exp（17.27*Tmin/（Tmin+237.3））;
????eS=（eTmax+eTmin）/2;
????
????%?********step?10:actual?vapor?pressure（eA）?derived?from?relative?hum