topmodel.cpp

这是一个GPS相关的程序

			for(n=0; n<gClass.Get_NCells(); n++)
			{
				iClass	= gClass.asInt(n);

				if( iClass >= 0 && iClass < nClasses )
				{
					pMoist->Set_Value(n, Vals.Get_Class(iClass)->S_);
				}
				else
				{
					pMoist->Set_NoData(n);
				}
			}

		//	DataObject_Update(pMoist);
			DataObject_Update(pMoist, 0, 0.35, true);
		}

		pRecord	= pTable->Add_Record();
		pRecord->Set_Value(0, Vals.Qt_[iTime]);		// QT
		pRecord->Set_Value(1, Vals.qt_Total);		// qt
		pRecord->Set_Value(2, Vals.qo_Total);		// q0
		pRecord->Set_Value(3, Vals.qs_Total);		// qs
		pRecord->Set_Value(4, Vals.qv_Total);		// qv
		pRecord->Set_Value(5, Vals.Sbar_);			// SBar
		pRecord->Set_Value(6, Infiltration);		// Infiltration
		pRecord->Set_Value(7, Infiltration_Excess);	// Infiltration Excess
		DataObject_Update(pTable);
	}

	//-----------------------------------------------------
	return( true );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CTOPMODEL::Run(double Evaporation, double Precipitation, double Infiltration_Excess)
{
	int				iClass;
	double			d, Excess;
	CTOPMODEL_Class	*pClass;

	Vals.qo_Total	= 0.0;
	Vals.qv_Total	= 0.0;
	Vals.qs_Total	= Vals._qs_ * exp(-Vals.Sbar_ / Vals.p_Model);

	for(iClass=0; iClass<Vals.Get_Count(); iClass++)
	{
		pClass			= Vals.Get_Class(iClass);


		//-------------------------------------------------
		//  CALCULATE LOCAL STORAGE DEFICIT

		pClass->S_		= Vals.Sbar_ + Vals.p_Model * (Vals.Get_Lambda() - pClass->AtanB);

		if( pClass->S_ < 0.0 )
		{
			pClass->S_		= 0.0;
		}


		//-------------------------------------------------
		//  ROOT ZONE CALCULATIONS

		pClass->Srz_	-= Precipitation;

		if( pClass->Srz_ < 0.0 )
		{
			pClass->Suz_	-= pClass->Srz_;
			pClass->Srz_	= 0.0;
		}


		//-------------------------------------------------
		//  UNSATURATED ZONE CALCULATIONS

		if( pClass->Suz_ > pClass->S_ )
		{
			Excess			= pClass->Suz_ - pClass->S_;
			pClass->Suz_	= pClass->S_;
		}
		else
		{
			Excess			= 0.0;
		}


		//-------------------------------------------------
		//  CALCULATE DRAINAGE FROM SUZ (Vertical Soil Water Flux (qv))...

		if( pClass->S_ > 0.0 )
		{
			if( Vals.p_Suz_TimeDelay > 0.0 )
			{	// Methode 1...
				d			= pClass->Suz_ / (pClass->S_ * Vals.p_Suz_TimeDelay) * dTime;	// GRASS
			}
			else
			{	// Methode 2...
				d			= -Vals.p_Suz_TimeDelay * Vals.p_K0 * exp(-pClass->S_ / Vals.p_Model);
			}

			if( d > pClass->Suz_ )
			{
				d			= pClass->Suz_;
			}

			pClass->Suz_	-= d;

			if( pClass->Suz_ < 0.0000001 )
			{
				pClass->Suz_	= 0.0;
			}

			pClass->qv_		= d * pClass->Area_Rel;
			Vals.qv_Total	+= pClass->qv_;
		}
		else
		{
			pClass->qv_		= 0.0;
		}


		//-------------------------------------------------
		//  CALCULATE EVAPOTRANSPIRATION FROM ROOT ZONE DEFICIT

		if( Evaporation > 0.0 )
		{
			d		= Evaporation * (1.0 - pClass->Srz_ / Vals.p_Srz_Max);

			if( d > Vals.p_Srz_Max - pClass->Srz_ )
			{
				d		= Vals.p_Srz_Max - pClass->Srz_;
			}

			pClass->Srz_	+= d;
		}


		//-------------------------------------------------
		pClass->qo_		= Excess * pClass->Area_Rel;
		Vals.qo_Total	+= pClass->qo_;

		pClass->qt_		= pClass->qo_ + Vals.qs_Total;
	}

	Vals.qo_Total	+= Infiltration_Excess;

	Vals.qt_Total	= Vals.qo_Total + Vals.qs_Total;

	Vals.Sbar_		+= Vals.qs_Total - Vals.qv_Total;
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CTOPMODEL::Get_Climate(int iTimeStep, double &Precipitation, double &Evaporation)
{
	CSG_Table_Record	*pRecord;

	if( pClimate && pClimate->Get_Field_Count() >= 2 && (pRecord = pClimate->Get_Record(iTimeStep)) != NULL )
	{
		Precipitation	= pRecord->asDouble(0);
		Evaporation		= pRecord->asDouble(1);

		return( true );
	}
	else
	{
		Precipitation	= 0.0;
		Evaporation		= 0.0;

		return( false );
	}
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#define NEWTON_EPSILON		0.001
#define NEWTON_MAXITER		100
#define	NEWTON_NTERMS		10

//---------------------------------------------------------
double CTOPMODEL::Get_Infiltration(double t, double R)
{
	int		i, j, factorial;

	double	f, f_, f1, f2, fc, R2, cnst, pt, psi_dtheta, sum;


	if( R <= 0.0 )
	{
		inf_cumf		= 0.0;
		inf_bPonding	= 0;

		return( 0.0 );
	}

	psi_dtheta	= Vals.p_Psi * Vals.p_dTheta;

	if( !inf_bPonding )
	{
		if( inf_cumf )
		{
			f1				= inf_cumf;
			R2				= -Vals.p_K0 / Vals.p_Model * (psi_dtheta + f1) / (1 - exp(f1 / Vals.p_Model));

			if( R2 < R )
			{
				f_				= inf_cumf;
				pt				= t - dTime;
				inf_bPonding	= 1;

				goto cont1;
			}
		}

		f2				= inf_cumf + R * dTime;
		R2				= -Vals.p_K0 / Vals.p_Model * (psi_dtheta + f2) / (1 - exp(f2 / Vals.p_Model));

		if( f2 == 0.0 || R2 > R )
		{
			f				= R;
			inf_cumf		+= f * dTime;
			inf_bPonding	= 0;

			return( f );
		}

		f_				= inf_cumf + R2 * dTime;

		for(i=0; i<NEWTON_MAXITER; i++)
		{
			R2				= -Vals.p_K0 / Vals.p_Model * (psi_dtheta + f_) / (1 - exp(f_ / Vals.p_Model));

			if( R2 > R )
			{
				f1				= f_;
				f_				= (f_ + f2) / 2.0;
				f				= f_ - f1;
			}
			else
			{
				f2				= f_;
				f_				= (f_ + f1) / 2.0;
				f				= f_ - f2;
			}

			if( fabs(f) < NEWTON_EPSILON )
				break;
		}

		if( i == NEWTON_MAXITER )
		{
			// G_set_d_null_value(&f, 1);
			return( 0.0 );
		}

		pt				= t - dTime + (f_ - inf_cumf) / R;

		if( pt > t )
		{
			f				= R;
			inf_cumf		+= f * dTime;
			inf_bPonding	= 0;

			return( f );
		}

cont1:
		cnst			= 0.0;
		factorial		= 1;
		fc				= (f_ + psi_dtheta);

		for(j=1; j<=NEWTON_NTERMS; j++)
		{
			factorial		*= j;
			cnst			+= pow(fc / Vals.p_Model, (double) j) / (double) (j * factorial);
		}

		cnst			= log(fc) - (log(fc) + cnst) / exp(psi_dtheta / Vals.p_Model);
		f_				+= R * (t - pt) / 2.0;
		inf_bPonding	= 1;
	}

	for(i=0; i<NEWTON_MAXITER; i++)
	{
		fc				= f_ + psi_dtheta;
		sum				= 0.0;
		factorial		= 1;

		for(j=1; j<=NEWTON_NTERMS; j++)
		{
			factorial		*= j;
			sum				+= pow(fc / Vals.p_Model, (double) j) / (double) (j * factorial);
		}

		f1				= - (log(fc) - (log(fc) + sum) / exp(psi_dtheta / Vals.p_Model) - cnst) / (Vals.p_K0 / Vals.p_Model) - (t - pt);
		f2				= (exp(f_ / Vals.p_Model) - 1.0) / (fc * Vals.p_K0 / Vals.p_Model);
		f				= - f1 / f2;
		f_				+= f;

		if( fabs(f) < NEWTON_EPSILON )
			break;
	}

	if( i == NEWTON_MAXITER )
	{
		// G_set_d_null_value(&f, 1);
		return( 0.0 );
	}

	if( f_ < inf_cumf + R )
	{
		f				= (f_ - inf_cumf) / dTime;
		inf_cumf		= f_;
		f_				+= f * dTime;
	}

	return( f );
}