WebSVN - public iLand - Diff - Rev 209 and 210


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

Subversion Repositories public iLand

(root)/src/core/speciesset.cpp @ 1222 - Rev 209 → 210