cropallocation.cpp 34 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963
  1. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  2. /// \file cropallocation.cpp
  3. /// \brief Crop allocation and growth
  4. /// \author Mats Lindeskog, Stefan Olin
  5. /// $Date: $
  6. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  7. #include "landcover.h"
  8. #include "cropallocation.h"
  9. // Seed carbon allocation to leaves and roots are done over a 10-day period
  10. const bool DELAYED_SEEDCARBON = false;
  11. /// Updates patch members fpc_total and fpc_rescale for crops (to be called after crop_phenology())
  12. void update_patch_fpc(Patch& patch) {
  13. if(patch.stand.landcover != CROPLAND) {
  14. return;
  15. }
  16. Vegetation& vegetation = patch.vegetation;
  17. patch.fpc_total = 0.0;
  18. vegetation.firstobj();
  19. while (vegetation.isobj) {
  20. Individual& indiv = vegetation.getobj();
  21. if(indiv.growingseason()) {
  22. patch.fpc_total += indiv.fpc;
  23. }
  24. vegetation.nextobj();
  25. }
  26. // Calculate rescaling factor to account for overlap between populations/
  27. // cohorts/individuals (i.e. total FPC > 1)
  28. // necessary to undate here after growingseason updated
  29. patch.fpc_rescale = 1.0 / max(patch.fpc_total, 1.0);
  30. }
  31. /// Updates lai_daily and fpc_daily from daily grs_cmass_leaf-value
  32. /** lai during senescence declines according the function senescence_curve()
  33. */
  34. void lai_crop(Patch& patch) {
  35. Vegetation& vegetation = patch.vegetation;
  36. vegetation.firstobj();
  37. while (vegetation.isobj) {
  38. Individual& indiv = vegetation.getobj();
  39. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  40. Patchpft& patchpft = patch.pft[indiv.pft.id];
  41. cropphen_struct& ppftcrop = *(patchpft.get_cropphen());
  42. if(indiv.pft.phenology == CROPGREEN) {
  43. if(ppftcrop.growingseason) {
  44. if(!ppftcrop.senescence || ifnlim)
  45. indiv.lai_daily = cropindiv.grs_cmass_leaf * indiv.pft.sla;
  46. else
  47. // Follow the senescence curve from leaf cmass at senescence (cmass_leaf_sen):
  48. indiv.lai_daily = cropindiv.cmass_leaf_sen * indiv.pft.sla * senescence_curve(indiv.pft, ppftcrop.fphu);
  49. if(indiv.lai_daily < 0.0)
  50. indiv.lai_daily = 0.0;
  51. indiv.fpc_daily = 1.0 - lambertbeer(indiv.lai_daily);
  52. indiv.lai_indiv_daily = indiv.lai_daily;
  53. }
  54. else if(date.day == ppftcrop.hdate) {
  55. indiv.lai_daily = 0.0;
  56. indiv.lai_indiv_daily = 0.0;
  57. indiv.fpc_daily = 0.0;
  58. }
  59. }
  60. vegetation.nextobj();
  61. }
  62. }
  63. /// Turnover function for continuous grass, to be called from any day of the year from allocation_crop_daily().
  64. void turnover_grass(Individual& indiv) {
  65. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  66. Patchpft& patchpft = indiv.patchpft();
  67. double cmass_leaf_inc = cropindiv.grs_cmass_leaf - indiv.cmass_leaf_post_turnover;
  68. double cmass_root_inc = cropindiv.grs_cmass_root - indiv.cmass_root_post_turnover;
  69. double grs_npp = cmass_leaf_inc + cmass_root_inc - CMASS_SEED;
  70. double cmass_leaf_pre_turnover = cropindiv.grs_cmass_leaf;
  71. double cmass_root_pre_turnover = cropindiv.grs_cmass_root;
  72. double cton_leaf_bg = indiv.cton_leaf(false);
  73. double cton_root_bg = indiv.cton_root(false);
  74. indiv.nstore_longterm += indiv.nstore_labile;
  75. indiv.nstore_labile = 0.0;
  76. turnover(indiv.pft.turnover_leaf, indiv.pft.turnover_root,
  77. indiv.pft.turnover_sap, indiv.pft.lifeform, indiv.pft.landcover,
  78. indiv.cropindiv->grs_cmass_leaf, indiv.cropindiv->grs_cmass_root, indiv.cmass_sap, indiv.cmass_heart,
  79. indiv.nmass_leaf, indiv.nmass_root, indiv.nmass_sap, indiv.nmass_heart,
  80. patchpft.litter_leaf,
  81. patchpft.litter_root,
  82. patchpft.nmass_litter_leaf,
  83. patchpft.nmass_litter_root,
  84. indiv.nstore_longterm, indiv.max_n_storage,
  85. true);
  86. indiv.cmass_leaf_post_turnover = cropindiv.grs_cmass_leaf;
  87. indiv.cmass_root_post_turnover = cropindiv.grs_cmass_root;
  88. // Nitrogen longtime storage
  89. // Nitrogen approx retranslocated next season
  90. double retransn_nextyear = cmass_leaf_pre_turnover * indiv.pft.turnover_leaf / cton_leaf_bg * nrelocfrac +
  91. cmass_root_pre_turnover * indiv.pft.turnover_root / cton_root_bg * nrelocfrac;
  92. // Max longterm nitrogen storage
  93. indiv.max_n_storage = max(0.0, min(cmass_root_pre_turnover * indiv.pft.fnstorage / cton_leaf_bg, retransn_nextyear));
  94. // Scale this year productivity to max storage
  95. if (grs_npp > 0.0) {
  96. indiv.scale_n_storage = max(indiv.max_n_storage * 0.1, indiv.max_n_storage - retransn_nextyear) * cton_leaf_bg / grs_npp;
  97. }
  98. indiv.nstore_labile = indiv.nstore_longterm;
  99. indiv.nstore_longterm = 0.0;
  100. }
  101. /// Help function used by allocation_crop_nlim(), described in Olin et al. 2015.
  102. void crop_allocation_devries(cropphen_struct& ppftcrop, Individual& indiv) {
  103. // Currently, the development stage (ds) calculations are done using a linear relationship
  104. // between ds and fphu to allow for dynamic variety selection (Lindeskog 2013).
  105. double t = 0.0;
  106. if(ppftcrop.fphu < 0.4367) {
  107. t = -0.07 + 2.45 * ppftcrop.fphu;
  108. }
  109. else {
  110. t = 0.2247 + 1.7753 * ppftcrop.fphu;
  111. }
  112. // Comment out this line if ds should be calculated according to Olin et al. 2015.
  113. ppftcrop.dev_stage = max(0.0,min(2.0,t));
  114. // Eq. 3-5, Olin 2015
  115. double f1 = min(1.0, max(0.0, richards_curve(indiv.pft.a1, indiv.pft.b1, indiv.pft.c1, indiv.pft.d1, ppftcrop.dev_stage)));
  116. double f2 = min(1.0, max(0.0, richards_curve(indiv.pft.a2, indiv.pft.b2, indiv.pft.c2, indiv.pft.d2, ppftcrop.dev_stage)));
  117. double f3 = min(1.0, max(0.0, richards_curve(indiv.pft.a3, indiv.pft.b3, indiv.pft.c3, indiv.pft.d3, ppftcrop.dev_stage)));
  118. // Eq. 15, Olin 2015
  119. if(indiv.daily_cmass_leafloss > 0.0)
  120. f2 *= f2;
  121. // Eq. 6, Olin 2015
  122. ppftcrop.f_alloc_root = f1 * (1-f3);
  123. ppftcrop.f_alloc_leaf = f2 * (1-f1)*(1-f3);
  124. ppftcrop.f_alloc_stem = (1.0 - f2)*(1.0 - f1)*(1.0 - f3);
  125. ppftcrop.f_alloc_horg = f3;
  126. }
  127. /// Daily allocation routine for crops with nitrogen limitation
  128. /** Allocates daily npp to leaf, roots and harvestable organs
  129. * according to Olin et al. 2015.
  130. */
  131. void allocation_crop_nlim(Individual& indiv, double cmass_seed, double nmass_seed) {
  132. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  133. Patch& patch = indiv.vegetation.patch;
  134. Patchpft& patchpft = patch.pft[indiv.pft.id];
  135. cropphen_struct& ppftcrop = *(patchpft.get_cropphen());
  136. if (!ppftcrop.growingseason) {
  137. return;
  138. }
  139. // report seed fluxes
  140. indiv.report_flux(Fluxes::SEEDC, -cmass_seed);
  141. indiv.report_flux(Fluxes::SEEDN, -nmass_seed);
  142. // add seed carbon
  143. double cmass_extra = cmass_seed;
  144. // add seed nitrogen, with the assumption that the seedling is initiated with 50% in roots and 50% in the leaves.
  145. indiv.nmass_leaf += nmass_seed / 2.0;
  146. indiv.nmass_root += nmass_seed / 2.0;
  147. // Get the daily allocation strategy.
  148. crop_allocation_devries(ppftcrop, indiv);
  149. // Use the fast C pool when NPP is negative.
  150. if(indiv.dnpp < 0.0){
  151. if(-indiv.dnpp < cropindiv.grs_cmass_agpool) {
  152. cropindiv.grs_cmass_agpool -= -indiv.dnpp;
  153. cropindiv.ycmass_agpool -= -indiv.dnpp;
  154. indiv.dnpp = 0.0;
  155. }
  156. else {
  157. indiv.report_flux(Fluxes::NPP, (-indiv.dnpp - cropindiv.grs_cmass_agpool));
  158. indiv.report_flux(Fluxes::GPP, (-indiv.dnpp - cropindiv.grs_cmass_agpool));
  159. cropindiv.ycmass_agpool -= cropindiv.grs_cmass_agpool;
  160. cropindiv.grs_cmass_agpool = 0.0;
  161. indiv.dnpp = 0.0;
  162. }
  163. }
  164. // Retranslocation from the fast C pool to the grains towards the end of the grainfilling period.
  165. // NB, only works for cereals (grasses), needs to be adjusted to work for herb and tuber crops.
  166. if (cropindiv.grs_cmass_agpool > 0.0 && patchpft.cropphen->f_alloc_horg > 0.95) {
  167. cmass_extra += 0.1 * cropindiv.grs_cmass_agpool;
  168. cropindiv.ycmass_agpool -= 0.1 * cropindiv.grs_cmass_agpool;
  169. cropindiv.grs_cmass_agpool *= 0.9;
  170. }
  171. indiv.daily_cmass_rootloss = 0.0;
  172. indiv.daily_nmass_rootloss = 0.0;
  173. // If senescense have occured this day.
  174. if (indiv.daily_cmass_leafloss > 0.0) {
  175. // Daily C mass leaf increment.
  176. cropindiv.dcmass_leaf = (indiv.dnpp + cmass_extra) * patchpft.cropphen->f_alloc_leaf - indiv.daily_cmass_leafloss;
  177. cropindiv.grs_cmass_dead_leaf += indiv.daily_cmass_leafloss;
  178. cropindiv.ycmass_dead_leaf += indiv.daily_cmass_leafloss;
  179. if (indiv.daily_cmass_leafloss / 100.0<indiv.nmass_leaf) {
  180. cropindiv.nmass_dead_leaf += indiv.daily_cmass_leafloss / 100.0; // Fixed low C:N in the dead leaf
  181. cropindiv.ynmass_dead_leaf += indiv.daily_cmass_leafloss / 100.0;
  182. indiv.nmass_leaf -= indiv.daily_cmass_leafloss / 100.0;
  183. }
  184. double new_CN = (cropindiv.grs_cmass_leaf + cropindiv.dcmass_leaf) / indiv.nmass_leaf;
  185. // If the result is smaller (higher [N]) than the min C:N then that N is
  186. // put in to the ag N pool
  187. if( new_CN < indiv.pft.cton_leaf_min ) {
  188. indiv.daily_nmass_leafloss = max(0.0, indiv.nmass_leaf - (cropindiv.grs_cmass_leaf + cropindiv.dcmass_leaf) / (1.33 * indiv.pft.cton_leaf_min));
  189. if(indiv.daily_nmass_leafloss > indiv.nmass_leaf) {
  190. indiv.daily_nmass_leafloss = 0.0;
  191. }
  192. } else {
  193. indiv.daily_nmass_leafloss = 0.0;
  194. }
  195. // Very experimental, root senescence
  196. // N and C loss when root senescence is allowed (f_HO > 0.5)
  197. // d3, the DS after which more than half of the daily assimilates are going to the grains.
  198. if(patchpft.cropphen->dev_stage > indiv.pft.d3) {
  199. // only have root senescence when leaf scenescence has occured
  200. if (indiv.daily_nmass_leafloss > 0.0) {
  201. double kC = 0.0;
  202. double kN = 0.0;
  203. // The root senescence is proportional to that of the leaves, Eq. 10 Olin 2015
  204. if(indiv.nmass_leaf > 0.0) {
  205. kN = indiv.daily_nmass_leafloss / indiv.nmass_leaf;
  206. }
  207. if (indiv.cmass_leaf_today() > 0.0) {
  208. kC = indiv.daily_cmass_leafloss / indiv.cmass_leaf_today();
  209. }
  210. indiv.daily_cmass_rootloss = indiv.cmass_root_today() * kC;
  211. indiv.daily_nmass_rootloss = indiv.nmass_root * kN;
  212. }
  213. }
  214. indiv.nmass_leaf -= indiv.daily_nmass_leafloss;
  215. cropindiv.nmass_agpool += indiv.daily_nmass_leafloss;
  216. }
  217. else {
  218. cropindiv.dcmass_leaf = (indiv.dnpp + cmass_extra) * patchpft.cropphen->f_alloc_leaf;
  219. }
  220. cropindiv.dcmass_stem = (indiv.dnpp + cmass_extra) * patchpft.cropphen->f_alloc_stem;
  221. if (indiv.daily_cmass_rootloss > indiv.cmass_root_today())
  222. indiv.daily_cmass_rootloss = 0.0;
  223. cropindiv.dcmass_root = (indiv.dnpp + cmass_extra) * patchpft.cropphen->f_alloc_root - indiv.daily_cmass_rootloss;
  224. // The lost root C is directly put into the litter.
  225. patch.soil.sompool[SOILMETA].cmass += indiv.daily_cmass_rootloss;
  226. patch.fluxes.report_flux(Fluxes::VEGLITTERC, indiv.daily_cmass_rootloss);
  227. if (indiv.daily_nmass_rootloss < indiv.nmass_root) {
  228. indiv.nmass_root -= indiv.daily_nmass_rootloss;
  229. // 50% of the N in the lost root is retranslocated, the rest is going to litter.
  230. cropindiv.nmass_agpool += indiv.daily_nmass_rootloss * 0.5;
  231. patch.soil.sompool[SOILMETA].nmass += indiv.daily_nmass_rootloss * 0.5;
  232. patch.fluxes.report_flux(Fluxes::VEGLITTERN, indiv.daily_nmass_rootloss * 0.5);
  233. }
  234. if (indiv.daily_cmass_rootloss > 0.0){
  235. patch.is_litter_day = true;
  236. }
  237. cropindiv.dcmass_ho = (indiv.dnpp + cmass_extra) * patchpft.cropphen->f_alloc_horg;
  238. cropindiv.dcmass_plant = cropindiv.dcmass_ho + cropindiv.dcmass_root + cropindiv.dcmass_stem + cropindiv.dcmass_leaf;
  239. // Add the daily increments to the annual and growing season variables.
  240. cropindiv.ycmass_leaf += cropindiv.dcmass_leaf;
  241. cropindiv.ycmass_root += cropindiv.dcmass_root;
  242. cropindiv.ycmass_ho += cropindiv.dcmass_ho;
  243. cropindiv.ycmass_plant += cropindiv.dcmass_plant;
  244. cropindiv.grs_cmass_leaf += cropindiv.dcmass_leaf;
  245. // 40% of the assimilates that goes to stem is put into the fast C pool (Sec. 2.1.1 Olin 2015)
  246. cropindiv.grs_cmass_stem += (1.0 - 0.4) * cropindiv.dcmass_stem;
  247. cropindiv.ycmass_stem += (1.0 - 0.4) * cropindiv.dcmass_stem;
  248. cropindiv.grs_cmass_agpool += 0.4 * cropindiv.dcmass_stem;
  249. cropindiv.ycmass_agpool += 0.4 * cropindiv.dcmass_stem;
  250. cropindiv.grs_cmass_root += cropindiv.dcmass_root;
  251. cropindiv.grs_cmass_ho += cropindiv.dcmass_ho;
  252. cropindiv.grs_cmass_plant += cropindiv.dcmass_plant;
  253. double ndemand_ho = 0.0;
  254. // The non-structural N that is potentially available for retranslocation in leaves, roots and stem.
  255. double avail_leaf_N = max(0.0, (1.0 / indiv.cton_leaf(false) - 1.0 / indiv.pft.cton_leaf_max) * indiv.cmass_leaf_today());
  256. double avail_root_N = max(0.0, (1.0 / indiv.cton_root(false) - 1.0 / indiv.pft.cton_root_max) * indiv.cmass_root_today());
  257. double avail_stem_N = max(0.0,cropindiv.nmass_agpool - 1.0 / indiv.pft.cton_stem_max * cropindiv.grs_cmass_stem);
  258. double avail_N = avail_leaf_N + avail_root_N + avail_stem_N;
  259. if (avail_N > 0.0 && cropindiv.dcmass_ho > 0.0) {
  260. ndemand_ho = cropindiv.dcmass_ho / indiv.pft.cton_leaf_avr;
  261. }
  262. // N mass to be translocated from leaves and roots
  263. double trans_leaf_N = 0.0;
  264. double trans_root_N = 0.0;
  265. if (ndemand_ho > 0.0) {
  266. if(avail_stem_N > 0.0) {
  267. if (ndemand_ho > avail_stem_N) {
  268. ndemand_ho -= avail_stem_N;
  269. cropindiv.dnmass_ho += avail_stem_N;
  270. cropindiv.nmass_agpool -= avail_stem_N;
  271. }
  272. else {
  273. cropindiv.nmass_agpool -= ndemand_ho;
  274. cropindiv.dnmass_ho += ndemand_ho;
  275. ndemand_ho = 0.0;
  276. }
  277. }
  278. // "willingness" to let go of the N in the organ to meet the demand from the storage organ, Eq. 17 Olin 2015
  279. double y0 = (1.0 / indiv.pft.cton_leaf_min + 1.0 / indiv.pft.cton_leaf_avr) / 2.0;
  280. double y = 1.0 / indiv.cton_leaf(false);
  281. double y2 = 1.0 / (1.0 * indiv.pft.cton_leaf_max);
  282. double z = (y0 - y)/(y0 - y2);
  283. // The leaves willingness to contribute to meet the storage organs N demand.
  284. double w_l = 1.0 - max(0.0, min(1.0, pow(1.0 - z, 2.0)));
  285. y0 = 1.0 / indiv.pft.cton_root_avr;
  286. y = 1.0 / indiv.cton_root(false);
  287. y2 = 1.0 / (1.0 * indiv.pft.cton_root_max);
  288. z = (y0 - y) / (y0 - y2);
  289. // The roots willingness to contribute to meet the storage organs N demand.
  290. double w_r = 1.0 - max(0.0, min(1.0, pow(1.0 - z, 2.0)));
  291. double w_s = w_r + w_l;
  292. // The total willingness to contribute to meet the storage organs N demand.
  293. double w = min(1.0, w_s);
  294. // If there is N availble for retranslocation in the roots or leaves
  295. if(w_s > 0.0) {
  296. trans_leaf_N = max(0.0, w_l * w * ndemand_ho / w_s);
  297. trans_root_N = max(0.0, w_r * w * ndemand_ho / w_s);
  298. // Compare the N the organ is willing to let go of to N available in the organ, based on the C:N.
  299. if(trans_leaf_N > avail_leaf_N) {
  300. trans_leaf_N = avail_leaf_N;
  301. }
  302. if(trans_root_N > avail_root_N) {
  303. trans_root_N = avail_root_N;
  304. }
  305. // Add the N to the havestable organ
  306. cropindiv.dnmass_ho += trans_leaf_N;
  307. cropindiv.dnmass_ho += trans_root_N;
  308. }
  309. }
  310. // Subtract the N that has been added to theharvestable organ from the donor organs
  311. indiv.nmass_leaf -= trans_leaf_N;
  312. indiv.nmass_root -= trans_root_N;
  313. cropindiv.nmass_ho += cropindiv.dnmass_ho;
  314. }
  315. /// Daily allocation routine for crops without nitrogen limitation
  316. /** Allocates daily npp to leaf, roots and harvestable organs
  317. * SWAT equations are from Neitsch et al. 2002.
  318. */
  319. void allocation_crop(Individual& indiv, double cmass_seed, double nmass_seed) {
  320. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  321. Patch& patch = indiv.vegetation.patch;
  322. Patchpft& patchpft = patch.pft[indiv.pft.id];
  323. cropphen_struct& ppftcrop = *(patchpft.get_cropphen());
  324. nmass_seed = 0.0; // for compatibility with previous code, doesn't affect crop growth
  325. // report seed flux
  326. indiv.report_flux(Fluxes::SEEDC, -cmass_seed);
  327. indiv.report_flux(Fluxes::SEEDN, -nmass_seed);
  328. // add seed carbon
  329. cropindiv.grs_cmass_plant += cmass_seed;
  330. cropindiv.ycmass_plant += cmass_seed;
  331. cropindiv.dcmass_plant += cmass_seed;
  332. // add seed nitrogen
  333. indiv.nmass_leaf += nmass_seed / 2.0;
  334. indiv.nmass_root += nmass_seed / 2.0;
  335. // add today's npp
  336. cropindiv.dcmass_plant += indiv.dnpp;
  337. cropindiv.grs_cmass_plant += indiv.dnpp;
  338. cropindiv.ycmass_plant += indiv.dnpp;
  339. // allocation to roots
  340. double froot = indiv.pft.frootstart -(indiv.pft.frootstart - indiv.pft.frootend) * ppftcrop.fphu; // SWAT 5:2,1,21
  341. double grs_cmass_root_old = cropindiv.grs_cmass_root;
  342. cropindiv.grs_cmass_root = froot * cropindiv.grs_cmass_plant;
  343. cropindiv.dcmass_root = cropindiv.grs_cmass_root - grs_cmass_root_old;
  344. cropindiv.ycmass_root += cropindiv.dcmass_root;
  345. // allocation to harvestable organs
  346. double grs_cmass_ag = (1.0-froot) * cropindiv.grs_cmass_plant;
  347. double grs_cmass_ho_old = cropindiv.grs_cmass_ho;
  348. if (indiv.pft.hiopt <= 1.0)
  349. cropindiv.grs_cmass_ho = ppftcrop.hi * grs_cmass_ag; // SWAT 5:2.4.2, 5:2.4.4
  350. else // below-ground harvestable organs
  351. cropindiv.grs_cmass_ho = (1.0 - 1.0 / (1.0 + ppftcrop.hi)) * cropindiv.grs_cmass_plant; // SWAT 5:2.4.3 8
  352. cropindiv.dcmass_ho = cropindiv.grs_cmass_ho - grs_cmass_ho_old;
  353. cropindiv.ycmass_ho += cropindiv.dcmass_ho;
  354. // allocation to leaves
  355. double grs_cmass_leaf_old = cropindiv.grs_cmass_leaf;
  356. cropindiv.grs_cmass_leaf = cropindiv.grs_cmass_plant - cropindiv.grs_cmass_root - cropindiv.grs_cmass_ho;
  357. cropindiv.dcmass_leaf = cropindiv.grs_cmass_leaf - grs_cmass_leaf_old;
  358. cropindiv.ycmass_leaf += cropindiv.dcmass_leaf;
  359. // allocation to above-ground pool (currently not used)
  360. cropindiv.dcmass_agpool = cropindiv.dcmass_plant - cropindiv.dcmass_root - cropindiv.dcmass_leaf - cropindiv.dcmass_ho;
  361. cropindiv.grs_cmass_agpool = cropindiv.grs_cmass_plant - cropindiv.grs_cmass_root - cropindiv.grs_cmass_leaf - cropindiv.grs_cmass_ho;
  362. cropindiv.ycmass_agpool = cropindiv.ycmass_plant - cropindiv.ycmass_root - cropindiv.ycmass_leaf - cropindiv.ycmass_ho;
  363. if (!largerthanzero(cropindiv.grs_cmass_agpool, -9))
  364. cropindiv.grs_cmass_agpool = 0.0;
  365. if (!largerthanzero(cropindiv.ycmass_agpool, -9))
  366. cropindiv.ycmass_agpool = 0.0;
  367. }
  368. /// Daily growth routine for crops
  369. /** Allocates daily npp to leaf, roots and harvestable organs by calling
  370. * allocation_crop_nlim() or allocation_crop()
  371. * Requires updated value of fphu and hi.
  372. */
  373. void growth_crop_daily(Patch& patch) {
  374. if(date.day == 0)
  375. patch.nharv = 0;
  376. patch.isharvestday = false;
  377. double nharv_today = 0;
  378. Landcover& lc = patch.stand.get_gridcell().landcover;
  379. Vegetation& vegetation = patch.vegetation;
  380. vegetation.firstobj();
  381. while (vegetation.isobj) {
  382. Individual& indiv = vegetation.getobj();
  383. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  384. Patchpft& patchpft = patch.pft[indiv.pft.id];
  385. cropphen_struct& ppftcrop = *(patchpft.get_cropphen());
  386. if(date.day == 0) {
  387. cropindiv.ycmass_plant = 0.0;
  388. cropindiv.ycmass_leaf = 0.0;
  389. cropindiv.ycmass_root = 0.0;
  390. cropindiv.ycmass_ho = 0.0;
  391. cropindiv.ycmass_agpool = 0.0;
  392. cropindiv.ycmass_dead_leaf = 0.0;
  393. cropindiv.ycmass_stem = 0.0;
  394. cropindiv.ynmass_leaf = 0.0;
  395. cropindiv.ynmass_root = 0.0;
  396. cropindiv.ynmass_ho = 0.0;
  397. cropindiv.ynmass_agpool = 0.0;
  398. cropindiv.ynmass_dead_leaf = 0.0;
  399. cropindiv.harv_cmass_plant = 0.0;
  400. cropindiv.harv_cmass_root = 0.0;
  401. cropindiv.harv_cmass_leaf = 0.0;
  402. cropindiv.harv_cmass_ho = 0.0;
  403. cropindiv.harv_cmass_agpool = 0.0;
  404. cropindiv.cmass_ho_harvest[0] = 0.0;
  405. cropindiv.cmass_ho_harvest[1] = 0.0;
  406. cropindiv.cmass_leaf_max = 0.0;
  407. if(indiv.pft.phenology == ANY && !cropindiv.isintercropgrass) { // zero of normal cc3g/cc4g cmass arbitrarily at new year
  408. cropindiv.grs_cmass_plant = 0.0;
  409. cropindiv.grs_cmass_root = 0.0;
  410. cropindiv.grs_cmass_ho = 0.0;
  411. cropindiv.grs_cmass_leaf = 0.0;
  412. cropindiv.grs_cmass_agpool = 0.0;
  413. }
  414. }
  415. cropindiv.dcmass_plant = 0.0;
  416. cropindiv.dcmass_leaf = 0.0;
  417. cropindiv.dcmass_root = 0.0;
  418. cropindiv.dcmass_ho = 0.0;
  419. cropindiv.dcmass_agpool = 0.0;
  420. cropindiv.dcmass_stem = 0.0;
  421. cropindiv.dnmass_ho = 0.0;
  422. // true crop allocation
  423. if(indiv.pft.phenology == CROPGREEN) {
  424. if(ppftcrop.growingseason) {
  425. double cmass_seed = 0.0;
  426. double nmass_seed = 0.0;
  427. if (DELAYED_SEEDCARBON) {
  428. // Seed carbon; portion the seed carbon over a 10-day period.
  429. if(dayinperiod(date.day, ppftcrop.sdate, stepfromdate(ppftcrop.sdate, 9))) {
  430. cmass_seed = 0.1 * CMASS_SEED;
  431. nmass_seed = 0.1 * CMASS_SEED / indiv.pft.cton_leaf_min;
  432. }
  433. } else {
  434. // add seed carbon on sowing date
  435. if(date.day == ppftcrop.sdate) {
  436. cmass_seed = CMASS_SEED;
  437. nmass_seed = CMASS_SEED / indiv.pft.cton_leaf_min;
  438. }
  439. }
  440. if(ifnlim)
  441. allocation_crop_nlim(indiv, cmass_seed, nmass_seed);
  442. else
  443. allocation_crop(indiv, cmass_seed, nmass_seed);
  444. // save this year's maximum leaf carbon mass
  445. if(cropindiv.grs_cmass_leaf > cropindiv.cmass_leaf_max)
  446. cropindiv.cmass_leaf_max = cropindiv.grs_cmass_leaf;
  447. // save leaf carbon mass at the beginning of senescence
  448. if(date.day == ppftcrop.sendate)
  449. cropindiv.cmass_leaf_sen = cropindiv.grs_cmass_leaf;
  450. // Check that no plant cmass or nmass is negative, if so, and correct fluxes
  451. double negative_cmass = indiv.check_C_mass();
  452. if(largerthanzero(negative_cmass, -14))
  453. dprintf("Year %d day %d Stand %d indiv %d: Negative main crop C mass in growth_crop_daily: %.15f\n", date.year, date.day, indiv.vegetation.patch.stand.id, indiv.id, -negative_cmass);
  454. double negative_nmass = indiv.check_N_mass();
  455. if(largerthanzero(negative_nmass, -14))
  456. dprintf("Year %d day %d Stand %d indiv %d: Negative main crop N mass in growth_crop_daily: %.15f\n", date.year, date.day, indiv.vegetation.patch.stand.id, indiv.id, -negative_nmass);
  457. }
  458. else if(date.day == ppftcrop.hdate) {
  459. patch.stand.isrotationday = true;
  460. cropindiv.harv_cmass_plant += cropindiv.grs_cmass_plant;
  461. cropindiv.harv_cmass_root += cropindiv.grs_cmass_root;
  462. cropindiv.harv_cmass_ho += cropindiv.grs_cmass_ho;
  463. cropindiv.harv_cmass_leaf += cropindiv.grs_cmass_leaf;
  464. cropindiv.harv_cmass_agpool += cropindiv.grs_cmass_agpool;
  465. cropindiv.harv_cmass_stem += cropindiv.grs_cmass_stem;
  466. cropindiv.harv_nmass_root += indiv.nmass_root;
  467. cropindiv.harv_nmass_ho += cropindiv.nmass_ho;
  468. cropindiv.harv_nmass_leaf += indiv.nmass_leaf;
  469. cropindiv.harv_nmass_agpool += cropindiv.nmass_agpool;
  470. // dead_leaf to be addad
  471. if(ppftcrop.nharv == 1)
  472. cropindiv.cmass_ho_harvest[0] = cropindiv.grs_cmass_ho;
  473. else if(ppftcrop.nharv == 2)
  474. cropindiv.cmass_ho_harvest[1] = cropindiv.grs_cmass_ho;
  475. patch.isharvestday = true;
  476. if(indiv.has_daily_turnover()) {
  477. if(lc.updated && patchpft.cropphen->nharv == 1)
  478. scale_indiv(indiv, true);
  479. harvest_crop(indiv, indiv.pft, indiv.alive, indiv.cropindiv->isintercropgrass, true);
  480. patch.is_litter_day = true;
  481. }
  482. patch.nharv++;
  483. cropindiv.grs_cmass_plant = 0.0;
  484. cropindiv.grs_cmass_root = 0.0;
  485. cropindiv.grs_cmass_ho = 0.0;
  486. cropindiv.grs_cmass_leaf = 0.0;
  487. cropindiv.grs_cmass_agpool = 0.0;
  488. cropindiv.cmass_leaf_sen = 0.0;
  489. cropindiv.grs_cmass_stem = 0.0;
  490. cropindiv.grs_cmass_dead_leaf = 0.0;
  491. cropindiv.nmass_dead_leaf = 0.0;
  492. cropindiv.nmass_agpool = 0.0;
  493. indiv.nmass_leaf = 0.0;
  494. indiv.nmass_root = 0.0;
  495. cropindiv.nmass_ho = 0.0;
  496. }
  497. }
  498. // crop grass allocation
  499. // NB: Only intercrop grass enters here ! normal cc3g/cc4g grass is treated just like natural grass
  500. else if(indiv.pft.phenology == ANY && indiv.pft.id != patch.stand.pftid) {
  501. if(ppftcrop.growingseason) {
  502. cropindiv.dcmass_plant = 0.0;
  503. // add seed carbon
  504. if(!indiv.continous_grass() && date.day == patch.pft[patch.stand.pftid].cropphen->bicdate
  505. || indiv.continous_grass() && date.day == stepfromdate(indiv.last_turnover_day, 1)) {
  506. double cmass_seed = CMASS_SEED;
  507. cropindiv.grs_cmass_plant += cmass_seed;
  508. cropindiv.ycmass_plant += cmass_seed;
  509. cropindiv.dcmass_plant += cmass_seed;
  510. double nmass_seed = cmass_seed / indiv.pft.cton_leaf_min;
  511. indiv.nmass_leaf += nmass_seed / 2.0;
  512. indiv.nmass_root += nmass_seed / 2.0;
  513. indiv.report_flux(Fluxes::SEEDC, -cmass_seed);
  514. indiv.report_flux(Fluxes::SEEDN, -nmass_seed);
  515. indiv.last_turnover_day = -1;
  516. }
  517. cropindiv.dcmass_plant += indiv.dnpp;
  518. cropindiv.grs_cmass_plant += indiv.dnpp;
  519. cropindiv.ycmass_plant += indiv.dnpp;
  520. indiv.ltor = indiv.wscal_mean() * indiv.pft.ltor_max;
  521. // allocation to roots
  522. double froot = 1.0 / (1.0 + indiv.ltor);
  523. double grs_cmass_root_old = cropindiv.grs_cmass_root;
  524. // Cumulative wscal-dependent root increase
  525. cropindiv.grs_cmass_root = froot * cropindiv.grs_cmass_plant;
  526. cropindiv.dcmass_root = cropindiv.grs_cmass_root - grs_cmass_root_old;
  527. cropindiv.ycmass_root += cropindiv.dcmass_root;
  528. // allocation to leaves
  529. double grs_cmass_leaf_old = cropindiv.grs_cmass_leaf;
  530. cropindiv.grs_cmass_leaf = cropindiv.grs_cmass_plant - cropindiv.grs_cmass_root;
  531. cropindiv.dcmass_leaf = cropindiv.grs_cmass_leaf - grs_cmass_leaf_old;
  532. cropindiv.ycmass_leaf += cropindiv.dcmass_leaf;
  533. // Check that no plant cmass is negative, if so, zero cmass and correct C fluxes
  534. double negative_cmass = indiv.check_C_mass();
  535. // save this year's maximum leaf carbon mass
  536. if(cropindiv.grs_cmass_leaf > cropindiv.cmass_leaf_max)
  537. cropindiv.cmass_leaf_max = cropindiv.grs_cmass_leaf;
  538. }
  539. else if(date.day == patch.pft[patch.stand.pftid].get_cropphen()->eicdate) {
  540. cropindiv.harv_cmass_plant += cropindiv.grs_cmass_plant;
  541. cropindiv.harv_cmass_root += cropindiv.grs_cmass_root;
  542. cropindiv.harv_cmass_leaf += cropindiv.grs_cmass_leaf;
  543. cropindiv.harv_cmass_ho += cropindiv.grs_cmass_ho;
  544. cropindiv.harv_cmass_agpool += cropindiv.grs_cmass_agpool;
  545. cropindiv.harv_nmass_root += indiv.nmass_root;
  546. cropindiv.harv_nmass_ho += cropindiv.nmass_ho;
  547. cropindiv.harv_nmass_leaf += indiv.nmass_leaf;
  548. cropindiv.harv_nmass_agpool += cropindiv.nmass_agpool;
  549. ppftcrop.nharv++;
  550. patch.isharvestday = true;
  551. nharv_today++;
  552. if(nharv_today > 1) // In case of both C3 and C4 growing
  553. patch.nharv--;
  554. if(indiv.has_daily_turnover()) {
  555. if(lc.updated && patchpft.cropphen->nharv == 1)
  556. scale_indiv(indiv, true);
  557. harvest_crop(indiv, indiv.pft, indiv.alive, indiv.cropindiv->isintercropgrass, true);
  558. patch.is_litter_day = true;
  559. }
  560. else {
  561. cropindiv.grs_cmass_root = 0.0;
  562. cropindiv.grs_cmass_ho = 0.0;
  563. cropindiv.grs_cmass_leaf = 0.0;
  564. cropindiv.grs_cmass_agpool = 0.0;
  565. }
  566. patch.nharv++;
  567. cropindiv.grs_cmass_plant = cropindiv.grs_cmass_root + cropindiv.grs_cmass_leaf;
  568. }
  569. if(indiv.continous_grass() && indiv.is_turnover_day()
  570. || ppftcrop.growingdays > date.year_length() && (indiv.is_turnover_day() || date.islastmonth && date.islastday)) {
  571. indiv.last_turnover_day = date.day;
  572. ppftcrop.growingdays = 0;
  573. cropindiv.harv_cmass_plant += cropindiv.grs_cmass_plant;
  574. cropindiv.harv_cmass_root += cropindiv.grs_cmass_root;
  575. cropindiv.harv_cmass_leaf += cropindiv.grs_cmass_leaf;
  576. cropindiv.harv_cmass_ho += cropindiv.grs_cmass_ho;
  577. cropindiv.harv_cmass_agpool += cropindiv.grs_cmass_agpool;
  578. cropindiv.harv_nmass_root += indiv.nmass_root;
  579. cropindiv.harv_nmass_ho += cropindiv.nmass_ho;
  580. cropindiv.harv_nmass_leaf += indiv.nmass_leaf;
  581. cropindiv.harv_nmass_agpool += cropindiv.nmass_agpool;
  582. ppftcrop.nharv++;
  583. patch.isharvestday = true;
  584. nharv_today++;
  585. if(nharv_today > 1) // In case of both C3 and C4 growing
  586. patch.nharv--;
  587. if(indiv.has_daily_turnover()) {
  588. if(lc.updated && patchpft.cropphen->nharv == 1)
  589. scale_indiv(indiv, true);
  590. turnover_grass(indiv);
  591. patch.is_litter_day = true;
  592. }
  593. else {
  594. cropindiv.grs_cmass_root = 0.0;
  595. cropindiv.grs_cmass_ho = 0.0;
  596. cropindiv.grs_cmass_leaf = 0.0;
  597. }
  598. patch.nharv++;
  599. cropindiv.grs_cmass_plant = cropindiv.grs_cmass_root + cropindiv.grs_cmass_leaf;
  600. cropindiv.grs_cmass_agpool = 0.0;
  601. }
  602. }
  603. vegetation.nextobj();
  604. }
  605. }
  606. /// Handles daily crop allocation and daily lai calculation
  607. /** Simple allocation based on heat unit accumulation.
  608. * LAI is set directly after allocation from leaf carbon mass.
  609. */
  610. void growth_daily(Patch& patch) {
  611. if(patch.stand.landcover == CROPLAND) {
  612. // allocate daily npp to leaf, roots and harvestable organs
  613. growth_crop_daily(patch);
  614. // update patchpft.lai_daily and fpc_daily
  615. lai_crop(patch);
  616. }
  617. }
  618. /// Handles yearly cropland lai and fpc calculation, called from allometry()
  619. /** Uses cmass_leaf_max for true crops and lai from whole-year grass stands for
  620. * cover-crop grass if present (NB. grass pft names must be as described in code !).
  621. * If grass pft not present or found in other stands, pft laimax is used.
  622. */
  623. void allometry_crop(Individual& indiv) {
  624. // crop grass compatible with natural grass
  625. if(indiv.pft.phenology == ANY) {
  626. indiv.lai_indiv = indiv.cmass_leaf * indiv.pft.sla;
  627. // For intercrop grass, use LAI of parent grass in its own stand.
  628. if(indiv.cropindiv->isintercropgrass) {
  629. bool done = false;
  630. Gridcell& gridcell = indiv.vegetation.patch.stand.get_gridcell();
  631. // First look in PASTURE.
  632. if(gridcell.landcover.frac[PASTURE] > 0.0) {
  633. for(unsigned int i = 0; i < gridcell.size() && !done; i++) {
  634. Stand& stand=gridcell[i];
  635. if(stand.landcover == PASTURE) {
  636. for(unsigned int j = 0; j < stand.nobj && !done; j++) {
  637. Patch& patch = stand[j];
  638. Vegetation& vegetation = patch.vegetation;
  639. for(unsigned int k = 0; k < vegetation.nobj && !done; k++) {
  640. Individual& grass_indiv = vegetation[k];
  641. if(grass_indiv.pft.phenology == ANY && grass_indiv.pft.pathway == indiv.pft.pathway) {
  642. indiv.lai_indiv = grass_indiv.lai_indiv;
  643. done = true;
  644. }
  645. }
  646. }
  647. }
  648. }
  649. }
  650. // If PASTURE landcover not used, look for crop stand with pasture grass.
  651. else {
  652. double highest_grass_lai = 0.0;
  653. double grass_cmass_leaf_sum = 0.0;
  654. // Get sum of intercrop grass cmass_leaf in this patch
  655. Vegetation& vegetation_self = indiv.vegetation;
  656. for(unsigned int k = 0; k < vegetation_self.nobj; k++) {
  657. Individual& indiv_veg = vegetation_self[k];
  658. if(indiv_veg.cropindiv->isintercropgrass)
  659. grass_cmass_leaf_sum += indiv_veg.cropindiv->cmass_leaf_max;
  660. }
  661. // Find highest lai in crop grass stands
  662. for(unsigned int i = 0; i < gridcell.size(); i++) {
  663. Stand& stand = gridcell[i];
  664. if(stand.landcover == CROPLAND && !stand.is_true_crop_stand()) {
  665. for(unsigned int j = 0; j < stand.nobj && !done; j++) {
  666. Patch& patch = stand[j];
  667. Vegetation& vegetation = patch.vegetation;
  668. for(unsigned int k = 0; k < vegetation.nobj && !done; k++) {
  669. Individual& grass_indiv = vegetation[k];
  670. if(grass_indiv.pft.phenology == ANY) {
  671. if(grass_indiv.lai_indiv > highest_grass_lai)
  672. highest_grass_lai = grass_indiv.lai_indiv;
  673. }
  674. }
  675. }
  676. }
  677. }
  678. if(grass_cmass_leaf_sum > 0.0)
  679. indiv.lai_indiv = indiv.cropindiv->cmass_leaf_max / grass_cmass_leaf_sum * highest_grass_lai;
  680. else
  681. indiv.lai_indiv = highest_grass_lai;
  682. if(highest_grass_lai)
  683. done = true;
  684. }
  685. // If no grass stand found in either cropland or pasture, use laimax value.
  686. if(!done) {
  687. double highest_grass_lai = 0.0;
  688. double grass_cmass_leaf_sum = 0.0;
  689. // Get sum of intercrop grass cmass_leaf and highest default laimax in this patch
  690. Vegetation& vegetation_self = indiv.vegetation;
  691. for(unsigned int k = 0; k < vegetation_self.nobj; k++) {
  692. Individual& indiv_veg = vegetation_self[k];
  693. if(indiv_veg.cropindiv->isintercropgrass) {
  694. grass_cmass_leaf_sum += indiv_veg.cropindiv->cmass_leaf_max;
  695. if(indiv_veg.pft.laimax > highest_grass_lai)
  696. highest_grass_lai = indiv_veg.pft.laimax;
  697. }
  698. }
  699. if(grass_cmass_leaf_sum > 0.0)
  700. indiv.lai_indiv = indiv.cropindiv->cmass_leaf_max / grass_cmass_leaf_sum * highest_grass_lai;
  701. else
  702. indiv.lai_indiv = highest_grass_lai;
  703. }
  704. }
  705. if(indiv.lai_indiv < 0.0)
  706. fail("lai_indiv negative for %s in stand %d year %d in growth: %f\n",
  707. (char*)indiv.pft.name, indiv.vegetation.patch.stand.id, date.year, indiv.lai_indiv);
  708. // FPC (Eqn 10)
  709. indiv.fpc = 1.0 - lambertbeer(indiv.lai_indiv);
  710. // Stand-level LAI
  711. indiv.lai = indiv.lai_indiv;
  712. }
  713. else { // cropgreen
  714. if (!negligible(indiv.cropindiv->cmass_leaf_max)) {
  715. // Grass "individual" LAI (Eqn 11)
  716. indiv.lai_indiv = indiv.cropindiv->cmass_leaf_max * indiv.pft.sla;
  717. // FPC (Eqn 10)
  718. // indiv.fpc = 1.0 - lambertbeer(indiv.lai_indiv);
  719. indiv.fpc = 1.0;
  720. // Stand-level LAI
  721. indiv.lai = indiv.lai_indiv;
  722. }
  723. }
  724. }
  725. /// Transfer of this year's growth (ycmass_xxx) to cmass_xxx_inc
  726. /** and pasture grass grown in cropland.
  727. * OUTPUT PARAMETERS
  728. * \param cmass_leaf_inc leaf C biomass (kgC/m2)
  729. * \param cmass_root_inc fine root C biomass (kgC/m2)
  730. * \param cmass_ho_inc harvestable organ C biomass (kgC/m2)
  731. * \param cmass_agpool_inc above-ground pool C biomass (kgC/m2)
  732. * \param cmass_stem_inc stem C biomass (kgC/m2)
  733. */
  734. void growth_crop_year(Individual& indiv, double& cmass_leaf_inc, double& cmass_root_inc, double& cmass_ho_inc, double& cmass_agpool_inc, double& cmass_stem_inc) {
  735. // true crop growth and grass intercrop growth; NB: bminit (cmass_repr & cmass_excess subtracted) not used !
  736. if(indiv.has_daily_turnover()) {
  737. indiv.cmass_leaf = 0.0;
  738. indiv.cmass_root = 0.0;
  739. indiv.cropindiv->cmass_ho = 0.0;
  740. indiv.cropindiv->cmass_agpool = 0.0;
  741. indiv.cropindiv->cmass_stem = 0.0;
  742. // Not completely accurate here when comparing this year's cmass after turnover with cmass increase (ycmass),
  743. // which could be from the preceding season, but probably OK, since values are not used for C balance.
  744. if(indiv.continous_grass()) {
  745. indiv.cmass_leaf = indiv.cmass_leaf_post_turnover;
  746. indiv.cmass_root = indiv.cmass_root_post_turnover;
  747. }
  748. }
  749. cmass_leaf_inc = indiv.cropindiv->ycmass_leaf + indiv.cropindiv->ycmass_dead_leaf;
  750. cmass_root_inc = indiv.cropindiv->ycmass_root;
  751. cmass_ho_inc = indiv.cropindiv->ycmass_ho;
  752. cmass_agpool_inc = indiv.cropindiv->ycmass_agpool;
  753. cmass_stem_inc = indiv.cropindiv->ycmass_stem;
  754. return;
  755. }
  756. //////////////////////////////////////////////////////////////////////////////////////////
  757. // REFERENCES
  758. //
  759. // Neitsch SL, Arnold JG, Kiniry JR et al.2002 Soil and Water Assessment Tool, Theorethical
  760. // Documentation + User's Manual. USDA_ARS-SR Grassland, Soil and Water Research Laboratory.
  761. // Agricultural Reasearch Service, Temple,Tx, US.
  762. // S. Olin, G. Schurgers, M. Lindeskog, D. W�rlind, B. Smith, P. Bodin, J. Holm�r, and A. Arneth. 2015
  763. // Biogeosciences 12, 2489-2515. Modelling the response of yields and tissue C:N to changes in
  764. // atmospheric CO2 and N management in the main wheat regions of western Europe