penman计算潜在蒸散 matlab

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）;  % 提取第一列站点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）;     % 经度
    Z=tmp（R4）;     % 海拔
    h=10;             % hight of the wind speed （风速 高度）
    
    % 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; % 风速 
%     ******************大气压*************
    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干湿常数********
    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