management.cpp 53 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340
  1. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  2. /// \file management.cpp
  3. /// \brief Harvest functions for cropland, managed forest and pasture
  4. /// \author Mats Lindeskog
  5. /// $Date: $
  6. ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  7. #include "landcover.h"
  8. #include "management.h"
  9. /// Harvest function used for managed forest and for clearing natural vegetation at land use change
  10. /** A fraction of trees is cut down (frac_cut)
  11. * A fraction of wood is harvested (pft.harv_eff) and returned as acflux_harvest
  12. * A fraction of harvested wood (pft.harvest_slow_frac) is returned as harvested_products_slow
  13. * The rest, including leaves and roots, is returned as litter, unless a fraction of twigs or roots removed.
  14. * Called from landcover_dynamics() first day of the year if any natural vegetation is transferred to another land use.
  15. * INPUT PARAMETER
  16. * \param frac_cut fraction of trees cut
  17. * \param harv_eff harvest efficiency
  18. * \param res_outtake_twig removed twig fraction
  19. * \param res_outtake_coarse_root removed course root fraction
  20. * INPUT/OUTPUT PARAMETERS
  21. * \param Harvest_CN& i struct containing the following indiv-specific public members:
  22. * - cmass_leaf leaf C biomass (kgC/m2)
  23. * - cmass_root fine root C biomass (kgC/m2)
  24. * - cmass_sap sapwood C biomass (kgC/m2)
  25. * - cmass_heart heartwood C biomass (kgC/m2)
  26. * - cmass_debt C "debt" (retrospective storage) (kgC/m2)
  27. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  28. * - nmass_root fine root nitrogen biomass (kgN/m2)
  29. * - nmass_sap sapwood nitrogen biomass (kgC/m2)
  30. * - nmass_heart heartwood nitrogen biomass (kgC/m2)
  31. * - nstore_labile labile nitrogen storage (kgC/m2)
  32. * - nstore_longterm longterm nitrogen storage (kgC/m2)
  33. * OUTPUT PARAMETERS
  34. * \param Harvest_CN& i struct containing the following patchpft-specific public members:
  35. * - litter_leaf new leaf C litter (kgC/m2)
  36. * - litter_root new root C litter (kgC/m2)
  37. * - litter_sap new sapwood C litter (kgC/m2)
  38. * - litter_heart new heartwood C litter (kgC/m2)
  39. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  40. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  41. * - nmass_litter_sap new sapwood nitrogen litter (kgN/m2)
  42. * - nmass_litter_heart new heartwood nitrogen litter (kgN/m2)
  43. * ,and the following patch-level public members:
  44. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  45. * - harvested_products_slow harvest products to slow pool (kgC/m2)
  46. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  47. * - harvested_products_slow_nmass harvest nitrogen products to slow pool (kgC/m2)
  48. */
  49. void harvest_wood(Harvest_CN& i, Pft& pft, bool alive, double frac_cut, double harv_eff, double res_outtake_twig, double res_outtake_coarse_root) {
  50. double harvest = 0.0;
  51. i.dcflux_harvest_wood = 0.0;
  52. double residue_outtake = 0.0;
  53. /// Fraction of wood cmass that are stems
  54. double stem_frac = 0.65; // Temporary values, should be pft-specific
  55. /// Fraction of wood cmass that are twigs
  56. double twig_frac = 0.13;
  57. /// Fraction of wood cmass that are coarse roots
  58. double coarse_root_frac = 1.0 - stem_frac - twig_frac; // 0.22 with default stem_frac and twig_frac values
  59. /// Fraction of leaves adhering to twigs at the time of removal
  60. double adhering_leaf_frac = 0.75;
  61. // only harvest trees
  62. if (pft.lifeform == GRASS)
  63. return;
  64. // all root carbon and nitrogen goes to litter
  65. if (alive) {
  66. i.litter_root += i.cmass_root * frac_cut;
  67. i.dcflux_harvest_wood_res += i.cmass_root * frac_cut;
  68. i.cmass_root *= (1.0 - frac_cut);
  69. }
  70. i.nmass_litter_root += i.nmass_root * frac_cut;
  71. i.nmass_litter_root += (i.nstore_labile + i.nstore_longterm) * frac_cut;
  72. i.nmass_root *= (1.0 - frac_cut);
  73. i.nstore_labile *= (1.0 - frac_cut);
  74. i.nstore_longterm *= (1.0 - frac_cut);
  75. if (alive) {
  76. // Carbon:
  77. if (i.cmass_debt <= i.cmass_sap + i.cmass_heart) {
  78. // harvested stem wood
  79. harvest += harv_eff * stem_frac * (i.cmass_sap + i.cmass_heart - i.cmass_debt) * frac_cut;
  80. // harvested products not consumed (oxidised) this year put into harvested_products_slow
  81. if (ifslowharvestpool) {
  82. i.harvested_products_slow += harvest * pft.harvest_slow_frac;
  83. if (harvest * pft.harvest_slow_frac > 0.0)
  84. i.dcflux_product_wood = harvest * pft.harvest_slow_frac;
  85. harvest = harvest * (1.0 - pft.harvest_slow_frac);
  86. }
  87. // harvested products consumed (oxidised) this year put into acflux_harvest
  88. i.acflux_harvest += harvest;
  89. i.dcflux_harvest_wood += harvest;
  90. // removed leaves adhering to twigs
  91. residue_outtake += res_outtake_twig * adhering_leaf_frac * i.cmass_leaf * frac_cut;
  92. // removed twigs
  93. residue_outtake += res_outtake_twig * twig_frac * (i.cmass_sap + i.cmass_heart - i.cmass_debt) * frac_cut;
  94. // removed coarse roots
  95. residue_outtake += res_outtake_coarse_root * coarse_root_frac * (i.cmass_sap + i.cmass_heart - i.cmass_debt) * frac_cut;
  96. // removed residues are oxidised
  97. i.acflux_harvest += residue_outtake;
  98. // not removed residues are put into litter
  99. i.litter_leaf += i.cmass_leaf * (1.0 - res_outtake_twig * adhering_leaf_frac) * frac_cut;
  100. i.dcflux_harvest_wood_res += i.cmass_leaf * (1.0 - res_outtake_twig * adhering_leaf_frac) * frac_cut;
  101. double to_partition_sap = 0.0;
  102. double to_partition_heart = 0.0;
  103. if (i.cmass_heart >= i.cmass_debt) {
  104. to_partition_sap = i.cmass_sap;
  105. to_partition_heart = i.cmass_heart - i.cmass_debt;
  106. }
  107. else {
  108. to_partition_sap = i.cmass_sap + i.cmass_heart - i.cmass_debt;
  109. // dprintf("ATTENTION: pft %s: cmass_debt > cmass_heart; difference=%f\n", (char*)pft.name, i.cmass_debt-i.cmass_heart);
  110. }
  111. i.litter_sap += to_partition_sap * (1.0 - res_outtake_twig * twig_frac - res_outtake_coarse_root * coarse_root_frac - harv_eff * stem_frac) * frac_cut;
  112. i.litter_heart += to_partition_heart * (1.0 - res_outtake_twig * twig_frac - res_outtake_coarse_root * coarse_root_frac - harv_eff * stem_frac) * frac_cut;
  113. i.dcflux_harvest_wood_res += (to_partition_sap + to_partition_heart) * (1.0 - res_outtake_twig * twig_frac - res_outtake_coarse_root * coarse_root_frac - harv_eff * stem_frac) * frac_cut;
  114. }
  115. // debt larger than existing wood biomass
  116. else {
  117. double debt_excess = i.cmass_debt - (i.cmass_sap + i.cmass_heart);
  118. dprintf("ATTENTION: cmass_debt > i.cmass_sap + i.cmass_heart; debt_excess=%f\n", debt_excess);
  119. // i.debt_excess += debt_excess * frac_cut; // debt_excess currently not dealt with during wood harvest
  120. }
  121. // unharvested trees:
  122. i.cmass_leaf *= (1.0 - frac_cut);
  123. i.cmass_sap *= (1.0 - frac_cut);
  124. i.cmass_heart *= (1.0 - frac_cut);
  125. i.cmass_debt *= (1.0 - frac_cut);
  126. //Nitrogen:
  127. harvest = 0.0;
  128. // harvested products
  129. harvest += harv_eff * stem_frac * (i.nmass_sap + i.nmass_heart) * frac_cut;
  130. // harvested products not consumed this year put into harvested_products_slow_nmass
  131. if (ifslowharvestpool) {
  132. i.harvested_products_slow_nmass += harvest * pft.harvest_slow_frac;
  133. if (harvest * pft.harvest_slow_frac > 0.0)
  134. i.dnflux_product += harvest * pft.harvest_slow_frac;
  135. harvest = harvest * (1.0 - pft.harvest_slow_frac);
  136. }
  137. // harvested products consumed this year put into anflux_harvest
  138. i.anflux_harvest += harvest;
  139. residue_outtake = 0.0;
  140. // removed leaves adhering to twigs
  141. residue_outtake += res_outtake_twig * adhering_leaf_frac * i.nmass_leaf * frac_cut;
  142. // removed twigs
  143. residue_outtake += res_outtake_twig * twig_frac * (i.nmass_sap + i.nmass_heart) * frac_cut;
  144. // removed coarse roots
  145. residue_outtake += res_outtake_coarse_root * coarse_root_frac * (i.nmass_sap + i.nmass_heart) * frac_cut;
  146. // removed residues are oxidised
  147. i.anflux_harvest += residue_outtake;
  148. // not removed residues are put into litter
  149. i.nmass_litter_leaf += i.nmass_leaf * (1.0 - res_outtake_twig * adhering_leaf_frac) * frac_cut;
  150. i.nmass_litter_sap += i.nmass_sap * (1.0 - res_outtake_twig * twig_frac - res_outtake_coarse_root * coarse_root_frac - harv_eff * stem_frac) * frac_cut;
  151. i.nmass_litter_heart += i.nmass_heart * (1.0 - res_outtake_twig * twig_frac - res_outtake_coarse_root * coarse_root_frac - harv_eff * stem_frac) * frac_cut;
  152. // unharvested trees:
  153. i.nmass_leaf *= (1.0 - frac_cut);
  154. i.nmass_sap *= (1.0 - frac_cut);
  155. i.nmass_heart *= (1.0 - frac_cut);
  156. }
  157. }
  158. /// Harvest function used for managed forest and for clearing natural vegetation at land use change
  159. /** A fraction of trees is cut down (frac_cut)
  160. * A fraction of wood is harvested (pft.harv_eff) and returned as acflux_harvest
  161. * A fraction of harvested wood (pft.harvest_slow_frac) is returned as harvested_products_slow
  162. * The rest, including leaves and roots, is returned as litter.
  163. * Called from landcover_dynamics() first day of the year if any natural vegetation is transferred to another land use.
  164. *
  165. * This function copies variables from an individual and it's associated patchpft and patch to
  166. * a Harvest_CN struct, which is then passed on to the main harvest_crop function.
  167. * After the execution of the main harvest_crop function, the output variables are copied
  168. * back to the individual and patchpft and the patch-level fluxes are updated.
  169. *
  170. * INPUT PARAMETER
  171. * \param frac_cut fraction of trees cut
  172. * \param harv_eff harvest efficiency
  173. * \param res_outtake_twig removed twig fraction
  174. * \param res_outtake_coarse_root removed course root fraction
  175. * \param lc_change whether to save harvest in gridcell-level luc variable
  176. * INPUT/OUTPUT PARAMETERS
  177. * \param indiv reference to an Individual containing the following indiv-specific public members:
  178. * - cmass_leaf leaf C biomass (kgC/m2)
  179. * - cmass_root fine root C biomass (kgC/m2)
  180. * - cmass_sap sapwood C biomass (kgC/m2)
  181. * - cmass_heart heartwood C biomass (kgC/m2)
  182. * - cmass_debt C "debt" (retrospective storage) (kgC/m2)
  183. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  184. * - nmass_root fine root nitrogen biomass (kgN/m2)
  185. * - nmass_sap sapwood nitrogen biomass (kgC/m2)
  186. * - nmass_heart heartwood nitrogen biomass (kgC/m2)
  187. * - nstore_labile labile nitrogen storage (kgC/m2)
  188. * - nstore_longterm longterm nitrogen storage (kgC/m2)
  189. * OUTPUT PARAMETERS
  190. * \param indiv reference to an Individual containing the following patchpft-specific public members:
  191. * - litter_leaf new leaf C litter (kgC/m2)
  192. * - litter_root new root C litter (kgC/m2)
  193. * - litter_sap new sapwood C litter (kgC/m2)
  194. * - litter_heart new heartwood C litter (kgC/m2)
  195. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  196. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  197. * - nmass_litter_sap new sapwood nitrogen litter (kgN/m2)
  198. * - nmass_litter_heart new heartwood nitrogen litter (kgN/m2)
  199. * ,and the following patch-level public members:
  200. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  201. * - harvested_products_slow harvest products to slow pool (kgC/m2)
  202. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  203. * - harvested_products_slow_nmass harvest nitrogen products to slow pool (kgC/m2)
  204. */
  205. void harvest_wood(Individual& indiv, double frac_cut, double harv_eff, double res_outtake_twig,
  206. double res_outtake_coarse_root, bool lc_change) {
  207. Harvest_CN indiv_cp;
  208. indiv_cp.copy_from_indiv(indiv);
  209. harvest_wood(indiv_cp, indiv.pft, indiv.alive, frac_cut, harv_eff, res_outtake_twig, res_outtake_coarse_root);
  210. indiv_cp.copy_to_indiv(indiv, false, lc_change);
  211. if (!lc_change) {
  212. return;
  213. }
  214. Stand& stand = indiv.vegetation.patch.stand;
  215. Landcover& lc = stand.get_gridcell().landcover;
  216. lc.harvest_product += stand.get_gridcell_fraction() * indiv_cp.harvested_products_slow / (double)stand.nobj;
  217. lc.harvest_product_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.harvested_products_slow / (double)stand.nobj;
  218. lc.harvest_product_nmass += stand.get_gridcell_fraction() * indiv_cp.harvested_products_slow_nmass / (double)stand.nobj;
  219. lc.harvest_product_nmass_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.harvested_products_slow_nmass / (double)stand.nobj;
  220. lc.acflux_harvest_wood_res += stand.get_gridcell_fraction() * indiv_cp.dcflux_harvest_wood_res / (double)stand.nobj;
  221. lc.dcflux_landuse_change += stand.get_gridcell_fraction() * indiv_cp.acflux_harvest / (double)stand.nobj; // ecev3 - reset to 0 each day
  222. lc.acflux_landuse_change += stand.get_gridcell_fraction() * indiv_cp.acflux_harvest / (double)stand.nobj;
  223. lc.acflux_harvest_wood_res_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.dcflux_harvest_wood_res / (double)stand.nobj;
  224. lc.acflux_landuse_change_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.acflux_harvest / (double)stand.nobj;
  225. lc.anflux_landuse_change += stand.get_gridcell_fraction() * indiv_cp.anflux_harvest / (double)stand.nobj;
  226. lc.anflux_landuse_change_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.anflux_harvest / (double)stand.nobj;
  227. }
  228. /// Use for normal forest management in calls from growth(). For clearcut during landcover change, use harvest_wood() and kill_remaining_vegetation()
  229. void clearcut(Individual& indiv, double anpp, bool& killed) {
  230. Patch& patch = indiv.vegetation.patch;
  231. Patchpft& ppft = patch.pft[indiv.pft.id];
  232. if (indiv.pft.lifeform == TREE) {
  233. if(indiv.alive)
  234. ppft.litter_sap += anpp;
  235. harvest_wood(indiv, 1.0, indiv.pft.harv_eff, indiv.pft.res_outtake); // frac_cut=1, harv_eff=pft.harv_eff, res_outtake_twig=pft.res_outtake, res_outtake_coarse_root=0
  236. indiv.kill();
  237. indiv.vegetation.killobj();
  238. killed = true;
  239. }
  240. patch.age = 0; //important for results
  241. patch.managed = true;
  242. patch.plant_this_year = true;
  243. }
  244. /// Set forest management intensity for all stands this year
  245. void cut_fractions(Gridcell& gridcell) {
  246. Gridcell::iterator gc_itr = gridcell.begin();
  247. while (gc_itr != gridcell.end()) {
  248. Stand& stand = *gc_itr;
  249. stand.firstobj();
  250. while (stand.isobj) {
  251. Patch& patch = stand.getobj();
  252. patch.man_strength = cut_fraction(patch);
  253. stand.nextobj();
  254. }
  255. ++gc_itr;
  256. }
  257. }
  258. // The following two functions are simplified adaptations (continous cutting) from Swedish forest management code by Fredrik Lagergren and should be
  259. // developed further. Specifically, the productivity values, which should ideally be observed values for each gridcell, are set to a static value.
  260. // Also, the calculated diameter limits and rotation times (which are dependent on productivity) are for Swedish forests.
  261. /// Determines whether this patch should be cut this year.
  262. double cut_fraction(Patch& patch) {
  263. if(!run_landcover)
  264. return 0.0;
  265. Stand& stand = patch.stand;
  266. xtring harvest_system = stlist[stand.stid].get_management(stand.current_rot).harvest_system;
  267. if(harvest_system == "")
  268. return 0.0;
  269. int first_cutyear = nyear_spinup; // Simulation year when forestry harvesting starts; default is directly after spinup.
  270. if(stlist[stand.stid].firstmanageyear < 100000) // Initialised to 100000; other values set in instruction file.
  271. first_cutyear = stlist[stand.stid].firstmanageyear - date.first_calendar_year;
  272. if(date.year < first_cutyear)
  273. return 0.0;
  274. const double minbon = 2.351; // The minimum average "bonitet" for a county in Sweden
  275. const double maxbon = 11.311; // The maximum average "bonitet" for a county in Sweden
  276. const double bonitet = 10.0; // Temporary static value (gives cut_int=17)
  277. double cut_fraction = 0.0;
  278. if(harvest_system == "CLEARCUT") {
  279. // First attempt to calculate optimum rotation age for clearcut
  280. if(patch.cmass_wood() / patch.age > patch.get_cmass_wood_inc_5() && patch.age > 20)
  281. cut_fraction = 1.0;
  282. }
  283. else if(harvest_system == "CONTINUOUS") {
  284. // Continuous forestry
  285. int cut_int; //Interval between cuttings
  286. int patch_order; //Which year in a cutting interval the patch belongs to
  287. // cut_int=30-(int)(15.0*(stand.bonitet-minbon)/(maxbon-minbon));
  288. cut_int=30-(int)(15.0*(bonitet-minbon)/(maxbon-minbon));
  289. patch_order = (int)(patch.id * cut_int * 1.0 / (1.0 * stand.npatch()));
  290. if (!((date.year - first_cutyear - patch_order) % cut_int)) // rule needs to be corrected
  291. cut_fraction = 0.40;
  292. }
  293. return cut_fraction;
  294. }
  295. /// Determines if and how much of this (average) individual is to be cut.
  296. /* Individual is cut by a fraction, determined by cut_fraction(), if diameter is above a calculated limit and
  297. * by 90 % if above a calculated maximum diameter.
  298. * If clearcut is selected (depending on result from cut_fraction()), individual is killed
  299. */
  300. void harvest_forest(Individual& indiv, Pft& pft, bool alive, double anpp, bool& killed) {
  301. Patch& patch = indiv.vegetation.patch;
  302. Patchpft& ppft = patch.pft[indiv.pft.id];
  303. const double minbon=2.351; // The minimum average "bonitet" for a county in Sweden
  304. const double maxbon=11.311; // The maximum average "bonitet" for a county in Sweden
  305. const double bonitet = 10.0; // Temporary static value
  306. int age_class = 0;
  307. double man_strength = patch.man_strength;
  308. if (pft.lifeform==TREE && man_strength > 0.00) {
  309. double diam = pow(indiv.height / indiv.pft.k_allom2, 1.0 / indiv.pft.k_allom3);
  310. if (man_strength == 1.00) {
  311. clearcut(indiv, anpp, killed);
  312. }
  313. else {
  314. double diam_limit=0.13+0.07*(bonitet-minbon)/(maxbon-minbon); // Harvest of trees > 19 cm
  315. double diam_max = diam_limit * 2.0;
  316. if (diam>diam_limit) {
  317. if (diam > diam_max)
  318. man_strength = 0.9;
  319. harvest_wood(indiv, man_strength, indiv.pft.harv_eff, indiv.pft.res_outtake); // frac_cut=man_strength, harv_eff=pft.harv_eff, res_outtake_twig=pft.res_outtake, res_outtake_coarse_root=0
  320. indiv.densindiv *= (1.0 - man_strength);
  321. }
  322. }
  323. // Will tell the program to skip establishment and mortality if management has been performed on this patch,
  324. patch.managed_this_year = true;
  325. patch.managed = true;
  326. }
  327. }
  328. /// Harvest function for pasture, representing grazing (previous year).
  329. /* Function for balancing carbon and nitrogen fluxes from last year's growth
  330. * A fraction of leaves is harvested (pft.harv_eff) and returned as acflux_harvest
  331. * This represents grazing minus return as manure.
  332. * The rest is handled like natural grass in turnover().
  333. * Called from growth() last day of the year for normal harvest/grazing.
  334. * Also called from landcover_dynamics() first day of the year if any natural vegetation
  335. * is transferred to another land use.
  336. * This calls for a scaling factor, when the pasture area has increased.
  337. *
  338. * INPUT/OUTPUT PARAMETERS
  339. * \param Harvest_CN& i struct containing the following indiv-specific public members:
  340. * - cmass_leaf leaf C biomass (kgC/m2)
  341. * - cmass_root fine root C biomass (kgC/m2)
  342. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  343. * - nmass_root fine root nitrogen biomass (kgN/m2)
  344. * OUTPUT PARAMETERS
  345. * \param Harvest_CN& i struct containing the following patchpft-specific public members:
  346. * - litter_leaf new leaf C litter (kgC/m2)
  347. * - litter_root new root C litter (kgC/m2)
  348. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  349. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  350. * ,and the following patch-level public members:
  351. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  352. * - harvested_products_slow harvest products to slow pool (kgC/m2)
  353. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  354. * - harvested_products_slow_nmass harvest nitrogen products to slow pool (kgC/m2)
  355. */
  356. void harvest_pasture(Harvest_CN& i, Pft& pft, bool alive) {
  357. double harvest;
  358. i.dcflux_harvest_pasture = 0.0;
  359. // harvest of leaves (grazing)
  360. // Carbon:
  361. harvest = pft.harv_eff * i.cmass_leaf;
  362. if (ifslowharvestpool) {
  363. i.harvested_products_slow += harvest * pft.harvest_slow_frac;
  364. if (harvest * pft.harvest_slow_frac > 0.0)
  365. i.dcflux_product_pasture = harvest * pft.harvest_slow_frac;
  366. harvest = harvest * (1 - pft.harvest_slow_frac);
  367. }
  368. if (alive) {
  369. i.acflux_harvest += harvest;
  370. i.dcflux_harvest_pasture += harvest;
  371. }
  372. i.cmass_leaf -= harvest;
  373. // Nitrogen
  374. // Reduced removal of N relative to C during grazing.
  375. double N_harvest_scale = 0.25; // Value that works. Needs to be verified in literature.
  376. harvest = pft.harv_eff * i.nmass_leaf * N_harvest_scale;
  377. if (ifslowharvestpool) {
  378. i.harvested_products_slow_nmass += harvest * pft.harvest_slow_frac;
  379. if (harvest * pft.harvest_slow_frac > 0.0)
  380. i.dnflux_product += harvest * pft.harvest_slow_frac;
  381. harvest = harvest * (1 - pft.harvest_slow_frac);
  382. }
  383. i.anflux_harvest += harvest;
  384. i.nmass_leaf -= harvest;
  385. if (grassforcrop && alive) {
  386. // Carbon:
  387. double residue_outtake = pft.res_outtake * i.cmass_leaf; // res_outtake currently set to 0.0,
  388. i.acflux_harvest += residue_outtake; // could be used for burning
  389. i.dcflux_harvest_pasture += residue_outtake;
  390. i.cmass_leaf -= residue_outtake;
  391. // Nitrogen:
  392. residue_outtake = pft.res_outtake * i.nmass_leaf;
  393. i.anflux_harvest += residue_outtake;
  394. i.nmass_leaf -= residue_outtake;
  395. }
  396. }
  397. /// Harvest function for pasture, representing grazing (previous year).
  398. /* Function for balancing carbon and nitrogen fluxes from last year's growth
  399. * A fraction of leaves is harvested (pft.harv_eff) and returned as acflux_harvest
  400. * This represents grazing minus return as manure.
  401. * The rest is handled like natural grass in turnover().
  402. * Called from growth() last day of the year for normal harvest/grazing.
  403. * Also called from landcover_dynamics() first day of the year if any natural vegetation
  404. * is transferred to another land use.
  405. * This calls for a scaling factor, when the pasture area has increased.
  406. *
  407. * This function copies variables from an individual and it's associated patchpft and patch to
  408. * a Harvest_CN struct, which is then passed on to the main harvest_crop function.
  409. * After the execution of the main harvest_crop function, the output variables are copied
  410. * back to the individual and patchpft and the patch-level fluxes are updated.
  411. *
  412. * INPUT/OUTPUT PARAMETERS
  413. * \param indiv reference to an Individual containing the following indiv-specific public members:
  414. * - cmass_leaf leaf C biomass (kgC/m2)
  415. * - cmass_root fine root C biomass (kgC/m2)
  416. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  417. * - nmass_root fine root nitrogen biomass (kgN/m2)
  418. * OUTPUT PARAMETERS
  419. * \param indiv reference to an Individual containing the following patchpft-specific public members:
  420. * - litter_leaf new leaf C litter (kgC/m2)
  421. * - litter_root new root C litter (kgC/m2)
  422. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  423. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  424. * ,and the following patch-level public members:
  425. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  426. * - harvested_products_slow harvest products to slow pool (kgC/m2)
  427. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  428. * - harvested_products_slow_nmass harvest nitrogen products to slow pool (kgC/m2)
  429. */
  430. void harvest_pasture(Individual& indiv, Pft& pft, bool alive) {
  431. Harvest_CN indiv_cp;
  432. indiv_cp.copy_from_indiv(indiv);
  433. harvest_pasture(indiv_cp, pft, alive);
  434. indiv_cp.copy_to_indiv(indiv);
  435. }
  436. /// Harvest function for cropland, including true crops, intercrop grass
  437. /** and pasture grass grown in cropland.
  438. * Function for balancing carbon and nitrogen fluxes from this year's harvested carbon and nitrogen.
  439. * A fraction of harvestable organs (grass:leaves) is harvested (pft.harv_eff) and returned as acflux_harvest.
  440. * A fraction of leaves is removed (pft.res_outtake) and returned as acflux_harvest
  441. * The rest, including roots, is returned as litter, leaving NO carbon or nitrogen in living tissue.
  442. * Called from growth() last day of the year for old-style harvest/grazing or, alternatively, from crop_growth_daily() at harvest day
  443. * (hdate) or last intercrop day (eicdate).
  444. * Also called from landcover_dynamics() first day of the year if any natural vegetation
  445. * is transferred to another land use.
  446. * This calls for a scaling factor, when the pasture area has increased.
  447. *
  448. * This function takes a Harvest_CN struct as an input parameter, copied from an individual and it's associated patchpft and patch.
  449. *
  450. * INPUT PARAMETERS
  451. * \param alive whether individual has survived the first year
  452. * \param isintercropgrass whether individual is cover crop grass
  453. *
  454. * INPUT/OUTPUT PARAMETERS
  455. * \param Harvest_CN& i struct containing the following indiv-specific public members:
  456. * - cmass_leaf leaf C biomass (kgC/m2)
  457. * - cmass_root fine root C biomass (kgC/m2)
  458. * - cmass_ho harvestable organ C biomass (kgC/m2)
  459. * - cmass_agpool above-ground pool C biomass (kgC/m2)
  460. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  461. * - nmass_root fine root nitrogen biomass (kgN/m2)
  462. * - param nmass_ho harvestable organ nitrogen biomass (kgC/m2)
  463. * - param nmass_agpool above-ground pool nitrogen biomass (kgC/m2)
  464. * - nstore_labile labile nitrogen storage (kgC/m2)
  465. * - nstore_longterm longterm nitrogen storage (kgC/m2)
  466. * OUTPUT PARAMETERS
  467. * \param Harvest_CN& i struct containing the following patchpft-specific public members:
  468. * - litter_leaf new leaf C litter (kgC/m2)
  469. * - litter_root new root C litter (kgC/m2)
  470. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  471. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  472. * ,and the following patch-level public members:
  473. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  474. * - harvested_products_slow harvest products to slow pool (kgC/m2)
  475. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  476. * - harvested_products_slow_nmass harvest nitrogen products to slow pool (kgC/m2)
  477. */
  478. void harvest_crop(Harvest_CN& i, Pft& pft, bool alive, bool isintercropgrass) {
  479. double residue_outtake, harvest;
  480. i.dcflux_harvest_crop = 0.0;
  481. if (pft.phenology==CROPGREEN) {
  482. // all root carbon and nitrogen goes to litter
  483. if (i.cmass_root > 0.0) {
  484. i.litter_root += i.cmass_root;
  485. i.dcflux_harvest_res += i.cmass_root;
  486. }
  487. i.cmass_root = 0.0;
  488. if (i.nmass_root > 0.0)
  489. i.nmass_litter_root += i.nmass_root;
  490. if (i.nstore_labile > 0.0)
  491. i.nmass_litter_root += i.nstore_labile;
  492. if (i.nstore_longterm > 0.0)
  493. i.nmass_litter_root += i.nstore_longterm;
  494. i.nmass_root = 0.0;
  495. i.nstore_labile = 0.0;
  496. i.nstore_longterm = 0.0;
  497. // harvest of harvestable organs
  498. // Carbon:
  499. if (i.cmass_ho > 0.0) {
  500. // harvested products
  501. harvest = pft.harv_eff * i.cmass_ho;
  502. // not removed harvestable organs are put into litter
  503. if (pft.aboveground_ho)
  504. i.litter_leaf += (i.cmass_ho - harvest);
  505. else
  506. i.litter_root += (i.cmass_ho - harvest);
  507. i.dcflux_harvest_res += (i.cmass_ho - harvest);
  508. // harvested products not consumed (oxidised) this year put into harvested_products_slow
  509. if (ifslowharvestpool) {
  510. i.harvested_products_slow += harvest * pft.harvest_slow_frac;
  511. if (harvest * pft.harvest_slow_frac > 0.0)
  512. i.dcflux_product_crop = harvest * pft.harvest_slow_frac;
  513. harvest = harvest * (1 - pft.harvest_slow_frac);
  514. }
  515. // harvested products consumed (oxidised) this year put into acflux_harvest
  516. i.acflux_harvest += harvest;
  517. i.dcflux_harvest_crop += harvest;
  518. }
  519. i.cmass_ho = 0.0;
  520. // Nitrogen:
  521. if (i.nmass_ho > 0.0) {
  522. // harvested products
  523. harvest = pft.harv_eff * i.nmass_ho;
  524. // not removed harvestable organs are put into litter
  525. if (pft.aboveground_ho)
  526. i.nmass_litter_leaf += (i.nmass_ho - harvest);
  527. else
  528. i.nmass_litter_root += (i.nmass_ho - harvest);
  529. // harvested products not consumed this year put into harvested_products_slow_nmass
  530. if (ifslowharvestpool) {
  531. i.harvested_products_slow_nmass += harvest * pft.harvest_slow_frac;
  532. if (harvest * pft.harvest_slow_frac > 0.0)
  533. i.dnflux_product += harvest * pft.harvest_slow_frac;
  534. harvest = harvest * (1 - pft.harvest_slow_frac);
  535. }
  536. // harvested products consumed this year put into anflux_harvest
  537. i.anflux_harvest += harvest;
  538. }
  539. i.nmass_ho = 0.0;
  540. // residues
  541. // Carbon
  542. if ((i.cmass_leaf + i.cmass_agpool + i.cmass_dead_leaf + i.cmass_stem) > 0.0) {
  543. // removed residues are oxidised
  544. residue_outtake = pft.res_outtake * (i.cmass_leaf + i.cmass_agpool + i.cmass_dead_leaf + i.cmass_stem);
  545. i.acflux_harvest += residue_outtake;
  546. i.dcflux_harvest_crop += residue_outtake;
  547. // not removed residues are put into litter
  548. i.litter_leaf += i.cmass_leaf + i.cmass_agpool + i.cmass_dead_leaf + i.cmass_stem - residue_outtake;
  549. i.dcflux_harvest_res += i.cmass_leaf + i.cmass_agpool + i.cmass_dead_leaf + i.cmass_stem - residue_outtake;
  550. }
  551. i.cmass_leaf = 0.0;
  552. i.cmass_agpool = 0.0;
  553. i.cmass_dead_leaf = 0.0;
  554. i.cmass_stem = 0.0;
  555. // Nitrogen:
  556. if ((i.nmass_leaf + i.nmass_agpool + i.nmass_dead_leaf) > 0.0) {
  557. // removed residues are oxidised
  558. residue_outtake = pft.res_outtake * (i.nmass_leaf + i.nmass_agpool + i.nmass_dead_leaf);
  559. i.nmass_litter_leaf += i.nmass_leaf + i.nmass_agpool + i.nmass_dead_leaf - residue_outtake;
  560. // not removed residues are put into litter
  561. i.anflux_harvest += residue_outtake;
  562. }
  563. i.nmass_leaf = 0.0;
  564. i.nmass_agpool = 0.0;
  565. i.nmass_dead_leaf = 0.0;
  566. }
  567. else if (pft.phenology == ANY) {
  568. // Intercrop grass
  569. if (isintercropgrass) {
  570. // roots
  571. // all root carbon and nitrogen goes to litter
  572. if (i.cmass_root > 0.0) {
  573. i.litter_root += i.cmass_root;
  574. i.dcflux_harvest_res += i.cmass_root;
  575. }
  576. if (i.nmass_root > 0.0)
  577. i.nmass_litter_root += i.nmass_root;
  578. if (i.nstore_labile > 0.0)
  579. i.nmass_litter_root += i.nstore_labile;
  580. if (i.nstore_longterm > 0.0)
  581. i.nmass_litter_root += i.nstore_longterm;
  582. i.cmass_root = 0.0;
  583. i.nmass_root = 0.0;
  584. i.nstore_labile = 0.0;
  585. i.nstore_longterm = 0.0;
  586. // leaves
  587. // Carbon:
  588. if (i.cmass_leaf > 0.0) {
  589. // Harvest/Grazing of leaves:
  590. harvest = pft.harv_eff_ic * i.cmass_leaf; // currently no harvest of intercrtop grass
  591. // not removed grass is put into litter
  592. i.litter_leaf += i.cmass_leaf - harvest;
  593. i.dcflux_harvest_res += i.cmass_leaf - harvest;
  594. if (ifslowharvestpool) {
  595. i.harvested_products_slow += harvest * pft.harvest_slow_frac; // no slow harvest for grass
  596. if (harvest * pft.harvest_slow_frac > 0.0)
  597. i.dcflux_product_any = harvest * pft.harvest_slow_frac;
  598. harvest = harvest * (1 - pft.harvest_slow_frac);
  599. }
  600. i.acflux_harvest += harvest;
  601. i.dcflux_harvest_crop += harvest;
  602. }
  603. i.cmass_leaf = 0.0;
  604. i.cmass_ho = 0.0;
  605. i.cmass_agpool = 0.0;
  606. // Nitrogen:
  607. if (i.nmass_leaf > 0.0) {
  608. // Harvest/Grazing of leaves:
  609. harvest = pft.harv_eff_ic * i.nmass_leaf; // currently no harvest of intercrtop grass
  610. // not removed grass is put into litter
  611. i.nmass_litter_leaf += i.nmass_leaf - harvest;
  612. if (ifslowharvestpool) {
  613. i.harvested_products_slow_nmass += harvest * pft.harvest_slow_frac; // no slow harvest for grass
  614. if (harvest * pft.harvest_slow_frac > 0.0)
  615. i.dnflux_product += harvest * pft.harvest_slow_frac;
  616. harvest = harvest * (1 - pft.harvest_slow_frac);
  617. }
  618. i.anflux_harvest += harvest;
  619. }
  620. i.nmass_leaf = 0.0;
  621. i.nmass_ho = 0.0;
  622. i.nmass_agpool = 0.0;
  623. }
  624. else { // pasture grass
  625. // harvest of leaves (grazing)
  626. // Carbon:
  627. harvest = pft.harv_eff * i.cmass_leaf;
  628. if (ifslowharvestpool) {
  629. i.harvested_products_slow += harvest * pft.harvest_slow_frac;
  630. if (harvest * pft.harvest_slow_frac > 0.0)
  631. i.dcflux_product_grass = harvest * pft.harvest_slow_frac;
  632. harvest = harvest * (1 - pft.harvest_slow_frac);
  633. }
  634. if (alive) {
  635. i.acflux_harvest += harvest;
  636. i.dcflux_harvest_crop += harvest;
  637. }
  638. i.cmass_leaf -= harvest;
  639. i.cmass_ho = 0.0;
  640. i.cmass_agpool = 0.0;
  641. // Nitrogen:
  642. // Reduced removal of N relative to C during grazing.
  643. double N_harvest_scale = 0.25; // Value that works. Needs to be verified in literature.
  644. harvest = pft.harv_eff * i.nmass_leaf * N_harvest_scale;
  645. if (ifslowharvestpool) {
  646. i.harvested_products_slow_nmass += harvest * pft.harvest_slow_frac;
  647. if (harvest * pft.harvest_slow_frac > 0.0)
  648. i.dnflux_product += harvest * pft.harvest_slow_frac;
  649. harvest = harvest * (1 - pft.harvest_slow_frac);
  650. }
  651. i.anflux_harvest += harvest;
  652. i.nmass_leaf -= harvest;
  653. i.nmass_ho=0.0;
  654. i.nmass_agpool=0.0;
  655. }
  656. }
  657. }
  658. /// Harvest function for cropland, including true crops, intercrop grass
  659. /** and pasture grass grown in cropland.
  660. * Function for balancing carbon and nitrogen fluxes from this year's harvested carbon and nitrogen.
  661. * A fraction of harvestable organs (grass:leaves) is harvested (pft.harv_eff) and returned as acflux_harvest.
  662. * A fraction of leaves is removed (pft.res_outtake) and returned as acflux_harvest
  663. * The rest, including roots, is returned as litter, leaving NO carbon or nitrogen in living tissue.
  664. * Called from growth() last day of the year for old-style harvest/grazing or, alternatively, from crop_growth_daily() at harvest day
  665. * (hdate) or last intercrop day (eicdate).
  666. * Also called from landcover_dynamics() first day of the year if any natural vegetation
  667. * is transferred to another land use.
  668. * This calls for a scaling factor, when the pasture area has increased.
  669. *
  670. * This function copies variables from an individual and it's associated patchpft and patch to
  671. * a Harvest_CN struct, which is then passed on to the main harvest_crop() function.
  672. * After the execution of the main harvest_crop function, the output variables are copied
  673. * back to the individual and patchpft and the patch-level fluxes are updated.
  674. *
  675. * INPUT PARAMETERS
  676. * \param alive whether individual has survived the first year
  677. * \param isintercropgrass whether individual is cover crop grass
  678. * \param harvest_grsC whether harvest daily carbon values are harvested
  679. *
  680. * INPUT/OUTPUT PARAMETERS
  681. * \param indiv reference to an Individual containing the following indiv-specific public members:
  682. * - cmass_leaf leaf C biomass (kgC/m2)
  683. * - cmass_root fine root C biomass (kgC/m2)
  684. * - cmass_ho harvestable organ C biomass (kgC/m2)
  685. * - cmass_agpool above-ground pool C biomass (kgC/m2)
  686. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  687. * - nmass_root fine root nitrogen biomass (kgN/m2)
  688. * - param nmass_ho harvestable organ nitrogen biomass (kgC/m2)
  689. * - param nmass_agpool above-ground pool nitrogen biomass (kgC/m2)
  690. * - nstore_labile labile nitrogen storage (kgC/m2)
  691. * - nstore_longterm longterm nitrogen storage (kgC/m2)
  692. * OUTPUT PARAMETERS
  693. * \param indiv reference to an Individual containing the following patchpft-specific public members:
  694. * - litter_leaf new leaf C litter (kgC/m2)
  695. * - litter_root new root C litter (kgC/m2)
  696. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  697. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  698. * ,and the following patch-level public members:
  699. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  700. * - harvested_products_slow harvest products to slow pool (kgC/m2)
  701. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  702. * - harvested_products_slow_nmass harvest nitrogen products to slow pool (kgC/m2)
  703. */
  704. void harvest_crop(Individual& indiv, Pft& pft, bool alive, bool isintercropgrass, bool harvest_grsC) {
  705. Harvest_CN indiv_cp;
  706. indiv_cp.copy_from_indiv(indiv, harvest_grsC);
  707. harvest_crop(indiv_cp, pft, alive, isintercropgrass);
  708. indiv_cp.copy_to_indiv(indiv, harvest_grsC);
  709. }
  710. /// Transfers all carbon and nitrogen from living tissue to litter
  711. /** Mainly used at land cover change when remaining vegetation after harvest (grass) is
  712. * killed by tillage, following an optional burning.
  713. *
  714. * This function takes a Harvest_CN struct as an input parameter, copied from an individual and it's associated patchpft and patch.
  715. *
  716. * INPUT PARAMETERS
  717. * \param alive whether individual has survived the first year
  718. * \param isintercropgrass whether individual is cover crop grass
  719. * \param burn whether above-ground vegetation C & N is sent to the atmosphere
  720. * rather than to litter
  721. * INPUT/OUTPUT PARAMETERS
  722. * \param Harvest_CN& i struct containing the following indiv-specific public members:
  723. * - cmass_leaf leaf C biomass (kgC/m2)
  724. * - cmass_root fine root C biomass (kgC/m2)
  725. * - cmass_ho harvestable organ C biomass (kgC/m2)
  726. * - cmass_agpool above-ground pool C biomass (kgC/m2)
  727. * - cmass_sap sapwood C biomass (kgC/m2)
  728. * - cmass_heart heartwood C biomass (kgC/m2)
  729. * - cmass_debt C "debt" (retrospective storage) (kgC/m2)
  730. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  731. * - nmass_root fine root nitrogen biomass (kgN/m2)
  732. * - nmass_sap sapwood nitrogen biomass (kgC/m2)
  733. * - nmass_heart heartwood nitrogen biomass (kgC/m2)
  734. * - param nmass_ho harvestable organ nitrogen biomass (kgC/m2)
  735. * - param nmass_agpool above-ground pool nitrogen biomass (kgC/m2)
  736. * - nstore_labile labile nitrogen storage (kgC/m2)
  737. * - nstore_longterm longterm nitrogen storage (kgC/m2)
  738. * OUTPUT PARAMETERS
  739. * \param Harvest_CN& i struct containing the following patchpft-specific public members:
  740. * - litter_leaf new leaf C litter (kgC/m2)
  741. * - litter_root new root C litter (kgC/m2)
  742. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  743. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  744. * ,and the following patch-level public members:
  745. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  746. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  747. */
  748. void kill_remaining_vegetation(Harvest_CN& cp, Pft& pft, bool alive, bool istruecrop_or_intercropgrass, bool burn) {
  749. if (alive || istruecrop_or_intercropgrass) {
  750. cp.litter_root += cp.cmass_root;
  751. if (burn) {
  752. cp.acflux_harvest += cp.cmass_leaf;
  753. cp.acflux_harvest += cp.cmass_sap;
  754. cp.acflux_harvest += cp.cmass_heart - cp.cmass_debt;
  755. }
  756. else {
  757. cp.litter_leaf += cp.cmass_leaf;
  758. cp.litter_sap += cp.cmass_sap;
  759. cp.litter_heart += cp.cmass_heart - cp.cmass_debt;
  760. }
  761. }
  762. cp.nmass_litter_root += cp.nmass_root;
  763. cp.nmass_litter_root += cp.nstore_longterm;
  764. cp.nmass_litter_root += cp.nstore_labile;
  765. if (burn) {
  766. cp.anflux_harvest += cp.nmass_leaf;
  767. cp.anflux_harvest += cp.nmass_sap;
  768. cp.anflux_harvest += cp.nmass_heart;
  769. }
  770. else {
  771. cp.nmass_litter_leaf += cp.nmass_leaf;
  772. cp.nmass_litter_sap += cp.nmass_sap;
  773. cp.nmass_litter_heart += cp.nmass_heart;
  774. }
  775. if (pft.landcover == CROPLAND) {
  776. if (pft.aboveground_ho) {
  777. if (burn) {
  778. cp.acflux_harvest += cp.cmass_ho;
  779. cp.anflux_harvest += cp.nmass_ho;
  780. }
  781. else {
  782. cp.litter_leaf += cp.cmass_ho;
  783. cp.nmass_litter_leaf += cp.nmass_ho;
  784. }
  785. }
  786. else {
  787. cp.litter_root += cp.cmass_ho;
  788. cp.nmass_litter_root += cp.nmass_ho;
  789. }
  790. if (burn) {
  791. cp.acflux_harvest += cp.cmass_agpool;
  792. cp.anflux_harvest += cp.nmass_agpool;
  793. }
  794. else {
  795. cp.litter_leaf += cp.cmass_agpool;
  796. cp.nmass_litter_leaf += cp.nmass_agpool;
  797. }
  798. }
  799. cp.cmass_leaf = 0.0;
  800. cp.cmass_root = 0.0;
  801. cp.cmass_sap = 0.0;
  802. cp.cmass_heart = 0.0;
  803. cp.cmass_debt = 0.0;
  804. cp.cmass_ho = 0.0;
  805. cp.cmass_agpool = 0.0;
  806. cp.nmass_leaf = 0.0;
  807. cp.nmass_root = 0.0;
  808. cp.nstore_longterm = 0.0;
  809. cp.nstore_labile = 0.0;
  810. cp.nmass_sap = 0.0;
  811. cp.nmass_heart = 0.0;
  812. cp.nmass_ho = 0.0;
  813. cp.nmass_agpool = 0.0;
  814. }
  815. /// Transfers all carbon and nitrogen from living tissue to litter
  816. /** Mainly used at land cover change when remaining vegetation after harvest (grass) is
  817. * killed by tillage, following an optional burning.
  818. *
  819. * This function copies variables from an individual and it's associated patchpft and patch to
  820. * a Harvest_CN struct, which is then passed on to the main harvest_crop() function.
  821. * After the execution of the main harvest_crop function, the output variables are copied
  822. * back to the individual and patchpft and the patch-level fluxes are updated.
  823. *
  824. * INPUT PARAMETERS
  825. * \param alive whether individual has survived the first year
  826. * \param isintercropgrass whether individual is cover crop grass
  827. * \param burn whether above-ground vegetation C & N is sent to the atmosphere
  828. * rather than to litter
  829. * INPUT/OUTPUT PARAMETERS
  830. * \param Harvest_CN& i struct containing the following indiv-specific public members:
  831. * - cmass_leaf leaf C biomass (kgC/m2)
  832. * - cmass_root fine root C biomass (kgC/m2)
  833. * - cmass_ho harvestable organ C biomass (kgC/m2)
  834. * - cmass_agpool above-ground pool C biomass (kgC/m2)
  835. * - cmass_sap sapwood C biomass (kgC/m2)
  836. * - cmass_heart heartwood C biomass (kgC/m2)
  837. * - cmass_debt C "debt" (retrospective storage) (kgC/m2)
  838. * - nmass_leaf leaf nitrogen biomass (kgN/m2)
  839. * - nmass_root fine root nitrogen biomass (kgN/m2)
  840. * - nmass_sap sapwood nitrogen biomass (kgC/m2)
  841. * - nmass_heart heartwood nitrogen biomass (kgC/m2)
  842. * - param nmass_ho harvestable organ nitrogen biomass (kgC/m2)
  843. * - param nmass_agpool above-ground pool nitrogen biomass (kgC/m2)
  844. * - nstore_labile labile nitrogen storage (kgC/m2)
  845. * - nstore_longterm longterm nitrogen storage (kgC/m2)
  846. * OUTPUT PARAMETERS
  847. * \param Harvest_CN& i struct containing the following patchpft-specific public members:
  848. * - litter_leaf new leaf C litter (kgC/m2)
  849. * - litter_root new root C litter (kgC/m2)
  850. * - nmass_litter_leaf new leaf nitrogen litter (kgN/m2)
  851. * - nmass_litter_root new root nitrogen litter (kgN/m2)
  852. * ,and the following patch-level public members:
  853. * - acflux_harvest harvest flux to atmosphere (kgC/m2)
  854. * - anflux_harvest harvest nitrogen flux out of system (kgC/m2)
  855. */
  856. void kill_remaining_vegetation(Individual& indiv, bool burn, bool lc_change) {
  857. Harvest_CN indiv_cp;
  858. indiv_cp.copy_from_indiv(indiv);
  859. kill_remaining_vegetation(indiv_cp, indiv.pft, indiv.alive, indiv.istruecrop_or_intercropgrass(), burn);
  860. indiv_cp.copy_to_indiv(indiv, false, lc_change);
  861. if (burn && lc_change) {
  862. Stand& stand = indiv.vegetation.patch.stand;
  863. Landcover& lc = stand.get_gridcell().landcover;
  864. lc.dcflux_landuse_change += stand.get_gridcell_fraction() * indiv_cp.acflux_harvest / (double)stand.nobj; // ecev3 - reset to 0 each day
  865. lc.acflux_landuse_change += stand.get_gridcell_fraction() * indiv_cp.acflux_harvest / (double)stand.nobj;
  866. lc.acflux_landuse_change_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.acflux_harvest / (double)stand.nobj;
  867. lc.anflux_landuse_change += stand.get_gridcell_fraction() * indiv_cp.anflux_harvest / (double)stand.nobj;
  868. lc.anflux_landuse_change_lc[stand.origin] += stand.get_gridcell_fraction() * indiv_cp.anflux_harvest / (double)stand.nobj;
  869. }
  870. }
  871. /// Scaling of last year's or harvest day individual carbon and nitrogen member values in stands that have increased their area fraction this year.
  872. /** Called immediately before harvest functions in growth() or allocation_crop_daily().
  873. */
  874. void scale_indiv(Individual& indiv, bool scale_grsC) {
  875. Stand& stand = indiv.vegetation.patch.stand;
  876. Gridcell& gridcell = stand.get_gridcell();
  877. if (stand.scale_LC_change >= 1.0) {
  878. return;
  879. }
  880. // Scale individual's C and N mass in stands that have increased in area
  881. // this year by (old area/new area):
  882. double scale = stand.scale_LC_change;
  883. if (scale_grsC) {
  884. if (indiv.pft.landcover == CROPLAND) {
  885. if (indiv.has_daily_turnover()) {
  886. indiv.cropindiv->grs_cmass_leaf -= indiv.cropindiv->grs_cmass_leaf_luc * (1.0 - scale);
  887. indiv.cropindiv->grs_cmass_root -= indiv.cropindiv->grs_cmass_root_luc * (1.0 - scale);
  888. indiv.cropindiv->grs_cmass_ho -= indiv.cropindiv->grs_cmass_ho_luc * (1.0 - scale);
  889. indiv.cropindiv->grs_cmass_agpool -= indiv.cropindiv->grs_cmass_agpool_luc * (1.0 - scale);
  890. indiv.cropindiv->grs_cmass_dead_leaf -= indiv.cropindiv->grs_cmass_dead_leaf_luc * (1.0 - scale);
  891. indiv.cropindiv->grs_cmass_stem -= indiv.cropindiv->grs_cmass_stem_luc * (1.0 - scale);
  892. indiv.check_C_mass();
  893. }
  894. else {
  895. indiv.cropindiv->grs_cmass_leaf *= scale;
  896. indiv.cropindiv->grs_cmass_root *= scale;
  897. indiv.cropindiv->grs_cmass_ho *= scale;
  898. indiv.cropindiv->grs_cmass_agpool *= scale;
  899. indiv.cropindiv->grs_cmass_plant *= scale; //grs_cmass_plant not used
  900. indiv.cropindiv->grs_cmass_dead_leaf *= scale;
  901. indiv.cropindiv->grs_cmass_stem *= scale;
  902. }
  903. }
  904. }
  905. else {
  906. indiv.cmass_root *= scale;
  907. indiv.cmass_leaf *= scale;
  908. indiv.cmass_heart *= scale;
  909. indiv.cmass_sap *= scale;
  910. indiv.cmass_debt *= scale;
  911. if (indiv.pft.landcover == CROPLAND) {
  912. indiv.cropindiv->cmass_agpool *= scale;
  913. indiv.cropindiv->cmass_ho *= scale;
  914. }
  915. }
  916. // Deduct individual N present day 0 this year in stands that have increased in area this year, scaled by (1 - old area/new area):
  917. indiv.nmass_root = indiv.nmass_root - indiv.nmass_root_luc * (1.0 - scale);
  918. indiv.nmass_leaf = indiv.nmass_leaf - indiv.nmass_leaf_luc * (1.0 - scale);
  919. indiv.nmass_heart = indiv.nmass_heart - indiv.nmass_heart_luc * (1.0 - scale);
  920. indiv.nmass_sap = indiv.nmass_sap - indiv.nmass_sap_luc * (1.0 - scale);
  921. if (indiv.pft.landcover == CROPLAND) {
  922. indiv.cropindiv->nmass_agpool = indiv.cropindiv->nmass_agpool - indiv.cropindiv->nmass_agpool_luc * (1.0 - scale);
  923. indiv.cropindiv->nmass_ho = indiv.cropindiv->nmass_ho - indiv.cropindiv->nmass_ho_luc * (1.0 - scale);
  924. indiv.cropindiv->nmass_dead_leaf =indiv.cropindiv->nmass_dead_leaf - indiv.cropindiv->nmass_dead_leaf_luc * (1.0 - scale);
  925. }
  926. if (indiv.nstore_labile > indiv.nstore_labile_luc * (1.0 - scale))
  927. indiv.nstore_labile -= indiv.nstore_labile_luc * (1.0 - scale);
  928. else
  929. indiv.nstore_longterm -= indiv.nstore_labile_luc * (1.0 - scale);
  930. indiv.nstore_longterm = indiv.nstore_longterm - indiv.nstore_longterm_luc * (1.0 - scale);
  931. indiv.check_N_mass();
  932. }
  933. /// Yearly function for harvest of all land covers that have yearly allocation, turnover and gridcell.expand_to_new_stand[lc] = false.
  934. /** Should only be called from growth().
  935. // Harvest functions are preceded by rescaling of living C.
  936. // Only affects natural stands if gridcell.expand_to_new_stand[NATURAL] is false.
  937. */
  938. bool harvest_year(Individual& indiv) {
  939. Stand& stand = indiv.vegetation.patch.stand;
  940. Landcover& landcover = stand.get_gridcell().landcover;
  941. bool killed = false;
  942. // Reduce individual's C and N mass in stands that have increased in area this year:
  943. if (landcover.updated && !indiv.has_daily_turnover()) {
  944. scale_indiv(indiv, false);
  945. }
  946. if (stand.landcover == CROPLAND) {
  947. if (!indiv.has_daily_turnover())
  948. harvest_crop(indiv, indiv.pft, indiv.alive, indiv.cropindiv->isintercropgrass, false);
  949. }
  950. else if (stand.landcover == PASTURE) {
  951. harvest_pasture(indiv, indiv.pft, indiv.alive);
  952. }
  953. else if(stand.landcover == FOREST || stand.landcover == NATURAL && run_landcover)
  954. harvest_forest(indiv, indiv.pft, indiv.alive, indiv.anpp, killed);
  955. return killed;
  956. }
  957. /// Yield function for true crops and intercrop grass.
  958. void yield_crop(Individual& indiv) {
  959. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  960. if (indiv.pft.phenology == ANY) { // grass intercrop yield
  961. // Yield dry wieght of allocated harvestable organs this year; NB independent from harvest calculation in harvest_crop (different years)
  962. if (cropindiv.ycmass_leaf > 0.0)
  963. cropindiv.yield = cropindiv.ycmass_leaf * indiv.pft.harv_eff_ic * 2.0;
  964. else
  965. cropindiv.yield = 0.0;
  966. // Yield dry wieght of actually harvest products this year; NB as above
  967. if (cropindiv.harv_cmass_leaf > 0.0)
  968. cropindiv.harv_yield = cropindiv.harv_cmass_leaf * indiv.pft.harv_eff_ic * 2.0;
  969. else
  970. cropindiv.harv_yield = 0.0;
  971. }
  972. else if (indiv.pft.phenology == CROPGREEN) { //true crop yield
  973. // Yield dry wieght of allocated harvestable organs this year; NB independent from harvest calculation in harvest_crop (different years)
  974. if (cropindiv.ycmass_ho > 0.0)
  975. cropindiv.yield = cropindiv.ycmass_ho * indiv.pft.harv_eff * 2.0;// Should be /0.446 instead
  976. else
  977. cropindiv.yield = 0.0;
  978. // Yield dry wieght of actually harvest products this year; NB as above
  979. if (cropindiv.harv_cmass_ho > 0.0)
  980. cropindiv.harv_yield=cropindiv.harv_cmass_ho * indiv.pft.harv_eff * 2.0;
  981. else
  982. cropindiv.harv_yield = 0.0;
  983. // Yield dry wieght of actually harvest products this year; NB as above
  984. for (int i=0;i<2;i++) {
  985. if (cropindiv.cmass_ho_harvest[i] > 0.0)
  986. cropindiv.yield_harvest[i] = cropindiv.cmass_ho_harvest[i] * indiv.pft.harv_eff * 2.0;
  987. else
  988. cropindiv.yield_harvest[i]=0.0;
  989. }
  990. }
  991. return;
  992. }
  993. /// Yield function for pasture grass grown in cropland landcover
  994. void yield_pasture(Individual& indiv, double cmass_leaf_inc) {
  995. cropindiv_struct& cropindiv = *(indiv.get_cropindiv());
  996. // Normal CC3G/CC4G stand growth (Pasture)
  997. // OK if turnover_leaf==1.0, else (cmass_leaf+cmass_leaf_inc)*indiv.pft.harv_eff*2.0
  998. cropindiv.yield = max(0.0, cmass_leaf_inc) * indiv.pft.harv_eff * 2.0;
  999. // Although no specified harvest date, harv_yield is set for compatibility.
  1000. cropindiv.harv_yield = cropindiv.yield;
  1001. }
  1002. /// Function that determines amount of nitrogen applied today. Crop-specific, pft-based.
  1003. void nfert_crop(Patch& patch) {
  1004. Gridcell& gridcell = patch.stand.get_gridcell();
  1005. patch.dnfert = 0.0;
  1006. pftlist.firstobj();
  1007. // Loop through PFTs
  1008. while(pftlist.isobj) {
  1009. Pft& pft = pftlist.getobj();
  1010. Patchpft& patchpft = patch.pft[pft.id];
  1011. Gridcellpft& gridcellpft = gridcell.pft[pft.id];
  1012. if (patch.stand.pft[pft.id].active && pft.phenology == CROPGREEN) {
  1013. cropphen_struct& ppftcrop = *(patchpft.get_cropphen());
  1014. if(!ppftcrop.growingseason) {
  1015. pftlist.nextobj();
  1016. continue;
  1017. }
  1018. double nfert = pft.N_appfert;
  1019. if (gridcellpft.Nfert_read >= 0.0) {
  1020. nfert = gridcellpft.Nfert_read;
  1021. }
  1022. if (!ppftcrop.fertilised[0] && ppftcrop.dev_stage > 0.0){
  1023. // Fertiliser application at dev_stage = 0, sowing.
  1024. patch.dnfert = nfert * (1.0 - pft.fertrate[0] - pft.fertrate[1]);
  1025. ppftcrop.fertilised[0] = true;
  1026. }
  1027. else if (!ppftcrop.fertilised[1] && ppftcrop.dev_stage > pft.fert_stages[0]){
  1028. patch.dnfert = nfert * pft.fertrate[0];
  1029. ppftcrop.fertilised[1] = true;
  1030. }
  1031. else if (!ppftcrop.fertilised[2] && ppftcrop.dev_stage > pft.fert_stages[1]){
  1032. patch.dnfert = nfert * (pft.fertrate[1]);
  1033. ppftcrop.fertilised[2] = true;
  1034. }
  1035. }
  1036. pftlist.nextobj();
  1037. }
  1038. patch.anfert += patch.dnfert;
  1039. patch.fluxes.report_flux(Fluxes::NFERT, patch.dnfert);
  1040. }
  1041. /// Function that determines amount of nitrogen applied today.
  1042. void nfert(Patch& patch) {
  1043. Stand& stand = patch.stand;
  1044. StandType& st = stlist[stand.stid];
  1045. Gridcell& gridcell = stand.get_gridcell();
  1046. if(stand.landcover == CROPLAND) {
  1047. nfert_crop(patch);
  1048. return;
  1049. }
  1050. // General code for applying nitrogen to other land cover types, an equal amount every day.
  1051. double nfert;
  1052. if(gridcell.st[st.id].nfert >= 0.0) { // todo: management type variable (mt.nfert)
  1053. nfert = gridcell.st[st.id].nfert;
  1054. }
  1055. else {
  1056. nfert = 0.0;
  1057. }
  1058. patch.dnfert = nfert / (double)date.year_length();
  1059. patch.anfert += patch.dnfert;
  1060. patch.fluxes.report_flux(Fluxes::NFERT, patch.dnfert);
  1061. }
  1062. /// Updates crop rotation status
  1063. /** Sets new crop management variables, typically on harvest day
  1064. */
  1065. void crop_rotation(Stand& stand) {
  1066. if (stand.landcover != CROPLAND) {
  1067. return;
  1068. }
  1069. CropRotation& rotation = stlist[stand.stid].rotation;
  1070. stand.ndays_inrotation++;
  1071. if (rotation.ncrops < 2 || !stand.isrotationday) {
  1072. return;
  1073. }
  1074. int firstrotyear = rotation.firstrotyear - date.first_calendar_year;
  1075. bool postpone_rotation = false;
  1076. // Alternative uses of firstrotyear:
  1077. // 1. Before firstrotyear, grow only crop1:
  1078. // if(date.year < firstrotyear)
  1079. // postpone_rotation = true;
  1080. // 2. Synchronise rotation with firstrotyear:
  1081. // A. At the creation of the stand:
  1082. if (date.year < stand.first_year + 3)
  1083. // B. At firstrotyear
  1084. // if(date.year == firstrotyear - 1)
  1085. // C. Continuously:
  1086. {
  1087. if ((abs(firstrotyear - date.year) % rotation.ncrops) != stand.current_rot)
  1088. postpone_rotation = true;
  1089. }
  1090. if (!postpone_rotation) {
  1091. if(stand.infallow) {
  1092. stand.infallow = false;
  1093. stand.get_gridcell().pft[stand.pftid].sowing_restriction = false;
  1094. }
  1095. int old_pftid = stand.pftid;
  1096. stand.rotate();
  1097. for (unsigned int p=0; p<stand.nobj; p++) {
  1098. cropphen_struct& previous = *(stand[p].pft[old_pftid].get_cropphen());
  1099. cropphen_struct& current = *(stand[p].pft[stand.pftid].get_cropphen());
  1100. previous.bicdate = -1;
  1101. if(!previous.intercropseason)
  1102. current.bicdate = stepfromdate(date.day, 15);
  1103. previous.eicdate = -1;
  1104. current.eicdate = -1;
  1105. previous.hdate = -1;
  1106. current.intercropseason = previous.intercropseason;
  1107. }
  1108. // Adds sowing and harvest dates for the second crop in a double cropping system
  1109. if (stlist[stand.stid].get_management(stand.current_rot).multicrop && rotation.ncrops == 2 && stand.current_rot == 1) {
  1110. if (stand.pft[stand.pftid].sdate_force < 0)
  1111. stand.pft[stand.pftid].sdate_force = stepfromdate(date.day, 10);
  1112. if (stand.pft[stand.pftid].hdate_force < 0) {
  1113. stand.pft[stand.pftid].hdate_force = stepfromdate(stand.pft[old_pftid].sdate_force, -10);
  1114. }
  1115. }
  1116. if(stlist[stand.stid].get_management(stand.current_rot).fallow) {
  1117. stand.infallow = true;
  1118. stand.get_gridcell().pft[stand.pftid].sowing_restriction = true;
  1119. }
  1120. }
  1121. stand.isrotationday = false;
  1122. }