新安江模型程序核心源代碼

1、% XAJ是新安江的運行程序Imp1=parameter.IMP ; Kc= parameter.Kc ;WMU=parameter.WMU ; WML=parameter.WML ;%新安江模型程序核心源代碼function Qr=XAJ_JUN(DAREA,DT,EM,WwU,WwL,WwD,PS0, FR0, Qrs0, Qrss0, Qrg0,parameter,Qm),用于單純形和遺傳算法調用 ,也用于新安江模型的預報%流域不透水面積比： 次洪%流域蒸散發(fā)折算系數：多年總徑流量決定 %流域上層蓄水容量%流域中層蓄水容量WMD = parameter.WMD ;%流域下層蓄水容量B =

2、 parameter.B ; C = parameter.C ;Ex = parameter.Ex;%流域蓄水容量分布曲線指數%流域深層蒸發(fā)系數%流域自由水分布曲線指數SM = parameter.SM ;Ki = parameter.Ki Kg = parameter.Kg Cs = parameter.CsCi = parameter.Ci Cg = parameter.Cg%流域自由水平均蓄水容量%自由水箱壤中流出流系數%自由水箱地下水出流系數%地面水線性水庫匯流系數%壤中流線性水庫匯流系數%地下水線性水庫匯流系數Ke = parameter.KeXe = parameter.Xe ;L

