Rev 1221 | Only display areas with differences | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 1221 | Rev 1222 | ||
---|---|---|---|
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 | /********************************************************************************************
|
2 | /********************************************************************************************
|
3 | ** iLand - an individual based forest landscape and disturbance model
|
3 | ** iLand - an individual based forest landscape and disturbance model
|
4 | ** http://iland.boku.ac.at
|
4 | ** http://iland.boku.ac.at
|
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl
|
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl
|
6 | **
|
6 | **
|
7 | ** This program is free software: you can redistribute it and/or modify
|
7 | ** This program is free software: you can redistribute it and/or modify
|
8 | ** it under the terms of the GNU General Public License as published by
|
8 | ** it under the terms of the GNU General Public License as published by
|
9 | ** the Free Software Foundation, either version 3 of the License, or
|
9 | ** the Free Software Foundation, either version 3 of the License, or
|
10 | ** (at your option) any later version.
|
10 | ** (at your option) any later version.
|
11 | **
|
11 | **
|
12 | ** This program is distributed in the hope that it will be useful,
|
12 | ** This program is distributed in the hope that it will be useful,
|
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of
|
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of
|
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
15 | ** GNU General Public License for more details.
|
15 | ** GNU General Public License for more details.
|
16 | **
|
16 | **
|
17 | ** You should have received a copy of the GNU General Public License
|
17 | ** You should have received a copy of the GNU General Public License
|
18 | ** along with this program. If not, see <http://www.gnu.org/licenses/>.
|
18 | ** along with this program. If not, see <http://www.gnu.org/licenses/>.
|
19 | ********************************************************************************************/
|
19 | ********************************************************************************************/
|
20 | 20 | ||
21 | #include <QtCore>
|
21 | #include <QtCore>
|
22 | #include <QtSql>
|
22 | #include <QtSql>
|
23 | #include <QVector>
|
23 | #include <QVector>
|
24 | #include "global.h"
|
24 | #include "global.h"
|
25 | #include "globalsettings.h"
|
25 | #include "globalsettings.h"
|
26 | #include "xmlhelper.h"
|
26 | #include "xmlhelper.h"
|
27 | #include "speciesset.h"
|
27 | #include "speciesset.h"
|
28 | #include "species.h"
|
28 | #include "species.h"
|
29 | #include "model.h"
|
29 | #include "model.h"
|
30 | #include "seeddispersal.h"
|
30 | #include "seeddispersal.h"
|
31 | #include "modelsettings.h"
|
31 | #include "modelsettings.h"
|
32 | #include "debugtimer.h"
|
32 | #include "debugtimer.h"
|
33 | 33 | ||
34 | /** @class SpeciesSet
|
34 | /** @class SpeciesSet
|
35 | A SpeciesSet acts as a container for individual Species objects. In iLand, theoretically,
|
35 | A SpeciesSet acts as a container for individual Species objects. In iLand, theoretically,
|
36 | multiple species sets can be used in parallel.
|
36 | multiple species sets can be used in parallel.
|
37 | */
|
37 | */
|
38 | 38 | ||
39 | SpeciesSet::SpeciesSet() |
39 | SpeciesSet::SpeciesSet() |
40 | {
|
40 | {
|
41 | mSetupQuery = 0; |
41 | mSetupQuery = 0; |
42 | }
|
42 | }
|
43 | 43 | ||
44 | SpeciesSet::~SpeciesSet() |
44 | SpeciesSet::~SpeciesSet() |
45 | {
|
45 | {
|
46 | clear(); |
46 | clear(); |
47 | }
|
47 | }
|
48 | 48 | ||
49 | void SpeciesSet::clear() |
49 | void SpeciesSet::clear() |
50 | {
|
50 | {
|
51 | qDeleteAll(mSpecies.values()); |
51 | qDeleteAll(mSpecies.values()); |
52 | qDeleteAll(mSeedDispersal); |
52 | qDeleteAll(mSeedDispersal); |
53 | mSpecies.clear(); |
53 | mSpecies.clear(); |
54 | mActiveSpecies.clear(); |
54 | mActiveSpecies.clear(); |
55 | }
|
55 | }
|
56 | 56 | ||
57 | const Species *SpeciesSet::species(const int &index) |
57 | const Species *SpeciesSet::species(const int &index) |
58 | {
|
58 | {
|
59 | foreach(Species *s, mSpecies) |
59 | foreach(Species *s, mSpecies) |
60 | if (s->index() == index) |
60 | if (s->index() == index) |
61 | return s; |
61 | return s; |
62 | return NULL; |
62 | return NULL; |
63 | }
|
63 | }
|
64 | 64 | ||
65 | /** loads active species from a database table and creates/setups the species.
|
65 | /** loads active species from a database table and creates/setups the species.
|
66 | The function uses the global database-connection.
|
66 | The function uses the global database-connection.
|
67 | */
|
67 | */
|
68 | int SpeciesSet::setup() |
68 | int SpeciesSet::setup() |
69 | {
|
69 | {
|
70 | const XmlHelper &xml = GlobalSettings::instance()->settings(); |
70 | const XmlHelper &xml = GlobalSettings::instance()->settings(); |
71 | QString tableName = xml.value("model.species.source", "species"); |
71 | QString tableName = xml.value("model.species.source", "species"); |
72 | mName = tableName; |
72 | mName = tableName; |
73 | QString readerFile = xml.value("model.species.reader", "reader.bin"); |
73 | QString readerFile = xml.value("model.species.reader", "reader.bin"); |
74 | readerFile = GlobalSettings::instance()->path(readerFile, "lip"); |
74 | readerFile = GlobalSettings::instance()->path(readerFile, "lip"); |
75 | mReaderStamp.load(readerFile); |
75 | mReaderStamp.load(readerFile); |
76 | if (GlobalSettings::instance()->settings().paramValueBool("debugDumpStamps", false) ) |
76 | if (GlobalSettings::instance()->settings().paramValueBool("debugDumpStamps", false) ) |
77 | qDebug() << mReaderStamp.dump(); |
77 | qDebug() << mReaderStamp.dump(); |
78 | 78 | ||
79 | 79 | ||
80 | QSqlQuery query(GlobalSettings::instance()->dbin()); |
80 | QSqlQuery query(GlobalSettings::instance()->dbin()); |
81 | mSetupQuery = &query; |
81 | mSetupQuery = &query; |
82 | QString sql = QString("select * from %1").arg(tableName); |
82 | QString sql = QString("select * from %1").arg(tableName); |
83 | query.exec(sql); |
83 | query.exec(sql); |
84 | if (query.lastError().isValid()){ |
84 | if (query.lastError().isValid()){ |
85 | throw IException(QString("Error loading species set: %1 \n %2").arg(sql, query.lastError().text()) ); |
85 | throw IException(QString("Error loading species set: %1 \n %2").arg(sql, query.lastError().text()) ); |
86 | }
|
86 | }
|
87 | 87 | ||
88 | clear(); |
88 | clear(); |
89 | qDebug() << "attempting to load a species set from" << tableName; |
89 | qDebug() << "attempting to load a species set from" << tableName; |
90 | while (query.next()) { |
90 | while (query.next()) { |
91 | if (var("active").toInt()==0) |
91 | if (var("active").toInt()==0) |
92 | continue; |
92 | continue; |
93 | 93 | ||
94 | Species *s = new Species(this); // create |
94 | Species *s = new Species(this); // create |
95 | // call setup routine (which calls SpeciesSet::var() to retrieve values
|
95 | // call setup routine (which calls SpeciesSet::var() to retrieve values
|
96 | s->setup(); |
96 | s->setup(); |
97 | 97 | ||
98 | mSpecies.insert(s->id(), s); // store |
98 | mSpecies.insert(s->id(), s); // store |
99 | if (s->active()) |
99 | if (s->active()) |
100 | mActiveSpecies.append(s); |
100 | mActiveSpecies.append(s); |
101 | 101 | ||
102 | Expression::addConstant(s->id(), s->index()); |
102 | Expression::addConstant(s->id(), s->index()); |
103 | } // while query.next() |
103 | } // while query.next() |
104 | qDebug() << "loaded" << mSpecies.count() << "active species:"; |
104 | qDebug() << "loaded" << mSpecies.count() << "active species:"; |
105 | qDebug() << "index, id, name"; |
105 | qDebug() << "index, id, name"; |
106 | foreach(const Species *s, mActiveSpecies) |
106 | foreach(const Species *s, mActiveSpecies) |
107 | qDebug() << s->index() << s->id() << s->name(); |
107 | qDebug() << s->index() << s->id() << s->name(); |
108 | 108 | ||
109 | mSetupQuery = 0; |
109 | mSetupQuery = 0; |
110 | 110 | ||
111 | // setup nitrogen response
|
111 | // setup nitrogen response
|
112 | XmlHelper resp(xml.node("model.species.nitrogenResponseClasses")); |
112 | XmlHelper resp(xml.node("model.species.nitrogenResponseClasses")); |
113 | if (!resp.isValid()) |
113 | if (!resp.isValid()) |
114 | throw IException("model.species.nitrogenResponseClasses not present!"); |
114 | throw IException("model.species.nitrogenResponseClasses not present!"); |
115 | mNitrogen_1a = resp.valueDouble("class_1_a"); |
115 | mNitrogen_1a = resp.valueDouble("class_1_a"); |
116 | mNitrogen_1b = resp.valueDouble("class_1_b"); |
116 | mNitrogen_1b = resp.valueDouble("class_1_b"); |
117 | mNitrogen_2a = resp.valueDouble("class_2_a"); |
117 | mNitrogen_2a = resp.valueDouble("class_2_a"); |
118 | mNitrogen_2b = resp.valueDouble("class_2_b"); |
118 | mNitrogen_2b = resp.valueDouble("class_2_b"); |
119 | mNitrogen_3a = resp.valueDouble("class_3_a"); |
119 | mNitrogen_3a = resp.valueDouble("class_3_a"); |
120 | mNitrogen_3b = resp.valueDouble("class_3_b"); |
120 | mNitrogen_3b = resp.valueDouble("class_3_b"); |
121 | if (mNitrogen_1a*mNitrogen_1b*mNitrogen_2a*mNitrogen_2b*mNitrogen_3a*mNitrogen_3b == 0) |
121 | if (mNitrogen_1a*mNitrogen_1b*mNitrogen_2a*mNitrogen_2b*mNitrogen_3a*mNitrogen_3b == 0) |
122 | throw IException("at least one parameter of model.species.nitrogenResponseClasses is not valid (value=0)!"); |
122 | throw IException("at least one parameter of model.species.nitrogenResponseClasses is not valid (value=0)!"); |
123 | 123 | ||
124 | // setup CO2 response
|
124 | // setup CO2 response
|
125 | XmlHelper co2(xml.node("model.species.CO2Response")); |
125 | XmlHelper co2(xml.node("model.species.CO2Response")); |
126 | mCO2base = co2.valueDouble("baseConcentration"); |
126 | mCO2base = co2.valueDouble("baseConcentration"); |
127 | mCO2comp = co2.valueDouble("compensationPoint"); |
127 | mCO2comp = co2.valueDouble("compensationPoint"); |
128 | mCO2beta0 = co2.valueDouble("beta0"); |
128 | mCO2beta0 = co2.valueDouble("beta0"); |
129 | mCO2p0 = co2.valueDouble("p0"); |
129 | mCO2p0 = co2.valueDouble("p0"); |
130 | if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0) |
130 | if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0) |
131 | throw IException("at least one parameter of model.species.CO2Response is not valid!"); |
131 | throw IException("at least one parameter of model.species.CO2Response is not valid!"); |
132 | 132 | ||
133 | // setup Light responses
|
133 | // setup Light responses
|
134 | XmlHelper light(xml.node("model.species.lightResponse")); |
134 | XmlHelper light(xml.node("model.species.lightResponse")); |
135 | mLightResponseTolerant.setAndParse(light.value("shadeTolerant")); |
135 | mLightResponseTolerant.setAndParse(light.value("shadeTolerant")); |
136 | mLightResponseIntolerant.setAndParse(light.value("shadeIntolerant")); |
136 | mLightResponseIntolerant.setAndParse(light.value("shadeIntolerant")); |
137 | mLightResponseTolerant.linearize(0., 1.); |
137 | mLightResponseTolerant.linearize(0., 1.); |
138 | mLightResponseIntolerant.linearize(0., 1.); |
138 | mLightResponseIntolerant.linearize(0., 1.); |
139 | if (mLightResponseTolerant.expression().isEmpty() || mLightResponseIntolerant.expression().isEmpty()) |
139 | if (mLightResponseTolerant.expression().isEmpty() || mLightResponseIntolerant.expression().isEmpty()) |
140 | throw IException("at least one parameter of model.species.lightResponse is empty!"); |
140 | throw IException("at least one parameter of model.species.lightResponse is empty!"); |
141 | // lri-correction
|
141 | // lri-correction
|
142 | mLRICorrection.setAndParse(light.value("LRImodifier","1")); |
142 | mLRICorrection.setAndParse(light.value("LRImodifier","1")); |
143 | // x: LRI, y: relative heigth
|
143 | // x: LRI, y: relative heigth
|
144 | mLRICorrection.linearize2d(0., 1., 0., 1.); |
144 | mLRICorrection.linearize2d(0., 1., 0., 1.); |
145 | 145 | ||
146 | createRandomSpeciesOrder(); |
146 | createRandomSpeciesOrder(); |
147 | return mSpecies.count(); |
147 | return mSpecies.count(); |
148 | 148 | ||
149 | }
|
149 | }
|
150 | 150 | ||
151 | void SpeciesSet::setupRegeneration() |
151 | void SpeciesSet::setupRegeneration() |
152 | {
|
152 | {
|
153 | SeedDispersal::setupExternalSeeds(); |
153 | SeedDispersal::setupExternalSeeds(); |
154 | foreach(Species *s, mActiveSpecies) { |
154 | foreach(Species *s, mActiveSpecies) { |
155 | SeedDispersal *sd = new SeedDispersal(s); |
155 | SeedDispersal *sd = new SeedDispersal(s); |
156 | sd->setup(); // setup memory for the seed map (grid) |
156 | sd->setup(); // setup memory for the seed map (grid) |
157 | s->setSeedDispersal(sd); // establish the link between species and the map |
157 | s->setSeedDispersal(sd); // establish the link between species and the map |
158 | }
|
158 | }
|
159 | SeedDispersal::finalizeExternalSeeds(); |
159 | SeedDispersal::finalizeExternalSeeds(); |
160 | qDebug() << "Setup of seed dispersal maps finished."; |
160 | qDebug() << "Setup of seed dispersal maps finished."; |
161 | }
|
161 | }
|
162 | 162 | ||
163 | void nc_seed_distribution(Species *species) |
163 | void nc_seed_distribution(Species *species) |
164 | {
|
164 | {
|
165 | species->seedDispersal()->execute(); |
165 | species->seedDispersal()->execute(); |
166 | }
|
166 | }
|
167 | 167 | ||
168 | void SpeciesSet::regeneration() |
168 | void SpeciesSet::regeneration() |
169 | {
|
169 | {
|
170 | if (!GlobalSettings::instance()->model()->settings().regenerationEnabled) |
170 | if (!GlobalSettings::instance()->model()->settings().regenerationEnabled) |
171 | return; |
171 | return; |
172 | DebugTimer t("seed dispersal (all species)"); |
172 | DebugTimer t("seed dispersal (all species)"); |
173 | 173 | ||
174 | ThreadRunner runner(mActiveSpecies); // initialize a thread runner object with all active species |
174 | ThreadRunner runner(mActiveSpecies); // initialize a thread runner object with all active species |
175 | runner.run(nc_seed_distribution); |
175 | runner.run(nc_seed_distribution); |
176 | 176 | ||
177 | if (logLevelDebug()) |
177 | if (logLevelDebug()) |
178 | qDebug() << "seed dispersal finished."; |
178 | qDebug() << "seed dispersal finished."; |
179 | }
|
179 | }
|
180 | 180 | ||
181 | /** newYear is called by Model::runYear at the beginning of a year before any growth occurs.
|
181 | /** newYear is called by Model::runYear at the beginning of a year before any growth occurs.
|
182 | This is used for various initializations, e.g. to clear seed dispersal maps
|
182 | This is used for various initializations, e.g. to clear seed dispersal maps
|
183 | */
|
183 | */
|
184 | void SpeciesSet::newYear() |
184 | void SpeciesSet::newYear() |
185 | {
|
185 | {
|
186 | if (!GlobalSettings::instance()->model()->settings().regenerationEnabled) |
186 | if (!GlobalSettings::instance()->model()->settings().regenerationEnabled) |
187 | return; |
187 | return; |
188 | foreach(Species *s, mActiveSpecies) { |
188 | foreach(Species *s, mActiveSpecies) { |
189 | s->newYear(); |
189 | s->newYear(); |
190 | }
|
190 | }
|
191 | }
|
191 | }
|
192 | 192 | ||
193 | /** retrieves variables from the datasource available during the setup of species.
|
193 | /** retrieves variables from the datasource available during the setup of species.
|
194 | */
|
194 | */
|
195 | QVariant SpeciesSet::var(const QString& varName) |
195 | QVariant SpeciesSet::var(const QString& varName) |
196 | {
|
196 | {
|
197 | Q_ASSERT(mSetupQuery!=0); |
197 | Q_ASSERT(mSetupQuery!=0); |
198 | 198 | ||
199 | int idx = mSetupQuery->record().indexOf(varName); |
199 | int idx = mSetupQuery->record().indexOf(varName); |
200 | if (idx>=0) |
200 | if (idx>=0) |
201 | return mSetupQuery->value(idx); |
201 | return mSetupQuery->value(idx); |
202 | throw IException(QString("SpeciesSet: variable not set: %1").arg(varName)); |
202 | throw IException(QString("SpeciesSet: variable not set: %1").arg(varName)); |
203 | //throw IException(QString("load species parameter: field %1 not found!").arg(varName));
|
203 | //throw IException(QString("load species parameter: field %1 not found!").arg(varName));
|
204 | // lookup in defaults
|
204 | // lookup in defaults
|
205 | //qDebug() << "variable" << varName << "not found - using default.";
|
205 | //qDebug() << "variable" << varName << "not found - using default.";
|
206 | //return GlobalSettings::instance()->settingDefaultValue(varName);
|
206 | //return GlobalSettings::instance()->settingDefaultValue(varName);
|
207 | }
|
207 | }
|
208 | 208 | ||
209 | void SpeciesSet::randomSpeciesOrder(QVector<int>::const_iterator &rBegin, QVector<int>::const_iterator &rEnd) |
209 | void SpeciesSet::randomSpeciesOrder(QVector<int>::const_iterator &rBegin, QVector<int>::const_iterator &rEnd) |
210 | {
|
210 | {
|
211 | int iset = irandom(0,mNRandomSets); |
211 | int iset = irandom(0,mNRandomSets); |
212 | rBegin=mRandomSpeciesOrder.begin() + iset * mActiveSpecies.size(); |
212 | rBegin=mRandomSpeciesOrder.begin() + iset * mActiveSpecies.size(); |
213 | rEnd=rBegin+mActiveSpecies.size(); |
213 | rEnd=rBegin+mActiveSpecies.size(); |
214 | }
|
214 | }
|
215 | 215 | ||
216 | //
|
216 | //
|
217 | void SpeciesSet::createRandomSpeciesOrder() |
217 | void SpeciesSet::createRandomSpeciesOrder() |
218 | {
|
218 | {
|
219 | 219 | ||
220 | mRandomSpeciesOrder.clear(); |
220 | mRandomSpeciesOrder.clear(); |
221 | mRandomSpeciesOrder.reserve(mActiveSpecies.size() * mNRandomSets); |
221 | mRandomSpeciesOrder.reserve(mActiveSpecies.size() * mNRandomSets); |
222 | for (int i=0;i<mNRandomSets;++i) { |
222 | for (int i=0;i<mNRandomSets;++i) { |
223 | QList<int> samples; |
223 | QList<int> samples; |
224 | // fill list
|
224 | // fill list
|
225 | foreach (const Species* s, mActiveSpecies) |
225 | foreach (const Species* s, mActiveSpecies) |
226 | samples.push_back(s->index()); |
226 | samples.push_back(s->index()); |
227 | // sample and reduce list
|
227 | // sample and reduce list
|
228 | while (!samples.isEmpty()) { |
228 | while (!samples.isEmpty()) { |
229 | mRandomSpeciesOrder.push_back( samples.takeAt(irandom(0, samples.size())) ); |
229 | mRandomSpeciesOrder.push_back( samples.takeAt(irandom(0, samples.size())) ); |
230 | }
|
230 | }
|
231 | }
|
231 | }
|
232 | }
|
232 | }
|
233 | 233 | ||
234 | inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const |
234 | inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const |
235 | {
|
235 | {
|
236 | if (availableNitrogen<=NB) |
236 | if (availableNitrogen<=NB) |
237 | return 0; |
237 | return 0; |
238 | double x = 1. - exp(NA * (availableNitrogen-NB)); |
238 | double x = 1. - exp(NA * (availableNitrogen-NB)); |
239 | return x; |
239 | return x; |
240 | }
|
240 | }
|
241 | 241 | ||
242 | 242 | ||
243 | 243 | ||
244 | /// calculate nitrogen response for a given amount of available nitrogen and a respone class
|
244 | /// calculate nitrogen response for a given amount of available nitrogen and a respone class
|
245 | /// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
|
245 | /// for fractional values, the response value is interpolated between the fixedly defined classes (1,2,3)
|
246 | double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const |
246 | double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const |
247 | {
|
247 | {
|
248 | double value1, value2, value3; |
248 | double value1, value2, value3; |
249 | if (responseClass>2.) { |
249 | if (responseClass>2.) { |
250 | if (responseClass==3.) |
250 | if (responseClass==3.) |
251 | return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
251 | return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
252 | else { |
252 | else { |
253 | // interpolate between 2 and 3
|
253 | // interpolate between 2 and 3
|
254 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
254 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
255 | value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
255 | value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
256 | return value2 + (responseClass-2)*(value3-value2); |
256 | return value2 + (responseClass-2)*(value3-value2); |
257 | }
|
257 | }
|
258 | }
|
258 | }
|
259 | if (responseClass==2.) |
259 | if (responseClass==2.) |
260 | return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
260 | return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
261 | if (responseClass==1.) |
261 | if (responseClass==1.) |
262 | return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
262 | return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
263 | // last ressort: interpolate between 1 and 2
|
263 | // last ressort: interpolate between 1 and 2
|
264 | value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
264 | value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
265 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
265 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
266 | return value1 + (responseClass-1)*(value2-value1); |
266 | return value1 + (responseClass-1)*(value2-value1); |
267 | }
|
267 | }
|
268 | 268 | ||
269 | /** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
|
269 | /** calculation for the CO2 response for the ambientCO2 for the water- and nitrogen responses given.
|
270 | The calculation follows Friedlingsstein 1995 (see also links to equations in code)
|
270 | The calculation follows Friedlingsstein 1995 (see also links to equations in code)
|
271 | see also: http://iland.boku.ac.at/CO2+response
|
271 | see also: http://iland.boku.ac.at/CO2+response
|
272 | @param ambientCO2 current CO2 concentration (ppm)
|
272 | @param ambientCO2 current CO2 concentration (ppm)
|
273 | @param nitrogenResponse (yearly) nitrogen response of the species
|
273 | @param nitrogenResponse (yearly) nitrogen response of the species
|
274 | @param soilWaterReponse soil water response (mean value for a month)
|
274 | @param soilWaterReponse soil water response (mean value for a month)
|
275 | */
|
275 | */
|
276 | double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const |
276 | double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const |
277 | {
|
277 | {
|
278 | if (nitrogenResponse==0.) |
278 | if (nitrogenResponse==0.) |
279 | return 0.; |
279 | return 0.; |
280 | 280 | ||
281 | double co2_water = 2. - soilWaterResponse; |
281 | double co2_water = 2. - soilWaterResponse; |
282 | double beta = mCO2beta0 * co2_water * nitrogenResponse; |
282 | double beta = mCO2beta0 * co2_water * nitrogenResponse; |
283 | 283 | ||
284 | double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17) |
284 | double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17) |
285 | 285 | ||
286 | // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
|
286 | // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions
|
287 | double deltaC = mCO2base - mCO2comp; |
287 | double deltaC = mCO2base - mCO2comp; |
288 | double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16 |
288 | double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16 |
289 | double K1 = (1. + K2*deltaC) / deltaC; |
289 | double K1 = (1. + K2*deltaC) / deltaC; |
290 | 290 | ||
291 | double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16 |
291 | double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16 |
292 | return response; |
292 | return response; |
293 | 293 | ||
294 | }
|
294 | }
|
295 | 295 | ||
296 | /** calculates the lightResponse based on a value for LRI and the species lightResponseClass.
|
296 | /** calculates the lightResponse based on a value for LRI and the species lightResponseClass.
|
297 | LightResponse is classified from 1 (very shade inolerant) and 5 (very shade tolerant) and interpolated for values between 1 and 5.
|
297 | LightResponse is classified from 1 (very shade inolerant) and 5 (very shade tolerant) and interpolated for values between 1 and 5.
|
298 | Returns a value between 0..1
|
298 | Returns a value between 0..1
|
299 | @sa http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth */
|
299 | @sa http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth */
|
300 | double SpeciesSet::lightResponse(const double lightResourceIndex, const double lightResponseClass) const |
300 | double SpeciesSet::lightResponse(const double lightResourceIndex, const double lightResponseClass) const |
301 | {
|
301 | {
|
302 | double low = mLightResponseIntolerant.calculate(lightResourceIndex); |
302 | double low = mLightResponseIntolerant.calculate(lightResourceIndex); |
303 | double high = mLightResponseTolerant.calculate(lightResourceIndex); |
303 | double high = mLightResponseTolerant.calculate(lightResourceIndex); |
304 | double result = low + 0.25*(lightResponseClass-1.)*(high-low); |
304 | double result = low + 0.25*(lightResponseClass-1.)*(high-low); |
305 | return limit(result, 0., 1.); |
305 | return limit(result, 0., 1.); |
306 | 306 | ||
307 | }
|
307 | }
|
308 | 308 | ||
309 | 309 | ||
310 | 310 | ||
311 | 311 |