dvwk_soilmoisture.cpp

这是一个GPS相关的程序

	Day		%= 365;

	if( Day < 0 )
	{
		Day 	+= 365;
	}

	for(iMonth=0, nDays=0; iMonth<12; iMonth++)
	{
		nDays	+= Days[iMonth];

		if( Day < nDays )
		{
			return( iMonth + 1 );
		}
	}

	return( 12 );
}

//---------------------------------------------------------
// VKR 4.10:	Bestandeskoeffizient z.Ber.d.pot.Verdunstung d.Pflbest. (kc)
//
//		Bestand		: Typ (als ENUM) des Pflanzenbestandes
//		Month		: 1=Jan, 2=Feb... 0=NA;
//
double CDVWK_SoilMoisture::Get_kc(int Bestand, int Day)
{
	if( Bestand >= 0 && Bestand < pCropCoeff->Get_Record_Count() )
	{
		return( pCropCoeff->Get_Record(Bestand)->asDouble(1 + Get_Month(Day)) );
	}

	return( 1.0 );
}


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

//---------------------------------------------------------
double CDVWK_SoilMoisture::Get_ETP_Haude(int Day)
{
	const double	f[]	=
	//	---   Jan   Feb   Mar   Apr   Mai   Jun   Jul   Aug   Sep   Okt   Nov   Dez
	{	0.00, 0.22, 0.22, 0.22, 0.29, 0.29, 0.28, 0.26, 0.25, 0.23, 0.22, 0.22, 0.22	};

	double	T14, U14, es, e, ETP_Haude;

	//-----------------------------------------------------
	CSG_Table_Record	*pRecord;

	if( (pRecord = pClimate->Get_Record(Day)) == NULL )
	{
		return( 0 );
	}

	T14		= pRecord->asDouble(1);
	U14		= pRecord->asDouble(2);


	//-----------------------------------------------------
	// VKR 4.6:		Berechnung des Saettigungsdampfdruckes (es [hPa])
	//
	//					               17,62 * T14
	//		es(T14)		= 6,11 * exp( -------------- );		bei T14 > 0
	//					               243,12 + T14
	//
	//					               22,46 * T14
	//		es(T14)		= 6,11 * exp( -------------- );		bei T14 < 0
	//					               272,62 + T14
	//
	//		T14			: Temperatur 14 Uhr [DegreeC]
	//

	if( T14 < 0 )
	{
		es	= 6.11 * exp( (22.46 * T14) / (272.62 + T14) );
	}
	else
	{
		es	= 6.11 * exp( (17.62 * T14) / (243.12 + T14) );
	}


	//-----------------------------------------------------
	// VKR 4.7:		Berechnung des aktuellen Dampfdruckes (e [hPa])
	//
	//		e		= es(T14) * U14 / 100
	//
	//		es(T14)		: Saettigungsdampfdruck 14 Uhr [hPa]
	//		U14			: relative Luftfeuchte 14 Uhr [%] (Dampfdruck in % des Saettigungsdampfdruckes)
	//

	e	= es * U14 / 100.0;


	//-----------------------------------------------------
	// VKR 4.8:		Berechnung der potentiellen Verdunstung nach Haude (ETPhaude [mm/d])
	//
	//		ETPhaude	= f(m) * (es - e)
	//
	//		f(m)		: Monatsfaktor [mm/mbar/d]
	//		e			: aktueller Dampfdruck [hPa]
	//		es			: Saettigungsdampfdruck
	//

	ETP_Haude	= f[Get_Month(Day)] * (es - e);

	return( ETP_Haude );
}

//---------------------------------------------------------
double CDVWK_SoilMoisture::Get_Pi(int Day)
{
	const double Messfehler[]	=
	//	---    Jan    Feb    Mar    Apr    Mai    Jun    Jul    Aug    Sep    Okt    Nov    Dez
	{	0.000, 0.228, 0.236, 0.200, 0.160, 0.120, 0.103, 0.105, 0.103, 0.115, 0.136, 0.162, 0.189	};

	double	P, Pi;

	//-----------------------------------------------------
	CSG_Table_Record	*pRecord;

	if( (pRecord = pClimate->Get_Record(Day)) == NULL )
	{
		return( 0 );
	}

	P		= pRecord->asDouble(0);


	//-----------------------------------------------------
	// VKR 4.11:	Berechnung der korrigierten Niederschlaege (Pi [mm/d])
	//
	//		Pi	= N * kn
	//
	//		N		: Niederschlag gemessen [mm/d]
	//		kn		: Korrekturfaktor ergibt sich aus den prozentualen Me遞ehlern...
	//		Month	: 1=Jan, 2=Feb... 0=NA;
	//

	Pi			= P + P * Messfehler[Get_Month(Day)];

	return( Pi );
}

