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

新安江模型程序C++代码以下是类的声明：

class XinanjiangModel

{

private:

// FORCING

double *m_pP; // 降水数据

double *m_pEm; // 水面蒸发数据

//

long m_nSteps; // 模型要运行的步长(一共m_nSteps步) long steps;

// OUTPUT

double *m_pR; // 流域内每一步长的产流量(径流深度) double *m_pRs; // 每一步长的地表径流深(毫米)

double *m_pRi; // 每一步长的壤中流深（毫米）

double *m_pRg; // 每一步长的地下径流深(毫米)

double *m_pE; // 每一步长的蒸发(毫米)

double *m_pQrs; // 流域出口地表径流量

double *m_pQri;// 流域出口壤中流径流流量

double *m_pQrg;// 流域出口地下径流量

double *m_pQ; // 流域出口的总流量

double m_U; // for 24h. U=A(km^2)/3.6/delta_t

// SOIL

double *m_pW; // 流域内土壤湿度

double *m_pWu; // 流域内上层土壤湿度

double *m_pWl; // 流域内下层土壤适度

double *m_pWd; // 流域内深层土壤湿度

double m_Wum; // 流域内上层土壤蓄水容量

double m_Wlm; // 流域内下层土壤蓄水容量

double m_Wdm; // 流域内深层土壤蓄水容量，WDM=WM-WUM-WLM // EVAPORATION

double *m_pEu; // 上层土壤蒸发量（毫米）

double *m_pEl; // 下层土壤蒸发量（毫米）

double *m_pEd; // 深层土壤蒸发量（毫米）

//runoff

double *RF;

// PARAMETER

double m_Kc; // 流域蒸散发能力与实测蒸散发值的比

double m_IM; // 不透水面积占全流域面积之比

double m_B; // 蓄水容量曲线的方次，小流域（几平方公里）B0.1左右// 中等面积（平方公里以内）.2~0.3，较大面积.3~0.4

double m_WM; // 流域平均蓄水容量（毫米）(WM=WUM+WLM+WDM) double m_C; // 流域内深层土壤蒸发系数，江南湿润地区：0.15-0.2，

//华北半湿润地区：.09-0.12

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:

XinanjiangModel(void);

~XinanjiangModel(void);

// 初始化模型

void InitModel(long nSteps, double Area,int DeltaT, int PD, char *ForcingFile);

// 设置模型参数

void SetParameters(double *Params);

// 运行新安江模型

void RunModel(void);

// 保存模拟结果到文件

void SaveResults(char *FileName);

// 记录出流数据，用以作图分析

void Runoff(char *runoff);

private:

// 进行汇流计算，将径流深度转换为流域出口的流量

void Routing(void);

};

以下是类的定义

#include"stdafx.h"

#include"xinanjiangmodel.h"

#include

#include

#include

usingnamespace std;

#include"math.h"

#include"stdio.h"

#include"conio.h"

XinanjiangModel::XinanjiangModel(void)

{

this->m_pP = NULL;

this->m_pEm = NULL;

this->m_pE = NULL;

this->m_pEd = 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;

}

XinanjiangModel::~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;

delete[] 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 = newdouble[this->steps];

this->m_pEm = newdouble[this->steps];

// 模型输出，蒸散发项

this->m_pE = newdouble[this->steps];

this->m_pEd = newdouble[this->steps];

this->m_pEl = newdouble[this->steps];

this->m_pEu = newdouble[this->steps];

// 模型输出，出流项，经过汇流的产流

this->m_pQrg = newdouble[this->steps];

this->m_pQrs = newdouble[this->steps];

this->m_pQri = newdouble[this->steps];

this->m_pQ = newdouble[this->steps];

// 模型输出，产流项

this->m_pR = newdouble[this->steps];

this->m_pRg= newdouble[this->steps];

this->m_pRi= newdouble[this->steps];

this->m_pRs= newdouble[this->steps];

// 模型状态量，土壤湿度

this->m_pW = newdouble[this->steps];

this->m_pWd = newdouble[this->steps];

this->m_pWl = newdouble[this->steps];

this->m_pWu = newdouble[this->steps];

//runoff值

this->RF = newdouble[this->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_pQrg[i] = 0.00;

this->m_pQrs[i] = 0.00;

this->m_pQri[i] = 0.00;

this->m_pQ[i] = 0.00;

// 模型输出，产流项

this->m_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_pWd[i] = 0.00;

this->m_pWl[i] = 0.00;

this->m_pWu[i] = 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文件格式：第一列：降水（单位毫米）空格第二列水面蒸发（毫米）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_pP[i]),&(this->m_pEm[i])); }

