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

新安江模型程序C+代码以下是类的声明：class XinanjiangModelprivate:/ FORCINGdouble *m_pP; / 降水数据double *m_pEm; / 水面蒸发数据/long m_nSteps; / 模型要运行的步长(一共m_nSteps步)long steps;/ OUTPUTdouble *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(km2)/3.6/delta_t/ SOILdouble *m_pW; / 流域内土壤湿度double *m_pWu; / 流域内上层土壤湿度double *m_pWl; / 流域内下层土壤适度double *m_pWd; / 流域内深层土壤湿度double m_Wum; / 流域内上层土壤蓄水容量double m_Wlm; / 流域内下层土壤蓄水容量double m_Wdm; / 流域内深层土壤蓄水容量，WDM=WM-WUM-WLM / EVAPORATIONdouble *m_pEu; / 上层土壤蒸发量（毫米）double *m_pEl; / 下层土壤蒸发量（毫米）double *m_pEd; / 深层土壤蒸发量（毫米）/runoffdouble *RF;/ PARAMETERdouble m_Kc; / 流域蒸散发能力与实测蒸散发值的比double m_IM; / 不透水面积占全流域面积之比double m_B; / 蓄水容量曲线的方次，小流域（几平方公里）B0.1左右 / 中等面积（平方公里以内）.20.3，较大面积.30.4 double m_WM; / 流域平均蓄水容量（毫米）(WM=WUM+WLM+WDM)double m_C; / 流域内深层土壤蒸发系数，江南湿润地区：0.15-0.2，/华北半湿润地区：.09-0.12double 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 using namespace std;#include math.h#include stdio.h#include conio.hXinanjiangModel: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 doublethis-steps;this-m_pEm = new doublethis-steps;/ 模型输出，蒸散发项this-m_pE = new doublethis-steps;this-m_pEd = new doublethis-steps;this-m_pEl = new doublethis-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;/ 模型输出，产流项this-m_pR = new doublethis-steps;this-m_pRg= new doublethis-steps;this-m_pRi= new doublethis-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值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_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文件格式：第一列：降水（单位毫米）空格第二列水面蒸发（毫米）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) 流域不透水面积占全流域面积之比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; / (9)自由水蓄水容量面积分布曲线指数this-m_KG = Params9; / (10)地下水日出流系数this-m_KI = Params10; / (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 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; / 深层土壤蒸散发量（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 = m_pPi * (1 - m_IM);Ep = m_Kc * m_pEmi;if (WU + P)= Ep) EU = Ep; EL = 0; ED = 0; else if(WU + P)= (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 (WLm_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; 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; /* 不透水面上产生的径流*/ 土壤湿度折算净雨量+降水后蒸发剩余雨量流域内最大含水容量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;if(WL - EL + ( WU + P - EU - R - m_Wum ) = m_Wlm )WL = 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 * 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_pEi = E; / 当前步长的蒸发（模型重要输出）/* 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; / 当前步长流域地下径流深度/* 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;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.3lfn,this-m_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;double 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.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+ )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,LFForci