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

1、新安江模型程序C+代碼以下是類的聲明：class XinanjiangModelprivate:/ FORCINGdouble *m_pP; / 降水?dāng)?shù)據(jù)double *m_pEm; / 水面蒸發(fā)數(shù)據(jù)/long m_nSteps; /模型要運(yùn)彳T的步長(zhǎng)（一共m_nStep涉）long steps;/ OUTPUTdouble *m_pR; / 流域內(nèi)每一步長(zhǎng)的產(chǎn)流量（徑流深度）double *m_pRs; / 每一步長(zhǎng)的地表徑流深（毫米 ）double *m_pRi; / 每一步長(zhǎng)的壤中流深（毫米）double *m_pRg; / 每一步長(zhǎng)的地下徑流深（毫米 ）double *m_pE; /

2、每一步長(zhǎng)的蒸發(fā) （毫米 ）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;/流域內(nèi)土壤濕度double *m_pWu; / 流域內(nèi)上層土壤濕度double *m_pWl; / 流域內(nèi)下層土壤適度double *m_pWd;/流域內(nèi)深層土壤濕度double m_Wum;/流域內(nèi)上層土壤蓄水容量double m_Wl

3、m; / 流域內(nèi)下層土壤蓄水容量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、/ 流域內(nèi)深層土壤蓄水容量， WDM=WM-WUM-WLM/ EVAPORATION / 上層土壤蒸發(fā)量（毫米）/ 下層土壤蒸發(fā)量（毫米）/ 深層土壤蒸發(fā)量（毫米）/ 流域蒸散發(fā)能力與實(shí)測(cè)蒸散發(fā)值的比/ 不透水面積占全流域面積之比/蓄水容量曲線的方次，小流域（幾平方公里）B0.1左右/ 中等面積（平方公里以內(nèi)）.20.3，較大面積.30.4/ 流域平均蓄水容量 （毫米） （WM=WUM+WLM+WDM）/ 流域內(nèi)深層土壤蒸發(fā)系數(shù)，江南濕潤(rùn)地區(qū)：0.15-0.2，/華北半濕潤(rùn)地區(qū)：.09-0.12/自由水蓄水容量自由水蓄水容量面積分布曲線指數(shù)/地下水日出流系數(shù)/壤中流日出流系數(shù)/地下水消退系數(shù)/

5、壤中流消退系數(shù)/ 單元流域上地面徑流的單位線/ 流域內(nèi)最大蓄水容量/ 流域面積/ 每一步長(zhǎng)的小時(shí)數(shù)/ 給定數(shù)據(jù)，用以判斷是否時(shí)行河道匯流計(jì)算XinanjiangModel(void);XinanjiangModel(void);/ 初始化模型void InitModel( long nSteps,double Area,int DeltaT, int PD, char *ForcingFile);/ 設(shè)置模型參數(shù)void SetParametersd(ouble *Params);/ 運(yùn)行新安江模型void RunModel(void);/ 保存模擬結(jié)果到文件void SaveResults(

6、char *FileName);/ 記錄出流數(shù)據(jù)，用以作圖分析void Runoff(char *runoff);private:/ 進(jìn)行匯流計(jì)算，將徑流深度轉(zhuǎn)換為流域出口的流量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;/ 驅(qū)動(dòng)數(shù)據(jù)this-m_pP = new doublethis-steps;this-m_pEm = new doublethis-steps;/

10、 模型輸出，蒸散發(fā)項(xiàng)this-m_pE = new doublethis-steps;this-m_pEd = new doublethis-steps;this-m_pEl = newdoublethis-steps;this-m_pEu = new doublethis-steps;/ 模型輸出，出流項(xiàng)，經(jīng)過(guò)匯流的產(chǎn)流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、/ 模型輸出，產(chǎn)流項(xiàng)this-m_pR = new doublethis-steps;this-m_pRg= newdoublethis-steps;this-m_pRi= newdoublethis-steps; this-m_pRs = new doublethis-steps; / 模型狀態(tài)量，土壤濕度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+ )/ 驅(qū)動(dòng)數(shù)據(jù)this-m_pP i = 0.00;this-m_pEm i = 0.00;/ 模型輸出，蒸散發(fā)項(xiàng)this-m_pE i = 0.00;this-m_pEd i = 0.00;this-m_pEl i = 0.00;this-m_pEu i = 0.00;/ 模型輸出，出流項(xiàng)，經(jīng)過(guò)匯流的產(chǎn)流this-m_pQrgi = 0.00;this-m_pQrsi = 0.00;this-m_pQrii = 0.00;this-m_pQi = 0.00;/ 模型輸出，產(chǎn)流項(xiàng)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;/ 模型狀態(tài)量，土壤濕度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文件格式：第一列：降水(單位毫米)空格第二列水面蒸發(fā)(毫

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);/ 設(shè)置模型參數(shù)void XinanjiangModel:SetParameters(double* Params)this-m_Kc = Params0; / (1) 流域蒸散發(fā)能力與實(shí)測(cè)水面蒸發(fā)之比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) 深層蒸散發(fā)系數(shù)this-m_SM = Params7; / (8)自由水蓄水容量this-m_EX = Params8; /自由水蓄水容量面積分布曲線指數(shù)this-m_KG = Params9;/ （10）地下水日出流系數(shù)this-m_KI = P

