新安江模型C++代码 VS2008平台下

// my xinanjiang.cpp : 定义控制台应用程序的入口点。
//新安江模型
#include “stdafx.h“
#include “Model parameters.h“
#include 
#include 
#include 
#include 
using namespace std;
void XAJ（）
{	ifstream ifile;
	ofstream ofile;
	ofile.open（“G:\\output.txt“）;
	ifile.open（“G:\\data.txt“）;
	ifile>>Date_Num;
	vector P;
	vector E;
	vector QR;
	P.resize（Date_Num）;
	E.resize（Date_Num）;
	QR.resize（Date_Num+3）;
	for （int  i=0;i	{
		ifile>>E[i]>>P[i];
		QR[i]=0;
		if （E[i]<0）
		{
			E[i]=0;
		}
		if （P[i]<0）
		{
			P[i]=0;
		}
	}
	ifile.close（）;
	double RRSRSSRGRIM;
	double EUELEDPE;
	double XQUDCICGKSSDKGD;
	double SSNNKSSDDKGDDSMMSMMFSMFAURGDRSDRSSDQRSQRSSQRGQTR;
	//初始化参数
	U= AREA / （DT * 3.6）;
 	if （DT <= 24） 
 	{
 		D = 24 / DT;
 		CI = pow（KKSS（1.0 / D））;
 		CG = pow（KKG（1.0 / D））;
 		KSSD = （1.0 -pow（（1.0-（KG+KSS））（1.0/D）））/（1.0+KG/KSS）;
 		KGD = KSSD * KG / KSS;
 	}
	
	//计算不透水面积的产流量且直接转化为地面径流
	
	//计算透水面积的产流量
	double WMM;
	double  A;
	for（int i=0;i	{
		WM=WLM+WUM+WDM;
		WO=WL+WU+WD;
		PE=P[i]-K*E[i];
		R=0;
		RIM=0;
	if （PE>0）//净雨大于0，产流，计算土壤含水量变化
		{
			//土壤的含水量WL，WU，WD会发生变化
			//流域土壤的WLM，WUM，WDM平均蓄水容量不会发生变化

			WMM=（（1+B）/（1-IMP））*WM;
			if （（WMM-WM）<=0.000001）
			{
				A=WMM;
			}
			else
			{
				A=WMM*（1-pow（（1-WO/WM）（1/（B+1））））;
			}
			if （PE+A			{
				R=PE-WM+WO+WM*pow（（1-（PE+A）/WMM）（B+1））;
			}
			else
			{
				R=PE-（WM-WO）;
			}
			RIM=PE*IMP;
		}
		//大于蒸发能力
		if（P[i]+WU>E[i]）
			{
				EU=K*E[i];
				ED=0;
				EL=0;
			}
		else//要求下层蒸发量与剩余蒸散发能力不小于深层蒸散发吸吮C
		{
			EU=WU+P[i];
			EL=（K*E[i]-EU）*WL/WLM;
			ED=0;
		if （WL<=C*WLM）
			{
				if （WL >= C * （K*E[i] - EU））
				{ 
					EL = C * （K*E[i] - EU）;
					ED=0;
				}
				else
				{ 
				EL = WL;
				ED = C * （K*E[i] - EU） - EL;
				}
			}
		}
	//土壤水计算
	WU = WU + P[i] - R - EU;
	WL = WL - EL;
	WD = WD - ED;
	if （WU < 0）  
	{
		WU = 0;
	}
	if （WU > WUM） 
	{
		WL = WU - WUM + WL;
		WU = WUM;
	if （WL > WLM） 
	{
		WD = WD + WL - WLM;
		WL = WLM;
	}
	}
	//三水源的划分
	X=FR;
	if （PE<=0）
	{
		RS=0;
		RSS=S*KSSD*FR;
		RGD=S*KGD*FR;
		S=S-（RSS+RG）/FR;
	}
	else
	{
		FR=R/PE;
		S=X*S/FR;
		SS=S;
		Q=R/FR;
		NN=int（Q/5.0）+1;
		Q=Q/NN;
		KSSDD = （1.0 -pow（（1.0-（KGD+KSSD））（1.0/NN）））/（1.0+KGD/KSSD）;
		KGDD = KSSDD * KGD / KSSD;
		RS=0;
		RSS=0;
		RG=0;
		SMM=（1.0+EX）*SM;
		if （EX<0.001）
		{
			SMMF=SMM;		
		}
		else
		{
			SMMF=SMM*（1.0 - pow（（1.0- FR）（1.0/ EX）））;
		}
	//}
	SMF=SMMF/（1.0+EX）;
	for （int k=1;k<=NN;k++）
	{
		if （S>SMF）
		{
			S=SMF;
		}
		AU=SMMF * （1 -pow（ （1 - S / SMF）  （1 / （1 + EX））））;
	if （（Q+AU）<=0）
	{
		RSD = 0;
		RSSD = 0;
		RGD = 0;
		S = 0;
	}
	else if （Q+AU>=SMMF）
	{
		RSD = （Q + S - SMF） * FR;
		RSSD = SMF * KSSDD * FR;
		RGD = SMF * KGDD * FR;
		S = SMF - （RSSD + RGD） / FR;
	}
	else if （Q+AU	{