3、 = parameter.L ;%次洪決定： WM， B，Imp%馬斯京根法河段傳播時間%馬斯京根法流量比重系數%滯后演算法參數%WwU(0)=WwU;WwL(0)=WwL;WwD(0)=WwD;%由于日模型與次洪模型的計算時段長不同，參數值不能全部通用，但K、 WM 、 WUM 、WLM、B、IMP、EX C與時段長無關，可以直接引用，%Kc SM KG KSS CS Cl、CG與時段長相關，不能直接引用，需要另外率定 %junjunzhu-XAJ-MODELU=DAREA/(DT*3.6); %單位轉換D=24/DT;KSSD = (1 - (1 - (Ki + Kg) A (1 / D)

4、 / (1 + Kg / Ki); % KSSD ki 出流系數 KGD消退系數 KGD = KSSD * Kg / Ki;%A_WM=A_WUM+A_WLM+A_WDM;%WMM=(1+B).*WM/(1-IMP);Epp=Kc*EM;% PE=P-K.*EM;for T=1:size(EM,1)% T以時段為單位計算%三層蒸散發(fā)計算if (WwU + P(T) = Epp(T)EU(T)Epp(T);% 上 層 蒸 %Epp 為 EMEL(T) = 0;ED(T) = 0;%中層%下層elseEU(T) = WwU + P(T) ; %Ww(1) + P 為 EUEL(T) = (Epp(

5、T) - EU(T) * WwL / WML;%要求計算的下層蒸發(fā)量與剩余蒸散發(fā)能力之比不小于深層蒸散發(fā)系數 cED(T) = 0;if WwL = C * (Epp(T) - EU(T)%要求計算的下層蒸發(fā)量與剩余蒸散發(fā)能力之比小于深層蒸散發(fā)系數 cEL(T) = C * (Epp(T) - EU(T) ;ED(T) = 0;elseEL(T) = WwL;ED(T) = C * (Epp(T) - EU(T) - EL(T) ; endendendPE(T) = P(T) - EU(T) - EL(T) - ED(T); %= %產流計算部分 %=Wm0 = WMU + WML + WMD

6、;% 平均蓄水容量W0 = WwU + WwL + WwD;%初始含水量R = 0;Rimp = 0;Wmm = (1 + B) * Wm0 / (1 - Imp1) ;% Imp1 不透水面積比， Wmm 為蓄水容量 極值if PE(T) 0 曲線的指數%Then GoTo 1000降雨小于蒸發(fā)，B 為蓄水容量if abs(Wm0 - W0) = 0.0001% Wmm為蓄水容量極值A =Wmm;elseA= Wmm * (1 - (1 - W0/ WmO) A (1 / (1 + B);%A 為與 W0 對應的在蓄水容量曲線的縱坐標% 部分產流endif (PE(T) + A) WmmR

7、= PE(T) - WmO + W0 + WmO * (1 - (PE(T) + A) / Wmm)人(1 + B); elseR = PE(T) - (WmO - WO) ;% 全部產流endif abs(R - PE(T) WMU% 由 Ww(1) = Ww(1) + P - RE(1)：E(1)兩斷 Epp 和 Ww1WwL = WwL + WwU - WMU;% 由 Ww(2) = Ww(2) + Ww(1) -WM(1) 檢查是否超標WwU = WMU;% 糾正if WwL WMLWwD = WwD + WwL - WML;WwL = WML;endendif WwU 0WwU =

8、0; end%= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = %匯流計算部分%= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =%水源劃分X = FRO ;% FRO產流面積if PE(T) = 0%認為單是地下自由水在產流面積上的深為Rs(T) = 0;Rss(T) = S0 * KSSD * FR0 ; 水出流)出流系數Rg(T) = S0 * KGD * FR0;S0 = S0 - (R

9、ss(T) + Rg(T) / FR0 ; 積上的平均蓄水深elseFR0 = R / PE(T);上的凈雨(可以理解為產流深)即得產流面積S0 = X * S0 / FR0 ;SS = S0;Q = R / FR0 ;NN = fix(Q / 5) + 1 ;整數 NN 為了消除前向差分誤差Q = Q / NN;%KSSD ki, KGD( KG地下% s 表示自由水在產流面% 用流量除以單位面積% 產流面積變化的影響% 為產流面積上的平均值% 每次入流按 5 毫米分成并取% 一天分為CSng(NN個時段Kssdd = (1 - (1 - (KGD + KSSD) A (1 / NN) /

10、(1 + KGD / KSSD); Kgdd = Kssdd * KGD / KSSD;Rs(T) = 0;Rss(T) = 0;Rg(T) = 0;Smm = (1 + Ex) * SM ;蓄水容量if Ex Smf%s 表示自由水在產流面積上的平均蓄水深S0 = Smf; endAU = Smmf * (1 - (1 - S0 / Smf) A (1 / (1 + Ex);% 當徑流與此時刻的平均蓄if Q + AU = Smmf最大平均蓄水深全面產壤中流% 當徑流與此時刻的平均蓄水深之和大于Rsd中d為分段的地面流Rsd 中 d 為分段的壤中流Rsd中d為分段的地下徑流 s 表示自由水在

11、產流面積上Rsd(T) = (Q + S0 - Smf) * FR0 ;%Rssd(T) = Smf * Kssdd * FR0 ;%Rgd(T) = Smf * Kgdd * FR0 ;%S0 = Smf - (Rssd(T) + Rgd(T) / FR0; %的平均蓄水深% 當徑流與此時刻的平均蓄水深之和大else if Q + AU Smmf 于最大平均蓄水深部分產壤中流Rsd(T) = (Q - Smf + SO + Smf * (1 - (Q + AU) / Smmf) A (1 + Ex) * FRO; Rssd(T) = (S0+ Q - Rsd(T) / FR0) * Kssd

12、d * FR0;Rgd(T) = (SO + Q - Rsd(T) / FRO) * Kgdd * FRO;SO = SO + Q - (Rsd(T) + Rssd(T) + Rgd(T) / FRO;endendendRs(T) = Rs(T) + Rsd(T) ;Rss(T) = Rss(T) + Rssd(T) ;Rg(T) = Rg(T) + Rgd(T); clear Rsd Rssd RgdendendOUT=Rs;Rss;Rg;%Rs=OUT(:,1); Rss=OUT(:,2);Rg=OUT(:,3);Rs(T) = Rs(T) * (1 - Imp1) ; Rss(T) =

13、Rss(T) * (1 - Imp1);Rg(T) = Rg(T) * (1 - Imp1);%累計三流% 累計% 扣除不透水面積%Qrs = (Rs + Rimp) * U%Qrss = Rss * U * (1 - Ci) + QrssO * Ci%Qrg = Rg * U * (1 - Cg) + QrgO * Cg%= = = = = = = = = = XX% 坡面匯流匯流% I #=CSCICgQrs(T) = (Rs(T) + Rimp) * U * (1 - Cs) + Qrs0 * Cs; % 地面水線性水庫匯流系數 Qrss(T) = Rss(T) * U * (1 - Ci) + Qrss0 * Ci ;% 壤中流線性水庫匯流系數Qrg(T) = Rg(T) * U * (1 - Cg) + Qrg0 * Cg ;% 地下水線性水庫匯流系數Qtr(T) = Qrs(T) + Qrss(T) + Qrg(T);