%?兩水源新安江模型程序核心源代碼

function?[objqqrr]=XAJ_fun（xajPara）
%obj為模擬的各目標函數(shù)值，qq為時段實際流量與模擬流量，rr為時段地面凈雨地下凈雨總凈雨
%?XAJ是新安江的運行程序用于參數(shù)優(yōu)化算法調(diào)用也用于新安江模型的預報
%?xajPara是調(diào)用的優(yōu)化參數(shù);
global?DATA
uh=load（‘.\data\uh.txt‘）;%導入地面徑流匯流單位線，為一列或一行數(shù)據(jù)
runoffPara=load（‘.\data\runoffPara.txt‘）;%讀入優(yōu)化好的產(chǎn)流參數(shù)
startn=100;%統(tǒng)計預報結(jié)果起算時段序號*****************
%?輸入起始值?WWUWLWDQG
WU=10;WL=61;WD=20;??%?WU-流域初始上層張力水蓄量;WL-流域初始下層張力水蓄量;WD-流域初始深層張力水蓄量；
AREA=12200;??%?流域面積，km2;
N=6;???%?一個時段長（h）；
U=AREA/3.6/N;???%?U-單位換算系數(shù)
W=WU+WL+WD;???%?流域初始蓄水量mm;

%輸入雨量P蒸散發(fā)能力EM實測流量QS
TIME=DATA（:1）;
P=DATA（:2）;???%?輸入雨量P
EM=DATA（:3）;??%?流域日蒸散發(fā)能力EM
QS=DATA（:4）;??%?實測流量QS

TRG0=0.4.*QS（1）;
%=================================
%產(chǎn)流參數(shù)
K=runoffPara（1）;%流域蒸發(fā)折算系數(shù)，是流域蒸散發(fā)能力與蒸發(fā)皿蒸發(fā)量之比
FC=runoffPara（2）;%?穩(wěn)滲率
WUM=runoffPara（3）;%?流域平均上層蓄水容量?mm?；
WLM=runoffPara（4）;%?流域平均下層蓄水容量?mm?；
WDM=runoffPara（5）;%?流域平均深層蓄水容量?mm
WM=WUM+WLM+WDM;
IMP=runoffPara（6）;??%?%?不透水面積占全流域面積的百分比，％，一般為0.01-0.02:[3664]p151倒數(shù)第一行；
B=runoffPara（7）;%?流域蓄水容量分布曲線指數(shù)；
C=runoffPara（8）;?%?流域蒸發(fā)擴散系數(shù)?OR?深層蒸散發(fā)系數(shù)???
%------------------------------
%匯流參數(shù)
KKG=xajPara（1）;?%?地下水消退系數(shù)；
CS=xajPara（2）;%?CS－河網(wǎng)蓄水消退系數(shù)；
L=xajPara（3）;%?L－滯后時間h；
L=round（L）;
KKS=xajPara（4）;%?地面徑流消退系數(shù)
%=================================
WMM=（1+B）.*WM/（1-IMP）;??%?是流域內(nèi)點最大蓄水容量?；
M=size（P1）;???%?M?-?降雨數(shù)據(jù)的數(shù)量
PE=P-K.*EM;??%?PE-扣除雨期蒸發(fā)后的降雨量mm;[3664]p145
for?T=1:M???????????????%%??T以時段為單位計算
???%===========================================
????%以下為產(chǎn)流計算;輸入數(shù)據(jù)：PE，W；?輸入?yún)?shù)：WM，WMM，B;?輸出數(shù)據(jù)：R
????if?PE（T）<0
????????R=0;?%?流域總徑流量，mm
????else
????????if??W>=WM
????????????A=WMM;??%?流域起始蓄水量為W0時，所對應的最大點蓄水容量mm或為對應的前期影響雨量[9343]p84;
????????else
????????????A=WMM*（1-（1-W/WM）.^（1/（1+B）））;??%?[3664]p145：式（5-9）
????????end
????????if?A+PE（T）>0
????????????if?A+PE（T）????????????????R=PE（T）-WM+W+WM.*（1-（PE（T）+A）./WMM）.^（1+B）;??%?[3664]p145：式（5-11）
????????????else
????????????????R=PE（T）+W-WM;???????????????????????????????%?[3664]p145：式（5-13）
????????????end
????????else
????????????R=0;
????????end
????end
????%===========================================
????%?以下為蒸發(fā)計算;;輸入數(shù)據(jù)：PE，W，R；?輸入?yún)?shù)：WM，WMM，B;?輸出數(shù)據(jù)：EW
????if?PE（T）<0???%?當降雨量小于蒸發(fā)量時時無產(chǎn)流R
????????if?WU+PE（T）>0
????????????EU=K*EM（T）;???%?EU等于蒸散發(fā)能力，即在充分供水的條件下，流域蒸散發(fā)達到最大;[3664]p143：式（1）
????????????ED=0;
????????????EL=0;
????????????WU=WU+PE（T）;
????????else??????????????%?[3664]p143：式（2）
????????????EU=WU+P（T）;
????????????WU=0;
????????????if?WL>C*WLM???%?[3664]p143：式（3）;C-深層蒸散發(fā)擴散系數(shù);此式的物理意義是如果下層