//---------------------------------------------------------
double CDVWK_SoilMoisture::Get_Wi(double Wi, double Pi, double ETP, double kc, double FK, double PWP)
{
	double	Ri, d, ETPi, ETAi;


	//-----------------------------------------------------
	// VKR 4.9:		Berechnung des Reduktionsfaktors f.Verdunstung (Ri)
	//
	//			        PWP
	//			   1 - -----
	//			        Wi               ETP		Ri = 1	bei Pi > ETPi
	//		Ri	= ------------- + 0,1 * ------;		Ri = 1	bei Ri > 1,0
	//			        PWP              ETPi		Ri = 0	bei Ri < 0,0
	//			   1 - -----
	//			        FK
	//
	//		Pi			: korrigierter Niederschlag [mm/d]
	//		Wi			: Bodenwassergehalt
	//		PWP			: Permanenter Welkepunkt
	//		FK			: Feldkapazitaet
	//		kc			: Koeffizient f.d. Pflz.Bestand
	//

	ETPi	= ETP * kc;

	if( Pi > ETPi )
	{
		Ri	= 1.0;
	}
	else if( ETPi <= 0.0 || Wi <= 0.0 || FK <= 0.0 )
	{
		Ri	= 0.0;
	}
	else
	{
		d	= 1.0 - PWP / FK;

		if( d == 0.0 )
		{
			Ri	= 0.0;
		}
		else
		{
			Ri	= (1.0 - PWP / Wi) / d + 0.1 * ETP / ETPi;

			if( Ri > 1.0 )
			{
				Ri	= 1.0;
			}
			else if( Ri < 0.0 )
			{
				Ri	= 0.0;
			}
		}
	}


	//-----------------------------------------------------
	// VKR 4.13:	Berechnung der tatsaechlichen Verdunstung (ETa [mm/d])
	//
	//		ETai	= Pi - (Pi - (ETPhaude * kc)) * Ri
	//
	//		Pi			: korrigierter Niederschlag [mm/d]
	//		ETPhaude	: potentielle Verdunstung nach Haude [mm/d]
	//		kc			: Bestandeskoeffizient
	//		Ri			: Reduktionsfaktor
	//

	ETAi	= Pi - (Pi - ETPi) * Ri;


	//-----------------------------------------------------
	// VKR 4.14:	Berechnung des Bodenwasservorrats (Wi [mm/d])
	//
	//		Wi+1	= Wi + Pi - ETai - SRi + KR
	//
	//		Pi		: korrigierter Niederschlag [mm/d]
	//		ETai	: tatsaechliche Verdunstung [mm/d]
	//		SRi		: taegliche Sickerwaserrate [mm/d]
	//					SRi	= 0					bei Wi+1 <= FK * ku
	//					SRi	= Wi+1 - FK * ku	bei Wi+1 >  FK * ku;
	//		FK		: Feldkapazitaet [mm]
	//		ku		: moegliche Uebersaettigung ueber FK	:= 1 (findet keine Beruecksichtigung)
	//		KR		: kapillare Aufstiegsrate			:= 0 (findet keine Beruecksichtigung)
	//

	Wi		+= Pi - ETAi;

	if( Wi > FK )
	{
		Wi	= FK;	// -= SRi...
	}
	else if( Wi < PWP )
	{
		Wi	= PWP;
	}


	return( Wi );
}


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

//---------------------------------------------------------
void CDVWK_SoilMoisture::Step_Day(int Day)
{
	int		x, y;
	double	Wi, Pi, ETP, kc, FK, PWP;

	ETP		= Get_ETP_Haude(Day);
	Pi		= Get_Pi(Day);

	for(y=0; y<Get_NY(); y++)
	{
		for(x=0; x<Get_NX(); x++)
		{
			kc		= Get_kc(pLandUse->asInt(x, y), Day);

			FK		= pFK_mm	? pFK_mm	->asDouble(x, y) : FK_mm_Def;
			PWP		= pPWP_mm	? pPWP_mm	->asDouble(x, y) : PWP_mm_Def;

			Wi		= pWi_mm->asDouble(x, y);
			Wi		= Get_Wi(Wi, Pi, ETP, kc, FK, PWP);
			pWi_mm->Set_Value(x, y, Wi);
		}
	}
}