growth.cpp 54 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604
  1. ///////////////////////////////////////////////////////////////////////////////////////
  2. /// \file growth.cpp
  3. /// \brief The growth module
  4. ///
  5. /// Vegetation C allocation, litter production, tissue turnover
  6. /// leaf phenology, allometry and growth
  7. ///
  8. /// (includes updated FPC formulation as required for "fast"
  9. /// cohort/individual mode - see canexch.cpp)
  10. ///
  11. /// \author Ben Smith
  12. /// $Date: 2013-10-21 14:55:54 +0200 (Mon, 21 Oct 2013) $
  13. ///
  14. ///////////////////////////////////////////////////////////////////////////////////////
  15. // WHAT SHOULD THIS FILE CONTAIN?
  16. // Module source code files should contain, in this order:
  17. // (1) a "#include" directive naming the framework header file. The framework header
  18. // file should define all classes used as arguments to functions in the present
  19. // module. It may also include declarations of global functions, constants and
  20. // types, accessible throughout the model code;
  21. // (2) other #includes, including header files for other modules accessed by the
  22. // present one;
  23. // (3) type definitions, constants and file scope global variables for use within
  24. // the present module only;
  25. // (4) declarations of functions defined in this file, if needed;
  26. // (5) definitions of all functions. Functions that are to be accessible to other
  27. // modules or to the calling framework should be declared in the module header
  28. // file.
  29. //
  30. // PORTING MODULES BETWEEN FRAMEWORKS:
  31. // Modules should be structured so as to be fully portable between models (frameworks).
  32. // When porting between frameworks, the only change required should normally be in the
  33. // "#include" directive referring to the framework header file.
  34. #include "config.h"
  35. #include "growth.h"
  36. #include "canexch.h"
  37. #include "landcover.h"
  38. #include <assert.h>
  39. ///////////////////////////////////////////////////////////////////////////////////////
  40. // FILE SCOPE GLOBAL CONSTANTS
  41. const double APHEN_MAX = 210.0;
  42. // Maximum number of equivalent days with full leaf cover per growing season
  43. // for summergreen PFTs
  44. ///////////////////////////////////////////////////////////////////////////////////////
  45. // LEAF PHENOLOGY
  46. // Call function leaf_phenology each simulation day prior to calculation of FPAR, to
  47. // calculate fractional leaf-out for each PFT and individual.
  48. // Function leaf_phenology_pft is not intended to be called directly by the framework,
  49. void leaf_phenology_pft(Pft& pft, Climate& climate, double wscal, double aphen,
  50. double& phen) {
  51. // DESCRIPTION
  52. // Calculates leaf phenological status (fractional leaf-out) for a individuals of
  53. // a given PFT, given current heat sum and length of chilling period (summergreen
  54. // PFTs) and water stress coefficient (raingreen PFTs)
  55. // INPUT PARAMETER
  56. // wscal = water stress coefficient (0-1; 1=minimum stress)
  57. // aphen = sum of daily fractional leaf cover (equivalent number of days with
  58. // full leaf cover) so far this growing season
  59. // OUTPUT PARAMETER
  60. // phen = fraction of full leaf cover for any individual of this PFT
  61. bool raingreen = pft.phenology == RAINGREEN || pft.phenology == ANY;
  62. bool summergreen = pft.phenology == SUMMERGREEN || pft.phenology == ANY;
  63. phen = 1.0;
  64. if (summergreen) {
  65. // Summergreen PFT - phenology based on GDD5 sum
  66. if (pft.lifeform == TREE) {
  67. // Calculate GDD base value for this PFT (if not already known) given
  68. // current length of chilling period (Sykes et al 1996, Eqn 1)
  69. if (pft.gdd0[climate.chilldays] < 0.0)
  70. pft.gdd0[climate.chilldays] = pft.k_chilla +
  71. pft.k_chillb * exp(-pft.k_chillk * (double)climate.chilldays);
  72. if (climate.gdd5 > pft.gdd0[climate.chilldays] && aphen < APHEN_MAX)
  73. phen = min(1.0,
  74. (climate.gdd5 - pft.gdd0[climate.chilldays]) / pft.phengdd5ramp);
  75. else
  76. phen = 0.0;
  77. }
  78. else if (pft.lifeform == GRASS) {
  79. // Summergreen grasses have no maximum number of leaf-on days per
  80. // growing season, and no chilling requirement
  81. phen = min(1.0, climate.gdd5 / pft.phengdd5ramp);
  82. }
  83. }
  84. if (raingreen && wscal < pft.wscal_min) {
  85. // Raingreen phenology based on water stress threshold
  86. phen = 0.0;
  87. }
  88. }
  89. void leaf_phenology(Patch& patch, Climate& climate) {
  90. // DESCRIPTION
  91. // Updates leaf phenological status (fractional leaf-out) for Patch PFT objects and
  92. // all individuals in a particular patch.
  93. // Updated by Ben Smith 2002-07-24 for compatability with "fast" canopy exchange
  94. // code (phenology assigned to patchpft for all vegetation modes)
  95. // Obtain reference to Vegetation object
  96. Vegetation& vegetation = patch.vegetation;
  97. // INDIVIDUAL AND COHORT MODES
  98. // Calculate phenology for each PFT at this patch
  99. // Loop through patch-PFTs
  100. patch.pft.firstobj();
  101. while (patch.pft.isobj) {
  102. Patchpft& ppft = patch.pft.getobj();
  103. // For this PFT ...
  104. if(patch.stand.pft[ppft.id].active) {
  105. if(ppft.pft.landcover == CROPLAND && patch.stand.landcover == CROPLAND)
  106. leaf_phenology_crop(ppft.pft, patch);
  107. else //natural, urban, pasture, forest and peatland stands/pft:s
  108. leaf_phenology_pft(ppft.pft, climate, ppft.wscal, ppft.aphen, ppft.phen);
  109. // Update annual leaf-on sum
  110. if ( (climate.lat >= 0.0 && date.day == COLDEST_DAY_NHEMISPHERE) ||
  111. (climate.lat < 0.0 && date.day == COLDEST_DAY_SHEMISPHERE) ) {
  112. ppft.aphen = 0.0;
  113. }
  114. ppft.aphen += ppft.phen;
  115. // ... on to next PFT
  116. }
  117. patch.pft.nextobj();
  118. }
  119. // Copy PFT-specific phenological status to individuals of each PFT
  120. // Loop through individuals
  121. vegetation.firstobj();
  122. while (vegetation.isobj) {
  123. Individual& indiv = vegetation.getobj();
  124. // For this individual ...
  125. indiv.phen = patch.pft[indiv.pft.id].phen;
  126. // Update annual leaf-day sum (raingreen PFTs)
  127. if (date.day == 0) indiv.aphen_raingreen = 0;
  128. indiv.aphen_raingreen += (indiv.phen != 0.0);
  129. // ... on to next individual
  130. vegetation.nextobj();
  131. }
  132. }
  133. /// Calculates nitrogen retranslocation fraction
  134. /* Calculates actual nitrogen retranslocation fraction so maximum
  135. * nitrogen storage capacity is not exceeded
  136. */
  137. double calc_nrelocfrac(lifeformtype lifeform, double turnover_leaf, double nmass_leaf,
  138. double turnover_root, double nmass_root, double turnover_sap, double nmass_sap, double max_n_storage, double longterm_nstore) {
  139. double turnover_nmass = turnover_leaf * nmass_leaf + turnover_root * nmass_root;
  140. if (lifeform == TREE)
  141. turnover_nmass += turnover_sap * nmass_sap;
  142. if (max_n_storage < longterm_nstore)
  143. return 0.0;
  144. else if (max_n_storage < longterm_nstore + turnover_nmass * nrelocfrac && !negligible(turnover_nmass))
  145. return (max_n_storage - longterm_nstore) / (turnover_nmass);
  146. else
  147. return nrelocfrac;
  148. }
  149. ///////////////////////////////////////////////////////////////////////////////////////
  150. // TURNOVER
  151. // Internal function (do not call directly from framework)
  152. void turnover(double turnover_leaf, double turnover_root, double turnover_sap,
  153. lifeformtype lifeform, landcovertype landcover, double& cmass_leaf, double& cmass_root, double& cmass_sap,
  154. double& cmass_heart, double& nmass_leaf, double& nmass_root, double& nmass_sap,
  155. double& nmass_heart, double& litter_leaf, double& litter_root,
  156. double& nmass_litter_leaf, double& nmass_litter_root,
  157. double& longterm_nstore, double &max_n_storage,
  158. bool alive) {
  159. // DESCRIPTION
  160. // Transfers carbon from leaves and roots to litter, and from sapwood to heartwood
  161. // Only turnover from 'alive' individuals is transferred to litter (Ben 2007-11-28)
  162. // INPUT PARAMETERS
  163. // turnover_leaf = leaf turnover per time period as a proportion of leaf C biomass
  164. // turnover_root = root turnover per time period as a proportion of root C biomass
  165. // turnover_sap = sapwood turnover to heartwood per time period as a proportion of
  166. // sapwood C biomass
  167. // lifeform = PFT life form class (TREE or GRASS)
  168. // alive = signifies new Individual object if false (see vegdynam.cpp)
  169. // INPUT AND OUTPUT PARAMETERS
  170. // cmass_leaf = leaf C biomass (kgC/m2)
  171. // cmass_root = fine root C biomass (kgC/m2)
  172. // cmass_sap = sapwood C biomass (kgC/m2)
  173. // nmass_leaf = leaf nitrogen biomass (kgN/m2)
  174. // nmass_root = fine root nitrogen biomass (kgN/m2)
  175. // nmass_sap = sapwood nitrogen biomass (kgN/m2)
  176. // OUTPUT PARAMETERS
  177. // litter_leaf = new leaf C litter (kgC/m2)
  178. // litter_root = new root C litter (kgC/m2)
  179. // nmass_litter_leaf = new leaf nitrogen litter (kgN/m2)
  180. // nmass_litter_root = new root nitrogen litter (kgN/m2)
  181. // cmass_heart = heartwood C biomass (kgC/m2)
  182. // nmass_heart = heartwood nitrogen biomass (kgC/m2)
  183. // longterm_nstore = longterm nitrogen storage (kgN/m2)
  184. double turnover = 0.0;
  185. // Calculate actual nitrogen retranslocation so maximum nitrogen storage capacity is not exceeded
  186. double actual_nrelocfrac = calc_nrelocfrac(lifeform, turnover_leaf, nmass_leaf, turnover_root, nmass_root,
  187. turnover_sap, nmass_sap, max_n_storage, longterm_nstore);
  188. // TREES AND GRASSES:
  189. // Leaf turnover
  190. turnover = turnover_leaf * cmass_leaf;
  191. cmass_leaf -= turnover;
  192. if (alive) litter_leaf += turnover;
  193. turnover = turnover_leaf * nmass_leaf;
  194. nmass_leaf -= turnover;
  195. nmass_litter_leaf += turnover * (1.0 - actual_nrelocfrac);
  196. longterm_nstore += turnover * actual_nrelocfrac;
  197. // Root turnover
  198. turnover = turnover_root * cmass_root;
  199. cmass_root -= turnover;
  200. if (alive) litter_root += turnover;
  201. turnover = turnover_root * nmass_root;
  202. nmass_root -= turnover;
  203. nmass_litter_root += turnover * (1.0 - actual_nrelocfrac);
  204. longterm_nstore += turnover * actual_nrelocfrac;
  205. if (lifeform == TREE) {
  206. // TREES ONLY:
  207. // Sapwood turnover by conversion to heartwood
  208. turnover = turnover_sap * cmass_sap;
  209. cmass_sap -= turnover;
  210. cmass_heart += turnover;
  211. // NB: assumes nitrogen is translocated from sapwood prior to conversion to
  212. // heartwood and that this is the same fraction that is conserved
  213. // in conjunction with leaf and root shedding
  214. turnover = turnover_sap * nmass_sap;
  215. nmass_sap -= turnover;
  216. nmass_heart += turnover * (1.0 - actual_nrelocfrac);
  217. longterm_nstore += turnover * actual_nrelocfrac;
  218. }
  219. }
  220. ///////////////////////////////////////////////////////////////////////////////////////
  221. // REPRODUCTION
  222. // Internal function (do not call directly from framework)
  223. void reproduction(double reprfrac, double npp, double& bminc, double& cmass_repr) {
  224. // DESCRIPTION
  225. // Allocation of net primary production (NPP) to reproduction and calculation of
  226. // assimilated carbon available for production of new biomass
  227. // INPUT PARAMETERS
  228. // reprfrac = fraction of NPP for this time period allocated to reproduction
  229. // npp = NPP (i.e. assimilation minus maintenance and growth respiration) for
  230. // this time period (kgC/m2)
  231. // OUTPUT PARAMETER
  232. // bminc = carbon biomass increment (component of NPP available for production
  233. // of new biomass) for this time period (kgC/m2)
  234. if (npp >= 0.0) {
  235. cmass_repr = npp * reprfrac;
  236. bminc = npp - cmass_repr;
  237. return;
  238. }
  239. // Negative NPP - no reproduction cost
  240. cmass_repr = 0.0;
  241. bminc = npp;
  242. }
  243. ///////////////////////////////////////////////////////////////////////////////////////
  244. // ALLOCATION
  245. // Function allocation is an internal function (do not call directly from framework);
  246. // function allocation_init may be called to distribute initial biomass among tissues
  247. // for a new individual.
  248. // File scope global variables: used by function f below (see function allocation)
  249. static double k1, k2, k3, b;
  250. static double ltor_g;
  251. static double cmass_heart_g;
  252. static double cmass_leaf_g;
  253. inline double f(double& cmass_leaf_inc) {
  254. // Returns value of f(cmass_leaf_inc), given by:
  255. //
  256. // f(cmass_leaf_inc) = 0 =
  257. // k1 * (b - cmass_leaf_inc - cmass_leaf_inc/ltor + cmass_heart) -
  258. // [ (b - cmass_leaf_inc - cmass_leaf_inc/ltor)
  259. // / (cmass_leaf + cmass_leaf_inc )*k3 ] ** k2
  260. //
  261. // See function allocation (below), Eqn (13)
  262. return k1 * (b - cmass_leaf_inc - cmass_leaf_inc / ltor_g + cmass_heart_g) -
  263. pow((b - cmass_leaf_inc - cmass_leaf_inc / ltor_g) / (cmass_leaf_g + cmass_leaf_inc) * k3,
  264. k2);
  265. }
  266. void allocation(double bminc,double cmass_leaf,double cmass_root,double cmass_sap,
  267. double cmass_debt,double cmass_heart,double ltor,double height,double sla,
  268. double wooddens,lifeformtype lifeform,double k_latosa,double k_allom2,
  269. double k_allom3,double& cmass_leaf_inc,double& cmass_root_inc,
  270. double& cmass_sap_inc,
  271. double& cmass_debt_inc,
  272. double& cmass_heart_inc,double& litter_leaf_inc,
  273. double& litter_root_inc,double& exceeds_cmass) {
  274. // DESCRIPTION
  275. // Calculates changes in C compartment sizes (leaves, roots, sapwood, heartwood)
  276. // and litter for a plant individual as a result of allocation of biomass increment.
  277. // Assumed allometric relationships are given in function allometry below.
  278. // INPUT PARAMETERS
  279. // bminc = biomass increment this time period on individual basis (kgC)
  280. // cmass_leaf = leaf C biomass for last time period on individual basis (kgC)
  281. // cmass_root = root C biomass for last time period on individual basis (kgC)
  282. // cmass_sap = sapwood C biomass for last time period on individual basis (kgC)
  283. // cmass_heart = heartwood C biomass for last time period on individual basis (kgC)
  284. // ltor = leaf to root mass ratio following allocation
  285. // height = individual height (m)
  286. // sla = specific leaf area (PFT-specific constant) (m2/kgC)
  287. // wooddens = wood density (PFT-specific constant) (kgC/m3)
  288. // lifeform = life form class (TREE or GRASS)
  289. // k_latosa = ratio of leaf area to sapwood cross-sectional area (PFT-specific
  290. // constant)
  291. // k_allom2 = constant in allometry equations
  292. // k_allom3 = constant in allometry equations
  293. // OUTPUT PARAMETERS
  294. // cmass_leaf_inc = increment (may be negative) in leaf C biomass following
  295. // allocation (kgC)
  296. // cmass_root_inc = increment (may be negative) in root C biomass following
  297. // allocation (kgC)
  298. // cmass_sap_inc = increment (may be negative) in sapwood C biomass following
  299. // allocation (kgC)
  300. // cmass_heart_inc = increment in heartwood C biomass following allocation (kgC)
  301. // litter_leaf_inc = increment in leaf litter following allocation, on individual
  302. // basis (kgC)
  303. // litter_root_inc = increment in root litter following allocation, on individual
  304. // basis (kgC)
  305. // exceeds_cmass = negative increment that exceeds existing biomass (kgC)
  306. // MATHEMATICAL DERIVATION FOR TREE ALLOCATION
  307. // Allocation attempts to distribute biomass increment (bminc) among the living
  308. // tissue compartments, i.e.
  309. // (1) bminc = cmass_leaf_inc + cmass_root_inc + cmass_sap_inc
  310. // while satisfying the allometric relationships (Shinozaki et al. 1964a,b; Waring
  311. // et al 1982, Huang et al 1992; see also function allometry, below) [** =
  312. // raised to the power of]:
  313. // (2) (leaf area) = k_latosa * (sapwood xs area)
  314. // (3) cmass_leaf = ltor * cmass_root
  315. // (4) height = k_allom2 * (stem diameter) ** k_allom3
  316. // From (1) and (3),
  317. // (5) cmass_sap_inc = bminc - cmass_leaf_inc -
  318. // (cmass_leaf + cmass_leaf_inc) / ltor + cmass_root
  319. // Let diam_new and height_new be stem diameter and height following allocation.
  320. // Then (see allometry),
  321. // (6) diam_new = 2 * [ ( cmass_sap + cmass_sap_inc + cmass_heart )
  322. // / wooddens / height_new / PI ]**(1/2)
  323. // From (4), (6) and (5),
  324. // (7) height_new**(1+2/k_allom3) =
  325. // k_allom2**(2/k_allom3) * 4 * [cmass_sap + bminc - cmass_leaf_inc
  326. // - (cmass_leaf + cmass_leaf_inc) / ltor + cmass_root + cmass_heart]
  327. // / wooddens / PI
  328. // Now,
  329. // (8) wooddens = cmass_sap / height / (sapwood xs area)
  330. // From (8) and (2),
  331. // (9) wooddens = cmass_sap / height / sla / cmass_leaf * k_latosa
  332. // From (9) and (1),
  333. // (10) wooddens = (cmass_sap + bminc - cmass_leaf_inc -
  334. // (cmass_leaf + cmass_leaf_inc) / ltor + cmass_root)
  335. // / height_new / sla / (cmass_leaf + cmass_leaf_inc) * k_latosa
  336. // From (10),
  337. // (11) height_new**(1+2/k_allom3) =
  338. // [ (cmass_sap + bminc - cmass_leaf_inc - (cmass_leaf + cmass_leaf_inc)
  339. // / ltor + cmass_root) / wooddens / sla
  340. // / (cmass_leaf + cmass_leaf_inc ) * k_latosa ] ** (1+2/k_allom3)
  341. //
  342. // Combining (7) and (11) gives a function of the unknown cmass_leaf_inc:
  343. //
  344. // (12) f(cmass_leaf_inc) = 0 =
  345. // k_allom2**(2/k_allom3) * 4/PI * [cmass_sap + bminc - cmass_leaf_inc
  346. // - (cmass_leaf + cmass_leaf_inc) / ltor + cmass_root + cmass_heart]
  347. // / wooddens -
  348. // [ (cmass_sap + bminc - cmass_leaf_inc - (cmass_leaf + cmass_leaf_inc)
  349. // / ltor + cmass_root) / (cmass_leaf + cmass_leaf_inc)
  350. // / wooddens / sla * k_latosa] ** (1+2/k_allom3)
  351. //
  352. // Let k1 = k_allom2**(2/k_allom3) * 4/PI / wooddens
  353. // k2 = 1+2/k_allom3
  354. // k3 = k_latosa / wooddens / sla
  355. // b = cmass_sap + bminc - cmass_leaf/ltor + cmass_root
  356. //
  357. // Then,
  358. // (13) f(cmass_leaf_inc) = 0 =
  359. // k1 * (b - cmass_leaf_inc - cmass_leaf_inc/ltor + cmass_heart) -
  360. // [ (b - cmass_leaf_inc - cmass_leaf_inc/ltor)
  361. // / (cmass_leaf + cmass_leaf_inc )*k3 ] ** k2
  362. //
  363. // Numerical methods are used to solve Eqn (13) for cmass_leaf_inc
  364. const int NSEG=20; // number of segments (parameter in numerical methods)
  365. const int JMAX=40; // maximum number of iterations (in numerical methods)
  366. const double XACC=0.0001; // threshold x-axis precision of allocation solution
  367. const double YACC=1.0e-10; // threshold y-axis precision of allocation solution
  368. const double CDEBT_MAXLOAN_DEFICIT=0.8; // maximum loan as a fraction of deficit
  369. const double CDEBT_MAXLOAN_MASS=0.2; // maximum loan as a fraction of (sapwood-cdebt)
  370. double cmass_leaf_inc_min;
  371. double cmass_root_inc_min;
  372. double x1,x2,dx,xmid,fx1,fmid,rtbis,sign;
  373. int j;
  374. double cmass_deficit,cmass_loan;
  375. // initialise
  376. litter_leaf_inc = 0.0;
  377. litter_root_inc = 0.0;
  378. exceeds_cmass = 0.0;
  379. cmass_leaf_inc = 0.0;
  380. cmass_root_inc = 0.0;
  381. cmass_sap_inc = 0.0;
  382. cmass_heart_inc = 0.0;
  383. cmass_debt_inc = 0.0;
  384. if (!largerthanzero(ltor, -10)) {
  385. // No leaf production possible - put all biomass into roots
  386. // (Individual will die next time period)
  387. cmass_leaf_inc = 0.0;
  388. // Make sure we don't end up with negative cmass_root
  389. if (bminc < -cmass_root) {
  390. exceeds_cmass = -(cmass_root + bminc);
  391. cmass_root_inc = -cmass_root;
  392. }
  393. else {
  394. cmass_root_inc=bminc;
  395. }
  396. if (lifeform==TREE) {
  397. cmass_sap_inc=-cmass_sap;
  398. cmass_heart_inc=-cmass_sap_inc;
  399. }
  400. }
  401. else if (lifeform==TREE) {
  402. // TREE ALLOCATION
  403. cmass_heart_inc=0.0;
  404. // Calculate minimum leaf increment to maintain current sapwood biomass
  405. // Given Eqn (2)
  406. if (height>0.0)
  407. cmass_leaf_inc_min=k_latosa*cmass_sap/(wooddens*height*sla)-cmass_leaf;
  408. else
  409. cmass_leaf_inc_min=0.0;
  410. // Calculate minimum root increment to support minimum resulting leaf biomass
  411. // Eqn (3)
  412. if (height>0.0)
  413. cmass_root_inc_min=k_latosa*cmass_sap/(wooddens*height*sla*ltor)-
  414. cmass_root;
  415. else
  416. cmass_root_inc_min=0.0;
  417. if (cmass_root_inc_min<0.0) { // some roots would have to be killed
  418. cmass_leaf_inc_min=cmass_root*ltor-cmass_leaf;
  419. cmass_root_inc_min=0.0;
  420. }
  421. // BLARP! C debt stuff
  422. if (ifcdebt) {
  423. cmass_deficit=cmass_leaf_inc_min+cmass_root_inc_min-bminc;
  424. if (cmass_deficit>0.0) {
  425. cmass_loan=max(min(cmass_deficit*CDEBT_MAXLOAN_DEFICIT,
  426. (cmass_sap-cmass_debt)*CDEBT_MAXLOAN_MASS),0.0);
  427. bminc+=cmass_loan;
  428. cmass_debt_inc=cmass_loan;
  429. }
  430. else cmass_debt_inc=0.0;
  431. }
  432. else cmass_debt_inc=0.0;
  433. if ( (cmass_root_inc_min >= 0.0 && cmass_leaf_inc_min >= 0.0 &&
  434. cmass_root_inc_min + cmass_leaf_inc_min <= bminc) || bminc<=0.0) {
  435. // Normal allocation (positive increment to all living C compartments)
  436. // Calculation of leaf mass increment (lminc_ind) satisfying Eqn (13)
  437. // using bisection method (Press et al 1986)
  438. // Set values for global variables for reuse by function f
  439. k1 = pow(k_allom2, 2.0 / k_allom3) * 4.0 / PI / wooddens;
  440. k2 = 1.0 + 2 / k_allom3;
  441. k3 = k_latosa / wooddens / sla;
  442. b = cmass_sap + bminc - cmass_leaf / ltor + cmass_root;
  443. ltor_g = ltor;
  444. cmass_leaf_g = cmass_leaf;
  445. cmass_heart_g = cmass_heart;
  446. x1 = 0.0;
  447. x2 = (bminc - (cmass_leaf / ltor - cmass_root)) / (1.0 + 1.0 / ltor);
  448. dx = (x2 - x1) / (double)NSEG;
  449. if (cmass_leaf < 1.0e-10) x1 += dx; // to avoid division by zero
  450. // Evaluate f(x1), i.e. Eqn (13) at cmass_leaf_inc = x1
  451. fx1 = f(x1);
  452. // Find approximate location of leftmost root on the interval
  453. // (x1,x2). Subdivide (x1,x2) into nseg equal segments seeking
  454. // change in sign of f(xmid) relative to f(x1).
  455. fmid = f(x1);
  456. xmid = x1;
  457. while (fmid * fx1 > 0.0 && xmid < x2) {
  458. xmid += dx;
  459. fmid = f(xmid);
  460. }
  461. x1 = xmid - dx;
  462. x2 = xmid;
  463. // Apply bisection to find root on new interval (x1,x2)
  464. if (f(x1) >= 0.0) sign = -1.0;
  465. else sign = 1.0;
  466. rtbis = x1;
  467. dx = x2 - x1;
  468. // Bisection loop
  469. // Search iterates on value of xmid until xmid lies within
  470. // xacc of the root, i.e. until |xmid-x|<xacc where f(x)=0
  471. fmid = 1.0; // dummy value to guarantee entry into loop
  472. j = 0; // number of iterations so far
  473. while (dx >= XACC && fabs(fmid) > YACC && j <= JMAX) {
  474. dx *= 0.5;
  475. xmid = rtbis + dx;
  476. fmid = f(xmid);
  477. if (fmid * sign <= 0.0) rtbis = xmid;
  478. j++;
  479. }
  480. // Now rtbis contains numerical solution for cmass_leaf_inc given Eqn (13)
  481. cmass_leaf_inc = rtbis;
  482. // Calculate increments in other compartments
  483. cmass_root_inc = (cmass_leaf_inc + cmass_leaf) / ltor - cmass_root; // Eqn (3)
  484. cmass_sap_inc = bminc - cmass_leaf_inc - cmass_root_inc; // Eqn (1)
  485. // guess2008 - extra check - abnormal allocation can still happen if ltor is very small
  486. if ((cmass_root_inc > 50 || cmass_root_inc < -50) && ltor < 0.0001) {
  487. cmass_leaf_inc = 0.0;
  488. cmass_root_inc = bminc;
  489. cmass_sap_inc = -cmass_sap;
  490. cmass_heart_inc = -cmass_sap_inc;
  491. }
  492. // Negative sapwood increment larger than existing sapwood or
  493. // if debt becomes larger than existing woody biomass
  494. if (cmass_sap < -cmass_sap_inc || cmass_sap + cmass_sap_inc + cmass_heart < cmass_debt + cmass_debt_inc) {
  495. // Abnormal allocation: reduction in some biomass compartment(s) to
  496. // satisfy allometry
  497. // Attempt to distribute this year's production among leaves and roots only
  498. // Eqn (3)
  499. cmass_leaf_inc = (bminc - cmass_leaf / ltor + cmass_root) / (1.0 + 1.0 / ltor);
  500. cmass_root_inc = bminc - cmass_leaf_inc;
  501. // Make sure we don't end up with negative cmass_leaf
  502. cmass_leaf_inc = max(-cmass_leaf, cmass_leaf_inc);
  503. // Make sure we don't end up with negative cmass_root
  504. cmass_root_inc = max(-cmass_root, cmass_root_inc);
  505. // If biomass of roots and leafs can't meet biomass decrease then
  506. // sapwood also needs to decrease
  507. cmass_sap_inc = bminc - cmass_leaf_inc - cmass_root_inc;
  508. // No sapwood turned into heartwood
  509. cmass_heart_inc = 0.0;
  510. // Make sure we don't end up with negative cmass_sap
  511. if (cmass_sap_inc < -cmass_sap) {
  512. exceeds_cmass = -(cmass_sap + cmass_sap_inc);
  513. cmass_sap_inc = -cmass_sap;
  514. }
  515. // Comment: Can happen that biomass decrease is larger than biomass in all compartments.
  516. // Then bminc is more negative than there is biomass to respire
  517. }
  518. }
  519. else {
  520. // Abnormal allocation: reduction in some biomass compartment(s) to
  521. // satisfy allometry
  522. // Attempt to distribute this year's production among leaves and roots only
  523. // Eqn (3)
  524. cmass_leaf_inc = (bminc - cmass_leaf / ltor + cmass_root) / (1.0 + 1.0 / ltor);
  525. if (cmass_leaf_inc > 0.0) {
  526. // Positive allocation to leaves
  527. cmass_root_inc = bminc - cmass_leaf_inc; // Eqn (1)
  528. // Add killed roots (if any) to litter
  529. // guess2008 - back to LPJF method in this case
  530. // if (cmass_root_inc<0.0) litter_root_inc=-cmass_root_inc;
  531. if (cmass_root_inc < 0.0) {
  532. cmass_leaf_inc = bminc;
  533. cmass_root_inc = (cmass_leaf_inc + cmass_leaf) / ltor - cmass_root; // Eqn (3)
  534. }
  535. }
  536. else {
  537. // Negative or zero allocation to leaves
  538. // Eqns (1), (3)
  539. cmass_root_inc = bminc;
  540. cmass_leaf_inc = (cmass_root + cmass_root_inc) * ltor - cmass_leaf;
  541. }
  542. // Make sure we don't end up with negative cmass_leaf
  543. if (cmass_leaf_inc < -cmass_leaf) {
  544. exceeds_cmass += -(cmass_leaf + cmass_leaf_inc);
  545. cmass_leaf_inc = -cmass_leaf;
  546. }
  547. // Make sure we don't end up with negative cmass_root
  548. if (cmass_root_inc < -cmass_root) {
  549. exceeds_cmass += -(cmass_root + cmass_root_inc);
  550. cmass_root_inc = -cmass_root;
  551. }
  552. // Add killed leaves to litter
  553. litter_leaf_inc = max(-cmass_leaf_inc, 0.0);
  554. // Add killed roots to litter
  555. litter_root_inc = max(-cmass_root_inc, 0.0);
  556. // Calculate increase in sapwood mass (which must be negative)
  557. // Eqn (2)
  558. cmass_sap_inc = (cmass_leaf_inc + cmass_leaf) * wooddens * height * sla / k_latosa -
  559. cmass_sap;
  560. // Make sure we don't end up with negative cmass_sap
  561. if (cmass_sap_inc < -cmass_sap) {
  562. exceeds_cmass += -(cmass_sap + cmass_sap_inc);
  563. cmass_sap_inc = -cmass_sap;
  564. }
  565. // Convert killed sapwood to heartwood
  566. cmass_heart_inc = -cmass_sap_inc;
  567. }
  568. }
  569. else if (lifeform == GRASS) {
  570. // GRASS ALLOCATION
  571. // Allocation attempts to distribute biomass increment (bminc) among leaf
  572. // and root compartments, i.e.
  573. // (14) bminc = cmass_leaf_inc + cmass_root_inc
  574. // while satisfying Eqn(3)
  575. cmass_leaf_inc = (bminc - cmass_leaf / ltor + cmass_root) / (1.0 + 1.0 / ltor);
  576. cmass_root_inc = bminc - cmass_leaf_inc;
  577. if (bminc >= 0.0) {
  578. // Positive biomass increment
  579. if (cmass_leaf_inc < 0.0) {
  580. // Positive bminc, but ltor causes negative allocation to leaves,
  581. // put all of bminc into roots
  582. cmass_root_inc = bminc;
  583. cmass_leaf_inc = (cmass_root + cmass_root_inc) * ltor - cmass_leaf; // Eqn (3)
  584. }
  585. else if (cmass_root_inc < 0.0) {
  586. // Positive bminc, but ltor causes negative allocation to roots,
  587. // put all of bminc into leaves
  588. cmass_leaf_inc = bminc;
  589. cmass_root_inc = (cmass_leaf + bminc) / ltor - cmass_root;
  590. }
  591. // Make sure we don't end up with negative cmass_leaf
  592. if (cmass_leaf_inc < -cmass_leaf) {
  593. exceeds_cmass += -(cmass_leaf + cmass_leaf_inc);
  594. cmass_leaf_inc = -cmass_leaf;
  595. }
  596. // Make sure we don't end up with negative cmass_root
  597. if (cmass_root_inc < -cmass_root) {
  598. exceeds_cmass += -(cmass_root + cmass_root_inc);
  599. cmass_root_inc = -cmass_root;
  600. }
  601. // Add killed leaves to litter
  602. litter_leaf_inc = max(-cmass_leaf_inc, 0.0);
  603. // Add killed roots to litter
  604. litter_root_inc = max(-cmass_root_inc, 0.0);
  605. }
  606. else if (bminc < 0) {
  607. // Abnormal allocation: negative biomass increment
  608. // Negative increment is respiration (neg bminc) or/and increment in other
  609. // compartments leading to no litter production
  610. if (bminc < -(cmass_leaf + cmass_root)) {
  611. // Biomass decrease is larger than existing biomass
  612. exceeds_cmass = -(bminc + cmass_leaf + cmass_root);
  613. cmass_leaf_inc = -cmass_leaf;
  614. cmass_root_inc = -cmass_root;
  615. }
  616. else if (cmass_root_inc < 0.0) {
  617. // Negative allocation to root
  618. // Make sure we don't end up with negative cmass_root
  619. if (cmass_root < -cmass_root_inc) {
  620. cmass_leaf_inc = bminc + cmass_root;
  621. cmass_root_inc = -cmass_root;
  622. }
  623. }
  624. else if (cmass_leaf_inc < 0.0) {
  625. // Negative allocation to leaf
  626. // Make sure we don't end up with negative cmass_leaf
  627. if (cmass_leaf < -cmass_leaf_inc) {
  628. cmass_root_inc = bminc + cmass_leaf;
  629. cmass_leaf_inc = -cmass_leaf;
  630. }
  631. }
  632. }
  633. }
  634. // Check C budget after allocation
  635. // maximum carbon mismatch
  636. double EPS = 1.0e-12;
  637. assert(fabs(bminc + exceeds_cmass - (cmass_leaf_inc + cmass_root_inc + cmass_sap_inc + cmass_heart_inc + litter_leaf_inc + litter_root_inc)) < EPS);
  638. }
  639. void allocation_init(double bminit, double ltor, Individual& indiv) {
  640. // DESCRIPTION
  641. // Allocates initial biomass among tissues for a new individual (tree or grass),
  642. // assuming standard LPJ allometry (see functions allocation, allometry).
  643. // INPUT PARAMETERS
  644. // bminit = initial total biomass (kgC)
  645. // ltor = initial leaf:root biomass ratio
  646. //
  647. // Note: indiv.densindiv (density of individuals across patch or modelled area)
  648. // should be set to a meaningful value before this function is called
  649. double dval;
  650. double cmass_leaf_ind;
  651. double cmass_root_ind;
  652. double cmass_sap_ind;
  653. allocation(bminit, 0.0, 0.0, 0.0, 0.0, 0.0, ltor, 0.0, indiv.pft.sla, indiv.pft.wooddens,
  654. indiv.pft.lifeform, indiv.pft.k_latosa, indiv.pft.k_allom2, indiv.pft.k_allom3,
  655. cmass_leaf_ind, cmass_root_ind, cmass_sap_ind, dval, dval, dval, dval, dval);
  656. indiv.cmass_leaf = cmass_leaf_ind * indiv.densindiv;
  657. indiv.cmass_root = cmass_root_ind * indiv.densindiv;
  658. indiv.nmass_leaf = 0.0;
  659. indiv.nmass_root = 0.0;
  660. if (indiv.pft.lifeform == TREE) {
  661. indiv.cmass_sap = cmass_sap_ind * indiv.densindiv;
  662. indiv.nmass_sap = 0.0;
  663. }
  664. }
  665. ///////////////////////////////////////////////////////////////////////////////////////
  666. // ALLOMETRY
  667. // Should be called to update allometry, FPC and FPC increment whenever biomass values
  668. // for a vegetation individual change.
  669. bool allometry(Individual& indiv) {
  670. // DESCRIPTION
  671. // Calculates tree allometry (height and crown area) and fractional projective
  672. // given carbon biomass in various compartments for an individual.
  673. // Returns true if the allometry is normal, otherwise false - guess2008
  674. // TREE ALLOMETRY
  675. // Trees aboveground allometry is modelled by a cylindrical stem comprising an
  676. // inner cylinder of heartwood surrounded by a zone of sapwood of constant radius,
  677. // and a crown (i.e. foliage) cylinder of known diameter. Sapwood and heartwood are
  678. // assumed to have the same, constant, density (wooddens). Tree height is related
  679. // to sapwood cross-sectional area by the relation:
  680. // (1) height = cmass_sap / (sapwood xs area)
  681. // Sapwood cross-sectional area is also assumed to be a constant proportion of
  682. // total leaf area (following the "pipe model"; Shinozaki et al. 1964a,b; Waring
  683. // et al 1982), i.e.
  684. // (2) (leaf area) = k_latosa * (sapwood xs area)
  685. // Leaf area is related to leaf biomass by specific leaf area:
  686. // (3) (leaf area) = sla * cmass_leaf
  687. // From (1), (2), (3),
  688. // (4) height = cmass_sap / wooddens / sla / cmass_leaf * k_latosa
  689. // Tree height is related to stem diameter by the relation (Huang et al 1992)
  690. // [** = raised to the power of]:
  691. // (5) height = k_allom2 * diam ** k_allom3
  692. // Crown area may be derived from stem diameter by the relation (Zeide 1993):
  693. // (6) crownarea = min ( k_allom1 * diam ** k_rp , crownarea_max )
  694. // Bole height (individual/cohort mode only; currently set to 0):
  695. // (7) boleht = 0
  696. // FOLIAR PROJECTIVE COVER (FPC)
  697. // The same formulation for FPC (Eqn 8 below) is now applied in all vegetation
  698. // modes (Ben Smith 2002-07-23). FPC is equivalent to fractional patch/grid cell
  699. // coverage for the purposes of canopy exchange calculations and, in population
  700. // mode, vegetation dynamics calculations.
  701. //
  702. // FPC on the modelled area (stand, patch, "grid-cell") basis is related to mean
  703. // individual leaf area index (LAI) by the Lambert-Beer law (Monsi & Saeki 1953,
  704. // Prentice et al 1993) based on the assumption that success of a PFT population
  705. // in competition for space will be proportional to competitive ability for light
  706. // in the vertical profile of the forest canopy:
  707. // (8) fpc = crownarea * densindiv * ( 1.0 - exp ( -0.5 * lai_ind ) )
  708. // (8a)ppc = crownarea * densindiv
  709. // where
  710. // (9) lai_ind = cmass_leaf/densindiv * sla / crownarea
  711. //
  712. // For grasses,
  713. // (10) fpc = ( 1.0 - exp ( -0.5 * lai_ind ) )
  714. // (11) lai_ind = cmass_leaf * sla
  715. double diam; // stem diameter (m)
  716. double fpc_new; // updated FPC
  717. double ppc_new;
  718. // guess2008 - max tree height allowed (metre).
  719. const double HEIGHT_MAX = 150.0;
  720. if (indiv.pft.lifeform == TREE) {
  721. // TREES
  722. // Height (Eqn 4)
  723. // guess2008 - new allometry check
  724. if (!negligible(indiv.cmass_leaf)) {
  725. indiv.height = indiv.cmass_sap / indiv.cmass_leaf / indiv.pft.sla * indiv.pft.k_latosa / indiv.pft.wooddens;
  726. // Stem diameter (Eqn 5)
  727. diam = pow(indiv.height / indiv.pft.k_allom2, 1.0 / indiv.pft.k_allom3);
  728. // Stem volume
  729. double vol = indiv.height * PI * diam * diam * 0.25;
  730. if (indiv.age && (indiv.cmass_heart + indiv.cmass_sap) / indiv.densindiv / vol < indiv.pft.wooddens * 0.9) {
  731. return false;
  732. }
  733. }
  734. else {
  735. indiv.height = 0.0;
  736. diam = 0.0;
  737. return false;
  738. }
  739. // guess2008 - extra height check
  740. if (indiv.height > HEIGHT_MAX) {
  741. indiv.height = 0.0;
  742. diam = 0.0;
  743. return false;
  744. }
  745. // Crown area (Eqn 6)
  746. indiv.crownarea = min(indiv.pft.k_allom1 * pow(diam, indiv.pft.k_rp),
  747. indiv.pft.crownarea_max);
  748. if (!negligible(indiv.crownarea)) {
  749. // Individual LAI (Eqn 9)
  750. indiv.lai_indiv = indiv.cmass_leaf / indiv.densindiv *
  751. indiv.pft.sla / indiv.crownarea;
  752. // FPC (Eqn 8)
  753. fpc_new = indiv.crownarea * indiv.densindiv *
  754. (1.0 - lambertbeer(indiv.lai_indiv));
  755. // Increment deltafpc
  756. indiv.deltafpc += fpc_new - indiv.fpc;
  757. indiv.fpc = fpc_new;
  758. }
  759. else {
  760. indiv.lai_indiv = 0.0;
  761. indiv.fpc = 0.0;
  762. }
  763. // Bole height (Eqn 7)
  764. indiv.boleht = 0.0;
  765. // Stand-level LAI
  766. indiv.lai = indiv.cmass_leaf * indiv.pft.sla;
  767. }
  768. else if (indiv.pft.lifeform == GRASS) {
  769. // GRASSES
  770. if(indiv.pft.landcover != CROPLAND) {
  771. // guess2008 - bugfix - added if
  772. if (!negligible(indiv.cmass_leaf)) {
  773. // Grass "individual" LAI (Eqn 11)
  774. indiv.lai_indiv = indiv.cmass_leaf * indiv.pft.sla;
  775. // FPC (Eqn 10)
  776. indiv.fpc = 1.0 - lambertbeer(indiv.lai_indiv);
  777. // Stand-level LAI
  778. indiv.lai = indiv.lai_indiv;
  779. }
  780. else {
  781. return false;
  782. }
  783. }
  784. else {
  785. // True crops use cmass_leaf_max, cover-crop grass uses lai of stands with whole-year grass growth
  786. allometry_crop(indiv);
  787. }
  788. }
  789. // guess2008 - new return value (was void)
  790. return true;
  791. }
  792. ///////////////////////////////////////////////////////////////////////////////////////
  793. // RELATIVE CHANGE IN BIOMASS
  794. // Call this function to calculate the change in biomass on a grid cell area basis
  795. // associated with a specified change in FPC
  796. double fracmass_lpj(double fpc_low,double fpc_high,Individual& indiv) {
  797. // DESCRIPTION
  798. // Calculates and returns new biomass as a fraction of old biomass given an FPC
  799. // reduction from fpc_high to fpc_low, assuming LPJ allometry (see function
  800. // allometry)
  801. // guess2008 - check
  802. if (fpc_high < fpc_low)
  803. fail("fracmass_lpj: fpc_high < fpc_low");
  804. if (indiv.pft.lifeform==TREE) {
  805. if (negligible(fpc_high)) return 1.0;
  806. // else
  807. return fpc_low/fpc_high;
  808. }
  809. else if (indiv.pft.lifeform==GRASS) { // grass
  810. if (fpc_high>=1.0 || fpc_low>=1.0 || negligible(indiv.cmass_leaf)) return 1.0;
  811. // else
  812. return 1.0+2.0/indiv.cmass_leaf/indiv.pft.sla*
  813. (log(1.0-fpc_high)-log(1.0-fpc_low));
  814. }
  815. else {
  816. fail("fracmass_lpj: unknown lifeform");
  817. return 0;
  818. }
  819. // This point will never be reached in practice, but to satisfy more pedantic
  820. // compilers ...
  821. return 1.0;
  822. }
  823. void flush_litter_repr(Patch& patch) {
  824. // Returns nitrogen-free reproduction "litter" to atmosphere
  825. patch.pft.firstobj();
  826. while (patch.pft.isobj) {
  827. Patchpft& pft = patch.pft.getobj();
  828. // Updated soil fluxes
  829. patch.fluxes.report_flux(Fluxes::REPRC, pft.litter_repr);
  830. #ifdef CRESCENDO
  831. patch.fluxes.report_flux(Fluxes::DREPRC, pft.litter_repr);
  832. #endif
  833. pft.litter_repr = 0.0;
  834. patch.pft.nextobj();
  835. }
  836. }
  837. /// GROWTH
  838. /** Tissue turnover and allocation of fixed carbon to reproduction and new biomass
  839. * Accumulated NPP (assimilation minus maintenance and growth respiration) on
  840. * patch or modelled area basis assumed to be given by 'anpp' member variable for
  841. * each individual.
  842. * Should be called by framework at the end of each simulation year for modelling
  843. * of turnover, allocation and growth, prior to vegetation dynamics and disturbance
  844. */
  845. void growth(Stand& stand, Patch& patch) {
  846. // minimum carbon mass allowed (kgC/m2)
  847. const double MINCMASS = 1.0e-8;
  848. // maximum carbon mass allowed (kgC/m2)
  849. const double MAXCMASS = 1.0e8;
  850. const double CDEBT_PAYBACK_RATE = 0.2;
  851. // carbon biomass increment (component of NPP available for production of
  852. // new biomass) for this time period on modelled area basis (kgC/m2)
  853. double bminc = 0.0;
  854. // C allocated to reproduction this time period on modelled area basis (kgC/m2)
  855. double cmass_repr = 0.0;
  856. // increment in leaf C biomass following allocation, on individual basis (kgC)
  857. double cmass_leaf_inc = 0.0;
  858. // increment in root C biomass following allocation, on individual basis (kgC)
  859. double cmass_root_inc = 0.0;
  860. // increment in sapwood C biomass following allocation, on individual basis (kgC)
  861. double cmass_sap_inc = 0.0;
  862. // increment in heartwood C biomass following allocation, on individual basis (kgC)
  863. double cmass_heart_inc = 0.0;
  864. // increment in heartwood C biomass following allocation, on individual basis (kgC)
  865. double cmass_debt_inc = 0.0;
  866. // increment in crop harvestable organ C biomass following allocation, on individual basis (kgC)
  867. double cmass_ho_inc = 0.0;
  868. // increment in crop above-ground pool C biomass following allocation, on individual basis (kgC)
  869. double cmass_agpool_inc = 0.0;
  870. // increment in crop stem C biomass following allocation, on individual basis (kgC)
  871. double cmass_stem_inc = 0.0;
  872. // increment in leaf litter following allocation, on individual basis (kgC)
  873. double litter_leaf_inc = 0.0;
  874. // increment in root litter following allocation, on individual basis (kgC)
  875. double litter_root_inc = 0.0;
  876. // negative increment that exceeds existing biomass following allocation, on individual basis (kgC)
  877. double exceeds_cmass;
  878. // C biomass of leaves in "excess" of set allocated last year to raingreen PFT last year (kgC/m2)
  879. double cmass_excess = 0.0;
  880. // Raingreen nitrogen demand for leaves dropped during the year
  881. double raingreen_ndemand;
  882. // Nitrogen stress scalar for leaf to root allocation
  883. double nscal;
  884. // Leaf C:N ratios before growth
  885. double cton_leaf_bg;
  886. // Root C:N ratios before growth
  887. double cton_root_bg;
  888. // Sap C:N ratios before growth
  889. double cton_sap_bg;
  890. double dval = 0.0;
  891. int p;
  892. #ifdef CRESCENDO_FACE
  893. // Reset CRESCENDO FACE outputs
  894. patch.gl = 0.0;
  895. patch.gw = 0.0;
  896. patch.gr = 0.0;
  897. #endif
  898. // Obtain reference to Vegetation object for this patch
  899. Vegetation& vegetation = patch.vegetation;
  900. Gridcell& gridcell = vegetation.patch.stand.get_gridcell();
  901. // On first call to function growth this year (patch #0), initialise stand-PFT
  902. // record of summed allocation to reproduction
  903. if (!patch.id)
  904. for (p=0; p<npft; p++)
  905. stand.pft[p].cmass_repr = 0.0;
  906. // Set forest management intensity for this year
  907. patch.man_strength = cut_fraction(patch);
  908. // Loop through individuals
  909. vegetation.firstobj();
  910. while (vegetation.isobj) {
  911. Individual& indiv = vegetation.getobj();
  912. // For this individual
  913. // Calculate vegetation carbon and nitrogen mass before growth to determine vegetation C:N ratios
  914. //indiv.cmass_veg = indiv.cmass_leaf + indiv.cmass_root + indiv.cmass_wood();
  915. //indiv.nmass_veg = indiv.nmass_leaf + indiv.nmass_root + indiv.nmass_wood();
  916. // Save real compartment C:N ratios before growth
  917. // phen is switched off for leaf and root
  918. // Leaf
  919. cton_leaf_bg = indiv.cton_leaf(false);
  920. // Root
  921. cton_root_bg = indiv.cton_root(false);
  922. // Sap
  923. cton_sap_bg = indiv.cton_sap();
  924. // Save leaf annual average C:N ratio
  925. if (!negligible(indiv.nday_leafon))
  926. indiv.cton_leaf_aavr /= indiv.nday_leafon;
  927. else
  928. indiv.cton_leaf_aavr = indiv.pft.cton_leaf_max;
  929. // Nitrogen stress scalar for leaf to root allocation (adopted from Zaehle and Friend 2010 SM eq 19)
  930. double cton_leaf_aopt = max(indiv.cton_leaf_aopt, indiv.pft.cton_leaf_avr);
  931. if (ifnlim)
  932. nscal = min(1.0, cton_leaf_aopt / indiv.cton_leaf_aavr);
  933. else
  934. nscal = 1.0;
  935. // Set leaf:root mass ratio based on water stress parameter
  936. // or nitrogen stress scalar
  937. indiv.ltor = min(indiv.wscal_mean(), nscal) * indiv.pft.ltor_max;
  938. // Move leftover compartment nitrogen storage to longterm storage
  939. if(!indiv.has_daily_turnover()) {
  940. indiv.nstore_longterm += indiv.nstore_labile;
  941. indiv.nstore_labile = 0.0;
  942. }
  943. indiv.deltafpc = 0.0;
  944. bool killed = false;
  945. if (negligible(indiv.densindiv))
  946. fail("growth: negligible densindiv for %s",(char*)indiv.pft.name);
  947. else {
  948. // Allocation to reproduction
  949. if(!indiv.istruecrop_or_intercropgrass())
  950. reproduction(indiv.pft.reprfrac,indiv.anpp,bminc,cmass_repr);
  951. raingreen_ndemand = 0.0;
  952. // added bminc check. Otherwise we get -ve litter_leaf for grasses when indiv.anpp < 0.
  953. //
  954. // cmass_excess-code inactivated for grass pft:s, due to frequent oscillations between high bminc and zero bminc in
  955. // certain grasslands using updated leaflong-value for grass-pft:s (0.5)
  956. if (bminc >= 0 && indiv.pft.phenology == RAINGREEN) {
  957. // Raingreen PFTs: reduce biomass increment to account for NPP
  958. // allocated to extra leaves during the past year.
  959. // Excess allocation to leaves given by:
  960. // aphen_raingreen / ( leaf_longevity * year_length) * cmass_leaf -
  961. // cmass_leaf
  962. // BLARP! excess allocation to roots now also included (assumes leaf longevity = root longevity)
  963. cmass_excess = max((double)indiv.aphen_raingreen /
  964. (indiv.pft.leaflong * (double)date.year_length()) * (indiv.cmass_leaf + indiv.cmass_root) -
  965. indiv.cmass_leaf - indiv.cmass_root, 0.0);
  966. if (cmass_excess > bminc) cmass_excess = bminc;
  967. // Transfer excess leaves to litter
  968. // only for 'alive' individuals
  969. if (indiv.alive) {
  970. patch.pft[indiv.pft.id].litter_leaf += cmass_excess;
  971. if (!negligible(cton_leaf_bg))
  972. raingreen_ndemand = min(indiv.nmass_leaf, cmass_excess / cton_leaf_bg);
  973. else
  974. raingreen_ndemand = 0.0;
  975. patch.pft[indiv.pft.id].nmass_litter_leaf += raingreen_ndemand * (1.0 - nrelocfrac);
  976. indiv.nstore_longterm += raingreen_ndemand * nrelocfrac;
  977. indiv.nmass_leaf -= raingreen_ndemand;
  978. }
  979. // Deduct from this year's C biomass increment
  980. // added alive check
  981. if (indiv.alive) bminc -= cmass_excess;
  982. }
  983. // All yearly harvest events
  984. killed = harvest_year(indiv);
  985. if (!killed) {
  986. if(!indiv.has_daily_turnover()) {
  987. // Tissue turnover and associated litter production
  988. turnover(indiv.pft.turnover_leaf, indiv.pft.turnover_root,
  989. indiv.pft.turnover_sap, indiv.pft.lifeform, indiv.pft.landcover,
  990. indiv.cmass_leaf, indiv.cmass_root, indiv.cmass_sap, indiv.cmass_heart,
  991. indiv.nmass_leaf, indiv.nmass_root, indiv.nmass_sap, indiv.nmass_heart,
  992. patch.pft[indiv.pft.id].litter_leaf,
  993. patch.pft[indiv.pft.id].litter_root,
  994. patch.pft[indiv.pft.id].nmass_litter_leaf,
  995. patch.pft[indiv.pft.id].nmass_litter_root,
  996. indiv.nstore_longterm,indiv.max_n_storage,
  997. indiv.alive);
  998. }
  999. // Update stand record of reproduction by this PFT
  1000. stand.pft[indiv.pft.id].cmass_repr += cmass_repr / (double)stand.npatch();
  1001. // Transfer reproduction straight to litter
  1002. // only for 'alive' individuals
  1003. if (indiv.alive) {
  1004. patch.pft[indiv.pft.id].litter_repr += cmass_repr;
  1005. }
  1006. if (indiv.pft.lifeform == TREE) {
  1007. // TREE GROWTH
  1008. // BLARP! Try and pay back part of cdebt
  1009. if (ifcdebt && bminc > 0.0) {
  1010. double cmass_payback = min(indiv.cmass_debt * CDEBT_PAYBACK_RATE, bminc);
  1011. bminc -= cmass_payback;
  1012. indiv.cmass_debt -= cmass_payback;
  1013. }
  1014. // Allocation: note conversion of mass values from grid cell area
  1015. // to individual basis
  1016. allocation(bminc / indiv.densindiv, indiv.cmass_leaf / indiv.densindiv,
  1017. indiv.cmass_root / indiv.densindiv, indiv.cmass_sap / indiv.densindiv,
  1018. indiv.cmass_debt / indiv.densindiv,
  1019. indiv.cmass_heart / indiv.densindiv, indiv.ltor,
  1020. indiv.height, indiv.pft.sla, indiv.pft.wooddens, TREE,
  1021. indiv.pft.k_latosa, indiv.pft.k_allom2, indiv.pft.k_allom3,
  1022. cmass_leaf_inc, cmass_root_inc, cmass_sap_inc, cmass_debt_inc,
  1023. cmass_heart_inc,
  1024. litter_leaf_inc, litter_root_inc, exceeds_cmass);
  1025. // Update carbon pools and litter (on area basis)
  1026. // (litter not accrued for not 'alive' individuals - Ben 2007-11-28)
  1027. // Leaves
  1028. indiv.cmass_leaf += cmass_leaf_inc * indiv.densindiv;
  1029. // Roots
  1030. indiv.cmass_root += cmass_root_inc * indiv.densindiv;
  1031. // Sapwood
  1032. indiv.cmass_sap += cmass_sap_inc * indiv.densindiv;
  1033. // Heartwood
  1034. indiv.cmass_heart += cmass_heart_inc * indiv.densindiv;
  1035. #ifdef CRESCENDO_FACE
  1036. patch.gl += cmass_leaf_inc * indiv.densindiv;
  1037. patch.gr += cmass_root_inc * indiv.densindiv;
  1038. patch.gw += (cmass_sap_inc + cmass_heart_inc) * indiv.densindiv;
  1039. #endif
  1040. indiv.cmass_wood_inc_5.add((cmass_sap_inc + cmass_heart_inc - cmass_debt_inc) * indiv.densindiv);
  1041. // If negative sap growth, then nrelocfrac of nitrogen will go to heart wood and
  1042. // (1.0 - nreloctrac) will go to storage
  1043. double nmass_sap_inc = cmass_sap_inc * indiv.densindiv / cton_sap_bg;
  1044. if (cmass_sap_inc < 0.0 && indiv.nmass_sap >= -nmass_sap_inc){
  1045. indiv.nmass_sap += nmass_sap_inc;
  1046. indiv.nmass_heart -= nmass_sap_inc * nrelocfrac;
  1047. indiv.nstore_longterm -= nmass_sap_inc * (1.0 - nrelocfrac);
  1048. // ecev3 - check
  1049. if (indiv.nmass_sap < 0) {
  1050. dprintf("ERROR! Year %d day %d Stand %d: Negative indiv.nmass_sap in growth: %.15f\n", date.year, date.day, patch.stand.id, indiv.nmass_sap);
  1051. }
  1052. assert(indiv.nmass_sap >= 0.0);
  1053. }
  1054. // C debt
  1055. indiv.cmass_debt += cmass_debt_inc * indiv.densindiv;
  1056. // alive check before ensuring C balance
  1057. if (indiv.alive) {
  1058. patch.pft[indiv.pft.id].litter_leaf += litter_leaf_inc * indiv.densindiv;
  1059. patch.pft[indiv.pft.id].litter_root += litter_root_inc * indiv.densindiv;
  1060. // C litter exceeding existing biomass
  1061. indiv.report_flux(Fluxes::DEBEXCC, exceeds_cmass * indiv.densindiv);
  1062. indiv.report_flux(Fluxes::RA, -exceeds_cmass * indiv.densindiv);
  1063. }
  1064. double leaf_inc = min(indiv.nmass_leaf, litter_leaf_inc * indiv.densindiv / cton_leaf_bg);
  1065. double root_inc = min(indiv.nmass_root, litter_root_inc * indiv.densindiv / cton_root_bg);
  1066. // Nitrogen litter always return to soil litter and storage
  1067. // Leaf
  1068. patch.pft[indiv.pft.id].nmass_litter_leaf += leaf_inc * (1.0 - nrelocfrac);
  1069. indiv.nstore_longterm += leaf_inc * nrelocfrac;
  1070. // Root
  1071. patch.pft[indiv.pft.id].nmass_litter_root += root_inc * (1.0 - nrelocfrac);
  1072. indiv.nstore_longterm += root_inc * nrelocfrac;
  1073. // Subtracting litter nitrogen from individuals
  1074. indiv.nmass_leaf -= leaf_inc;
  1075. indiv.nmass_root -= root_inc;
  1076. // Update individual age
  1077. indiv.age++;
  1078. // Kill individual and transfer biomass to litter if any biomass
  1079. // compartment negative
  1080. if (indiv.cmass_leaf < MINCMASS || indiv.cmass_root < MINCMASS ||
  1081. indiv.cmass_sap < MINCMASS) {
  1082. indiv.kill();
  1083. vegetation.killobj();
  1084. killed = true;
  1085. }
  1086. }
  1087. else if (indiv.pft.lifeform == GRASS) {
  1088. // GRASS GROWTH
  1089. //True crops do not use bminc.or cmass_leaf etc.
  1090. if(indiv.istruecrop_or_intercropgrass()) {
  1091. // transfer crop cmass increase values to common variables
  1092. growth_crop_year(indiv, cmass_leaf_inc, cmass_root_inc, cmass_ho_inc, cmass_agpool_inc, cmass_stem_inc);
  1093. exceeds_cmass = 0.0; //exceeds_cmass not used for true crops
  1094. }
  1095. else {
  1096. allocation(bminc, indiv.cmass_leaf, indiv.cmass_root,
  1097. 0.0, 0.0, 0.0, indiv.ltor, 0.0, 0.0, 0.0, GRASS, 0.0,
  1098. 0.0, 0.0, cmass_leaf_inc, cmass_root_inc, dval, dval, dval,
  1099. litter_leaf_inc, litter_root_inc, exceeds_cmass);
  1100. }
  1101. if(indiv.pft.landcover == CROPLAND) {
  1102. if(indiv.istruecrop_or_intercropgrass())
  1103. yield_crop(indiv);
  1104. else
  1105. yield_pasture(indiv, cmass_leaf_inc);
  1106. }
  1107. // Update carbon pools and litter (on area basis)
  1108. // only litter in the case of 'alive' individuals
  1109. // Leaves
  1110. indiv.cmass_leaf += cmass_leaf_inc;
  1111. // Roots
  1112. indiv.cmass_root += cmass_root_inc;
  1113. #ifdef CRESCENDO_FACE
  1114. patch.gl += cmass_leaf_inc;
  1115. patch.gr += cmass_root_inc;
  1116. #endif
  1117. if(indiv.pft.landcover == CROPLAND) {
  1118. indiv.cropindiv->cmass_ho += cmass_ho_inc;
  1119. indiv.cropindiv->cmass_agpool += cmass_agpool_inc;
  1120. indiv.cropindiv->cmass_stem += cmass_stem_inc;
  1121. }
  1122. if(indiv.pft.phenology != CROPGREEN && !(indiv.has_daily_turnover() && indiv.continous_grass())) {
  1123. // alive check before ensuring C balance
  1124. if (indiv.alive && !indiv.istruecrop_or_intercropgrass()) {
  1125. patch.pft[indiv.pft.id].litter_leaf += litter_leaf_inc;
  1126. patch.pft[indiv.pft.id].litter_root += litter_root_inc;
  1127. // C litter exceeding existing biomass
  1128. indiv.report_flux(Fluxes::DEBEXCC, exceeds_cmass * indiv.densindiv);
  1129. indiv.report_flux(Fluxes::RA, -exceeds_cmass * indiv.densindiv);
  1130. }
  1131. // Nitrogen always return to soil litter and storage
  1132. // Leaf
  1133. if (indiv.nmass_leaf > 0.0){
  1134. patch.pft[indiv.pft.id].nmass_litter_leaf += litter_leaf_inc * indiv.densindiv /
  1135. cton_leaf_bg * (1.0 - nrelocfrac);
  1136. indiv.nstore_longterm += litter_leaf_inc * indiv.densindiv / cton_leaf_bg * nrelocfrac;
  1137. // Subtracting litter nitrogen from individuals
  1138. indiv.nmass_leaf -= min(indiv.nmass_leaf, litter_leaf_inc * indiv.densindiv / cton_leaf_bg);
  1139. }
  1140. // Root
  1141. if (indiv.nmass_root > 0.0){
  1142. patch.pft[indiv.pft.id].nmass_litter_root += litter_root_inc * indiv.densindiv /
  1143. cton_root_bg * (1.0 - nrelocfrac);
  1144. indiv.nstore_longterm += litter_root_inc / cton_root_bg * nrelocfrac;
  1145. // Subtracting litter nitrogen from individuals
  1146. indiv.nmass_root -= min(indiv.nmass_root, litter_root_inc * indiv.densindiv / cton_root_bg);
  1147. }
  1148. }
  1149. // Kill individual and transfer biomass to litter if either biomass
  1150. // compartment negative
  1151. if ((indiv.cmass_leaf < MINCMASS || indiv.cmass_root < MINCMASS) && !indiv.istruecrop_or_intercropgrass()) {
  1152. indiv.kill();
  1153. vegetation.killobj();
  1154. killed = true;
  1155. }
  1156. }
  1157. if (!killed && indiv.pft.phenology != CROPGREEN && !(indiv.has_daily_turnover() && indiv.continous_grass())) {
  1158. // Update nitrogen longtime storage
  1159. // Nitrogen approx retranslocated next year
  1160. double retransn_nextyear = indiv.cmass_leaf * indiv.pft.turnover_leaf / cton_leaf_bg * nrelocfrac +
  1161. indiv.cmass_root * indiv.pft.turnover_root / cton_root_bg * nrelocfrac;
  1162. if (indiv.pft.lifeform == TREE)
  1163. retransn_nextyear += indiv.cmass_sap * indiv.pft.turnover_sap / cton_sap_bg * nrelocfrac;
  1164. // Assume that raingreen will lose same amount of N through extra leaves next year
  1165. if (indiv.alive && bminc >= 0 && (indiv.pft.phenology == RAINGREEN || indiv.pft.phenology == ANY))
  1166. retransn_nextyear -= min(raingreen_ndemand, retransn_nextyear);
  1167. // Max longterm nitrogen storage
  1168. if (indiv.pft.lifeform == TREE)
  1169. indiv.max_n_storage = max(0.0, min(indiv.cmass_sap * indiv.pft.fnstorage / cton_leaf_bg, retransn_nextyear));
  1170. else // GRASS
  1171. indiv.max_n_storage = max(0.0, min(indiv.cmass_root * indiv.pft.fnstorage / cton_leaf_bg, retransn_nextyear));
  1172. // Scale this year productivity to max storage
  1173. if (indiv.anpp > 0.0) {
  1174. indiv.scale_n_storage = max(indiv.max_n_storage * 0.1, indiv.max_n_storage - retransn_nextyear) * cton_leaf_bg / indiv.anpp;
  1175. } // else use last years scaling factor
  1176. }
  1177. }
  1178. }
  1179. if (!killed) {
  1180. if (!allometry(indiv)) {
  1181. indiv.kill();
  1182. vegetation.killobj();
  1183. killed = true;
  1184. }
  1185. if (!killed) {
  1186. if (!indiv.alive) {
  1187. // The individual has survived its first year...
  1188. indiv.alive = true;
  1189. // ...now we can start counting its fluxes,
  1190. // debit current biomass as establishment flux
  1191. if (!indiv.istruecrop_or_intercropgrass()) {
  1192. indiv.report_flux(Fluxes::ESTC,
  1193. - (indiv.cmass_leaf + indiv.cmass_root + indiv.cmass_sap +
  1194. indiv.cmass_heart - indiv.cmass_debt));
  1195. #ifdef CRESCENDO_FACE
  1196. indiv.report_flux(Fluxes::DESTC,
  1197. -(indiv.cmass_leaf + indiv.cmass_root + indiv.cmass_sap +
  1198. indiv.cmass_heart - indiv.cmass_debt));
  1199. #endif
  1200. }
  1201. }
  1202. // Move long-term nitrogen storage pool to labile storage pool for usage next year
  1203. if(!indiv.has_daily_turnover()) {
  1204. indiv.nstore_labile = indiv.nstore_longterm;
  1205. indiv.nstore_longterm = 0.0;
  1206. }
  1207. // ... on to next individual
  1208. vegetation.nextobj();
  1209. }
  1210. }
  1211. }
  1212. // Flush nitrogen free litter from reproduction straight to atmosphere
  1213. flush_litter_repr(patch);
  1214. }
  1215. ///////////////////////////////////////////////////////////////////////////////////////
  1216. // REFERENCES
  1217. //
  1218. // LPJF refers to the original FORTRAN implementation of LPJ as described by Sitch
  1219. // et al 2001
  1220. // Huang, S, Titus, SJ & Wiens, DP (1992) Comparison of nonlinear height-diameter
  1221. // functions for major Alberta tree species. Canadian Journal of Forest Research 22:
  1222. // 1297-1304
  1223. // Monsi M & Saeki T 1953 Ueber den Lichtfaktor in den Pflanzengesellschaften und
  1224. // seine Bedeutung fuer die Stoffproduktion. Japanese Journal of Botany 14: 22-52
  1225. // Prentice, IC, Sykes, MT & Cramer W (1993) A simulation model for the transient
  1226. // effects of climate change on forest landscapes. Ecological Modelling 65: 51-70.
  1227. // Press, WH, Teukolsky, SA, Vetterling, WT & Flannery, BT. (1986) Numerical
  1228. // Recipes in FORTRAN, 2nd ed. Cambridge University Press, Cambridge
  1229. // Sitch, S, Prentice IC, Smith, B & Other LPJ Consortium Members (2000) LPJ - a
  1230. // coupled model of vegetation dynamics and the terrestrial carbon cycle. In:
  1231. // Sitch, S. The Role of Vegetation Dynamics in the Control of Atmospheric CO2
  1232. // Content, PhD Thesis, Lund University, Lund, Sweden.
  1233. // Shinozaki, K, Yoda, K, Hozumi, K & Kira, T (1964) A quantitative analysis of
  1234. // plant form - the pipe model theory. I. basic analyses. Japanese Journal of
  1235. // Ecology 14: 97-105
  1236. // Shinozaki, K, Yoda, K, Hozumi, K & Kira, T (1964) A quantitative analysis of
  1237. // plant form - the pipe model theory. II. further evidence of the theory and
  1238. // its application in forest ecology. Japanese Journal of Ecology 14: 133-139
  1239. // Sykes, MT, Prentice IC & Cramer W 1996 A bioclimatic model for the potential
  1240. // distributions of north European tree species under present and future climates.
  1241. // Journal of Biogeography 23: 209-233.
  1242. // Waring, RH Schroeder, PE & Oren, R (1982) Application of the pipe model theory
  1243. // to predict canopy leaf area. Canadian Journal of Forest Research 12:
  1244. // 556-560
  1245. // Zaehle, S. & Friend, A. D. 2010. Carbon and nitrogen cycle dynamics in the O-CN
  1246. // land surface model: 1. Model description, site-scale evaluation, and sensitivity
  1247. // to parameter estimates. Global Biogeochemical Cycles, 24.
  1248. // Zeide, B (1993) Primary unit of the tree crown. Ecology 74: 1598-1602.