Rev 211 | Rev 214 | 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 | { |
211 | werner | 39 | // setup phenology |
40 | setupPhenology(); |
||
41 | |||
102 | Werner | 42 | const XmlHelper &xml = GlobalSettings::instance()->settings(); |
191 | werner | 43 | QString tableName = xml.value("model.species.source", "species"); |
44 | QString readerFile = xml.value("model.species.reader", "reader.bin"); |
||
102 | Werner | 45 | readerFile = GlobalSettings::instance()->path(readerFile, "lip"); |
46 | mReaderStamp.load(readerFile); |
||
47 | |||
91 | Werner | 48 | QSqlQuery query(GlobalSettings::instance()->dbin()); |
49 | mSetupQuery = &query; |
||
50 | QString sql = QString("select * from %1").arg(tableName); |
||
51 | query.exec(sql); |
||
52 | clear(); |
||
53 | qDebug() << "attempting to load a species set from" << tableName; |
||
54 | while (query.next()) { |
||
55 | if (var("active").toInt()==0) |
||
56 | continue; |
||
57 | |||
58 | Species *s = new Species(this); // create |
||
99 | Werner | 59 | // call setup routine (which calls SpeciesSet::var() to retrieve values |
91 | Werner | 60 | s->setup(); |
61 | |||
62 | mSpecies.insert(s->id(), s); // store |
||
179 | werner | 63 | if (s->active()) |
64 | mActiveSpecies.append(s); |
||
91 | Werner | 65 | } // while query.next() |
66 | qDebug() << "loaded" << mSpecies.count() << "active species:"; |
||
67 | qDebug() << mSpecies.keys(); |
||
68 | |||
69 | mSetupQuery = 0; |
||
209 | werner | 70 | |
71 | // setup nitrogen response |
||
72 | XmlHelper resp(xml.node("model.species.nitrogenResponseClasses")); |
||
73 | if (!resp.isValid()) |
||
74 | throw IException("model.species.nitrogenResponseClasses not present!"); |
||
75 | mNitrogen_1a = resp.valueDouble("class_1_a"); |
||
76 | mNitrogen_1b = resp.valueDouble("class_1_b"); |
||
77 | mNitrogen_2a = resp.valueDouble("class_2_a"); |
||
78 | mNitrogen_2b = resp.valueDouble("class_2_b"); |
||
79 | mNitrogen_3a = resp.valueDouble("class_3_a"); |
||
80 | mNitrogen_3b = resp.valueDouble("class_3_b"); |
||
81 | 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)!"); |
||
83 | |||
84 | // setup CO2 response |
||
85 | XmlHelper co2(xml.node("model.species.CO2Response")); |
||
86 | mCO2base = co2.valueDouble("baseConcentration"); |
||
87 | mCO2comp = co2.valueDouble("compensationPoint"); |
||
88 | mCO2beta0 = co2.valueDouble("beta0"); |
||
89 | mCO2p0 = co2.valueDouble("p0"); |
||
90 | if (mCO2base*mCO2comp*(mCO2base-mCO2comp)*mCO2beta0*mCO2p0==0) |
||
91 | throw IException("at least one parameter of model.species.CO2Response is not valid!"); |
||
92 | |||
91 | Werner | 93 | return mSpecies.count(); |
94 | |||
95 | } |
||
211 | werner | 96 | |
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(""))); |
||
213 | werner | 111 | xml.setCurrentNode(n); |
211 | werner | 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 | |||
91 | Werner | 125 | /** retrieves variables from the datasource available during the setup of species. |
126 | */ |
||
127 | QVariant SpeciesSet::var(const QString& varName) |
||
128 | { |
||
94 | Werner | 129 | Q_ASSERT(mSetupQuery!=0); |
91 | Werner | 130 | |
131 | int idx = mSetupQuery->record().indexOf(varName); |
||
132 | if (idx>=0) |
||
133 | return mSetupQuery->value(idx); |
||
125 | Werner | 134 | throw IException(QString("SpeciesSet: variable not set: %1").arg(varName)); |
120 | Werner | 135 | //throw IException(QString("load species parameter: field %1 not found!").arg(varName)); |
91 | Werner | 136 | // lookup in defaults |
119 | Werner | 137 | //qDebug() << "variable" << varName << "not found - using default."; |
138 | //return GlobalSettings::instance()->settingDefaultValue(varName); |
||
91 | Werner | 139 | } |
209 | werner | 140 | |
141 | inline double SpeciesSet::nitrogenResponse(const double &availableNitrogen, const double &NA, const double &NB) const |
||
142 | { |
||
143 | if (availableNitrogen<=NB) |
||
144 | return 0; |
||
145 | double x = 1. - exp(NA * (availableNitrogen-NB)); |
||
146 | return x; |
||
147 | } |
||
148 | |||
149 | /// 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) |
||
151 | double SpeciesSet::nitrogenResponse(const double availableNitrogen, const double &responseClass) const |
||
152 | { |
||
153 | double value1, value2, value3; |
||
154 | if (responseClass>2.) { |
||
155 | if (responseClass==3.) |
||
156 | return nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
||
157 | else { |
||
158 | // interpolate between 2 and 3 |
||
159 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
||
160 | value3 = nitrogenResponse(availableNitrogen, mNitrogen_3a, mNitrogen_3b); |
||
161 | return value2 + (responseClass-2)*(value3-value2); |
||
162 | } |
||
163 | } |
||
164 | if (responseClass==2) |
||
165 | return nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
||
166 | if (responseClass==1) |
||
167 | return nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
||
168 | // last ressort: interpolate between 1 and 2 |
||
169 | value1 = nitrogenResponse(availableNitrogen, mNitrogen_1a, mNitrogen_1b); |
||
170 | value2 = nitrogenResponse(availableNitrogen, mNitrogen_2a, mNitrogen_2b); |
||
171 | return value1 + (responseClass-1)*(value2-value1); |
||
172 | } |
||
173 | |||
174 | /** 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) |
||
176 | */ |
||
177 | double SpeciesSet::co2Response(const double ambientCO2, const double nitrogenResponse, const double soilWaterResponse) const |
||
178 | { |
||
210 | werner | 179 | if (nitrogenResponse==0) |
180 | return 0.; |
||
181 | |||
209 | werner | 182 | double co2_water = 2. - soilWaterResponse; |
183 | double beta = mCO2beta0 * co2_water * nitrogenResponse; |
||
184 | |||
185 | double r =1. + M_LN2 * beta; // NPP increase for a doubling of atmospheric CO2 (Eq. 17) |
||
186 | |||
187 | // fertilization function (cf. Farquhar, 1980) based on Michaelis-Menten expressions |
||
188 | double deltaC = mCO2base - mCO2comp; |
||
189 | double K2 = ((2*mCO2base - mCO2comp) - r*deltaC ) / ((r-1.)*deltaC*(2*mCO2base - mCO2comp)); // Eq. 16 |
||
190 | double K1 = (1. + K2*deltaC) / deltaC; |
||
191 | |||
192 | double response = mCO2p0 * K1*(ambientCO2 - mCO2comp) / (1 + K2*(ambientCO2-mCO2comp)); // Eq. 16 |
||
193 | return response; |
||
194 | |||
195 | } |
||
211 | werner | 196 | |
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 | } |