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

新安江模型程序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 <iostream>
#include <fstream>
#include <iomanip>
using namespace 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 = new double[this->steps];
this->m_pEm = new double[this->steps];
// 模型输出，蒸散发项
this->m_pE = new double[this->steps];
this->m_pEd = new double[this->steps];
this->m_pEl = new double[this->steps];
this->m_pEu = new double[this->steps];
// 模型输出，出流项，经过汇流的产流
this->m_pQrg = new double[this->steps];
this->m_pQrs = new double[this->steps];
this->m_pQri = new double[this->steps];
this->m_pQ = new double[this->steps];
// 模型输出，产流项
this->m_pR = new double[this->steps];
this->m_pRg= new double[this->steps];
this->m_pRi= new double[this->steps];
this->m_pRs = new double[this->steps];
// 模型状态量，土壤湿度
this->m_pW = new double[this->steps];
this->m_pWd = new double[this->steps];
this->m_pWl = new double[this->steps];
this->m_pWu = new double[this->steps];
//runoff值
this->RF = new double[this->steps];
for(i = 0;i<this->steps;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;i<this->m_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;i<this->m_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; }
else if((WU + P)<Ep)
{
EU = WU + P;
if(WL>= (m_C * m_Wlm))
{ EL = (Ep - EU) * WL/m_Wlm; ED = 0; }
else if ((m_C * (Ep - EU))<= WL&&WL<(m_C * m_Wlm))
{ EL = m_C * (Ep - EU); ED = 0; }
else if (WL<m_C * (Ep - EU))
{ 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)<this->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<MS)
RS = FR * (PE + So * FRo/FR - m_SM + m_SM * pow((1 - (PE
+ AU) / MS),m_EX + 1));
else if(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;i<this->steps;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;i<this->m_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;i<this->m_nSteps;i++ )
{
for(j = 0;j<18;j++ )
{ B[i + j] += this->m_pRs[i] * UH[j]/10.0; }
}
}
for(i = 0;i<this->steps;i++ )
this->m_pQrs[i] = B[i];
///// 壤中流汇流计算:线性水库
for(i = 1;i<this->steps;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;i<this->steps;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;i<this->steps;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;i<this->steps;i++ ) { outfile<<setprecision(3)<<setiosflags(ios::fixed)<<this->m_pQ[i]<<endl; }
}
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<<setprecision(3)<<setiosflags(ios::fixed)<<Q[i]<<endl;
}
fclose(fp1);
fclose(fp2);
}
outfile.close() ;
return 0;
}