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

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

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

TRG0=0.4.*QS（1）;
%=================================
%產流參數
K=runoffPara（1）;%流域蒸發折算系數，是流域蒸散發能力與蒸發皿蒸發量之比
FC=runoffPara（2）;%?穩滲率
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倒數第一行；
B=runoffPara（7）;%?流域蓄水容量分布曲線指數；
C=runoffPara（8）;?%?流域蒸發擴散系數?OR?深層蒸散發系數???
%------------------------------
%匯流參數
KKG=xajPara（1）;?%?地下水消退系數；
CS=xajPara（2）;%?CS－河網蓄水消退系數；
L=xajPara（3）;%?L－滯后時間h；
L=round（L）;
KKS=xajPara（4）;%?地面徑流消退系數
%=================================
WMM=（1+B）.*WM/（1-IMP）;??%?是流域內點最大蓄水容量?；
M=size（P1）;???%?M?-?降雨數據的數量
PE=P-K.*EM;??%?PE-扣除雨期蒸發后的降雨量mm;[3664]p145
for?T=1:M???????????????%%??T以時段為單位計算
???%===========================================
????%以下為產流計算;輸入數據：PE，W；?輸入參數：WM，WMM，B;?輸出數據：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
????%===========================================
????%?以下為蒸發計算;;輸入數據：PE，W，R；?輸入參數：WM，WMM，B;?輸出數據：EW
????if?PE（T）<0???%?當降雨量小于蒸發量時時無產流R
????????if?WU+PE（T）>0
????????????EU=K*EM（T）;???%?EU等于蒸散發能力，即在充分供水的條件下，流域蒸散發達到最大;[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-深層蒸散發擴散系數;此式的物理意義是如果下層