新安江模型程序C++代码

1、新安江模型程序C+代码以下是类的声明：class XinanjiangModelprivate:/ FORCINGdouble *m_pP; / 降水数据double *m_pEm; / 水面蒸发数据/long m_nSteps; /模型要运彳T的步长（一共m_nStep涉）long steps;/ OUTPUTdouble *m_pR; / 流域内每一步长的产流量（径流深度）double *m_pRs; / 每一步长的地表径流深（毫米 ）double *m_pRi; / 每一步长的壤中流深（毫米）double *m_pRg; / 每一步长的地下径流深（毫米 ）double *m_pE; /

2、每一步长的蒸发 （毫米 ）double *m_pQrs; / 流域出口地表径流量double *m_pQri; / 流域出口壤中流径流流量double *m_pQrg; / 流域出口地下径流量double *m_pQ; / 流域出口的总流量double m_U;/for24h. U=A（kmA2）/3.6/delta_t/ SOILdouble *m_pW;/流域内土壤湿度double *m_pWu; / 流域内上层土壤湿度double *m_pWl; / 流域内下层土壤适度double *m_pWd;/流域内深层土壤湿度double m_Wum;/流域内上层土壤蓄水容量double m_Wl

3、m; / 流域内下层土壤蓄水容量double *m_pEu;double *m_pEl;double *m_pEd;/runoffdouble *RF;/ PARAMETERdouble m_Kc;double m_IM;double m_B;double m_WM;double m_C;double m_SM;double m_EX;double m_KG;double m_KI;double m_CG;double m_CI;double *m_UH;double m_WMM;double m_Area;int m_DeltaT;int m_PD;public:double m_Wdm;

4、/ 流域内深层土壤蓄水容量， WDM=WM-WUM-WLM/ EVAPORATION / 上层土壤蒸发量（毫米）/ 下层土壤蒸发量（毫米）/ 深层土壤蒸发量（毫米）/ 流域蒸散发能力与实测蒸散发值的比/ 不透水面积占全流域面积之比/蓄水容量曲线的方次，小流域（几平方公里）B0.1左右/ 中等面积（平方公里以内）.20.3，较大面积.30.4/ 流域平均蓄水容量 （毫米） （WM=WUM+WLM+WDM）/ 流域内深层土壤蒸发系数，江南湿润地区：0.15-0.2，/华北半湿润地区：.09-0.12/自由水蓄水容量自由水蓄水容量面积分布曲线指数/地下水日出流系数/壤中流日出流系数/地下水消退系数/

5、壤中流消退系数/ 单元流域上地面径流的单位线/ 流域内最大蓄水容量/ 流域面积/ 每一步长的小时数/ 给定数据，用以判断是否时行河道汇流计算XinanjiangModel(void);XinanjiangModel(void);/ 初始化模型void InitModel( long nSteps,double Area,int DeltaT, int PD, char *ForcingFile);/ 设置模型参数void SetParametersd(ouble *Params);/ 运行新安江模型void RunModel(void);/ 保存模拟结果到文件void SaveResults(

6、char *FileName);/ 记录出流数据，用以作图分析void Runoff(char *runoff);private:/ 进行汇流计算，将径流深度转换为流域出口的流量void Routing(void);以下是类的定义#include stdafx.h#include xinanjiangmodel.h#include #include #include using namespacestd;#include math.h#include stdio.h#include conio.hXinanjiangModel:XinanjiangModel( void)this-m_pP =

7、 NULL;this-m_pEm = NULL;this-m_pE = NULL;this-m_pEl = NULL;this-m_pEu = NULL;this-m_pW = NULL;this-m_pWd = NULL;this-m_pWl = NULL;this-m_pWu = NULL;this-m_pR = NULL;this-m_pRg = NULL;this-m_pRi = NULL;this-m_pRs = NULL;this-m_pQ = NULL;this-m_pQrg = NULL;this-m_pQri = NULL;this-m_pQrs = NULL;Xinanji

8、angModel:XinanjiangModel( void) delete this-m_pP;delete this-m_pEm;delete this-m_pE;delete this-m_pEd;delete this-m_pEl;delete this-m_pEu;delete this-m_pW;delete this-m_pWd;delete this-m_pWl;delete this-m_pWu;delete this-m_pR;delete this-m_pRg;delete this-m_pRi;delete this-m_pRs;delete this-m_pQ;del

