% evap 三层蒸发模型
% 输入：降雨量p，水面蒸发e0，初始土壤含水量wu,wl,wd
% 参数：蒸散发折算系数kc,下层张力水容量lm,深层蒸散发折算系数c
% 输出：各层蒸发eu,el,ed，总蒸发e，净雨量pe
function [eu,el,ed,e,pe]=evap(p,e0,wu,wl,wd,kc,lm,c)
ep=kc*e0;
if wu+p>ep
    eu=ep;el=0;ed=0;
elseif wl>c*lm
    eu=wu+p;el=(ep-eu)*wl/lm;ed=0;
    % 只有当(ep-eu)>lm时会出现el>wl的情况
elseif wl>c*(ep-wu-p)  % 这里不能是(ep-eu)
    eu=wu+p;el=c*(ep-eu);ed=0;    
else
    eu=wu+p;el=wl;ed=c*(ep-eu)-el;
    if ed>wd
        ed=wd;
    end
end
e=eu+el+ed;
pe=p-e;

二、产流量计算

采用蓄满产流模型计算产流量，需注意的是：

（1）只有当PE>0时才会有产流量，否则产流量R为0。

（2）有人会如流程图中在不透水面积比IM上计算产流量RB再加到R中，不需要。

% runoff 蓄满产流模型
% 输入：净雨量pe，初始土壤含水量w
% 参数：流域平均张力水容量wm，张力水蓄水容量曲线方次b，不透水面积比例im
% 输出：产流量r
function r=runoff(pe,w,wm,b,im)
if pe>0  % 只有pe>0时才产流
    mm=wm*(1+b)/(1-im);
    a=mm*(1-(1-w/wm)^(1/(1+b)));
    if pe+a<mm
        r=pe-wm+w+wm*(1-(pe+a)/mm)^(1+b);
    else
        r=pe-wm+w;
    end
else
    r=0;
end

三、土壤含水量更新

降雨量、蒸发量、产流量以及土壤含水量的变化量之间存在水量平衡关系。利用前面计算得到的蒸发量及产流量，进行各层土壤含水量的更新。多余的水量先补充上层，再补充下层及深层。

这里单独用一个模块来计算土壤含水量的变化与更新，使程序更加简明易读。

% update_w 土壤水更新
% 输入：初始土壤含水量，各层蒸发，降雨量，产流量
% 参数：各层含水能力um,lm,dm
% 输出：各层土壤含水量和总量
function [wu,wl,wd,w]=update_w(wu0,wl0,wd0,eu,el,ed,p,r,um,lm,dm)
wu=wu0+p-eu-r;
wl=wl0-el;
wd=wd0-ed;
% 检查土壤水超量
if wu>um
    wl=wl+(wu-um);wu=um;
    if wl>lm
        wd=wd+(wl-lm);wl=lm;
        if wd>dm
            wd=dm;
        end
    end
end
w=wu+wl+wd;

四、分水源计算

产流量R进入自由水蓄水库，根据出流情况划分为地面径流RS、壤中流RI和地下径流RG。三分水源计算是新安江模型中较难理解和编程易出错的部分，我总结了有以下要点：

（1）底宽的变化

自由水蓄水库的底宽为产流面积比FR，并且FR是随时段而变化的。可以理解为S*FR为实际蓄量，或者说RS、RI、RG都是发生在产流面积比FR之上。

（2）差分误差的处理

为了减小因向前差分所造成的误差，每个计算时段的入流量R，按5mm为一段划分为N段进入水库进行水量平衡计算，计算时段DT也分为了N段。

（3）参数的时段转换

参数KI，KG，包括后面的CI、CG、CR，都是以日（24h）为时段长定义的，它们需要根据实际计算时段长而改变。