16、arams10;/ （11）壤中流日出流系數(shù)this-m_CG = Params11;/ （12）地下水消退系數(shù)this-m_CI = Params12;/ （13）壤中流消退系數(shù)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）;/ 運(yùn)行新安江模型void XinanjiangModel:RunModel（void）long i;/ 模型的狀態(tài)變量double PE; / 0 時(shí)為凈雨量; 0 為蒸發(fā)不足量（ mm）double

17、 Ep; /m_Kc * m_pEmidouble P;double R; / 產(chǎn)流深度，包括地表徑流、壤中流和地下徑流（ mm）double RB; / 不透水面上產(chǎn)生的徑流深度（ mm）double RG; / 地下徑流深度（mm）double RI; / 壤中流深度（mm）double RS; / 地表徑流深（mm）double A;/土壤濕度為W時(shí)土壤含水量折算成的徑流深度（mm）double E = 0.0;/ 蒸散發(fā) （mm）double EU = 0.0;/ 上層土壤蒸散發(fā)量（mm）double EL = 0.0;/下層土壤蒸散發(fā)量（mm）double ED =0.0;/深層土壤

18、蒸散發(fā)量（mm）double S;double FRo;double FR;double MS;double AU;double WU = 5.0;/ 流域內(nèi)上層土壤濕度double WL = 55.0; / 流域內(nèi)下層土壤適度double WD = 40.0; / 流域內(nèi)深層土壤濕度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+ )/ 蒸散發(fā)計(jì)算 /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;/* 蒸散發(fā)計(jì)算結(jié)束 */子流域產(chǎn)流量計(jì)算 /if(PE= 0) R = 0.00; W = W + PE; elseA = m_WMM * (1 - pow( (1.0 - W/m_WM), 1.0/(1 + m_B) ) );/ 土壤濕度折算凈雨量+降水后蒸發(fā)剩余雨量流域內(nèi)最大含水容量if(A + PE)m_WMM)/ 流域內(nèi)的產(chǎn)流深度計(jì)算R = PE/* 降水蒸發(fā)后的剩余量*/+ W/* 流域內(nèi)土壤濕度*/+ m_WM * pow(1 - (PE + A)/m_WMM),(1 + m_B)- m_WM /* 減去流域平均蓄水容量( m_WM: 參數(shù)) */+ RB; /*

21、 不透水面上產(chǎn)生的徑流*/ 土壤濕度折算凈雨量+降水后蒸發(fā)剩余雨量流域內(nèi)最大含水容量else/ 流域內(nèi)的產(chǎn)流深度計(jì)算R = PE/* 降水蒸發(fā)后的剩余量+ W /* 流域內(nèi)土壤濕度*/- m_WM /* 減去流域平均蓄水容量*/+ RB; /* 不透水面上產(chǎn)生的徑流*/三層蓄水量的計(jì)算: 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;/三水源劃分計(jì)算結(jié)束*/* 以下部分是狀態(tài)量：總蒸發(fā)量、上、下和深層土壤的蒸發(fā)的保存/* 1 */this-m_pEi = E;/ 當(dāng)前步長(zhǎng)的蒸發(fā)(模型重要輸

24、出)/* 2 */this-m_pEui = EU;/ 當(dāng)前步長(zhǎng)上層土壤蒸發(fā)/* 3 */this-m_pEli = EL;/ 當(dāng)前步長(zhǎng)下層土壤蒸發(fā)/* 4 */this-m_pEdi = ED;/ 當(dāng)前步長(zhǎng)深層土壤蒸發(fā)/* 5 */this-m_pWi = W;/ 當(dāng)前步長(zhǎng)流域平均土壤含水量/* 6 */ this-m_pWui = WU; / 當(dāng)前步長(zhǎng)流域上層土壤含水量/* 7 */ this-m_pWli = WL; / 當(dāng)前步長(zhǎng)流域下層土壤含水量/* 8 */ this-m_pWdi = WD; / 當(dāng)前步長(zhǎng)流域深層土壤含水量/* 9 */ this-m_pRgi = RG; / 當(dāng)前

25、步長(zhǎng)流域地下徑流深度/* 10 */ this-m_pRii = RI; / 當(dāng)前步長(zhǎng)流域壤中流深度/* 11 */ this-m_pRsi = RS; / 當(dāng)前步長(zhǎng)流域地表徑流徑流深度/* 12 */ this-m_pRi = R;/ 當(dāng)前步長(zhǎng)的總產(chǎn)流徑流深度this-Routing();/ 保存模擬結(jié)果到文件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);/ 進(jìn)行匯流計(jì)算，將徑流深度轉(zhuǎn)換為流域出口的流量void XinanjiangModel:Routing(void)/ 地面徑流匯流計(jì)算：?jiǎn)挝痪€法 /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;/ 壤中流匯流計(jì)算:線性水庫(kù)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; / 地下徑流匯流計(jì)算：線性水庫(kù)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; /單元面積總?cè)肓饔?jì)算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()函數(shù)語(yǔ)句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 );/模型參數(shù) /*Kc,