9、ete this-m_pQrg;delete this-m_pQrs;delete this-m_pQri;/ 初始化模型void XinanjiangModel:InitModel( long nSteps,double Area, int DeltaT,int PD, char* ForcingFile)FILE * fp;int i;this-m_nSteps = nSteps;this-steps = this-m_nSteps + 18;/ 驱动数据this-m_pP = new doublethis-steps;this-m_pEm = new doublethis-steps;/

10、 模型输出，蒸散发项this-m_pE = new doublethis-steps;this-m_pEd = new doublethis-steps;this-m_pEl = newdoublethis-steps;this-m_pEu = new doublethis-steps;/ 模型输出，出流项，经过汇流的产流this-m_pQrg = new doublethis-steps;this-m_pQrs = new doublethis-steps;this-m_pQri = new doublethis-steps;this-m_pQ = new doublethis-steps;

11、/ 模型输出，产流项this-m_pR = new doublethis-steps;this-m_pRg= newdoublethis-steps;this-m_pRi= newdoublethis-steps; this-m_pRs = new doublethis-steps; / 模型状态量，土壤湿度this-m_pW = new doublethis-steps;this-m_pWd = new doublethis-steps;this-m_pWl = new doublethis-steps; this-m_pWu = new doublethis-steps; /runoff

12、值this-RF = new doublethis-steps;for(i = 0;isteps;i+ )/ 驱动数据this-m_pP i = 0.00;this-m_pEm i = 0.00;/ 模型输出，蒸散发项this-m_pE i = 0.00;this-m_pEd i = 0.00;this-m_pEl i = 0.00;this-m_pEu i = 0.00;/ 模型输出，出流项，经过汇流的产流this-m_pQrgi = 0.00;this-m_pQrsi = 0.00;this-m_pQrii = 0.00;this-m_pQi = 0.00;/ 模型输出，产流项this-m

13、_pR i = 0.00;this-m_pRg i = 0.00;this-m_pRi i = 0.00;this-m_pRs i = 0.00;/ 模型状态量，土壤湿度this-m_pW i = 0.00;this-m_pWdi = 0.00;this-m_pWli = 0.00;this-m_pWui = 0.00;this-m_Area = Area;this-m_DeltaT = DeltaT;this-m_PD = PD;this-m_U = this-m_Area/(3.6 * this-m_DeltaT);/ Forcing文件格式：第一列：降水(单位毫米)空格第二列水面蒸发(毫

14、米) if(fp = fopen(ForcingFile, r) = NULL)printf( Can not open forcing file!n);return; for(i = 0;im_nSteps;i+ ) fscanf(fp,%lf%lf ,&(this-m_pPi),&( this-m_pEmi); fclose(fp);/ 设置模型参数void XinanjiangModel:SetParameters(double* Params)this-m_Kc = Params0; / (1) 流域蒸散发能力与实测水面蒸发之比this-m_IM = Params1; / (2) 流域

15、不透水面积占全流域面积之比this-m_B = Params2; / (3) 蓄水容量曲线的方次this-m_Wum = Params3; / (4) 上层蓄水容量this-m_Wlm = Params4; / (5) 下层蓄水容量this-m_Wdm = Params5; / (6) 深层蓄水容量this-m_C = Params6; / (7) 深层蒸散发系数this-m_SM = Params7; / (8)自由水蓄水容量this-m_EX = Params8; /自由水蓄水容量面积分布曲线指数this-m_KG = Params9;/ （10）地下水日出流系数this-m_KI = P

16、arams10;/ （11）壤中流日出流系数this-m_CG = Params11;/ （12）地下水消退系数this-m_CI = Params12;/ （13）壤中流消退系数this-m_WM = this-m_Wum + this-m_Wlm + this-m_Wdm;this-m_WMM = this-m_WM * （1.0 + this-m_B）/（1.0 - this-m_IM）;/ 运行新安江模型void XinanjiangModel:RunModel（void）long i;/ 模型的状态变量double PE; / 0 时为净雨量; 0 为蒸发不足量（ mm）double

17、 Ep; /m_Kc * m_pEmidouble P;double R; / 产流深度，包括地表径流、壤中流和地下径流（ mm）double RB; / 不透水面上产生的径流深度（ mm）double RG; / 地下径流深度（mm）double RI; / 壤中流深度（mm）double RS; / 地表径流深（mm）double A;/土壤湿度为W时土壤含水量折算成的径流深度（mm）double E = 0.0;/ 蒸散发 （mm）double EU = 0.0;/ 上层土壤蒸散发量（mm）double EL = 0.0;/下层土壤蒸散发量（mm）double ED =0.0;/深层土壤

18、蒸散发量（mm）double S;double FRo;double FR;double MS;double AU;double WU = 5.0;/ 流域内上层土壤湿度double WL = 55.0; / 流域内下层土壤适度double WD = 40.0; / 流域内深层土壤湿度double W = 100.0;double So = 5.0;MS = m_SM * (1 + m_EX);FRo = 1 - pow(1 - So/MS),m_EX);for(i = 0;im_nSteps;i+ )/ 蒸散发计算 /RB = m_pPi * m_IM;/ RB是降在不透水面的降雨量P =