% separate_r 三分水源
% 输入：净雨pe，产流r，时段初自由水蓄量s，产流面积比fr
% 参数：不透水面积比im，自由水蓄水容量sm，自由水蓄水容量曲线方次ex，
%       壤中流出流系数ki，地下水出流系数kg
% 输出：三水源rs,ri,rg,时段末s,fr
function [rs,ri,rg,fr,s]=separate_r(pe,r,fr,s,im,sm,ex,ki,kg)
ms=sm*(1+ex);
if pe>0
    x=fr;
    fr=(r-im*pe)/pe;  %扣除不透水面积比来计算产流面积比
    s=s*x/fr;      %s有变化底宽fr的影响
    ss=s;          %保存初始s用于后面计算
    q=r/fr;        %产流面积比上的产流深
    % 处理向前差分误差
    n=fix(q/5)+1;  %按5mm一段划分n段
    q=q/n;
    kid=(1-(1-(ki+kg))^(1/n))/(1+kg/ki);  %出流系数按时段长改变
    kgd=kid*kg/ki;
    rs=0;ri=0;rg=0;
    for i=1:n
        if s>sm;s=sm;end;
        au=ms*(1-(1-s/sm)^(1/(1+ex)));
        if q+au<ms
            rs=fr*(q+s-sm+sm*(1-(q+au)/ms)^(1+ex))+rs;
        else
            rs=fr*(q+s-sm)+rs;
        end
        s=s+q*i-rs/fr;
        ri=kid*s*fr+ri;
        rg=kgd*s*fr+rg;
        s=ss+q*i-(rs+ri+rg)/fr;
    end   
else
    rs=0;
    ri=s*ki*fr;
    rg=s*kg*fr;
    s=s*(1-ki-kg);    
end

说明：

（1）在计算产流面积比FR时，需要考虑不透水面积比IM的影响，所以仍需要引入此参数

（2）输入和输出的S采用同一个变量，FR也是如此，在差分误差的处理计算中，S、RS、RI、RG是不断更新变化的

（3）这里的KI、KG是已经根据计算时段DT转换后的值，在后面调用此模块的程序中会有所体现。在这个程序段中，可以看到KI、KG仍需要根据DT分段后的新时段进行再一次转换。

五、汇流计算

（1）坡地汇流

采用线性水库法，将RS、RI、RG转变为QS、QI、QG。相应的消退系数参数分别为CS、CI、CG，一般CS取0，所以略去参数CS。

（2）单元面积河网汇流

单元面积河网总入流QT=QS+QI+QG，采用滞后演算法计算单元面积河网汇流QTR。

（3）河道汇流

采用马斯京根分段连续演算法，计算QTR由单元面积出口至流域总出口的河道汇流过程。

% conflu 三水源汇流计算
% 输入：三水源产流rs,ri,rg，流域面积F
% 参数：壤中地下消退系数ci,cg,河网水流消退系数cr,滞时L,马斯京根参数k,x,n
% 输出：流量q
function [qt,q]=conflu(rs,ri,rg,F,cs,ci,cg,cr,L,k,x,n,qt0)
dt=k;
U=F/(3.6*dt);  %单位换算系数
%坡地汇流
qs(1)=rs(1)*(1-cs)*U;  %初始流量计算
qi(1)=ri(1)*(1-ci)*U;
qg(1)=rg(1)*(1-cg)*U;
qt(1)=qs(1)+qi(1)+qg(1);
for i=2:length(rs);
    qs(i)=qs(i-1)*cs+rs(i)*(1-cs)*U;
    qi(i)=qi(i-1)*ci+ri(i)*(1-ci)*U;  %线性水库法
    qg(i)=qg(i-1)*cg+rg(i)*(1-cg)*U;
    qt(i)=qs(i)+qi(i)+qg(i);
end
%单元面积河网汇流
exn=10;   %扩展
qtt=[ones(1,exn)*qt0,qt];
qtr=zeros(size(qtt));
for i=(exn+1):length(qtr)    
    qtr(i)=cr*qtr(i-1)+(1-cr)*qtt(i-L); %滞后演算法
end
qtr=qtr(exn+1:end);

if n>0
    %单元面积以下河道汇流
    c0=(0.5*dt-k*x)/(0.5*dt+k-k*x);
    c1=(0.5*dt+k*x)/(0.5*dt+k-k*x);
    c2=(k-0.5*dt-k*x)/(0.5*dt+k-k*x);
    qq=zeros(length(qtr)+1,n+1);
    qq(2:end,1)=qtr';
    %用分段马斯京根法
    for j=2:size(qq,2)
        for i=2:size(qq,1)
            qq(i,j)=c0*qq(i,j-1)+c1*qq(i-1,j-1)+c2*qq(i-1,j); 
        end
    end
    q=qq(2:end,n+1);
else
    q=qtr';
end

qt=qtt(exn:end-1);

六、单元的综合计算

新安江模型根据雨量站的分布，采用泰森多边形法，将整个流域分为若干单元。每个单元将上面的各个模块串联起来计算，便可以得到该单元对流域出口的流量过程。

