【MATLAB】高塔观测数据通量梯度法计算感热及潜热通量

1. 数据导入
   利用MATLAB菜单栏的“导入数据”功能，将高塔观测数据导入
   要点

Loggernet将原始二进制数据转为ASCII
【data】【card convert】按文件夹选择

MATLAB导入数据 输出类型中选择【列向量】

导入65080 30min数据和9345flux的amb_press_Avg（气压）
因为二者观测时间并不完全一致，因此在计算中应将时间取二者并集

2. 运行代码

main.m

%flux gradient method 
%author:Iris Lee
M=96;%number of time series

%% define and input
kappa=0.4;%Karman constant
H1=20;
H2=60;%height
    
z1=H1-(2/3)*2;
z2=H2-(2/3)*2;%consider zero-plane displacement

for i=1:M
    u1 =WS_20m_Avg(i) ;
    u2 =WS_60m_Avg(i) ;%wind speed
    
    T1 =Ta_20m_Avg (i);
    T2 =Ta_60m_Avg (i);%tempurature(C)
    Tmean =mean([T1 ,T2 ]);
    
    theta1 =T1 +273;
    theta2 =T2 +273;%potential temperature(K)
    thetamean =mean([theta1 ,theta2 ]);
    
    p =amb_press_Avg(i) *10;%air pressure(kPa->hPa)
    e1 =e_20m_Avg(i) *10;
    e2 =e_60m_Avg (i)*10;%vapor preesure(kPa->hPa)
    q1 =0.622.*(e1 /(p -0.378*e1 ));
    q2 =0.622*(e2 /(p -0.378*e2 ));%Specific humidity(kg/kg)
    qmean =mean([q1 ,q2 ]);
    
    Cp=1004.67*(1+0.84*qmean );%constant pressure specific heat capacity
    lambda =(2.501-0.00237*Tmean )*10^6;%latent heat of vaporization
    rho =1.293*(p /1000).*(273.15/(Tmean +273.15));%air density
    %% iteration
    
    L0 =inf;
    ustar =kappa*(u2 -u1 )/flux_profile_u(z1,z2,L0);
    thetastar =kappa*(theta2 -theta1 )/flux_profile_theta(z1,z2,L0);
    qstar =kappa*(q2 -q1 )/flux_profile_q(z1,z2,L0);
    while 1
        
        L =M_O_Length(ustar ,thetastar ,thetamean );
        
        
        if abs((L -L0 )/L0)<0.005
            break;
        else
        	ustar =kappa*(u2 -u1 )/flux_profile_u(z1,z2,L);
        	thetastar =kappa*(theta2 -theta1 )/flux_profile_theta(z1,z2,L0);
        	qstar =kappa*(q2 -q1 )/flux_profile_q(z1,z2,L);
        
            L0 =L ;
        end
    end
    
    
    %% output
%     tau(i)=rho *ustar *ustar ;%momentum flux
    Heat(i) =-rho *Cp*ustar *thetastar ;%sensible heat flux
    LHeat (i)=-rho *lambda *ustar *qstar ;%latent heat flux
    %% loop
     i=i+1;
end

Cal_PsiM.m

function PsiM = Cal_PsiM(z,L)
%M_O_LENGTH author: Iris Lee
%   Calculate PsiM
x=(1-16*z/L)^0.25;
if (z/L)<=0
    PsiM=-2*log((1+x)/2)-log((1+x^2)/2)+2*atan(x)-pi/2;
else
    PsiM=5*z/L;
end

Cal_PsiH_q.m

function PsiHq = Cal_PsiH_q(z,L)
%PSIHQ author: Iris Lee
%   Calculate PsiH and Psiq
y=(1-16*z/L)^0.5;
if (z/L)<=0
    PsiHq=-2*log((1+y)/2);
else
    PsiHq=5*z/L;
end

flux_profile_u.m

function f_p_u= flux_profile_u(z1,z2,L)
%FLUX_GRADIENT_U author:Iris Lee
%   通量-廓线公式动量公式的右项
f_p_u=log(z2/z1)+Cal_PsiM(z2,L)-Cal_PsiM(z1,L);
end

flux_profile_theta.m

function f_p_theta= flux_profile_theta(z1,z2,L)
%FLUX_GRADIENT_U author: Iris Lee
%   通量-廓线公式位温公式的右项
f_p_theta=log(z2/z1)+Cal_PsiH_q(z2,L)-Cal_PsiH_q(z1,L);
end

flux_profile_q.m

function f_p_q= flux_profile_q(z1,z2,L)
%FLUX_GRADIENT_U author: Iris Lee
%   通量-廓线公式比湿公式的右项
f_p_q=log(z2/z1)+Cal_PsiH_q(z2,L)-Cal_PsiH_q(z1,L);
end

M_O_Length

function L = M_O_Length(ustar,thetastar,thetamean)
%M_O_LENGTH author: Iris Lee
%Calculate Monin–Obukhov length
kappa=0.4;
g=9.8;
L=thetamean.*ustar.*ustar./(kappa*g*thetastar);
end

output处隐去了动量通量tau，是因为计算的出来的tau精度太差

M是输入的向量行数

输出的感热和潜热通量是1xM的矩阵

3. 要点

注意原始数据的单位，比如水汽压和大气压强都是kpa而不是hpa,
相对湿度RH自带百分号%等

因为存在迭代，检验L，因此只能使用按数字算外套循环的办法，需自己设定循环次数M，
否则针对时间序列的数据，可以在MATLAB中采用按列向量计算的方法，但是需注意矩阵按元素的乘除法需要把*变成.*,/变成./
还需注意矩阵与数字之间计算时，如果出现数字分别除矩阵中每个元素，需要变换形式
如

计算括号内项时，Ta是带有时间序列的列向量，则

(273.15/(Ta +273.15))

这种写法得到的ans是一个行向量，即ans*（Ta+273.15）的矩阵相乘为273.15
因此，需要得到数字与矩阵的每一个元素分别相除，然后按元素输出列向量的话，需要

（Ta+273.15）.\273.15

此处参考matlab中的乘除法