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