这里输入、输出的降雨、蒸发、产流等，都是时间序列的数组。前面模块函数中定义的是单个值，所以在各计算时段调用这些模块得到时间序列值。该程序段的流程包括：参数的读取，参数的时段转换，初值的设置，产流计算，汇流计算。

% 新安江模型主模块
% 输入：降雨，水面蒸发，参数，初值，面积权重，至出口断面河段数
% 输出：流域蒸发，产流，土壤含水量，流量
function [e,r,w,wu,wl,wd,fr,s,qt,q] = xaj_mdl(p,e0,Para,S0,FW,NE)
    
    %disp('Input parameters ...');
    KC = Para(1) ; UM = Para(2) ; LM = Para(3) ; C  = Para(4);
    WM = Para(5) ; B  = Para(6) ; IM = Para(7) ;
    SM = Para(8) ; EX = Para(9) ; KI = Para(10); KG = Para(11);
    CI = Para(12); CG = Para(13); CR = Para(14); LT = Para(15);
    XE = Para(16); KE = Para(17); F  = Para(18)*FW;

    DM=WM-UM-LM; CS=0; DT=KE;     % time step (h)
    KIt=(1-(1-KI-KG)^(DT/24))/(1+KG/KI);
    KGt=KIt*KG/KI;
    CIt=CI^(DT/24);
    CGt=CG^(DT/24);
    CRt=CR^(DT/24);
    
    %disp('Calculate initial state ...');
    % initial value
    [wu(1),wl(1),wd(1),w(1),fr(1),s(1),qt0]=ini_state(S0,WM,UM,LM,DM,SM);
        
    %disp('Calculate runoff product ...');
    for i=1:length(p)
        [eu,el,ed,e(i),pe]=evap(p(i),e0(i),wu(i),wl(i),wd(i),KC,LM,C);    
        r(i)=runoff(pe,w(i),WM,B,IM);
        [wu(i+1),wl(i+1),wd(i+1),w(i+1)]=update_w(wu(i),wl(i),wd(i),eu,el,ed,p(i),r(i),UM,LM,DM);    
        [rs(i),ri(i),rg(i),fr(i+1),s(i+1)]=separate_r(pe,r(i),fr(i),s(i),IM,SM,EX,KIt,KGt);    
    end

    %disp('Calculate flow confluence ...');
    [qt,q]=conflu(rs,ri,rg,F,CS,CIt,CGt,CRt,LT,KE,XE,NE,qt0);

end

七、总流域的综合计算

前面得到的是计算单元的结果，将各个计算单元结果综合得到总流域的结果。

function [re_s,re_dt,re_val] = xaj_cal(dh,ST,ET,Pdat,Edat,Qdat,Para,Zdat,Sdat,Sdat0)

DT=Para(17); F=Para(18); ZN=size(Zdat,1);  % 时段，面积，分块数
%disp('Extract data ...');
[T,TL,p,e0,qob,Tq,TLq,Qobq]=extra_dat(dh,ST,ET,DT,Pdat,Edat,Qdat,ZN);
N=length(e0);

% 分区域计算结果------------------------------------------------------------
e_z=[]; r_z=[]; w_z=[]; q_z=[]; re_s=[];
for zn=1:ZN
    if isempty(Sdat)
        disp(['zone ',num2str(zn),' has no initial state value, use default']);
        S0=Sdat0;
    else        
        idx=find(Sdat(:,1)==str2num(datestr(ST,'yyyymmdd'))&Sdat(:,2)==zn);
        if isempty(idx)
            disp(['zone ',num2str(zn),' have no state value, use default']);
            S0=Sdat0;
        else
            S0=Sdat(idx,3:end);
        end   
    end
    
    [e,r,w,wu,wl,~,fr,s,qt,q] = xaj_mdl(p(:,zn),e0,Para,S0,Zdat(zn,1),Zdat(zn,2));
    
    e_z=[e_z,e'];r_z=[r_z,r'];w_z=[w_z,w(2:end)'];q_z=[q_z,q]; %保存各分区结果
    
    if (dh=='d')
        re_s=[re_s;TL',zn*ones(size(TL')),w(1:end-1)',wu(1:end-1)', ...
            wl(1:end-1)',fr(1:end-1)',s(1:end-1)',qt'];  %日模型保存初始状态值
    end
    
end

% 计算流域平均结果----------------------------------------------------------
zn_w=Zdat(:,1);  % 面积权重矩阵（zn*1）
p_m=p*zn_w; e_m=e_z*zn_w; r_m=r_z*zn_w; w_m=w_z*zn_w; qcal=sum(q_z,2);

