根据太阳辐射来计算日照时数DLc 和 根据比湿计算相对湿度RHu

% 张旭东更新20240724
% DLc——平均日照时数[hour] 
% Rn——冠层表面净辐射[MJ/m^2/day]
% Rns——净短波辐射[MJ/m^2/day]
% Rnl——净长波辐射[MJ/m^2/day]
% Rs——入射的太阳辐射[MJ/m^2/day]
% Ra24——碧空太阳总辐射[MJ/m^2/day]
% latitude——地理纬度转换为[度] 
% elevation——海拔高度，[m]
% RHmean——平均相对湿度，%
% uh——高度为h（m）的风速m/s
% DLc——日照时数 h
clear
clc
load('D:\matlab_work2022\Climate\Read_nc\Hebei_Point_CliaPara_6.mat')
% 计算日序数；
T_num=datevec(T_day);
DateVector1=[T_num(:,1:3)]; % 某某年某月某日
DateVector0=[T_num(:,1),ones(size(T_num,1),2)]; % 某某年1月1日
jday=datenum(DateVector1)-datenum(DateVector0)+1;
% 计算日序数；

elevation=51;
latitude=38.7025; % 目标站点的纬度

Gsc =0.082; % Gsc——太阳常数，[MJ*m^(-2/)*min^(-1)]
as=0.25; % as——地区修正系数 FAO 式（35）的修正值
bs=0.5; % bs——地区修正系数
%% ====转换风速Conversion of wind speed 
% h——实测风速处距地面的距离[m] “数据集”中一般为10m
% u2=4.87*uh./(log(67.8*h-5.42)); % uh——测量点平均风速[m*s^(-1)] 
% uh=u2.*(log(67.8*h-5.42))/4.87;
u2= WIN(:,27); % u2——2m处风速
%% ====比湿/平均温度/
SHum_p=SHum(:,27);
TEM_avg_p=TEM_avg(:,27);
TEM_max_p=TEM_max(:,27);
TEM_min_p=TEM_min(:,27);
PRS_p=PRS(:,27);
SSRA_p=SSRA(:,27);
PRE_p= PRE(:,27);

%% ====参数Parameters
P=PRS_p/1000; % P——大气压强[Kpa];
Rs = 0.0864.*SSRA_p; % 短波辐射 [MJ/m^2/day]
fai=latitude*pi./180; % fai——地理纬度,[rad]

% Slope of saturation vapour pressure curve FAO 式（13）
Tmean=TEM_avg_p-273.15;
derta=(4098.*(0.6108*exp(17.27*Tmean./(Tmean+237.3))))./(Tmean+237.3).^2; % FAO 式（13）
% derta= FAOTable2_4(Tmean); % 采用FAO附录2中TABLE 2.4. Slope of vapour pressure curve for different temperatures

% Psychrometric constant——garma
% lanmta=2.501-2.361*10^(-3).*Tmean; % lanmta——水的汽化潜热 % 修正
lanmta=2.45; % lanmta——水的汽化潜热 2.45 [MJ kg-1]
Cp=1.013/1000; % Cp——标准大气压下的显热常数[MJ/Kg/度C]
efsong=0.622; % efsong——水蒸气与干空气分子的比率常数
garma=Cp/(efsong.*lanmta).*P; % FAO 式（8） 
%% Radiation
% FAO 式（23）
dr=1+0.033*cos(2*pi*jday./365); % dr——地球和太阳的相对距离
% FAO 式（24）
solar_angle=0.409*sin(2*pi*jday./365-1.39); % solar_angle——太阳角[rad]
% FAO 式（25）
ws=acos(-tan(fai).*tan(solar_angle)); % ws——日落角[rad]
DL=24*ws./pi; % DL——最大可能日照时间[hour]
% Extraterrestrial radiation for daily periods (Ra24) FAO 式(21)
Ra24=24*60/pi*Gsc*dr.*(ws.*sin(fai).*sin(solar_angle)+cos(fai).*cos(solar_angle).*sin(ws));
DLc=(Rs./Ra24-as).*(DL./bs);
DLc(DLc<0)=0;
for i=1:size(DL,1)
if DLc(i)>DL(i)
DLc(i)=DL(i);
end
end
%% 计算相对湿度RHu 
% SHum——比湿
% 参考:https://blog.csdn.net/Soul_taker/article/details/126534580
% 参考:https://earthscience.stackexchange.com/questions/2360/how-do-i-convert-specific-humidity-to-relative-humidity?newreg=db443585846b4ee9b1fc57be895988f6
% TEM_avg——平均温度，K。

RHu = 0.263 .* PRS_p .* SHum_p .* exp((17.67 * (TEM_avg_p- 273.15))./(TEM_avg_p - 29.65)).^(-1);
RHu(RHu>99)=99;
h = 10; % 高度10m
uh = u2.*(log(67.8*h-5.42))/4.87; % uh——10m处风速