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/snag.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/snag.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 "snag.h"
|
21 | #include "snag.h"
|
22 | #include "tree.h"
|
22 | #include "tree.h"
|
23 | #include "species.h"
|
23 | #include "species.h"
|
24 | #include "globalsettings.h"
|
24 | #include "globalsettings.h"
|
25 | #include "expression.h"
|
25 | #include "expression.h"
|
26 | // for calculation of climate decomposition
|
26 | // for calculation of climate decomposition
|
27 | #include "resourceunit.h"
|
27 | #include "resourceunit.h"
|
28 | #include "watercycle.h"
|
28 | #include "watercycle.h"
|
29 | #include "climate.h"
|
29 | #include "climate.h"
|
30 | #include "model.h"
|
30 | #include "model.h"
|
31 | 31 | ||
32 | /** @class Snag
|
32 | /** @class Snag
|
33 | @ingroup core
|
33 | @ingroup core
|
34 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
|
34 | Snag deals with carbon / nitrogen fluxes from the forest until the reach soil pools.
|
35 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
|
35 | Snag lives on the level of the ResourceUnit; carbon fluxes from trees enter Snag, and parts of the biomass of snags
|
36 | is subsequently forwarded to the soil sub model.
|
36 | is subsequently forwarded to the soil sub model.
|
37 | Carbon is stored in three classes (depending on the size)
|
37 | Carbon is stored in three classes (depending on the size)
|
38 | The Snag dynamics class uses the following species parameter:
|
38 | The Snag dynamics class uses the following species parameter:
|
39 | cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
|
39 | cnFoliage, cnFineroot, cnWood, snagHalflife, snagKSW
|
40 | 40 | ||
41 | */
|
41 | */
|
42 | // static variables
|
42 | // static variables
|
43 | double Snag::mDBHLower = -1.; |
43 | double Snag::mDBHLower = -1.; |
44 | double Snag::mDBHHigher = 0.; |
44 | double Snag::mDBHHigher = 0.; |
45 | double Snag::mCarbonThreshold[] = {0., 0., 0.}; |
45 | double Snag::mCarbonThreshold[] = {0., 0., 0.}; |
46 | 46 | ||
47 | double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
47 | double CNPair::biomassCFraction = biomassCFraction; // get global from globalsettings.h |
48 | 48 | ||
49 | /// add biomass and weigh the parameter_value with the current C-content of the pool
|
49 | /// add biomass and weigh the parameter_value with the current C-content of the pool
|
50 | void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value) |
50 | void CNPool::addBiomass(const double biomass, const double CNratio, const double parameter_value) |
51 | {
|
51 | {
|
52 | if (biomass==0.) |
52 | if (biomass==0.) |
53 | return; |
53 | return; |
54 | double new_c = biomass*biomassCFraction; |
54 | double new_c = biomass*biomassCFraction; |
55 | double p_old = C / (new_c + C); |
55 | double p_old = C / (new_c + C); |
56 | mParameter = mParameter*p_old + parameter_value*(1.-p_old); |
56 | mParameter = mParameter*p_old + parameter_value*(1.-p_old); |
57 | CNPair::addBiomass(biomass, CNratio); |
57 | CNPair::addBiomass(biomass, CNratio); |
58 | }
|
58 | }
|
59 | 59 | ||
60 | // increase pool (and weigh the value)
|
60 | // increase pool (and weigh the value)
|
61 | void CNPool::operator+=(const CNPool &s) |
61 | void CNPool::operator+=(const CNPool &s) |
62 | {
|
62 | {
|
63 | if (s.C==0.) |
63 | if (s.C==0.) |
64 | return; |
64 | return; |
65 | mParameter = parameter(s); // calculate weighted parameter |
65 | mParameter = parameter(s); // calculate weighted parameter |
66 | C+=s.C; |
66 | C+=s.C; |
67 | N+=s.N; |
67 | N+=s.N; |
68 | }
|
68 | }
|
69 | 69 | ||
70 | double CNPool::parameter(const CNPool &s) const |
70 | double CNPool::parameter(const CNPool &s) const |
71 | {
|
71 | {
|
72 | if (s.C==0.) |
72 | if (s.C==0.) |
73 | return parameter(); |
73 | return parameter(); |
74 | double p_old = C / (s.C + C); |
74 | double p_old = C / (s.C + C); |
75 | double result = mParameter*p_old + s.parameter()*(1.-p_old); |
75 | double result = mParameter*p_old + s.parameter()*(1.-p_old); |
76 | return result; |
76 | return result; |
77 | }
|
77 | }
|
78 | 78 | ||
79 | 79 | ||
80 | void Snag::setupThresholds(const double lower, const double upper) |
80 | void Snag::setupThresholds(const double lower, const double upper) |
81 | {
|
81 | {
|
82 | if (mDBHLower == lower) |
82 | if (mDBHLower == lower) |
83 | return; |
83 | return; |
84 | mDBHLower = lower; |
84 | mDBHLower = lower; |
85 | mDBHHigher = upper; |
85 | mDBHHigher = upper; |
86 | mCarbonThreshold[0] = lower / 2.; |
86 | mCarbonThreshold[0] = lower / 2.; |
87 | mCarbonThreshold[1] = lower + (upper - lower)/2.; |
87 | mCarbonThreshold[1] = lower + (upper - lower)/2.; |
88 | mCarbonThreshold[2] = upper + (upper - lower)/2.; |
88 | mCarbonThreshold[2] = upper + (upper - lower)/2.; |
89 | //# threshold levels for emptying out the dbh-snag-classes
|
89 | //# threshold levels for emptying out the dbh-snag-classes
|
90 | //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
|
90 | //# derived from Psme woody allometry, converted to C, with a threshold level set to 10%
|
91 | //# values in kg!
|
91 | //# values in kg!
|
92 | for (int i=0;i<3;i++) |
92 | for (int i=0;i<3;i++) |
93 | mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1; |
93 | mCarbonThreshold[i] = 0.10568*pow(mCarbonThreshold[i],2.4247)*0.5*0.1; |
94 | }
|
94 | }
|
95 | 95 | ||
96 | 96 | ||
97 | Snag::Snag() |
97 | Snag::Snag() |
98 | {
|
98 | {
|
99 | mRU = 0; |
99 | mRU = 0; |
100 | CNPair::setCFraction(biomassCFraction); |
100 | CNPair::setCFraction(biomassCFraction); |
101 | }
|
101 | }
|
102 | 102 | ||
103 | void Snag::setup( const ResourceUnit *ru) |
103 | void Snag::setup( const ResourceUnit *ru) |
104 | {
|
104 | {
|
105 | mRU = ru; |
105 | mRU = ru; |
106 | mClimateFactor = 0.; |
106 | mClimateFactor = 0.; |
107 | // branches
|
107 | // branches
|
108 | mBranchCounter=0; |
108 | mBranchCounter=0; |
109 | for (int i=0;i<3;i++) { |
109 | for (int i=0;i<3;i++) { |
110 | mTimeSinceDeath[i] = 0.; |
110 | mTimeSinceDeath[i] = 0.; |
111 | mNumberOfSnags[i] = 0.; |
111 | mNumberOfSnags[i] = 0.; |
112 | mAvgDbh[i] = 0.; |
112 | mAvgDbh[i] = 0.; |
113 | mAvgHeight[i] = 0.; |
113 | mAvgHeight[i] = 0.; |
114 | mAvgVolume[i] = 0.; |
114 | mAvgVolume[i] = 0.; |
115 | mKSW[i] = 0.; |
115 | mKSW[i] = 0.; |
116 | mCurrentKSW[i] = 0.; |
116 | mCurrentKSW[i] = 0.; |
117 | mHalfLife[i] = 0.; |
117 | mHalfLife[i] = 0.; |
118 | }
|
118 | }
|
119 | mTotalSnagCarbon = 0.; |
119 | mTotalSnagCarbon = 0.; |
120 | if (mDBHLower<=0) |
120 | if (mDBHLower<=0) |
121 | throw IException("Snag::setupThresholds() not called or called with invalid parameters."); |
121 | throw IException("Snag::setupThresholds() not called or called with invalid parameters."); |
122 | 122 | ||
123 | // Inital values from XML file
|
123 | // Inital values from XML file
|
124 | XmlHelper xml=GlobalSettings::instance()->settings(); |
124 | XmlHelper xml=GlobalSettings::instance()->settings(); |
125 | double kyr = xml.valueDouble("model.site.youngRefractoryDecompRate", -1); |
125 | double kyr = xml.valueDouble("model.site.youngRefractoryDecompRate", -1); |
126 | // put carbon of snags to the middle size class
|
126 | // put carbon of snags to the middle size class
|
127 | xml.setCurrentNode("model.initialization.snags"); |
127 | xml.setCurrentNode("model.initialization.snags"); |
128 | mSWD[1].C = xml.valueDouble(".swdC"); |
128 | mSWD[1].C = xml.valueDouble(".swdC"); |
129 | mSWD[1].N = mSWD[1].C / xml.valueDouble(".swdCN", 50.); |
129 | mSWD[1].N = mSWD[1].C / xml.valueDouble(".swdCN", 50.); |
130 | mSWD[1].setParameter(kyr); |
130 | mSWD[1].setParameter(kyr); |
131 | mKSW[1] = xml.valueDouble(".swdDecompRate"); |
131 | mKSW[1] = xml.valueDouble(".swdDecompRate"); |
132 | mNumberOfSnags[1] = xml.valueDouble(".swdCount"); |
132 | mNumberOfSnags[1] = xml.valueDouble(".swdCount"); |
133 | mHalfLife[1] = xml.valueDouble(".swdHalfLife"); |
133 | mHalfLife[1] = xml.valueDouble(".swdHalfLife"); |
134 | // and for the Branch/coarse root pools: split the init value into five chunks
|
134 | // and for the Branch/coarse root pools: split the init value into five chunks
|
135 | CNPool other(xml.valueDouble(".otherC"), xml.valueDouble(".otherC")/xml.valueDouble(".otherCN", 50.), kyr ); |
135 | CNPool other(xml.valueDouble(".otherC"), xml.valueDouble(".otherC")/xml.valueDouble(".otherCN", 50.), kyr ); |
136 | 136 | ||
137 | mTotalSnagCarbon = other.C + mSWD[1].C; |
137 | mTotalSnagCarbon = other.C + mSWD[1].C; |
138 | 138 | ||
139 | other *= 0.2; |
139 | other *= 0.2; |
140 | for (int i=0;i<5;i++) |
140 | for (int i=0;i<5;i++) |
141 | mOtherWood[i] = other; |
141 | mOtherWood[i] = other; |
142 | }
|
142 | }
|
143 | 143 | ||
144 | void Snag::scaleInitialState() |
144 | void Snag::scaleInitialState() |
145 | {
|
145 | {
|
146 | double area_factor = mRU->stockableArea() / cRUArea; // fraction stockable area |
146 | double area_factor = mRU->stockableArea() / cRUArea; // fraction stockable area |
147 | // avoid huge snag pools on very small resource units (see also soil.cpp)
|
147 | // avoid huge snag pools on very small resource units (see also soil.cpp)
|
148 | // area_factor = std::max(area_factor, 0.1);
|
148 | // area_factor = std::max(area_factor, 0.1);
|
149 | mSWD[1] *= area_factor; |
149 | mSWD[1] *= area_factor; |
150 | mNumberOfSnags[1] *= area_factor; |
150 | mNumberOfSnags[1] *= area_factor; |
151 | for (int i=0;i<5;i++) |
151 | for (int i=0;i<5;i++) |
152 | mOtherWood[i]*= area_factor; |
152 | mOtherWood[i]*= area_factor; |
153 | mTotalSnagCarbon *= area_factor; |
153 | mTotalSnagCarbon *= area_factor; |
154 | 154 | ||
155 | }
|
155 | }
|
156 | 156 | ||
157 | // debug outputs
|
157 | // debug outputs
|
158 | QList<QVariant> Snag::debugList() |
158 | QList<QVariant> Snag::debugList() |
159 | {
|
159 | {
|
160 | // list columns
|
160 | // list columns
|
161 | // for three pools
|
161 | // for three pools
|
162 | QList<QVariant> list; |
162 | QList<QVariant> list; |
163 | 163 | ||
164 | // totals
|
164 | // totals
|
165 | list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N; |
165 | list << mTotalSnagCarbon << mTotalIn.C << mTotalToAtm.C << mSWDtoSoil.C << mSWDtoSoil.N; |
166 | // fluxes to labile soil pool and to refractory soil pool
|
166 | // fluxes to labile soil pool and to refractory soil pool
|
167 | list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N; |
167 | list << mLabileFlux.C << mLabileFlux.N << mRefractoryFlux.C << mRefractoryFlux.N; |
168 | 168 | ||
169 | for (int i=0;i<3;i++) { |
169 | for (int i=0;i<3;i++) { |
170 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
|
170 | // pools "swdx_c", "swdx_n", "swdx_count", "swdx_tsd", "toswdx_c", "toswdx_n"
|
171 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
171 | list << mSWD[i].C << mSWD[i].N << mNumberOfSnags[i] << mTimeSinceDeath[i] << mToSWD[i].C << mToSWD[i].N; |
172 | list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i]; |
172 | list << mAvgDbh[i] << mAvgHeight[i] << mAvgVolume[i]; |
173 | }
|
173 | }
|
174 | 174 | ||
175 | // branch/coarse wood pools (5 yrs)
|
175 | // branch/coarse wood pools (5 yrs)
|
176 | for (int i=0;i<5;i++) { |
176 | for (int i=0;i<5;i++) { |
177 | list << mOtherWood[i].C << mOtherWood[i].N; |
177 | list << mOtherWood[i].C << mOtherWood[i].N; |
178 | }
|
178 | }
|
179 | // list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
|
179 | // list << mOtherWood[mBranchCounter].C << mOtherWood[mBranchCounter].N
|
180 | // << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
|
180 | // << mOtherWood[(mBranchCounter+1)%5].C << mOtherWood[(mBranchCounter+1)%5].N
|
181 | // << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
|
181 | // << mOtherWood[(mBranchCounter+2)%5].C << mOtherWood[(mBranchCounter+2)%5].N
|
182 | // << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
|
182 | // << mOtherWood[(mBranchCounter+3)%5].C << mOtherWood[(mBranchCounter+3)%5].N
|
183 | // << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
|
183 | // << mOtherWood[(mBranchCounter+4)%5].C << mOtherWood[(mBranchCounter+4)%5].N;
|
184 | return list; |
184 | return list; |
185 | }
|
185 | }
|
186 | 186 | ||
187 | 187 | ||
188 | void Snag::newYear() |
188 | void Snag::newYear() |
189 | {
|
189 | {
|
190 | for (int i=0;i<3;i++) { |
190 | for (int i=0;i<3;i++) { |
191 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
191 | mToSWD[i].clear(); // clear transfer pools to standing-woody-debris |
192 | mCurrentKSW[i] = 0.; |
192 | mCurrentKSW[i] = 0.; |
193 | }
|
193 | }
|
194 | mLabileFlux.clear(); |
194 | mLabileFlux.clear(); |
195 | mRefractoryFlux.clear(); |
195 | mRefractoryFlux.clear(); |
196 | mTotalToAtm.clear(); |
196 | mTotalToAtm.clear(); |
197 | mTotalToExtern.clear(); |
197 | mTotalToExtern.clear(); |
198 | mTotalToDisturbance.clear(); |
198 | mTotalToDisturbance.clear(); |
199 | mTotalIn.clear(); |
199 | mTotalIn.clear(); |
200 | mSWDtoSoil.clear(); |
200 | mSWDtoSoil.clear(); |
201 | }
|
201 | }
|
202 | 202 | ||
203 | /// calculate the dynamic climate modifier for decomposition 're'
|
203 | /// calculate the dynamic climate modifier for decomposition 're'
|
204 | /// calculation is done on the level of ResourceUnit because the water content per day is needed.
|
204 | /// calculation is done on the level of ResourceUnit because the water content per day is needed.
|
205 | double Snag::calculateClimateFactors() |
205 | double Snag::calculateClimateFactors() |
206 | {
|
206 | {
|
207 | // the calculation of climate factors requires calculated evapotranspiration. In cases without vegetation (trees or saplings)
|
207 | // the calculation of climate factors requires calculated evapotranspiration. In cases without vegetation (trees or saplings)
|
208 | // we have to trigger the water cycle calculation for ourselves [ the waterCycle checks if it has already been run in a year and doesn't run twice in that case ]
|
208 | // we have to trigger the water cycle calculation for ourselves [ the waterCycle checks if it has already been run in a year and doesn't run twice in that case ]
|
209 | const_cast<WaterCycle*>(mRU->waterCycle())->run(); |
209 | const_cast<WaterCycle*>(mRU->waterCycle())->run(); |
210 | double ft, fw; |
210 | double ft, fw; |
211 | double f_sum = 0.; |
211 | double f_sum = 0.; |
212 | int iday=0; |
212 | int iday=0; |
213 | // calculate the water-factor for each month (see Adair et al 2008)
|
213 | // calculate the water-factor for each month (see Adair et al 2008)
|
214 | double fw_month[12]; |
214 | double fw_month[12]; |
215 | double ratio; |
215 | double ratio; |
216 | for (int m=0;m<12;m++) { |
216 | for (int m=0;m<12;m++) { |
217 | if (mRU->waterCycle()->referenceEvapotranspiration()[m]>0.) |
217 | if (mRU->waterCycle()->referenceEvapotranspiration()[m]>0.) |
218 | ratio = mRU->climate()->precipitationMonth()[m] / mRU->waterCycle()->referenceEvapotranspiration()[m]; |
218 | ratio = mRU->climate()->precipitationMonth()[m] / mRU->waterCycle()->referenceEvapotranspiration()[m]; |
219 | else
|
219 | else
|
220 | ratio = 0; |
220 | ratio = 0; |
221 | fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio)); |
221 | fw_month[m] = 1. / (1. + 30.*exp(-8.5*ratio)); |
222 | if (logLevelDebug()) qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->referenceEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m]; |
222 | if (logLevelDebug()) qDebug() <<"month"<< m << "PET" << mRU->waterCycle()->referenceEvapotranspiration()[m] << "prec" <<mRU->climate()->precipitationMonth()[m]; |
223 | }
|
223 | }
|
224 | 224 | ||
225 | for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday) |
225 | for (const ClimateDay *day=mRU->climate()->begin(); day!=mRU->climate()->end(); ++day, ++iday) |
226 | {
|
226 | {
|
227 | ft = exp(308.56*(1./56.02-1./((273.+day->temperature)-227.13))); // empirical variable Q10 model of Lloyd and Taylor (1994), see also Adair et al. (2008) |
227 | ft = exp(308.56*(1./56.02-1./((273.+day->temperature)-227.13))); // empirical variable Q10 model of Lloyd and Taylor (1994), see also Adair et al. (2008) |
228 | fw = fw_month[day->month-1]; |
228 | fw = fw_month[day->month-1]; |
229 | 229 | ||
230 | f_sum += ft*fw; |
230 | f_sum += ft*fw; |
231 | }
|
231 | }
|
232 | // the climate factor is defined as the arithmentic annual mean value
|
232 | // the climate factor is defined as the arithmentic annual mean value
|
233 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
233 | mClimateFactor = f_sum / double(mRU->climate()->daysOfYear()); |
234 | return mClimateFactor; |
234 | return mClimateFactor; |
235 | }
|
235 | }
|
236 | 236 | ||
237 | /// do the yearly calculation
|
237 | /// do the yearly calculation
|
238 | /// see http://iland.boku.ac.at/snag+dynamics
|
238 | /// see http://iland.boku.ac.at/snag+dynamics
|
239 | void Snag::calculateYear() |
239 | void Snag::calculateYear() |
240 | {
|
240 | {
|
241 | mSWDtoSoil.clear(); |
241 | mSWDtoSoil.clear(); |
242 | 242 | ||
243 | // calculate anyway, because also the soil module needs it (and currently one can have Snag and Soil only as a couple)
|
243 | // calculate anyway, because also the soil module needs it (and currently one can have Snag and Soil only as a couple)
|
244 | calculateClimateFactors(); |
244 | calculateClimateFactors(); |
245 | const double climate_factor_re = mClimateFactor; |
245 | const double climate_factor_re = mClimateFactor; |
246 | 246 | ||
247 | if (isEmpty()) // nothing to do |
247 | if (isEmpty()) // nothing to do |
248 | return; |
248 | return; |
249 | 249 | ||
250 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
|
250 | // process branches: every year one of the five baskets is emptied and transfered to the refractory soil pool
|
251 | mRefractoryFlux+=mOtherWood[mBranchCounter]; |
251 | mRefractoryFlux+=mOtherWood[mBranchCounter]; |
252 | mOtherWood[mBranchCounter].clear(); |
252 | mOtherWood[mBranchCounter].clear(); |
253 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
253 | mBranchCounter= (mBranchCounter+1) % 5; // increase index, roll over to 0. |
254 | 254 | ||
255 | // decay of branches/coarse roots
|
255 | // decay of branches/coarse roots
|
256 | for (int i=0;i<5;i++) { |
256 | for (int i=0;i<5;i++) { |
257 | if (mOtherWood[i].C>0.) { |
257 | if (mOtherWood[i].C>0.) { |
258 | double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value... |
258 | double survive_rate = exp(- climate_factor_re * mOtherWood[i].parameter() ); // parameter: the "kyr" value... |
259 | mTotalToAtm.C += mOtherWood[i].C * (1. - survive_rate); // flux to atmosphere (decayed carbon) |
259 | mTotalToAtm.C += mOtherWood[i].C * (1. - survive_rate); // flux to atmosphere (decayed carbon) |
260 | mOtherWood[i].C *= survive_rate; |
260 | mOtherWood[i].C *= survive_rate; |
261 | }
|
261 | }
|
262 | }
|
262 | }
|
263 | 263 | ||
264 | // process standing snags.
|
264 | // process standing snags.
|
265 | // the input of the current year is in the mToSWD-Pools
|
265 | // the input of the current year is in the mToSWD-Pools
|
266 | for (int i=0;i<3;i++) { |
266 | for (int i=0;i<3;i++) { |
267 | 267 | ||
268 | // update the swd-pool with this years' input
|
268 | // update the swd-pool with this years' input
|
269 | if (!mToSWD[i].isEmpty()) { |
269 | if (!mToSWD[i].isEmpty()) { |
270 | // update decay rate (apply average yearly input to the state parameters)
|
270 | // update decay rate (apply average yearly input to the state parameters)
|
271 | mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C)); |
271 | mKSW[i] = mKSW[i]*(mSWD[i].C/(mSWD[i].C+mToSWD[i].C)) + mCurrentKSW[i]*(mToSWD[i].C/(mSWD[i].C+mToSWD[i].C)); |
272 | //move content to the SWD pool
|
272 | //move content to the SWD pool
|
273 | mSWD[i] += mToSWD[i]; |
273 | mSWD[i] += mToSWD[i]; |
274 | }
|
274 | }
|
275 | 275 | ||
276 | if (mSWD[i].C > 0) { |
276 | if (mSWD[i].C > 0) { |
277 | // reduce the Carbon (note: the N stays, thus the CN ratio changes)
|
277 | // reduce the Carbon (note: the N stays, thus the CN ratio changes)
|
278 | // use the decay rate that is derived as a weighted average of all standing woody debris
|
278 | // use the decay rate that is derived as a weighted average of all standing woody debris
|
279 | double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep |
279 | double survive_rate = exp(-mKSW[i] *climate_factor_re * 1. ); // 1: timestep |
280 | mTotalToAtm.C += mSWD[i].C * (1. - survive_rate); |
280 | mTotalToAtm.C += mSWD[i].C * (1. - survive_rate); |
281 | mSWD[i].C *= survive_rate; |
281 | mSWD[i].C *= survive_rate; |
282 | 282 | ||
283 | // transition to downed woody debris
|
283 | // transition to downed woody debris
|
284 | // update: use negative exponential decay, species parameter: half-life
|
284 | // update: use negative exponential decay, species parameter: half-life
|
285 | // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
|
285 | // modified for the climatic effect on decomposition, i.e. if decomp is slower, snags stand longer and vice versa
|
286 | // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
|
286 | // this is loosely oriented on Standcarb2 (http://andrewsforest.oregonstate.edu/pubs/webdocs/models/standcarb2.htm),
|
287 | // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
|
287 | // where lag times for cohort transitions are linearly modified with re although here individual good or bad years have
|
288 | // an immediate effect, the average climatic influence should come through (and it is inherently transient)
|
288 | // an immediate effect, the average climatic influence should come through (and it is inherently transient)
|
289 | // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
|
289 | // note that swd.hl is species-specific, and thus a weighted average over the species in the input (=mortality)
|
290 | // needs to be calculated, followed by a weighted update of the previous swd.hl.
|
290 | // needs to be calculated, followed by a weighted update of the previous swd.hl.
|
291 | // As weights here we use stem number, as the processes here pertain individual snags
|
291 | // As weights here we use stem number, as the processes here pertain individual snags
|
292 | // calculate the transition probability of SWD to downed dead wood
|
292 | // calculate the transition probability of SWD to downed dead wood
|
293 | 293 | ||
294 | double half_life = mHalfLife[i] / climate_factor_re; |
294 | double half_life = mHalfLife[i] / climate_factor_re; |
295 | double rate = -M_LN2 / half_life; // M_LN2: math. constant |
295 | double rate = -M_LN2 / half_life; // M_LN2: math. constant |
296 | 296 | ||
297 | // higher decay rate for the class with smallest diameters
|
297 | // higher decay rate for the class with smallest diameters
|
298 | if (i==0) |
298 | if (i==0) |
299 | rate*=2.; |
299 | rate*=2.; |
300 | 300 | ||
301 | double transfer = 1. - exp(rate); |
301 | double transfer = 1. - exp(rate); |
302 | 302 | ||
303 | // calculate flow to soil pool...
|
303 | // calculate flow to soil pool...
|
304 | mSWDtoSoil += mSWD[i] * transfer; |
304 | mSWDtoSoil += mSWD[i] * transfer; |
305 | mRefractoryFlux += mSWD[i] * transfer; |
305 | mRefractoryFlux += mSWD[i] * transfer; |
306 | mSWD[i] *= (1.-transfer); // reduce pool |
306 | mSWD[i] *= (1.-transfer); // reduce pool |
307 | // calculate the stem number of remaining snags
|
307 | // calculate the stem number of remaining snags
|
308 | mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer); |
308 | mNumberOfSnags[i] = mNumberOfSnags[i] * (1. - transfer); |
309 | 309 | ||
310 | mTimeSinceDeath[i] += 1.; |
310 | mTimeSinceDeath[i] += 1.; |
311 | // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
|
311 | // if stems<0.5, empty the whole cohort into DWD, i.e. release the last bit of C and N and clear the stats
|
312 | // also, if the Carbon of an average snag is less than 10% of the original average tree
|
312 | // also, if the Carbon of an average snag is less than 10% of the original average tree
|
313 | // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
|
313 | // (derived from allometries for the three diameter classes), the whole cohort is emptied out to DWD
|
314 | if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) { |
314 | if (mNumberOfSnags[i] < 0.5 || mSWD[i].C / mNumberOfSnags[i] < mCarbonThreshold[i]) { |
315 | // clear the pool: add the rest to the soil, clear statistics of the pool
|
315 | // clear the pool: add the rest to the soil, clear statistics of the pool
|
316 | mRefractoryFlux += mSWD[i]; |
316 | mRefractoryFlux += mSWD[i]; |
317 | mSWDtoSoil += mSWD[i]; |
317 | mSWDtoSoil += mSWD[i]; |
318 | mSWD[i].clear(); |
318 | mSWD[i].clear(); |
319 | mAvgDbh[i] = 0.; |
319 | mAvgDbh[i] = 0.; |
320 | mAvgHeight[i] = 0.; |
320 | mAvgHeight[i] = 0.; |
321 | mAvgVolume[i] = 0.; |
321 | mAvgVolume[i] = 0.; |
322 | mKSW[i] = 0.; |
322 | mKSW[i] = 0.; |
323 | mCurrentKSW[i] = 0.; |
323 | mCurrentKSW[i] = 0.; |
324 | mHalfLife[i] = 0.; |
324 | mHalfLife[i] = 0.; |
325 | mTimeSinceDeath[i] = 0.; |
325 | mTimeSinceDeath[i] = 0.; |
326 | }
|
326 | }
|
327 | 327 | ||
328 | }
|
328 | }
|
329 | 329 | ||
330 | }
|
330 | }
|
331 | // total carbon in the snag-container on the RU *after* processing is the content of the
|
331 | // total carbon in the snag-container on the RU *after* processing is the content of the
|
332 | // standing woody debris pools + the branches
|
332 | // standing woody debris pools + the branches
|
333 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
333 | mTotalSnagCarbon = mSWD[0].C + mSWD[1].C + mSWD[2].C + |
334 | mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C; |
334 | mOtherWood[0].C + mOtherWood[1].C + mOtherWood[2].C + mOtherWood[3].C + mOtherWood[4].C; |
335 | mTotalSWD = mSWD[0] + mSWD[1] + mSWD[2]; |
335 | mTotalSWD = mSWD[0] + mSWD[1] + mSWD[2]; |
336 | mTotalOther = mOtherWood[0] + mOtherWood[1] + mOtherWood[2] + mOtherWood[3] + mOtherWood[4]; |
336 | mTotalOther = mOtherWood[0] + mOtherWood[1] + mOtherWood[2] + mOtherWood[3] + mOtherWood[4]; |
337 | }
|
337 | }
|
338 | 338 | ||
339 | /// foliage and fineroot litter is transferred during tree growth.
|
339 | /// foliage and fineroot litter is transferred during tree growth.
|
340 | void Snag::addTurnoverLitter(const Species *species, const double litter_foliage, const double litter_fineroot) |
340 | void Snag::addTurnoverLitter(const Species *species, const double litter_foliage, const double litter_fineroot) |
341 | {
|
341 | {
|
342 | mLabileFlux.addBiomass(litter_foliage, species->cnFoliage(), species->snagKyl()); |
342 | mLabileFlux.addBiomass(litter_foliage, species->cnFoliage(), species->snagKyl()); |
343 | mLabileFlux.addBiomass(litter_fineroot, species->cnFineroot(), species->snagKyl()); |
343 | mLabileFlux.addBiomass(litter_fineroot, species->cnFineroot(), species->snagKyl()); |
344 | DBGMODE(
|
344 | DBGMODE(
|
345 | if (isnan(mLabileFlux.C)) |
345 | if (isnan(mLabileFlux.C)) |
346 | qDebug("Snag::addTurnoverLitter: NaN"); |
346 | qDebug("Snag::addTurnoverLitter: NaN"); |
347 | ); |
347 | ); |
348 | }
|
348 | }
|
349 | 349 | ||
350 | void Snag::addTurnoverWood(const Species *species, const double woody_biomass) |
350 | void Snag::addTurnoverWood(const Species *species, const double woody_biomass) |
351 | {
|
351 | {
|
352 | mRefractoryFlux.addBiomass(woody_biomass, species->cnWood(), species->snagKyr()); |
352 | mRefractoryFlux.addBiomass(woody_biomass, species->cnWood(), species->snagKyr()); |
353 | DBGMODE(
|
353 | DBGMODE(
|
354 | if (isnan(mRefractoryFlux.C)) |
354 | if (isnan(mRefractoryFlux.C)) |
355 | qDebug("Snag::addTurnoverWood: NaN"); |
355 | qDebug("Snag::addTurnoverWood: NaN"); |
356 | ); |
356 | ); |
357 | 357 | ||
358 | }
|
358 | }
|
359 | 359 | ||
360 | 360 | ||
361 | /** process the remnants of a single tree.
|
361 | /** process the remnants of a single tree.
|
362 | The part of the stem / branch not covered by snag/soil fraction is removed from the system (e.g. harvest, fire)
|
362 | The part of the stem / branch not covered by snag/soil fraction is removed from the system (e.g. harvest, fire)
|
363 | @param tree the tree to process
|
363 | @param tree the tree to process
|
364 | @param stem_to_snag fraction (0..1) of the stem biomass that should be moved to a standing snag
|
364 | @param stem_to_snag fraction (0..1) of the stem biomass that should be moved to a standing snag
|
365 | @param stem_to_soil fraction (0..1) of the stem biomass that should go directly to the soil
|
365 | @param stem_to_soil fraction (0..1) of the stem biomass that should go directly to the soil
|
366 | @param branch_to_snag fraction (0..1) of the branch biomass that should be moved to a standing snag
|
366 | @param branch_to_snag fraction (0..1) of the branch biomass that should be moved to a standing snag
|
367 | @param branch_to_soil fraction (0..1) of the branch biomass that should go directly to the soil
|
367 | @param branch_to_soil fraction (0..1) of the branch biomass that should go directly to the soil
|
368 | @param foliage_to_soil fraction (0..1) of the foliage biomass that should go directly to the soil
|
368 | @param foliage_to_soil fraction (0..1) of the foliage biomass that should go directly to the soil
|
369 | 369 | ||
370 | */
|
370 | */
|
371 | void Snag::addBiomassPools(const Tree *tree, |
371 | void Snag::addBiomassPools(const Tree *tree, |
372 | const double stem_to_snag, const double stem_to_soil, |
372 | const double stem_to_snag, const double stem_to_soil, |
373 | const double branch_to_snag, const double branch_to_soil, |
373 | const double branch_to_snag, const double branch_to_soil, |
374 | const double foliage_to_soil) |
374 | const double foliage_to_soil) |
375 | {
|
375 | {
|
376 | const Species *species = tree->species(); |
376 | const Species *species = tree->species(); |
377 | 377 | ||
378 | double branch_biomass = tree->biomassBranch(); |
378 | double branch_biomass = tree->biomassBranch(); |
379 | // fine roots go to the labile pool
|
379 | // fine roots go to the labile pool
|
380 | mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), species->snagKyl()); |
380 | mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), species->snagKyl()); |
381 | 381 | ||
382 | // a part of the foliage goes to the soil
|
382 | // a part of the foliage goes to the soil
|
383 | mLabileFlux.addBiomass(tree->biomassFoliage() * foliage_to_soil, species->cnFoliage(), species->snagKyl()); |
383 | mLabileFlux.addBiomass(tree->biomassFoliage() * foliage_to_soil, species->cnFoliage(), species->snagKyl()); |
384 | 384 | ||
385 | //coarse roots and a part of branches are equally distributed over five years:
|
385 | //coarse roots and a part of branches are equally distributed over five years:
|
386 | double biomass_rest = (tree->biomassCoarseRoot() + branch_to_snag*branch_biomass) * 0.2; |
386 | double biomass_rest = (tree->biomassCoarseRoot() + branch_to_snag*branch_biomass) * 0.2; |
387 | for (int i=0;i<5; i++) |
387 | for (int i=0;i<5; i++) |
388 | mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), species->snagKyr()); |
388 | mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), species->snagKyr()); |
389 | 389 | ||
390 | // the other part of the branches goes directly to the soil
|
390 | // the other part of the branches goes directly to the soil
|
391 | mRefractoryFlux.addBiomass(branch_biomass*branch_to_soil, species->cnWood(), species->snagKyr() ); |
391 | mRefractoryFlux.addBiomass(branch_biomass*branch_to_soil, species->cnWood(), species->snagKyr() ); |
392 | // a part of the stem wood goes directly to the soil
|
392 | // a part of the stem wood goes directly to the soil
|
393 | mRefractoryFlux.addBiomass(tree->biomassStem()*stem_to_soil, species->cnWood(), species->snagKyr() ); |
393 | mRefractoryFlux.addBiomass(tree->biomassStem()*stem_to_soil, species->cnWood(), species->snagKyr() ); |
394 | 394 | ||
395 | // just for book-keeping: keep track of all inputs of branches / roots / swd into the "snag" pools
|
395 | // just for book-keeping: keep track of all inputs of branches / roots / swd into the "snag" pools
|
396 | mTotalIn.addBiomass(tree->biomassBranch()*branch_to_snag + tree->biomassCoarseRoot() + tree->biomassStem()*stem_to_snag, species->cnWood()); |
396 | mTotalIn.addBiomass(tree->biomassBranch()*branch_to_snag + tree->biomassCoarseRoot() + tree->biomassStem()*stem_to_snag, species->cnWood()); |
397 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
397 | // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
398 | int pi = poolIndex(tree->dbh()); // get right transfer pool |
398 | int pi = poolIndex(tree->dbh()); // get right transfer pool |
399 | 399 | ||
400 | if (stem_to_snag>0.) { |
400 | if (stem_to_snag>0.) { |
401 | // update statistics - stemnumber-weighted averages
|
401 | // update statistics - stemnumber-weighted averages
|
402 | // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
402 | // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
403 | double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers) |
403 | double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers) |
404 | double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1). |
404 | double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1). |
405 | mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new; |
405 | mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new; |
406 | mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new; |
406 | mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new; |
407 | mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new; |
407 | mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new; |
408 | mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new; |
408 | mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new; |
409 | mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new; |
409 | mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new; |
410 | 410 | ||
411 | // average the decay rate (ksw); this is done based on the carbon content
|
411 | // average the decay rate (ksw); this is done based on the carbon content
|
412 | // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
412 | // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
413 | p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
413 | p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
414 | p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
414 | p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction); |
415 | mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new; |
415 | mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new; |
416 | mNumberOfSnags[pi]++; |
416 | mNumberOfSnags[pi]++; |
417 | }
|
417 | }
|
418 | 418 | ||
419 | // finally add the biomass of the stem to the standing snag pool
|
419 | // finally add the biomass of the stem to the standing snag pool
|
420 | CNPool &to_swd = mToSWD[pi]; |
420 | CNPool &to_swd = mToSWD[pi]; |
421 | to_swd.addBiomass(tree->biomassStem()*stem_to_snag, species->cnWood(), species->snagKyr()); |
421 | to_swd.addBiomass(tree->biomassStem()*stem_to_snag, species->cnWood(), species->snagKyr()); |
422 | 422 | ||
423 | // the biomass that is not routed to snags or directly to the soil
|
423 | // the biomass that is not routed to snags or directly to the soil
|
424 | // is removed from the system (to atmosphere or harvested)
|
424 | // is removed from the system (to atmosphere or harvested)
|
425 | mTotalToExtern.addBiomass(tree->biomassFoliage()* (1. - foliage_to_soil) + |
425 | mTotalToExtern.addBiomass(tree->biomassFoliage()* (1. - foliage_to_soil) + |
426 | branch_biomass*(1. - branch_to_snag - branch_to_soil) + |
426 | branch_biomass*(1. - branch_to_snag - branch_to_soil) + |
427 | tree->biomassStem()*(1. - stem_to_snag - stem_to_soil), species->cnWood()); |
427 | tree->biomassStem()*(1. - stem_to_snag - stem_to_soil), species->cnWood()); |
428 | 428 | ||
429 | }
|
429 | }
|
430 | 430 | ||
431 | 431 | ||
432 | /// after the death of the tree the five biomass compartments are processed.
|
432 | /// after the death of the tree the five biomass compartments are processed.
|
433 | void Snag::addMortality(const Tree *tree) |
433 | void Snag::addMortality(const Tree *tree) |
434 | {
|
434 | {
|
435 | addBiomassPools(tree, 1., 0., // all stem biomass goes to snag |
435 | addBiomassPools(tree, 1., 0., // all stem biomass goes to snag |
436 | 1., 0., // all branch biomass to snag |
436 | 1., 0., // all branch biomass to snag |
437 | 1.); // all foliage to soil |
437 | 1.); // all foliage to soil |
438 | 438 | ||
439 | // const Species *species = tree->species();
|
439 | // const Species *species = tree->species();
|
440 | 440 | ||
441 | // // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
|
441 | // // immediate flows: 100% of foliage, 100% of fine roots: they go to the labile pool
|
442 | // mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());
|
442 | // mLabileFlux.addBiomass(tree->biomassFoliage(), species->cnFoliage(), tree->species()->snagKyl());
|
443 | // mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
|
443 | // mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
|
444 | 444 | ||
445 | // // branches and coarse roots are equally distributed over five years:
|
445 | // // branches and coarse roots are equally distributed over five years:
|
446 | // double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;
|
446 | // double biomass_rest = (tree->biomassBranch()+tree->biomassCoarseRoot()) * 0.2;
|
447 | // for (int i=0;i<5; i++)
|
447 | // for (int i=0;i<5; i++)
|
448 | // mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
|
448 | // mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
|
449 | 449 | ||
450 | // // just for book-keeping: keep track of all inputs into branches / roots / swd
|
450 | // // just for book-keeping: keep track of all inputs into branches / roots / swd
|
451 | // mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood());
|
451 | // mTotalIn.addBiomass(tree->biomassBranch() + tree->biomassCoarseRoot() + tree->biomassStem(), species->cnWood());
|
452 | // // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
452 | // // stem biomass is transferred to the standing woody debris pool (SWD), increase stem number of pool
|
453 | // int pi = poolIndex(tree->dbh()); // get right transfer pool
|
453 | // int pi = poolIndex(tree->dbh()); // get right transfer pool
|
454 | 454 | ||
455 | // // update statistics - stemnumber-weighted averages
|
455 | // // update statistics - stemnumber-weighted averages
|
456 | // // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
456 | // // note: here the calculations are repeated for every died trees (i.e. consecutive weighting ... but delivers the same results)
|
457 | // double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)
|
457 | // double p_old = mNumberOfSnags[pi] / (mNumberOfSnags[pi] + 1); // weighting factor for state vars (based on stem numbers)
|
458 | // double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).
|
458 | // double p_new = 1. / (mNumberOfSnags[pi] + 1); // weighting factor for added tree (p_old + p_new = 1).
|
459 | // mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;
|
459 | // mAvgDbh[pi] = mAvgDbh[pi]*p_old + tree->dbh()*p_new;
|
460 | // mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;
|
460 | // mAvgHeight[pi] = mAvgHeight[pi]*p_old + tree->height()*p_new;
|
461 | // mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;
|
461 | // mAvgVolume[pi] = mAvgVolume[pi]*p_old + tree->volume()*p_new;
|
462 | // mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;
|
462 | // mTimeSinceDeath[pi] = mTimeSinceDeath[pi]*p_old + 1.*p_new;
|
463 | // mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new;
|
463 | // mHalfLife[pi] = mHalfLife[pi]*p_old + species->snagHalflife()* p_new;
|
464 | 464 | ||
465 | // // average the decay rate (ksw); this is done based on the carbon content
|
465 | // // average the decay rate (ksw); this is done based on the carbon content
|
466 | // // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
466 | // // aggregate all trees that die in the current year (and save weighted decay rates to CurrentKSW)
|
467 | // if (tree->biomassStem()==0)
|
467 | // if (tree->biomassStem()==0)
|
468 | // throw IException("Snag::addMortality: tree without stem biomass!!");
|
468 | // throw IException("Snag::addMortality: tree without stem biomass!!");
|
469 | // p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
|
469 | // p_old = mToSWD[pi].C / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
|
470 | // p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
|
470 | // p_new =tree->biomassStem()* biomassCFraction / (mToSWD[pi].C + tree->biomassStem()* biomassCFraction);
|
471 | // mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new;
|
471 | // mCurrentKSW[pi] = mCurrentKSW[pi]*p_old + species->snagKsw() * p_new;
|
472 | // mNumberOfSnags[pi]++;
|
472 | // mNumberOfSnags[pi]++;
|
473 | 473 | ||
474 | // // finally add the biomass
|
474 | // // finally add the biomass
|
475 | // CNPool &to_swd = mToSWD[pi];
|
475 | // CNPool &to_swd = mToSWD[pi];
|
476 | // to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());
|
476 | // to_swd.addBiomass(tree->biomassStem(), species->cnWood(), tree->species()->snagKyr());
|
477 | }
|
477 | }
|
478 | 478 | ||
479 | /// add residual biomass of 'tree' after harvesting.
|
479 | /// add residual biomass of 'tree' after harvesting.
|
480 | /// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation [0..1] (i.e.: not to stay in the system)
|
480 | /// remove_{stem, branch, foliage}_fraction: percentage of biomass compartment that is *removed* by the harvest operation [0..1] (i.e.: not to stay in the system)
|
481 | /// records on harvested biomass is collected (mTotalToExtern-pool).
|
481 | /// records on harvested biomass is collected (mTotalToExtern-pool).
|
482 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
482 | void Snag::addHarvest(const Tree* tree, const double remove_stem_fraction, const double remove_branch_fraction, const double remove_foliage_fraction ) |
483 | {
|
483 | {
|
484 | addBiomassPools(tree,
|
484 | addBiomassPools(tree,
|
485 | 0., 1.-remove_stem_fraction, // "remove_stem_fraction" is removed -> the rest goes to soil |
485 | 0., 1.-remove_stem_fraction, // "remove_stem_fraction" is removed -> the rest goes to soil |
486 | 0., 1.-remove_branch_fraction, // "remove_branch_fraction" is removed -> the rest goes directly to the soil |
486 | 0., 1.-remove_branch_fraction, // "remove_branch_fraction" is removed -> the rest goes directly to the soil |
487 | 1.-remove_foliage_fraction); // the rest of foliage is routed to the soil |
487 | 1.-remove_foliage_fraction); // the rest of foliage is routed to the soil |
488 | // const Species *species = tree->species();
|
488 | // const Species *species = tree->species();
|
489 | 489 | ||
490 | // // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
|
490 | // // immediate flows: 100% of residual foliage, 100% of fine roots: they go to the labile pool
|
491 | // mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());
|
491 | // mLabileFlux.addBiomass(tree->biomassFoliage() * (1. - remove_foliage_fraction), species->cnFoliage(), tree->species()->snagKyl());
|
492 | // mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
|
492 | // mLabileFlux.addBiomass(tree->biomassFineRoot(), species->cnFineroot(), tree->species()->snagKyl());
|
493 | 493 | ||
494 | // // for branches, add all biomass that remains in the forest to the soil
|
494 | // // for branches, add all biomass that remains in the forest to the soil
|
495 | // mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());
|
495 | // mRefractoryFlux.addBiomass(tree->biomassBranch()*(1.-remove_branch_fraction), species->cnWood(), tree->species()->snagKyr());
|
496 | // // the same treatment for stem residuals
|
496 | // // the same treatment for stem residuals
|
497 | // mRefractoryFlux.addBiomass(tree->biomassStem() * (1. - remove_stem_fraction), species->cnWood(), tree->species()->snagKyr());
|
497 | // mRefractoryFlux.addBiomass(tree->biomassStem() * (1. - remove_stem_fraction), species->cnWood(), tree->species()->snagKyr());
|
498 | 498 | ||
499 | // // split the corase wood biomass into parts (slower decay)
|
499 | // // split the corase wood biomass into parts (slower decay)
|
500 | // double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;
|
500 | // double biomass_rest = (tree->biomassCoarseRoot()) * 0.2;
|
501 | // for (int i=0;i<5; i++)
|
501 | // for (int i=0;i<5; i++)
|
502 | // mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
|
502 | // mOtherWood[i].addBiomass(biomass_rest, species->cnWood(), tree->species()->snagKyr());
|
503 | 503 | ||
504 | 504 | ||
505 | // // for book-keeping...
|
505 | // // for book-keeping...
|
506 | // mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +
|
506 | // mTotalToExtern.addBiomass(tree->biomassFoliage()*remove_foliage_fraction +
|
507 | // tree->biomassBranch()*remove_branch_fraction +
|
507 | // tree->biomassBranch()*remove_branch_fraction +
|
508 | // tree->biomassStem()*remove_stem_fraction, species->cnWood());
|
508 | // tree->biomassStem()*remove_stem_fraction, species->cnWood());
|
509 | }
|
509 | }
|
510 | 510 | ||
511 | // add flow from regeneration layer (dead trees) to soil
|
511 | // add flow from regeneration layer (dead trees) to soil
|
512 | void Snag::addToSoil(const Species *species, const CNPair &woody_pool, const CNPair &litter_pool) |
512 | void Snag::addToSoil(const Species *species, const CNPair &woody_pool, const CNPair &litter_pool) |
513 | {
|
513 | {
|
514 | mLabileFlux.add(litter_pool, species->snagKyl()); |
514 | mLabileFlux.add(litter_pool, species->snagKyl()); |
515 | mRefractoryFlux.add(woody_pool, species->snagKyr()); |
515 | mRefractoryFlux.add(woody_pool, species->snagKyr()); |
516 | DBGMODE(
|
516 | DBGMODE(
|
517 | if (isnan(mLabileFlux.C) || isnan(mRefractoryFlux.C)) |
517 | if (isnan(mLabileFlux.C) || isnan(mRefractoryFlux.C)) |
518 | qDebug("Snag::addToSoil: NaN in C Pool"); |
518 | qDebug("Snag::addToSoil: NaN in C Pool"); |
519 | ); |
519 | ); |
520 | }
|
520 | }
|
521 | 521 | ||
522 | /// disturbance function: remove the fraction of 'factor' of biomass from the SWD pools; 0: remove nothing, 1: remove all
|
522 | /// disturbance function: remove the fraction of 'factor' of biomass from the SWD pools; 0: remove nothing, 1: remove all
|
523 | /// biomass removed by this function goes to the atmosphere
|
523 | /// biomass removed by this function goes to the atmosphere
|
524 | void Snag::removeCarbon(const double factor) |
524 | void Snag::removeCarbon(const double factor) |
525 | {
|
525 | {
|
526 | // reduce pools of currently standing dead wood and also of pools that are added
|
526 | // reduce pools of currently standing dead wood and also of pools that are added
|
527 | // during (previous) management operations of the current year
|
527 | // during (previous) management operations of the current year
|
528 | for (int i=0;i<3;i++) { |
528 | for (int i=0;i<3;i++) { |
529 | mTotalToDisturbance += (mSWD[i] + mToSWD[i]) * factor; |
529 | mTotalToDisturbance += (mSWD[i] + mToSWD[i]) * factor; |
530 | mSWD[i] *= (1. - factor); |
530 | mSWD[i] *= (1. - factor); |
531 | mToSWD[i] *= (1. - factor); |
531 | mToSWD[i] *= (1. - factor); |
532 | }
|
532 | }
|
533 | 533 | ||
534 | for (int i=0;i<5;i++) { |
534 | for (int i=0;i<5;i++) { |
535 | mTotalToDisturbance += mOtherWood[i]*factor; |
535 | mTotalToDisturbance += mOtherWood[i]*factor; |
536 | mOtherWood[i] *= (1. - factor); |
536 | mOtherWood[i] *= (1. - factor); |
537 | }
|
537 | }
|
538 | }
|
538 | }
|
539 | 539 | ||
540 | 540 | ||
541 | /// cut down swd (and branches) and move to soil pools
|
541 | /// cut down swd (and branches) and move to soil pools
|
542 | /// @param factor 0: cut 0%, 1: cut and slash 100% of the wood
|
542 | /// @param factor 0: cut 0%, 1: cut and slash 100% of the wood
|
543 | void Snag::management(const double factor) |
543 | void Snag::management(const double factor) |
544 | {
|
544 | {
|
545 | if (factor<0. || factor>1.) |
545 | if (factor<0. || factor>1.) |
546 | throw IException(QString("Invalid factor in Snag::management: '%1'").arg(factor)); |
546 | throw IException(QString("Invalid factor in Snag::management: '%1'").arg(factor)); |
547 | // swd pools
|
547 | // swd pools
|
548 | for (int i=0;i<3;i++) { |
548 | for (int i=0;i<3;i++) { |
549 | mSWDtoSoil += mSWD[i] * factor; |
549 | mSWDtoSoil += mSWD[i] * factor; |
550 | mRefractoryFlux += mSWD[i] * factor; |
550 | mRefractoryFlux += mSWD[i] * factor; |
551 | mSWD[i] *= (1. - factor); |
551 | mSWD[i] *= (1. - factor); |
552 | //mSWDtoSoil += mToSWD[i] * factor;
|
552 | //mSWDtoSoil += mToSWD[i] * factor;
|
553 | //mToSWD[i] *= (1. - factor);
|
553 | //mToSWD[i] *= (1. - factor);
|
554 | }
|
554 | }
|
555 | // what to do with the branches: now move also all wood to soil (note: this is note
|
555 | // what to do with the branches: now move also all wood to soil (note: this is note
|
556 | // very good w.r.t the coarse roots...
|
556 | // very good w.r.t the coarse roots...
|
557 | for (int i=0;i<5;i++) { |
557 | for (int i=0;i<5;i++) { |
558 | mRefractoryFlux+=mOtherWood[i]*factor; |
558 | mRefractoryFlux+=mOtherWood[i]*factor; |
559 | mOtherWood[i]*=(1. - factor); |
559 | mOtherWood[i]*=(1. - factor); |
560 | }
|
560 | }
|
561 | 561 | ||
562 | }
|
562 | }
|
563 | 563 | ||
564 | 564 | ||
565 | 565 |