WebSVN - public iLand - Diff - Rev 270 and 274


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/speciesset.cpp':
#include <QtCore>
#include <QtSql>
#include "global.h"
#include "xmlhelper.h"
#include "speciesset.h"
#include "species.h"

SpeciesSet::SpeciesSet()
{
    mSetupQuery = 0;
}

SpeciesSet::~SpeciesSet()
{
   clear();
}

void SpeciesSet::clear()
{
    qDeleteAll(mSpecies.values());
    mSpecies.clear();
    mActiveSpecies.clear();
}

const Species *SpeciesSet::species(const int &index)
{
    foreach(Species *s, mSpecies)
        if (s->index() == index)
            return s;
    return NULL;
}

/** loads active species from a database table and creates/setups the species.
    The function uses the global database-connection.
  */
int SpeciesSet::setup()
{
    const XmlHelper &xml = GlobalSettings::instance()->settings();
    QString tableName = xml.value("model.species.source", "species");
    QString readerFile = xml.value("model.species.reader", "reader.bin");
    readerFile = GlobalSettings::instance()->path(readerFile, "lip");
    mReaderStamp.load(readerFile);

    QSqlQuery query(GlobalSettings::instance()->dbin());
    mSetupQuery = &query;
    QString sql = QString("select * from %1").arg(tableName);
    query.exec(sql);
    if (query.lastError().isValid()){
        throw IException(QString("Error loading species set: %1 \n %2").arg(sql, query.lastError().text()) );
    }

    clear();
    qDebug() << "attempting to load a species set from" << tableName;
    while (query.next()) {
        if (var("active").toInt()==0)
            continue;

        Species *s = new Species(this); // create
        // call setup routine (which calls SpeciesSet::var() to retrieve values
        s->setup();

        mSpecies.insert(s->id(), s); // store
        if (s->active())
            mActiveSpecies.append(s);
    } // while query.next()
    qDebug() << "loaded" << mSpecies.count() << "active species:";
    qDebug() << mSpecies.keys();

    mSetupQuery = 0;

    // setup nitrogen response
    XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));
    if (!resp.isValid())
        throw IException("model.species.nitrogenResponseClasses not present!");
    mNitrogen_1a = resp.valueDouble("class_1_a");
    mNitrogen_1b = resp.valueDouble("class_1_b");
    mNitrogen_2a = resp.valueDouble("class_2_a");
    mNitrogen_2b = resp.valueDouble("class_2_b");
    mNitrogen_3a = resp.valueDouble("class_3_a");
    mNitrogen_3b = resp.valueDouble("class_3_b");
    if (mNitrogen_1a*mNitrogen_1b*mNitrogen_2a*mNitrogen_2b*mNitrogen_3a*mNitrogen_3b == 0)
        throw IException("at least one parameter of model.species.nitrogenResponseClasses is not valid (value=0)!");

    // setup CO2 response
    XmlHelper co2(xml.node("model.species.CO2Response"));
    mCO2base = co2.valueDouble("baseConcentration");
    mCO2comp = co2.valueDouble("compensationPoint");
    mCO2beta0 = co2.valueDouble("beta0");
    mCO2p0 = co2.valueDouble("p0");
    if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0)
        throw IException("at least one parameter of model.species.CO2Response is not valid!");

    // setup Light responses
    XmlHelper light(xml.node("model.species.lightResponse"));
    mLightResponseTolerant.setAndParse(light.value("shadeTolerant"));
    mLightResponseIntolerant.setAndParse(light.value("shadeIntolerant"));
    if (mLightResponseTolerant.expression().isEmpty() || mLightResponseIntolerant.expression().isEmpty())
        throw IException("at least one parameter of model.species.lightResponse is empty!");

    return mSpecies.count();

}


/** retrieves variables from the datasource available during the setup of species.
  */
