WebSVN - public iLand - Diff - Rev 214 and 270


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!");

    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;

}


 

Rev 214	Rev 270
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()){
-		50	throw IException(QString("Error loading species set: %1 \n %2").arg(sql, query.lastError().text()) );
-		51	}
-		52
49	clear();	53	clear();
50	qDebug() << "attempting to load a species set from" << tableName;	54	qDebug() << "attempting to load a species set from" << tableName;
51	while (query.next()) {	55	while (query.next()) {
52	if (var("active").toInt()==0)	56	if (var("active").toInt()==0)
53	continue;	57	continue;
54		58
55	Species *s = new Species(this); // create	59	Species *s = new Species(this); // create
56	// call setup routine (which calls SpeciesSet::var() to retrieve values	60	// call setup routine (which calls SpeciesSet::var() to retrieve values
57	s->setup();	61	s->setup();
58		62
59	mSpecies.insert(s->id(), s); // store	63	mSpecies.insert(s->id(), s); // store
60	if (s->active())	64	if (s->active())
61	mActiveSpecies.append(s);	65	mActiveSpecies.append(s);
62	} // while query.next()	66	} // while query.next()
63	qDebug() << "loaded" << mSpecies.count() << "active species:";	67	qDebug() << "loaded" << mSpecies.count() << "active species:";
64	qDebug() << mSpecies.keys();	68	qDebug() << mSpecies.keys();
65		69
66	mSetupQuery = 0;	70	mSetupQuery = 0;
67		71
68	// setup nitrogen response	72	// setup nitrogen response
69	XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));	73	XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));
70	if (!resp.isValid())	74	if (!resp.isValid())
71	throw IException("model.species.nitrogenResponseClasses not present!");	75	throw IException("model.species.nitrogenResponseClasses not present!");
72	mNitrogen_1a = resp.valueDouble("class_1_a");	76	mNitrogen_1a = resp.valueDouble("class_1_a");
73	mNitrogen_1b = resp.valueDouble("class_1_b");	77	mNitrogen_1b = resp.valueDouble("class_1_b");
74	mNitrogen_2a = resp.valueDouble("class_2_a");	78	mNitrogen_2a = resp.valueDouble("class_2_a");
75	mNitrogen_2b = resp.valueDouble("class_2_b");	79	mNitrogen_2b = resp.valueDouble("class_2_b");
76	mNitrogen_3a = resp.valueDouble("class_3_a");	80	mNitrogen_3a = resp.valueDouble("class_3_a");
77	mNitrogen_3b = resp.valueDouble("class_3_b");	81	mNitrogen_3b = resp.valueDouble("class_3_b");
78	if (mNitrogen_1amNitrogen_1bmNitrogen_2amNitrogen_2bmNitrogen_3a*mNitrogen_3b == 0)	82	if (mNitrogen_1amNitrogen_1bmNitrogen_2amNitrogen_2bmNitrogen_3a*mNitrogen_3b == 0)
79	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)!");
80		84
81	// setup CO2 response	85	// setup CO2 response
82	XmlHelper co2(xml.node("model.species.CO2Response"));	86	XmlHelper co2(xml.node("model.species.CO2Response"));
83	mCO2base = co2.valueDouble("baseConcentration");	87	mCO2base = co2.valueDouble("baseConcentration");
84	mCO2comp = co2.valueDouble("compensationPoint");	88	mCO2comp = co2.valueDouble("compensationPoint");
85	mCO2beta0 = co2.valueDouble("beta0");	89	mCO2beta0 = co2.valueDouble("beta0");
86	mCO2p0 = co2.valueDouble("p0");	90	mCO2p0 = co2.valueDouble("p0");
87	if (mCO2basemCO2comp(mCO2base-mCO2comp)mCO2beta0mCO2p0==0)	91	if (mCO2basemCO2comp(mCO2base-mCO2comp)mCO2beta0mCO2p0==0)
88	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!");
89		93
90	return mSpecies.count();	94	return mSpecies.count();
91		95
92	}	96	}
93		97
94		98
95	/** retrieves variables from the datasource available during the setup of species.	99	/** retrieves variables from the datasource available during the setup of species.
96	*/	100	*/
97	QVariant SpeciesSet::var(const QString& varName)	101	QVariant SpeciesSet::var(const QString& varName)
98	{	102	{
99	Q_ASSERT(mSetupQuery!=0);	103	Q_ASSERT(mSetupQuery!=0);
100		104
101	int idx = mSetupQuery->record().indexOf(varName);	105	int idx = mSetupQuery->record().indexOf(varName);
102	if (idx>=0)	106	if (idx>=0)
103	return mSetupQuery->value(idx);	107	return mSetupQuery->value(idx);
104	throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));	108	throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));
105	//throw IException(QString("load species parameter: field %1 not found!").arg(varName));	109	//throw IException(QString("load species parameter: field %1 not found!").arg(varName));
106	// lookup in defaults	110	// lookup in defaults
107	//qDebug() << "variable" << varName << "not found - using default.";	111	//qDebug() << "variable" << varName << "not found - using default.";
108	//return GlobalSettings::instance()->settingDefaultValue(varName);	112	//return GlobalSettings::instance()->settingDefaultValue(varName);
109	}	113	}
110		114
111	inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const	115	inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const
112	{	116	{
113	if (availableNitrogen<=NB)	117	if (availableNitrogen<=NB)
114	return 0;	118	return 0;
115	double x = 1. - exp(NA * (availableNitrogen-NB));	119	double x = 1. - exp(NA * (availableNitrogen-NB));
116	return x;	120	return x;
117	}	121	}
118		122
119	/// calculate nitrogen response for a given amount of available nitrogen and a respone class	123	/// calculate nitrogen response for a given amount of available nitrogen and a respone class
120	/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)	124	/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
121	double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const	125	double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const
122	{	126	{
123	double value1, value2, value3;	127	double value1, value2, value3;
124	if (responseClass>2.) {	128	if (responseClass>2.) {
125	if (responseClass==3.)	129	if (responseClass==3.)
126	return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);	130	return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
127	else {	131	else {
128	// interpolate between 2 and 3	132	// interpolate between 2 and 3
129	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);	133	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
130	value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);	134	value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
131	return value2 + (responseClass-2)*(value3-value2);	135	return value2 + (responseClass-2)*(value3-value2);
132	}	136	}
133	}	137	}
134	if (responseClass==2)	138	if (responseClass==2)
135	return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);	139	return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
136	if (responseClass==1)	140	if (responseClass==1)
137	return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);	141	return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
138	// last ressort: interpolate between 1 and 2	142	// last ressort: interpolate between 1 and 2
139	value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);	143	value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
140	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);	144	value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
141	return value1 + (responseClass-1)*(value2-value1);	145	return value1 + (responseClass-1)*(value2-value1);
142	}	146	}
143		147
144	/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.	148	/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
145	The calculation follows Friedlingsstein 1995 (see also links to equations in code)	149	The calculation follows Friedlingsstein 1995 (see also links to equations in code)
146	*/	150	*/
147	double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const	151	double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const
148	{	152	{
149	if (nitrogenResponse==0)	153	if (nitrogenResponse==0)
150	return 0.;	154	return 0.;
151		155
152	double co2_water = 2. - soilWaterResponse;	156	double co2_water = 2. - soilWaterResponse;
153	double beta = mCO2beta0 * co2_water * nitrogenResponse;	157	double beta = mCO2beta0 * co2_water * nitrogenResponse;
154		158
155	double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)	159	double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)
156		160
157	// fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions	161	// fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
158	double deltaC = mCO2base - mCO2comp;	162	double deltaC = mCO2base - mCO2comp;
159	double K2 = ((2mCO2base - mCO2comp) - rdeltaC ) / ((r-1.)deltaC(2*mCO2base - mCO2comp)); // Eq. 16	163	double K2 = ((2mCO2base - mCO2comp) - rdeltaC ) / ((r-1.)deltaC(2*mCO2base - mCO2comp)); // Eq. 16
160	double K1 = (1. + K2*deltaC) / deltaC;	164	double K1 = (1. + K2*deltaC) / deltaC;
161		165
162	double response = mCO2p0 * K1(ambientCO2 - mCO2comp) / (1 + K2(ambientCO2-mCO2comp)); // Eq. 16	166	double response = mCO2p0 * K1(ambientCO2 - mCO2comp) / (1 + K2(ambientCO2-mCO2comp)); // Eq. 16
163	return response;	167	return response;
164		168
165	}	169	}
166		170
167		171
168		172

Subversion Repositories public iLand

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