Rev 209 | Rev 211 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 1 | |||
| 91 | Werner | 2 | #include <QtCore> |
| 3 | #include <QtSql> |
||
| 4 | #include "global.h" |
||
| 102 | Werner | 5 | #include "xmlhelper.h" |
| 90 | Werner | 6 | #include "speciesset.h" |
| 91 | Werner | 7 | #include "species.h" |
| 90 | Werner | 8 | |
| 9 | SpeciesSet::SpeciesSet() |
||
| 10 | { |
||
| 91 | Werner | 11 | mSetupQuery = 0; |
| 90 | Werner | 12 | } |
| 91 | Werner | 13 | |
| 14 | SpeciesSet::~SpeciesSet() |
||
| 15 | { |
||
| 16 | clear(); |
||
| 17 | } |
||
| 18 | |||
| 19 | void SpeciesSet::clear() |
||
| 20 | { |
||
| 21 | qDeleteAll(mSpecies.values()); |
||
| 22 | mSpecies.clear(); |
||
| 179 | werner | 23 | mActiveSpecies.clear(); |
| 91 | Werner | 24 | } |
| 25 | |||
| 111 | Werner | 26 | const Species *SpeciesSet::species(const int &index) |
| 27 | { |
||
| 28 | foreach(Species *s, mSpecies) |
||
| 29 | if (s->index() == index) |
||
| 30 | return s; |
||
| 31 | return NULL; |
||
| 32 | } |
||
| 91 | Werner | 33 | |
| 34 | /** loads active species from a database table and creates/setups the species. |
||
| 35 | The function uses the global database-connection. |
||
| 36 | */ |
||
| 102 | Werner | 37 | int SpeciesSet::setup() |
| 91 | Werner | 38 | { |
| 102 | Werner | 39 | const XmlHelper &xml = GlobalSettings::instance()->settings(); |
| 191 | werner | 40 | QString tableName = xml.value("model.species.source", "species"); |
| 41 | QString readerFile = xml.value("model.species.reader", "reader.bin"); |
||
| 102 | Werner | 42 | readerFile = GlobalSettings::instance()->path(readerFile, "lip"); |
| 43 | mReaderStamp.load(readerFile); |
||
| 44 | |||
| 91 | Werner | 45 | QSqlQuery query(GlobalSettings::instance()->dbin()); |
| 46 | mSetupQuery = &query; |
||
| 47 | QString sql = QString("select * from %1").arg(tableName); |
||
| 48 | query.exec(sql); |
||
| 49 | clear(); |
||
| 50 | qDebug() << "attempting to load a species set from" << tableName; |
||
| 51 | while (query.next()) { |
||
| 52 | if (var("active").toInt()==0) |
||
| 53 | continue; |
||
| 54 | |||
| 55 | Species *s = new Species(this); // create |
||
| 99 | Werner | 56 | // call setup routine (which calls SpeciesSet::var() to retrieve values |
| 91 | Werner | 57 | s->setup(); |
| 58 | |||
| 59 | mSpecies.insert(s->id(), s); // store |
||
| 179 | werner | 60 | if (s->active()) |
| 61 | mActiveSpecies.append(s); |
||
| 91 | Werner | 62 | } // while query.next() |
| 63 | qDebug() << "loaded" << mSpecies.count() << "active species:"; |
||
| 64 | qDebug() << mSpecies.keys(); |
||
| 65 | |||
| 66 | mSetupQuery = 0; |
||
| 209 | werner | 67 | |
| 68 | // setup nitrogen response |
||
| 69 | XmlHelper resp(xml.node("model.species.nitrogenResponseClasses")); |
||
| 70 | if (!resp.isValid()) |
||
| 71 | throw IException("model.species.nitrogenResponseClasses not present!"); |
||
| 72 | mNitrogen_1a = resp.valueDouble("class_1_a"); |
||
| 73 | mNitrogen_1b = resp.valueDouble("class_1_b"); |
||
| 74 | mNitrogen_2a = resp.valueDouble("class_2_a"); |
||
| 75 | mNitrogen_2b = resp.valueDouble("class_2_b"); |
||
| 76 | mNitrogen_3a = resp.valueDouble("class_3_a"); |
||
| 77 | mNitrogen_3b = resp.valueDouble("class_3_b"); |
||
| 78 | if (mNitrogen_1a*mNitrogen_1b*mNitrogen_2a*mNitrogen_2b*mNitrogen_3a*mNitrogen_3b == 0) |
||
| 79 | throw IException("at least one parameter of model.species.nitrogenResponseClasses is not valid (value=0)!"); |
||
| 80 | |||
| 81 | // setup CO2 response |
||
| 82 | XmlHelper co2(xml.node("model.species.CO2Response")); |
||
| 83 | mCO2base = co2.valueDouble("baseConcentration"); |
||
| 84 | mCO2comp = co2.valueDouble("compensationPoint"); |
||
| 85 | mCO2beta0 = co2.valueDouble("beta0"); |
||
| 86 | mCO2p0 = co2.valueDouble("p0"); |
||
| 87 | if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0) |
||
| 88 | throw IException("at least one parameter of model.species.CO2Response is not valid!"); |
||
| 89 | |||
| 91 | Werner | 90 | return mSpecies.count(); |
| 91 | |||
| 92 | } |
||
| 93 | /** retrieves variables from the datasource available during the setup of species. |
||
| 94 | */ |
||
| 95 | QVariant SpeciesSet::var(const QString& varName) |
||
| 96 | { |
||
| 94 | Werner | 97 | Q_ASSERT(mSetupQuery!=0); |
| 91 | Werner | 98 | |
| 99 | int idx = mSetupQuery->record().indexOf(varName); |
||
| 100 | if (idx>=0) |
||
| 101 | return mSetupQuery->value(idx); |
||
| 125 | Werner | 102 | throw IException(QString("SpeciesSet: variable not set: %1").arg(varName)); |
| 120 | Werner | 103 | //throw IException(QString("load species parameter: field %1 not found!").arg(varName)); |
| 91 | Werner | 104 | // lookup in defaults |
| 119 | Werner | 105 | //qDebug() << "variable" << varName << "not found - using default."; |
| 106 | //return GlobalSettings::instance()->settingDefaultValue(varName); |
||
| 91 | Werner | 107 | } |
| 209 | werner | 108 | |
| 109 | inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const |
||
| 110 | { |
||
| 111 | if (availableNitrogen<=NB) |
||
| 112 | return 0; |
||
| 113 | double x = 1. - exp(NA * (availableNitrogen-NB)); |
||
| 114 | return x; |
||
| 115 | } |
||
| 116 | |||
| 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) |
||
| 119 | double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const |
||
| 120 | { |
||
| 121 | double value1, value2, value3; |
||
| 122 | if (responseClass>2.) { |
||
| 123 | if (responseClass==3.) |
||
| 124 | return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
||
| 125 | else { |
||
| 126 | // interpolate between 2 and 3 |
||
| 127 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
||
| 128 | value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
||
| 129 | return value2 + (responseClass-2)*(value3-value2); |
||
| 130 | } |
||
| 131 | } |
||
| 132 | if (responseClass==2) |
||
| 133 | return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
||
| 134 | if (responseClass==1) |
||
| 135 | return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
||
| 136 | // last ressort: interpolate between 1 and 2 |
||
| 137 | value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
||
| 138 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
||
| 139 | return value1 + (responseClass-1)*(value2-value1); |
||
| 140 | } |
||
| 141 | |||
| 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) |
||
| 144 | */ |
||
| 145 | double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const |
||
| 146 | { |
||
| 210 | werner | 147 | if (nitrogenResponse==0) |
| 148 | return 0.; |
||
| 149 | |||
| 209 | werner | 150 | double co2_water = 2. - soilWaterResponse; |
| 151 | double beta = mCO2beta0 * co2_water * nitrogenResponse; |
||
| 152 | |||
| 153 | double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17) |
||
| 154 | |||
| 155 | // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions |
||
| 156 | double deltaC = mCO2base - mCO2comp; |
||
| 157 | double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16 |
||
| 158 | double K1 = (1. + K2*deltaC) / deltaC; |
||
| 159 | |||
| 160 | double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16 |
||
| 161 | return response; |
||
| 162 | |||
| 163 | } |