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

1、%新安江模型程序核心源代码function qr=xaj_jun(darea,dt,em,wwu,wwl,wwd,p,s0, fr0, qrs0, qrss0, qrg0,parameter,qm)% xaj是新安江的运行程序,用于单纯形和遗传算法调用,也用于新安江模型的预报 imp1=parameter.imp ; %流域不透水面积比： 次洪 kc= parameter.kc ; %流域蒸散发折算系数：多年总径流量决定 wmu=parameter.wmu ; %流域上层蓄水容量 wml=parameter.wml ; %流域中层蓄水容量 wmd = parameter.wmd ; %流域下层

2、蓄水容量 b = parameter.b ; %流域蓄水容量分布曲线指数 c = parameter.c ; %流域深层蒸发系数 ex = parameter.ex; %流域自由水分布曲线指数 sm = parameter.sm ; %流域自由水平均蓄水容量 ki = parameter.ki ; %自由水箱壤中流出流系数 kg = parameter.kg ; %自由水箱地下水出流系数 cs = parameter.cs ; %地面水线性水库汇流系数 ci = parameter.ci ; %壤中流线性水库汇流系数 cg = parameter.cg ; %地下水线性水库汇流系数 ke =

3、parameter.ke ; %马斯京根法河段传播时间 xe = parameter.xe ; %马斯京根法流量比重系数 l = 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、ci、cg与时段长相关，不能直接引用，需要另外率定 %junjunzhu-xaj-model u=darea/(dt*3.6); %单位转换d=24/d

4、t;kssd = (1 - (1 - (ki + kg) (1 / d) / (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以时段为单位计算% %三层蒸散发计算 if (wwu + p(t) = epp(t) eu(t) = epp(t); %上层蒸发 %epp为em el(t) = 0; %中层 ed(t) = 0; %下层 else eu(

5、t) = wwu + p(t) ; %ww(1) + p为eu el(t) = (epp(t) - eu(t) * wwl / wml; %要求计算的下层蒸发量与剩余蒸散发能力之比不小于深层蒸散发系数c ed(t) = 0; if wwl = c * (epp(t) - eu(t) %要求计算的下层蒸发量与剩余蒸散发能力之比小于深层蒸散发系数c el(t) = c * (epp(t) - eu(t) ; ed(t) = 0; else el(t) = wwl; ed(t) = c * (epp(t) - eu(t) - el(t) ; end end end pe(t) = p(t) - eu

6、(t) - el(t) - ed(t);%= %产流计算部分%= wm0 = wmu + wml + wmd; % 平均蓄水容量 w0 = wwu + wwl + wwd; %初始含水量 r = 0; rimp = 0; wmm = (1 + b) * wm0 / (1 - imp1) ; % imp1不透水面积比，wmm为蓄水容量 极值 if pe(t) 0 %then goto 1000 降雨小于蒸发，b为蓄水容量曲线的指数 if abs(wm0 - w0) = 0.0001 % wmm为蓄水容量极值 a = wmm; else a = wmm * (1 - (1 - w0 / wm0)

7、(1 / (1 + b); %a为与w0对应的在蓄水容量曲线的纵坐标 end if (pe(t) + a) wmm % 部分产流 r = pe(t) - wm0 + w0 + wm0 * (1 - (pe(t) + a) / wmm) (1 + b); else r = pe(t) - (wm0 - w0) ; % 全部产流 end if abs(r - pe(t) wmu % 由ww(1) = ww(1) + p - re(1)：e(1)两断epp和ww1 wwl = wwl + wwu - wmu; % 由ww(2) = ww(2) + ww(1) - wm(1)检查是否超标 wwu =

8、wmu; % 纠正 if wwl wml wwd = wwd + wwl - wml; wwl = wml; end end if wwu 0 wwu = 0; end%汇流计算部分%水源划分 x = fr0 ; % fr0产流面积 if pe(t) = 0 %认为单是地下自由水在产流面积上的深为 rs(t) = 0; rss(t) = s0 * kssd * fr0 ; %kssd，ki，kgd（kg地下水出流）出流系数 rg(t) = s0 * kgd * fr0; s0 = s0 - (rss(t) + rg(t) / fr0 ; % s表示自由水在产流面积上的平均蓄水深 else fr

9、0 = r / pe(t); % 用流量除以单位面积上的净雨（可以理解为产流深）即得产流面积 s0 = x * s0 / fr0 ; % 产流面积变化的影响 ss = s0; q = r / fr0 ; % 为产流面积上的平均值 nn = fix(q / 5) + 1 ; % 每次入流按5毫米分成并取整数nn为了消除前向差分误差 q = q / nn; % 一天分为csng(nn)个时段 kssdd = (1 - (1 - (kgd + kssd) (1 / nn) / (1 + kgd / kssd); kgdd = kssdd * kgd / kssd; rs(t) = 0; rss(t)

10、 = 0; rg(t) = 0; % sm流域的平均自由水容量 smm = (1 + ex) * sm ; % smm全流域最大的自由水蓄水容量 if ex smf %s 表示自由水在产流面积上的平均蓄水深 s0 = smf; end au = smmf * (1 - (1 - s0 / smf) (1 / (1 + ex); if q + au = smmf % 当径流与此时刻的平均蓄水深之和大于最大平均蓄水深全面产壤中流 rsd(t) = (q + s0 - smf) * fr0 ; % rsd中d为分段的地面流 rssd(t) = smf * kssdd * fr0 ; % rsd中d为

11、分段的壤中流 rgd(t) = smf * kgdd * fr0 ; % rsd中d为分段的地下径流 s0 = smf - (rssd(t) + rgd(t) / fr0; % s表示自由水在产流面积上的平均蓄水深 else if q + au smmf % 当径流与此时刻的平均蓄水深之和大于最大平均蓄水深部分产壤中流 rsd(t) = (q - smf + s0 + smf * (1 - (q + au) / smmf) (1 + ex) * fr0; rssd(t) = (s0+ q - rsd(t) / fr0) * kssdd * fr0; rgd(t) = (s0 + q - rsd

12、(t) / fr0) * kgdd * fr0; s0 = s0 + q - (rsd(t) + rssd(t) + rgd(t) / fr0; end end end rs(t) = rs(t) + rsd(t) ; % 累计三流 rss(t) = rss(t) + rssd(t) ; % 累计 rg(t) = rg(t) + rgd(t); clear rsd rssd rgd end end out=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) + qrss0 * ci%qrg = rg * u * (1 - cg) + qrg0 * cg% 坡面汇流汇流%！￥ qrs(t) = (rs(t) + rimp) * u * (1 - cs) + qrs0 * cs; % 地面水线性水库汇流系数cs qrss(t) = rss(t) * u * (1 - ci) + qrss0 * ci ; % 壤中流线性水库汇流系数ci qrg(t) = rg(t) * u * (1 - cg) + qrg0 * cg ; % 地下水线性水库汇流系数cg qtr(t) = qrs(t) + qrss(t) + qrg(t); qsn(t) = (rs(t) + rim