Subversion Repositories public iLand

Rev

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