Subversion Repositories public iLand

Rev

Rev 213 | Rev 270 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 213 Rev 214
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
    // setup phenology
-
 
40
    setupPhenology();
-
 
41
-
 
42
    const XmlHelper &xml = GlobalSettings::instance()->settings();
39
    const XmlHelper &xml = GlobalSettings::instance()->settings();
43
    QString tableName = xml.value("model.species.source", "species");
40
    QString tableName = xml.value("model.species.source", "species");
44
    QString readerFile = xml.value("model.species.reader", "reader.bin");
41
    QString readerFile = xml.value("model.species.reader", "reader.bin");
45
    readerFile = GlobalSettings::instance()->path(readerFile, "lip");
42
    readerFile = GlobalSettings::instance()->path(readerFile, "lip");
46
    mReaderStamp.load(readerFile);
43
    mReaderStamp.load(readerFile);
47
44
48
    QSqlQuery query(GlobalSettings::instance()->dbin());
45
    QSqlQuery query(GlobalSettings::instance()->dbin());
49
    mSetupQuery = &query;
46
    mSetupQuery = &query;
50
    QString sql = QString("select * from %1").arg(tableName);
47
    QString sql = QString("select * from %1").arg(tableName);
51
    query.exec(sql);
48
    query.exec(sql);
52
    clear();
49
    clear();
53
    qDebug() << "attempting to load a species set from" << tableName;
50
    qDebug() << "attempting to load a species set from" << tableName;
54
    while (query.next()) {
51
    while (query.next()) {
55
        if (var("active").toInt()==0)
52
        if (var("active").toInt()==0)
56
            continue;
53
            continue;
57
54
58
        Species *s = new Species(this); // create
55
        Species *s = new Species(this); // create
59
        // call setup routine (which calls SpeciesSet::var() to retrieve values
56
        // call setup routine (which calls SpeciesSet::var() to retrieve values
60
        s->setup();
57
        s->setup();
61
58
62
        mSpecies.insert(s->id(), s); // store
59
        mSpecies.insert(s->id(), s); // store
63
        if (s->active())
60
        if (s->active())
64
            mActiveSpecies.append(s);
61
            mActiveSpecies.append(s);
65
    } // while query.next()
62
    } // while query.next()
66
    qDebug() << "loaded" << mSpecies.count() << "active species:";
63
    qDebug() << "loaded" << mSpecies.count() << "active species:";
67
    qDebug() << mSpecies.keys();
64
    qDebug() << mSpecies.keys();
68
65
69
    mSetupQuery = 0;
66
    mSetupQuery = 0;
70
67
71
    // setup nitrogen response
68
    // setup nitrogen response
72
    XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));
69
    XmlHelper resp(xml.node("model.species.nitrogenResponseClasses"));
73
    if (!resp.isValid())
70
    if (!resp.isValid())
74
        throw IException("model.species.nitrogenResponseClasses not present!");
71
        throw IException("model.species.nitrogenResponseClasses not present!");
75
    mNitrogen_1a = resp.valueDouble("class_1_a");
72
    mNitrogen_1a = resp.valueDouble("class_1_a");
76
    mNitrogen_1b = resp.valueDouble("class_1_b");
73
    mNitrogen_1b = resp.valueDouble("class_1_b");
77
    mNitrogen_2a = resp.valueDouble("class_2_a");
74
    mNitrogen_2a = resp.valueDouble("class_2_a");
78
    mNitrogen_2b = resp.valueDouble("class_2_b");
75
    mNitrogen_2b = resp.valueDouble("class_2_b");
79
    mNitrogen_3a = resp.valueDouble("class_3_a");
76
    mNitrogen_3a = resp.valueDouble("class_3_a");
80
    mNitrogen_3b = resp.valueDouble("class_3_b");
77
    mNitrogen_3b = resp.valueDouble("class_3_b");
81
    if (mNitrogen_1a*mNitrogen_1b*mNitrogen_2a*mNitrogen_2b*mNitrogen_3a*mNitrogen_3b == 0)
78
    if (mNitrogen_1a*mNitrogen_1b*mNitrogen_2a*mNitrogen_2b*mNitrogen_3a*mNitrogen_3b == 0)
82
        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)!");
83
80
84
    // setup CO2 response
81
    // setup CO2 response
85
    XmlHelper co2(xml.node("model.species.CO2Response"));
82
    XmlHelper co2(xml.node("model.species.CO2Response"));
86
    mCO2base = co2.valueDouble("baseConcentration");
83
    mCO2base = co2.valueDouble("baseConcentration");
87
    mCO2comp = co2.valueDouble("compensationPoint");
84
    mCO2comp = co2.valueDouble("compensationPoint");
88
    mCO2beta0 = co2.valueDouble("beta0");
85
    mCO2beta0 = co2.valueDouble("beta0");
89
    mCO2p0 = co2.valueDouble("p0");
86
    mCO2p0 = co2.valueDouble("p0");
90
    if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0)
87
    if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0)
91
        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!");
92
89
93
    return mSpecies.count();
90
    return mSpecies.count();
94
91
95
}
92
}
96
93
97
void SpeciesSet::setupPhenology()
-
 