% 计算特征值并显示----------------------------------------------------------
[r_cal,r_obs,r_er,qm_cal,qm_obs,qm_er,Dc]=obj_fun(TL,qcal,TLq,Qobq,DT,F);
p_t=sum(p_m); e_t=sum(e_m); r_t=3.6*sum(qcal)*DT/F; qm_t=max(qcal); %计算期总量

disp(strcat('Total-----------------------Rainfall (mm): ',num2str(p_t)));
disp(strcat('                         Evaporation (mm): ',num2str(e_t)));
disp(strcat('                              Runoff (mm): ',num2str(r_t)));
disp(strcat('                              Peak (m3/s): ',num2str(qm_t)));
disp(strcat('Measure Time---------Observed Runoff (mm): ',num2str(r_obs)));
disp(strcat('                   Calculated Runoff (mm): ',num2str(r_cal)));
disp(strcat('                        Runoff Difference: ',num2str(r_er),'%'));
disp(strcat('                     Observed Peak (m3/s): ',num2str(qm_obs)));
disp(strcat('                   Calculated Peak (m3/s): ',num2str(qm_cal)));
disp(strcat('                          Peak Difference: ',num2str(qm_er),'%'));
disp(strcat('--------------------------N-S Coefficient: ',num2str(Dc)));
disp(' ');

re_val=[TL(1),TL(end),p_t,e_t,r_t,qm_t,r_obs,r_cal,r_er,qm_obs,qm_cal,qm_er,Dc];
re_dt=[TL',p_m,e_m,r_m,w_m,qcal,qob'];

% 画图---------------------------------------------------------------------
figure;
if isempty(Tq)
    subplot(3,1,1); bar(1:N,p_m); ylabel('P（mm）');
    set(gca,'XTick',[]); set(gca,'YDir','reverse'); %对Y方向反转
    title([num2str(TL(1)),' - ',num2str(TL(end)),' hydrograph']);
    subplot(3,1,2:3); plot(T,qcal,'b-'); ylabel('Q（m^3/s）');
    legend('Calculated');    
else    
    subplot(3,1,1); bar(1:N,p_m); ylabel('P（mm）');
    set(gca,'XTick',[]); set(gca,'YDir','reverse'); %对Y方向反转
    title([num2str(TL(1)),' - ',num2str(TL(end)),' hydrograph']);
    subplot(3,1,2:3); plot(Tq,Qobq,'r.',T,qcal,'b-'); ylabel('Q（m^3/s）');
    legend('Observed','Calculated');    
end

end

说明：

（1）这里输入的降雨是二维数组，即各雨量站（计算单元）、各时段的值。

（2）该程序段的流程可概括为：读取参数，提取所选时间的资料，状态初值的处理，新安江模型主模块的运行，流域结果的综合，特征值计算，画图。

（3）该程序段中还用到其他模块程序，后文说明。

八、其他模块

以上所述的是新安江模型的核心原理及模块程序。对于完整的程序系统，还有其他功能需实现，所以还有其他模块程序编制。

上图为各程序文件，简述如下：

（1）输入和输出数据库

输入：input.xlsx，存放日模型和次洪模型的参数、流域分区面积权重以及到流域出口的河段数、洪水场次、降雨、蒸发、实测流量资料

输出：result.xlsx，存放日模型状态结果（作为初值调用），以及日模型库中没有结果时的默认初始值，日模型和次洪模型的时段结果（包括时间、降雨、蒸发、产流、土壤含水量、计算流量、实测流量），日模型和次洪模型的结果特征值（包括计算起始时间、降雨量、蒸发量、径流量、洪峰值、径流误差、洪峰误差、峰现时差、确定性系数）

（2）核心模块

核心模块程序在上文中已讲述，包括evap（蒸散发计算），runoff（产流量计算），update_w（土壤含水量更新），separate_r（分水源），conflu（汇流计算），xaj_mdl（单元的综合），xaj_cal（流域的综合）

（3）其他模块

extra_dat：根据设置的起止时间从数据库提取降雨、蒸发及实测流量资料，并进行缺测值的识别与处理

ini_state：初值的读取与判定

obj_fun：计算相对误差、确定性系数等特征值

XAJ_DP：主运行程序，包括日模型或次洪模型的选择，数据的读取，场次洪水的选择，新安江模型的运行，数据的保存