19、m_pPi * (1 - m_IM);Ep = m_Kc * m_pEmi;if (WU + P)= Ep)EU = Ep; EL = 0; ED = 0; elseif(WU + P)= (m_C * m_Wlm) EL = (Ep - EU) * WL/m_Wlm; ED = 0; elseif (m_C * (Ep - EU)= WL&WL(m_C * m_Wlm) EL = m_C * (Ep - EU); ED = 0; elseif (WLm_C * (Ep - EU) EL = WL; ED = m_C * (Ep - EU) - EL; E = EU + EL + ED;PE

20、= P - E;/* 蒸散发计算结束 */子流域产流量计算 /if(PE= 0) R = 0.00; W = W + PE; elseA = m_WMM * (1 - pow( (1.0 - W/m_WM), 1.0/(1 + m_B) ) );/ 土壤湿度折算净雨量+降水后蒸发剩余雨量流域内最大含水容量if(A + PE)m_WMM)/ 流域内的产流深度计算R = PE/* 降水蒸发后的剩余量*/+ W/* 流域内土壤湿度*/+ m_WM * pow(1 - (PE + A)/m_WMM),(1 + m_B)- m_WM /* 减去流域平均蓄水容量( m_WM: 参数) */+ RB; /*

21、 不透水面上产生的径流*/ 土壤湿度折算净雨量+降水后蒸发剩余雨量流域内最大含水容量else/ 流域内的产流深度计算R = PE/* 降水蒸发后的剩余量+ W /* 流域内土壤湿度*/- m_WM /* 减去流域平均蓄水容量*/+ RB; /* 不透水面上产生的径流*/三层蓄水量的计算: WU, WL, WDif(WU + P - EU - R = m_Wum) WU = WU + P - EU - R; WL = WL - EL; WD = WD - ED;elseWU = m_Wum;m_Wlm )if(WL - EL + ( WU + P - EU - R - m_Wum ) = WL

22、= WL - EL + ( WU + P - EU - R - m_Wum );WD = WD - ED;elseWL = m_Wlm;if(WD - ED + WL - EL + ( WU + P - EU - R - m_Wum )- m_Wlm 0)FR = (R - RB) / PE;AU = MS * (1 - pow(1 - So * FRo/FR/m_SM),1/(1 + m_EX);if(PE + AU= MS)RS = FR * ( PE + So * Fro / FR - m_SM );S = So * Fro / FR + ( R -RS ) / FR;RI = m_KI

23、 * S * FR;RG = m_KG * S * FR;RS += RB;R = RS + RI + RG;So = S * ( 1 - m_KI - m_KG );FRo = FR;elseS = So;FR = 1 - pow(1 - S / MS) , m_EX );RI = 0.00;RG = 0.00;So = S * ( 1 - m_KI - m_KG );RS = RB;R = RS + RI + RG;FRo = FR;/三水源划分计算结束*/* 以下部分是状态量：总蒸发量、上、下和深层土壤的蒸发的保存/* 1 */this-m_pEi = E;/ 当前步长的蒸发(模型重要输

24、出)/* 2 */this-m_pEui = EU;/ 当前步长上层土壤蒸发/* 3 */this-m_pEli = EL;/ 当前步长下层土壤蒸发/* 4 */this-m_pEdi = ED;/ 当前步长深层土壤蒸发/* 5 */this-m_pWi = W;/ 当前步长流域平均土壤含水量/* 6 */ this-m_pWui = WU; / 当前步长流域上层土壤含水量/* 7 */ this-m_pWli = WL; / 当前步长流域下层土壤含水量/* 8 */ this-m_pWdi = WD; / 当前步长流域深层土壤含水量/* 9 */ this-m_pRgi = RG; / 当前

25、步长流域地下径流深度/* 10 */ this-m_pRii = RI; / 当前步长流域壤中流深度/* 11 */ this-m_pRsi = RS; / 当前步长流域地表径流径流深度/* 12 */ this-m_pRi = R;/ 当前步长的总产流径流深度this-Routing();/ 保存模拟结果到文件void XinanjiangModel:SaveResults(char* FileName)int i;FILE * fp;if(fp = fopen(FileName, w) = NULL)printf(Can not create output file!n);return;f

26、printf(fp, -n);fprintf(fp, E(mm) EU(mm) EL(mm) ED(mm) W(mm) WU(mm) WL(mm) WD(mm) R(mm) RS(mm) RI(mm) RG(mm) Q(m3/d) QS(m3/d) QI(m3/d) QG(m3/d)n);fprintf(fp, n );for(i = 0;isteps;i+ )fprintf(fp, %9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lf% 9.3lf%9.3lf%9.3lf%9.3lf%9.3lf%9.3lfn ,this-m

27、_pEi, this-m_pEui, this-m_pEli, this-m_pEdi, this-m_pWi, this-m_pWui, this-m_pWli, this-m_pWdi, this-m_pRi, this-m_pRsi, this-m_pRii, this-m_pRgi, this-m_pQi, this-m_pQrsi, this-m_pQrii, this-m_pQrgi);fclose(fp);/ 进行汇流计算，将径流深度转换为流域出口的流量void XinanjiangModel:Routing(void)/ 地面径流汇流计算：单位线法 /int i,j;doubl

28、e B10000 = 0.00;if (this-m_PD = 1)double UH =3.71,12.99,38.96,94.63,131.74,154.00,166.99,176.27,178.12, 172.55,146.58, 90.91,53.80, 31.54,18.55, 9.27, 3.71,0.00;for(i = 0;im_nSteps;i+ )for(j = 0;jm_pRsi * UHj/10.0;elsedouble UH = 7.18,23.38,63.20,143.10,221.75,365.18,447.40,491.29,506.93,504.82,468.

29、46,388.56,309.91,166.49,84.26,40.37,17.56,3.46;for(i = 0;im_nSteps;i+ )for(j = 0;jm_pRsi * UHj/10.0;for(i = 0;isteps;i+ )this-m_pQrsi = Bi;/ 壤中流汇流计算:线性水库for(i = 1;isteps;i+ ) this-m_pQrii = this-m_CI * this-m_pQrii - 1+ (1.0 - this-m_CI) * this-m_pRii * this-m_U; / 地下径流汇流计算：线性水库for(i = 1;isteps;i+ )

30、 this-m_pQrgi = this-m_pQrgi - 1 * this-m_CG+ this-m_pRgi * (1.0 - this-m_CG) * this-m_U; /单元面积总入流计算for(i = 0;isteps;i+ ) this-m_pQi = this-m_pQrsi + this-m_pQrii + this-m_pQrgi; void XinanjiangModel:Runoff( char * runoff)int i;ofstream outfile;outfile.open(runoff);if(outfile.is_open()for(i = 0;isteps;i+ ) outfilesetprecision(3)setiosflags(ios:fixed)m_pQiendl; outfile.close();以下是main()函数语句int _tmain(int argc, _TCHAR * argv) long nSteps = 942;int DeltaT = 24;double Area1 = 1603;XinanjiangModel Model1;Model1.InitModel(nSteps, Area1,DeltaT,1,LFForcingfile.txt );/模型参数 /*Kc,