QVariant SpeciesSet::var(const QString& varName)
{
    Q_ASSERT(mSetupQuery!=0);

    int idx = mSetupQuery->record().indexOf(varName);
    if (idx>=0)
        return mSetupQuery->value(idx);
    throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));
    //throw IException(QString("load species parameter: field %1 not found!").arg(varName));
    // lookup in defaults
    //qDebug() << "variable" << varName << "not found - using default.";
    //return GlobalSettings::instance()->settingDefaultValue(varName);
}

inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const
{
    if (availableNitrogen<=NB)
        return 0;
    double x = 1. - exp(NA * (availableNitrogen-NB));
    return x;
}

/// calculate nitrogen response for a given amount of available nitrogen and a respone class
/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const
{
    double value1, value2, value3;
    if (responseClass>2.) {
        if (responseClass==3.)
            return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
        else {
            // interpolate between 2 and 3
            value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
            value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
            return value2 + (responseClass-2)*(value3-value2);
        }
    }
    if (responseClass==2)
        return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
    if (responseClass==1)
        return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
    // last ressort: interpolate between 1 and 2
    value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
    value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
    return value1 + (responseClass-1)*(value2-value1);
}

/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
    The calculation follows Friedlingsstein 1995 (see also links to equations in code)
*/
double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const
{
    if (nitrogenResponse==0)
        return 0.;

    double co2_water = 2. - soilWaterResponse;
    double beta = mCO2beta0 * co2_water * nitrogenResponse;

    double r =1. +  M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)

    // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
    double deltaC = mCO2base - mCO2comp;
    double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16
    double K1 = (1. + K2*deltaC) / deltaC;

    double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16
    return response;

}

/** calculates the lightResponse based on a value for LRI and the species lightResponseClass.
    LightResponse is classified from 1 (very shade inolerant) and 5 (very shade tolerant) and interpolated for values between 1 and 5.
    Returns a value between 0..1 */
double SpeciesSet::lightResponse(const double lightResourceIndex, const double lightResponseClass)
{
    QMutexLocker l(&mMutex); // serialize access to calculations
    double low = mLightResponseIntolerant.calculate(lightResourceIndex);
    double high = mLightResponseTolerant.calculate(lightResourceIndex);
    double result = low + 0.25*(lightResponseClass-1.)*(high-low);
    return limit(result, 0., 1.);

}



 