98
{
-
 
99
    mPhenology.clear();
-
 
100
    mPhenology.push_back(Phenology()); // id=0
-
 
101
    XmlHelper xml(GlobalSettings::instance()->settings().node("model.species.phenology"));
-
 
102
    int i=0;
-
 
103
    do {
-
 
104
        QDomElement n = xml.node(QString("type[%1]").arg(i));
-
 
105
        if (n.isNull())
-
 
106
            break;
-
 
107
        i++;
-
 
108
        int id;
-
 
109
        id = n.attribute("id", "-1").toInt();
-
 
110
        if (id<0) throw IException(QString("Error setting up phenology: id invalid\ndump: %1").arg(xml.dump("")));
-
 
111
        xml.setCurrentNode(n);
-
 
112
        Phenology item( id,
-
 
113
                        xml.valueDouble(".vpdMin",0.5), // use relative access to node (".x")
-
 
114
                        xml.valueDouble(".vpdMax", 5),
-
 
115
                        xml.valueDouble(".minDayLength",10),
-
 
116
                        xml.valueDouble(".maxDayLength",11),
-
 
117
                        xml.valueDouble(".tempMin", 2),
-
 
118
                        xml.valueDouble(".tempMax", 9) );
-
 
119
-
 
120
        mPhenology.push_back(item);
-
 
121
    } while(true);
-
 
122
-
 
123
}
-
 
124
94
125
/** retrieves variables from the datasource available during the setup of species.
95
/** retrieves variables from the datasource available during the setup of species.
126
  */
96
  */
127
QVariant SpeciesSet::var(const QString& varName)
97
QVariant SpeciesSet::var(const QString& varName)
128
{
98
{
129
    Q_ASSERT(mSetupQuery!=0);
99
    Q_ASSERT(mSetupQuery!=0);
130
100
131
    int idx = mSetupQuery->record().indexOf(varName);
101
    int idx = mSetupQuery->record().indexOf(varName);
132
    if (idx>=0)
102
    if (idx>=0)
133
        return mSetupQuery->value(idx);
103
        return mSetupQuery->value(idx);
134
    throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));
104
    throw IException(QString("SpeciesSet: variable not set: %1").arg(varName));
135
    //throw IException(QString("load species parameter: field %1 not found!").arg(varName));
105
    //throw IException(QString("load species parameter: field %1 not found!").arg(varName));
136
    // lookup in defaults
106
    // lookup in defaults
137
    //qDebug() << "variable" << varName << "not found - using default.";
107
    //qDebug() << "variable" << varName << "not found - using default.";
138
    //return GlobalSettings::instance()->settingDefaultValue(varName);
108
    //return GlobalSettings::instance()->settingDefaultValue(varName);
139
}
109
}
140
110
141
inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const
111
inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const
142
{
112
{
143
    if (availableNitrogen<=NB)
113
    if (availableNitrogen<=NB)
144
        return 0;
114
        return 0;
145
    double x = 1. - exp(NA * (availableNitrogen-NB));
115
    double x = 1. - exp(NA * (availableNitrogen-NB));
146
    return x;
116
    return x;
147
}
117
}
148
118
149
/// calculate nitrogen response for a given amount of available nitrogen and a respone class
119
/// calculate nitrogen response for a given amount of available nitrogen and a respone class
150
/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
120
/// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
151
double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const
121
double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const
152
{
122
{
153
    double value1, value2, value3;
123
    double value1, value2, value3;
154
    if (responseClass>2.) {
124
    if (responseClass>2.) {
155
        if (responseClass==3.)
125
        if (responseClass==3.)
156
            return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
126
            return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
157
        else {
127
        else {
158
            // interpolate between 2 and 3
128
            // interpolate between 2 and 3
159
            value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
129
            value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
160
            value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
130
            value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b);
161
            return value2 + (responseClass-2)*(value3-value2);
131
            return value2 + (responseClass-2)*(value3-value2);
162
        }
132
        }
163
    }
133
    }
164
    if (responseClass==2)
134
    if (responseClass==2)
165
        return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
135
        return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
166
    if (responseClass==1)
136
    if (responseClass==1)
167
        return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
137
        return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
168
    // last ressort: interpolate between 1 and 2
138
    // last ressort: interpolate between 1 and 2
169
    value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
139
    value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b);
170
    value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
140
    value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b);
171
    return value1 + (responseClass-1)*(value2-value1);
141
    return value1 + (responseClass-1)*(value2-value1);
172
}
142
}
173
143
174
/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
144
/** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
175
    The calculation follows Friedlingsstein 1995 (see also links to equations in code)
145
    The calculation follows Friedlingsstein 1995 (see also links to equations in code)
176
*/
146
*/
177
double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const
147
double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const
178
{
148
{
179
    if (nitrogenResponse==0)
149
    if (nitrogenResponse==0)
180
        return 0.;
150
        return 0.;
181
151
182
    double co2_water = 2. - soilWaterResponse;
152
    double co2_water = 2. - soilWaterResponse;
183
    double beta = mCO2beta0 * co2_water * nitrogenResponse;
153
    double beta = mCO2beta0 * co2_water * nitrogenResponse;
184
154
185
    double r =1. +  M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)
155
    double r =1. +  M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17)
186
156
187
    // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
157
    // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
188
    double deltaC = mCO2base - mCO2comp;
158
    double deltaC = mCO2base - mCO2comp;
189
    double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16
159
    double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16
190
    double K1 = (1. + K2*deltaC) / deltaC;
160
    double K1 = (1. + K2*deltaC) / deltaC;
191
161
192
    double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16
162
    double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16
193
    return response;
163
    return response;
194
164
195
}
165
}
196
166
197
/** return the phenology of the group... */
-
 
198
const Phenology &SpeciesSet::phenology(const int phenologyGroup)
-
 
199
{
-
 
200
    const Phenology &p = mPhenology.at(phenologyGroup);
-
 
201
    if (p.id() == phenologyGroup)
-
 
202
        return p;
-
 
203
    // search...
-
 
204
    for (int i=0;i<mPhenology.count();i++)
-
 
205
        if (mPhenology[i].id()==phenologyGroup)
-
 
206
            return mPhenology[i];
-
 
207
    throw IException(QString("Error at SpeciesSEt::phenology(): invalid group: %1").arg(phenologyGroup));
-
 
208
}
-
 
-
 
167
209
 
168