标签：name -- output mes 记录 mil ++ man model
以下是类的声明：
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;
}
标签：name -- output mes 记录 mil ++ man model
原文地址：https://www.cnblogs.com/yellowhh/p/12038310.html