Rev 270	Rev 274
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/speciesset.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/speciesset.cpp':
2	#include <QtCore>	2	#include <QtCore>
3	#include <QtSql>	3	#include <QtSql>
4	#include "global.h"	4	#include "global.h"
5	#include "xmlhelper.h"	5	#include "xmlhelper.h"
6	#include "speciesset.h"	6	#include "speciesset.h"
7	#include "species.h"	7	#include "species.h"
8		8
9	SpeciesSet::SpeciesSet()	9	SpeciesSet::SpeciesSet()
10	{	10	{
11	mSetupQuery = 0;	11	mSetupQuery = 0;
12	}	12	}
13		13
14	SpeciesSet::~SpeciesSet()	14	SpeciesSet::~SpeciesSet()
15	{	15	{
16	clear();	16	clear();
17	}	17	}
18		18
19	void SpeciesSet::clear()	19	void SpeciesSet::clear()
20	{	20	{
21	qDeleteAll(mSpecies.values());	21	qDeleteAll(mSpecies.values());
22	mSpecies.clear();	22	mSpecies.clear();
23	mActiveSpecies.clear();	23	mActiveSpecies.clear();
24	}	24	}
25		25
26	const Species *SpeciesSet::species(const int &index)	26	const Species *SpeciesSet::species(const int &index)
27	{	27	{
28	foreach(Species *s, mSpecies)	28	foreach(Species *s, mSpecies)
29	if (s->index() == index)	29	if (s->index() == index)
30	return s;	30	return s;
31	return NULL;	31	return NULL;
32	}	32	}
33		33
34	/** loads active species from a database table and creates/setups the species.	34	/** loads active species from a database table and creates/setups the species.
35	The function uses the global database-connection.	35	The function uses the global database-connection.
36	*/	36	*/
37	int SpeciesSet::setup()	37	int SpeciesSet::setup()
38	{	38	{
39	const XmlHelper &xml = GlobalSettings::instance()->settings();	39	const XmlHelper &xml = GlobalSettings::instance()->settings();
40	QString tableName = xml.value("model.species.source", "species");	40	QString tableName = xml.value("model.species.source", "species");
41	QString readerFile = xml.value("model.species.reader", "reader.bin");	41	QString readerFile = xml.value("model.species.reader", "reader.bin");
42	readerFile = GlobalSettings::instance()->path(readerFile, "lip");	42	readerFile = GlobalSettings::instance()->path(readerFile, "lip");
43	mReaderStamp.load(readerFile);	43	mReaderStamp.load(readerFile);
44		44
45	QSqlQuery query(GlobalSettings::instance()->dbin());	45	QSqlQuery query(GlobalSettings::instance()->dbin());
46	mSetupQuery = &query;	46	mSetupQuery = &query;
47	QString sql = QString("select * from %1").arg(tableName);	47	QString sql = QString("select * from %1").arg(tableName);
48	query.exec(sql);	48	query.exec(sql);
49	if (query.lastError().isValid()){	49	if (query.lastError().isValid()){
50	throw IException(QString("Error loading species set: %1 \n %2").arg(sql, query.lastError().text()) );	50	throw IException(QString("Error loading species set: %1 \n %2").arg(sql, query.lastError().text()) );
51	}	51	}
52		52
53	clear();	53	clear();
54	qDebug() << "attempting to load a species set from" << tableName;	54	qDebug() << "attempting to load a species set from" << tableName;
55	while (query.next()) {	55	while (query.next()) {
56	if (var("active").toInt()==0)	56	if (var("active").toInt()==0)
57	continue;	57	continue;
58		58
59	Species *s = new Species(this); // create	59	Species *s = new Species(this); // create
60	// call setup routine (which calls SpeciesSet::var() to retrieve values	60	// call setup routine (which calls SpeciesSet::var() to retrieve values
61	s->setup();	61	s->setup();
62		62
63	mSpecies.insert(s->id(), s); // store	63	mSpecies.insert(s->id(), s); // store
64	if (s->active())	64	if (s->active())
65	mActiveSpecies.append(s);	65	mActiveSpecies.append(s);
66	} // while query.next()	66	} // while query.next()
67	qDebug() << "loaded" << mSpecies.count() << "active species:";	67	qDebug() << "loaded" << mSpecies.count() << "active species:";
68	qDebug() << mSpecies.keys();	68	qDebug() << mSpecies.keys();
69		69
70	mSetupQuery = 0;	70	mSetupQuery = 0;
71		71
72	// setup nitrogen response	72	// setup nitrogen response
73	XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));	73	XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));
74	if (!resp.isValid())	74	if (!resp.isValid())
75	throw IException("model.species.nitrogenResponseClasses not present!");	75	throw IException("model.species.nitrogenResponseClasses not present!");
76	mNitrogen_1a = resp.valueDouble("class_1_a");	76	mNitrogen_1a = resp.valueDouble("class_1_a");
77	mNitrogen_1b = resp.valueDouble("class_1_b");	77	mNitrogen_1b = resp.valueDouble("class_1_b");
78	mNitrogen_2a = resp.valueDouble("class_2_a");	78	mNitrogen_2a = resp.valueDouble("class_2_a");
79	mNitrogen_2b = resp.valueDouble("class_2_b");	79	mNitrogen_2b = resp.valueDouble("class_2_b");
80	mNitrogen_3a = resp.valueDouble("class_3_a");	80	mNitrogen_3a = resp.valueDouble("class_3_a");
81	mNitrogen_3b = resp.valueDouble("class_3_b");	81	mNitrogen_3b = resp.valueDouble("class_3_b");
82	if (mNitrogen_1amNitrogen_1bmNitrogen_2amNitrogen_2bmNitrogen_3a*mNitrogen_3b == 0)	82	if (mNitrogen_1amNitrogen_1bmNitrogen_2amNitrogen_2bmNitrogen_3a*mNitrogen_3b == 0)
83	throw IException("at least one parameter of model.species.nitrogenResponseClasses is not valid (value=0)!");	83	throw IException("at least one parameter of model.species.nitrogenResponseClasses is not valid (value=0)!");
84		84
85	// setup CO2 response	85	// setup CO2 response
86	XmlHelper co2(xml.node("model.species.CO2Response"));	86	XmlHelper co2(xml.node("model.species.CO2Response"));
87	mCO2base = co2.valueDouble("baseConcentration");	87	mCO2base = co2.valueDouble("baseConcentration");
88	mCO2comp = co2.valueDouble("compensationPoint");	88	mCO2comp = co2.valueDouble("compensationPoint");
89	mCO2beta0 = co2.valueDouble("beta0");	89	mCO2beta0 = co2.valueDouble("beta0");
90	mCO2p0 = co2.valueDouble("p0");	90	mCO2p0 = co2.valueDouble("p0");
91	if (mCO2basemCO2comp(mCO2base-mCO2comp)mCO2beta0mCO2p0==0)	91	if (mCO2basemCO2comp(mCO2base-mCO2comp)mCO2beta0mCO2p0==0)
92	throw IException("at least one parameter of model.species.CO2Response is not valid!");	92	throw IException("at least one parameter of model.species.CO2Response is not valid!");
-		93
-		94	// setup Light responses
-		95	XmlHelper light(xml.node("model.species.lightResponse"));
-		96	mLightResponseTolerant.setAndParse(light.value("shadeTolerant"));
-		97	mLightResponseIntolerant.setAndParse(light.value("shadeIntolerant"));
-		98	if (mLightResponseTolerant.expression().isEmpty() \|\| mLightResponseIntolerant.expression().isEmpty())
-		99	throw IException("at least one parameter of model.species.lightResponse is empty!");
93		100
94	return mSpecies.count();	101	return mSpecies.count();
95		102
96	}	103	}
97		104
98		105
99	/** retrieves variables from the datasource available during the setup of species.	106	/** retrieves variables from the datasource available during the setup of species.
100	*/	107	*/
101	QVariant SpeciesSet::var(const QString& varName)	108	QVariant SpeciesSet::var(const QString& varName)
102	{	109	{
103	Q_ASSERT(mSetupQuery!=0);	110	Q_ASSERT(mSetupQuery!=0);
104		111
105	int idx = mSetupQuery->record().indexOf(varName);	112	int idx = mSetupQuery->record().indexOf(varName);
106	if (idx>=0)	113	if (idx>=0)
107	return mSetupQuery->value(idx);	114	return mSetupQuery->value(idx);
108	throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));	115	throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));
109	//throw IException(QString("load species parameter: field %1 not found!").arg(varName));	116	//throw IException(QString("load species parameter: field %1 not found!").arg(varName));
110	// lookup in defaults	117	// lookup in defaults
111	//qDebug() << "variable" << varName << "not found - using default.";	118	//qDebug() << "variable" << varName << "not found - using default.";
112	//return GlobalSettings::instance()->settingDefaultValue(varName);	119	//return GlobalSettings::instance()->settingDefaultValue(varName);
113	}	120	}
114		121
115	inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const	122	inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const
116	{	123	{
117	if (availableNitrogen<=NB)	124	if (availableNitrogen<=NB)
118	return 0;	125	return 0;
119	double x = 1. - exp(NA * (availableNitrogen-NB));	126	double x = 1. - exp(NA * (availableNitrogen-NB));
120	return x;	127	return x;
121	}	128	}
122		129
123	/// calculate nitrogen response for a given amount of available nitrogen and a respone class	130	/// calculate nitrogen response for a given amount of available nitrogen and a respone class
124	/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)	131	/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
125	double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const	132	double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const
126	{	133	{
127	double value1, value2, value3;	134	double value1, value2, value3;
128	if (responseClass>2.) {	135	if (responseClass>2.) {
129	if (responseClass==3.)	136	if (responseClass==3.)
130	return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);	137	return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
131	else {	138	else {
132	// interpolate between 2 and 3	139	// interpolate between 2 and 3
133	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);	140	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
134	value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);	141	value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
135	return value2 + (responseClass-2)*(value3-value2);	142	return value2 + (responseClass-2)*(value3-value2);
136	}	143	}
137	}	144	}
138	if (responseClass==2)	145	if (responseClass==2)
139	return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);	146	return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
140	if (responseClass==1)	147	if (responseClass==1)
141	return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);	148	return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
142	// last ressort: interpolate between 1 and 2	149	// last ressort: interpolate between 1 and 2
143	value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);	150	value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
144	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);	151	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
145	return value1 + (responseClass-1)*(value2-value1);	152	return value1 + (responseClass-1)*(value2-value1);
146	}	153	}
147		154
148	/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.	155	/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
149	The calculation follows Friedlingsstein 1995 (see also links to equations in code)	156	The calculation follows Friedlingsstein 1995 (see also links to equations in code)
150	*/	157	*/
151	double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const	158	double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const
152	{	159	{
153	if (nitrogenResponse==0)	160	if (nitrogenResponse==0)
154	return 0.;	161	return 0.;
155		162
156	double co2_water = 2. - soilWaterResponse;	163	double co2_water = 2. - soilWaterResponse;
157	double beta = mCO2beta0 * co2_water * nitrogenResponse;	164	double beta = mCO2beta0 * co2_water * nitrogenResponse;
158		165
159	double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)	166	double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)
160		167
161	// fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions	168	// fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
162	double deltaC = mCO2base - mCO2comp;	169	double deltaC = mCO2base - mCO2comp;
163	double K2 = ((2mCO2base - mCO2comp) - rdeltaC ) / ((r-1.)deltaC(2*mCO2base - mCO2comp)); // Eq. 16	170	double K2 = ((2mCO2base - mCO2comp) - rdeltaC ) / ((r-1.)deltaC(2*mCO2base - mCO2comp)); // Eq. 16
164	double K1 = (1. + K2*deltaC) / deltaC;	171	double K1 = (1. + K2*deltaC) / deltaC;
165		172
166	double response = mCO2p0 * K1(ambientCO2 - mCO2comp) / (1 + K2(ambientCO2-mCO2comp)); // Eq. 16	173	double response = mCO2p0 * K1(ambientCO2 - mCO2comp) / (1 + K2(ambientCO2-mCO2comp)); // Eq. 16
167	return response;	174	return response;
168		175
169	}	176	}
-		177
-		178	/** calculates the lightResponse based on a value for LRI and the species lightResponseClass.
-		179	LightResponse is classified from 1 (very shade inolerant) and 5 (very shade tolerant) and interpolated for values between 1 and 5.
-		180	Returns a value between 0..1 */
-		181	double SpeciesSet::lightResponse(const double lightResourceIndex, const double lightResponseClass)
-		182	{
-		183	QMutexLocker l(&mMutex); // serialize access to calculations
-		184	double low = mLightResponseIntolerant.calculate(lightResourceIndex);
-		185	double high = mLightResponseTolerant.calculate(lightResourceIndex);
-		186	double result = low + 0.25(lightResponseClass-1.)(high-low);
-		187	return limit(result, 0., 1.);
-		188
-		189	}
-		190
170		191
171		192
172		193

Subversion Repositories public iLand

(root)/src/core/speciesset.cpp @ 1222 - Rev 270 → 274