fclose(fp);

}

// 设置模型参数

void XinanjiangModel::SetParameters(double* Params)

{

this->m_Kc = Params[0]; // (1) 流域蒸散发能力与实测水面蒸发之比

this->m_IM = Params[1]; // (2) 流域不透水面积占全流域面积之比

this->m_B = Params[2]; // (3) 蓄水容量曲线的方次

this->m_Wum = Params[3]; // (4) 上层蓄水容量

this->m_Wlm = Params[4]; // (5) 下层蓄水容量

this->m_Wdm = Params[5]; // (6) 深层蓄水容量

this->m_C = Params[6]; // (7) 深层蒸散发系数

this->m_SM = Params[7]; // (8)自由水蓄水容量

this->m_EX = Params[8]; // (9)自由水蓄水容量～面积分布曲线指数

this->m_KG = Params[9]; // (10)地下水日出流系数

this->m_KI = Params[10]; // (11)壤中流日出流系数

this->m_CG = Params[11]; // (12)地下水消退系数

this->m_CI = Params[12]; // (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 Ep; //m_Kc * m_pEm[i]

double 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; // 深层土壤蒸散发量（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_pP[i] * m_IM; // RB是降在不透水面的降雨量

P = m_pP[i] * (1 - m_IM);

Ep = m_Kc * m_pEm[i];

if ((WU + P)>= Ep)

{EU = Ep; EL = 0; ED = 0; }

elseif((WU + P)

{

EU = WU + P;

if(WL>= (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 (WL

{EL = WL; ED = m_C * (Ep - EU) - EL; }

}

E = EU + EL + ED;

PE = P - E;

/* ———————蒸散发计算结束——————————— */ //——————子流域产流量计算————————————// if(PE<= 0)

{R = 0.00;W = W + PE; }

else

{

A = 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; /* 不透水面上产生的径流*/

}

// 土壤湿度折算净雨量+降水后蒸发剩余雨量<流域内最大含水容量else

{

// 流域内的产流深度计算

R = PE /* 降水蒸发后的剩余量+ W /* 流域内土壤湿度*/

- m_WM /* 减去流域平均蓄水容量*/

+ RB; /* 不透水面上产生的径流*/

}

}

//三层蓄水量的计算: WU, WL, WD

if(WU + P - EU - R <= m_Wum)

{WU = WU + P - EU - R;WL = WL - EL; WD = WD – ED; }

else

{

WU = m_Wum;

if(WL - EL + ( WU + P - EU - R - m_Wum ) <= m_Wlm )

{

WL = WL – EL + ( WU + P - EU - R - m_Wum );

WD = WD - ED;

}

else

{

WL = m_Wlm;

if(WD - ED + WL - EL + ( WU + P - EU - R - m_Wum )

- m_Wlm <= m_Wdm )

WD = WD - ED + WL - EL

+ (WU + P - EU - R - m_Wum ) - m_Wlm;

else

WD = m_Wdm;

}

}

W = WU + WL + WD;

////三水源划分汇流计算

if(PE>0)

{

FR = (R - RB) / PE;

AU = MS * (1 - pow((1 - So * FRo/FR/m_SM),1/(1 + m_EX)));

if(PE + AU

RS = FR * (PE + So * FRo/FR - m_SM + m_SM * pow((1 - (PE + AU) / MS),m_EX + 1));

elseif(PE + AU>= MS)

RS = FR * ( PE + So * Fro / FR - m_SM );

S = So * Fro / FR + ( R – RS ) / FR;

RI = m_KI * S * FR;

RG = m_KG * S * FR;

RS += RB;

R = RS + RI + RG;

So = S * ( 1 - m_KI - m_KG );

FRo = FR;

}

else

{

S = 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_pE[i] = E; // 当前步长的蒸发（模型重要输出）/* 2 */this->m_pEu[i] = EU; // 当前步长上层土壤蒸发

/* 3 */this->m_pEl[i] = EL; // 当前步长下层土壤蒸发

/* 4 */this->m_pEd[i] = ED; // 当前步长深层土壤蒸发

/* 5 */this->m_pW[i] = W; // 当前步长流域平均土壤含水量

/* 6 */this->m_pWu[i] = WU; // 当前步长流域上层土壤含水量

/* 7 */this->m_pWl[i] = WL; // 当前步长流域下层土壤含水量

/* 8 */this->m_pWd[i] = WD; // 当前步长流域深层土壤含水量

/* 9 */this->m_pRg[i] = RG; // 当前步长流域地下径流深度

/* 10 */this->m_pRi[i] = RI; // 当前步长流域壤中流深度

/* 11 */this->m_pRs[i] = RS; // 当前步长流域地表径流径流深度

/* 12 */this->m_pR[i] = 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;

}

fprintf(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.3lf\n",

this->m_pE[i],this->m_pEu[i],this->m_pEl[i],this->m_pEd[i],

this->m_pW[i],this->m_pWu[i],this->m_pWl[i],this->m_pWd[i],

this->m_pR[i],this->m_pRs[i],this->m_pRi[i],this->m_pRg[i],

this->m_pQ[i],this->m_pQrs[i],this->m_pQri[i],this->m_pQrg[i]);

}

fclose(fp);

}

// 进行汇流计算，将径流深度转换为流域出口的流量

void XinanjiangModel::Routing(void)

{

///////////// 地面径流汇流计算：单位线法///////////////////////

int i,j;

double B[10000] = {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;j<18;j++ )

{B[i + j] += this->m_pRs[i] * UH[j]/10.0;}

}

}

else

{

double UH[] ={7.18,23.38,63.20,143.10,221.75,365.18,447.40,491.29,

506.93,504.82,468.46,388.56,309.91,166.49,84.26,40.37,17.56,3.46};

for(i = 0;im_nSteps;i++ )

{

for(j = 0;j<18;j++ )

{B[i + j] += this->m_pRs[i] * UH[j]/10.0; }

}

}

for(i = 0;isteps;i++ )

this->m_pQrs[i] = B[i];

/////壤中流汇流计算:线性水库

for(i = 1;isteps;i++ ) {this->m_pQri[i] = this->m_CI * this->m_pQri[i - 1]

+ (1.0 - this->m_CI) * this->m_pRi[i] * this->m_U;} /////地下径流汇流计算：线性水库

for(i = 1;isteps;i++ ) {this->m_pQrg[i] = this->m_pQrg[i - 1] *

this->m_CG

+ this->m_pRg[i] * (1.0 - this->m_CG) * this->m_U; } //////单元面积总入流计算

for(i = 0;isteps;i++ )

{this->m_pQ[i] = this->m_pQrs[i] + this->m_pQri[i] + this->m_pQrg[i]; }

}

void XinanjiangModel::Runoff(char* runoff)

{

int i;

ofstream outfile;

outfile.open(runoff);

if(outfile.is_open())

{

for(i = 0;isteps;i++ )

{outfile<m_pQ[i]<

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,IM,B,m_Wum,Wlm,Wdm,C,SM,EX,KG,KI,CG,CI

double Params1[] = {0.50,0.01,0.30,10,60,40,0.18,32,1.2,0.075,0.072,0.94,0.7}; Model1.SetParameters(Params1);

Model1.RunModel();

Model1.SaveResults("来凤站日模型计算结果.txt");

Model1.Runoff("LF_Q.txt");

Model1.~XinanjiangModel();

double Area2 = 2991;

XinanjiangModel Model2;

Model2.InitModel(nSteps, Area2,DeltaT,0,"YCForcingfile.txt");

//模型参数/*Kc,IM,B,m_Wum,Wlm,Wdm,C,SM,EX,KG,KI,CG,CI

double Params2[] = {0.75,0.01,0.32,10,60,40,0.18,27,1.2,0.065,0.067,0.96,0.8}; Model2.SetParameters(Params2);

Model2.RunModel();

Model2.SaveResults("file.txt");

Model2.Runoff("YC_Q.txt");

Model2.~XinanjiangModel();

FILE *fp1,*fp2;

double Q1[1000],Q2[1000],Q[1000]={0.00};

ofstream outfile;

outfile.open("Q.txt");

if((fp1=fopen("LF_Q.txt","r"))==NULL)

{printf("Can not open the file!\n");return 0;}

if((fp2=fopen("YC_Q.txt","r"))==NULL)

{printf("Can not open the file!\n");return 0; }

if(outfile.is_open())

{

for(int i=0;i<960;i++)

{

fscanf(fp1,"%lf",&Q1[i]);

fscanf(fp2,"%lf",&Q2[i]);

Q[i]=Q1[i]+Q2[i];

outfile<

}

fclose(fp1);

fclose(fp2);

}

outfile.close() ;

return 0;

}