landcover.cpp 103 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887
  1. ///////////////////////////////////////////////////////////////////////////////////////
  2. /// \file landcover.cpp
  3. /// \brief Functions handling landcover aspects, such as creating or resizing Stands
  4. ///
  5. /// Landcover change.
  6. ///
  7. /// \author Mats Lindeskog,
  8. /// \Part of code in this file as well as in cropphenology.cpp, cropallocation.cpp and
  9. /// \management.cpp based on LPJ-mL C++ code received from Alberte Bondeau in 2008.
  10. /// $Date$
  11. ///
  12. ///////////////////////////////////////////////////////////////////////////////////////
  13. #include "config.h"
  14. #include "landcover.h"
  15. #include "guessmath.h"
  16. #define CHECK_NLIM -13
  17. //#define CHECK_NLIM -8
  18. //#define CHECK_NLIM -5
  19. // prevent a check error to terminate the run, this should only be used in a test run.
  20. #define FAIL dprintf
  21. //#define FAIL fail
  22. /// Query whether a date is within a period spanned by two dates.
  23. bool dayinperiod(int day, int start, int end) {
  24. bool acrossnewyear = false;
  25. if(day < 0 || start < 0 || end < 0) // a negative value should not be a valid day
  26. return false;
  27. if(start > end)
  28. acrossnewyear = true;
  29. if(day >= start && day <= end && !acrossnewyear || (day >= start || day <= end) && acrossnewyear)
  30. return true;
  31. else
  32. return false;
  33. }
  34. /// Step n days from a date.
  35. /** Local proxy to the Date class function in order to simplify call syntax
  36. from the modules that #includes landcover.h anyway.
  37. */
  38. int stepfromdate(int day, int step) {
  39. return Date::stepfromdate(day, step);
  40. }
  41. /// Help function to access two-dimentional arrays that have been created dynamically
  42. int index(int from, int to, int ncols = nst) {
  43. return from * ncols + to;
  44. }
  45. /// Initial creation of stands when run_landcover==true based on gridcell stand type area fractions.
  46. void landcover_init(Gridcell& gridcell, InputModule* input_module) {
  47. // Set CFT-specific members of gridcellpft:
  48. for(unsigned int p = 0; p < gridcell.pft.nobj; p++) {
  49. Gridcellpft& gcpft = gridcell.pft[p];
  50. if (gridcell.get_lat() >= 0.0) {
  51. gcpft.sdate_default = gcpft.pft.sdatenh;
  52. gcpft.hlimitdate_default = gcpft.pft.hlimitdatenh;
  53. }
  54. else {
  55. gcpft.sdate_default = gcpft.pft.sdatesh;
  56. gcpft.hlimitdate_default = gcpft.pft.hlimitdatesh;
  57. }
  58. }
  59. // get landcover and crop area fractions from landcover input file(s) or ins-file.
  60. input_module->getlandcover(gridcell);
  61. stlist.firstobj();
  62. while (stlist.isobj) {
  63. StandType& st = stlist.getobj();
  64. Gridcellst& gcst = gridcell.st[st.id];
  65. gcst.frac_old = gcst.frac;
  66. if(gcst.frac > 0.0) {
  67. gridcell.create_stand_lu(st, gcst.frac);
  68. }
  69. stlist.nextobj();
  70. }
  71. }
  72. enum {SECONDARY_MATURE, SECONDARY_YOUNG, PRIMARY, NFORESTCLASSES};
  73. /// identifies which stands to reduce in area and sets standtype.nstands
  74. /** Updates frac, frac_change, frac_old and gross_frac_decrease for reduced stands
  75. * Updates frac_old and frac_temp for all stands
  76. *
  77. * INTPUT PARAMETERS
  78. *
  79. * \param st_frac_transfer array with this year's transferred area fractions between the different stand types
  80. */
  81. void reduce_stands(Gridcell& gridcell, double* st_frac_transfer, forest_st_frac_transfer& forest_st_frac_transfer_s) {
  82. for(unsigned int i = 0; i < gridcell.st.nobj; i++) {
  83. Gridcellst& gcst = gridcell.st[i];
  84. gcst.nstands = 0;
  85. }
  86. for(unsigned int i = 0; i < gridcell.nbr_stands(); i++) {
  87. Stand& stand = gridcell[i];
  88. stand.frac_old = stand.get_gridcell_fraction();
  89. stand.frac_temp = stand.frac_old;
  90. stand.frac_change = 0.0;
  91. stand.gross_frac_increase = 0.0;
  92. stand.gross_frac_decrease = 0.0;
  93. stand.cloned_fraction = 0.0;
  94. if(stand.transfer_area_st) {
  95. for(int i=0;i<nst;i++)
  96. stand.transfer_area_st[i] = 0.0;
  97. }
  98. gridcell.st[stand.stid].nstands++;
  99. }
  100. for(int from=0; from<nst; from++) {
  101. StandType& st = stlist[from];
  102. Gridcellst& gcst = gridcell.st[from];
  103. if(gcst.gross_frac_decrease > 0.0) {
  104. if(gcst.nstands > 1) {
  105. /* Rules for selecting stands:
  106. natural to cropland: young_stands_first=true
  107. ifprimary_lc_transfer:
  108. secondary mature lap: young_stands_first=true, secondary_harvest=true (will use young -> old secondary stands over age limit)
  109. restart: young_stands_first=false, secondary_harvest=false, ignore age_limit_reduce (will use primary stand before age-limited stands)
  110. (secondary young lap): young_stands_first=true, secondary_harvest=true (will use young -> old secondary stands over age limit)
  111. restart: young_stands_first=true, secondary_harvest=false, ignore age_limit_reduce (will use primary stand before age-limited stands)
  112. primary lap: young_stands_first=false, secondary_harvest=false (will use old -> young secondary stands if no more primary stand)
  113. restart: young_stands_first=false, secondary_harvest=false, ignore age_limit_reduce (will use age-limited stands)
  114. natural to natural: young_stands_first=false
  115. ifprimary_lc_transfer:
  116. secondary mature lap: young_stands_first=false, secondary_harvest=true (will use old -> young secondary stands over age limit)
  117. restart: young_stands_first=false, secondary_harvest=false, ignore age_limit_reduce (will use primary stand before age-limited stands)
  118. secondary young lap: young_stands_first=true, secondary_harvest=true (will use young -> old secondary stands over age limit)
  119. restart: young_stands_first=true, secondary_harvest=false, ignore age_limit_reduce (will use age-limited stands before primary stand)
  120. primary lap: young_stands_first=false, secondary_harvest=false (will use old -> young secondary stands if no more primary stand)
  121. restart: young_stands_first=false, secondary_harvest=false, ignore age_limit_reduce (will use age-limited stands)
  122. */
  123. int nlaps = 1;
  124. if(use_primary_lc_transfer)
  125. nlaps = NFORESTCLASSES;
  126. for(int n=0; n<nlaps; n++) {
  127. for(int to=0; to<nst; to++) {
  128. if(st_frac_transfer[index(from, to)] > 0.0) {
  129. double st_change_remain = -st_frac_transfer[index(from, to)];
  130. bool secondary_harvest = false;
  131. bool young_stands_first = true; //convert area from youngest stands first
  132. if(st.landcover == NATURAL || st.landcover == FOREST) {
  133. if(stlist[to].landcover == NATURAL || stlist[to].landcover == FOREST)
  134. young_stands_first = false;
  135. else
  136. young_stands_first = true;
  137. }
  138. if(use_primary_lc_transfer) {
  139. if(n == SECONDARY_MATURE) {
  140. // First reduce secondary (mature) stands (fraction not in primary_st_frac_transfer array + secondary_young_st_frac_transfer[index(from, to))
  141. st_change_remain += forest_st_frac_transfer_s.primary[index(from, to)] + forest_st_frac_transfer_s.secondary_young[index(from, to)];
  142. secondary_harvest = true;
  143. } else if(n == SECONDARY_YOUNG) {
  144. // Then reduce young secondary stands
  145. st_change_remain = -forest_st_frac_transfer_s.secondary_young[index(from, to)];
  146. young_stands_first = true;
  147. secondary_harvest = true;
  148. }
  149. else if(n == PRIMARY) {
  150. // Finally, reduce primary stands (fraction in primary_st_frac_transfer array)
  151. st_change_remain = -forest_st_frac_transfer_s.primary[index(from, to)];
  152. young_stands_first = false;
  153. }
  154. }
  155. bool restart = false;
  156. while(st_change_remain < 0.0) {
  157. int count_st = 0;
  158. for(unsigned int i = 0; i < gridcell.size(); i++) {
  159. int index;
  160. double stands_frac_sum = 0.0;
  161. for(unsigned int j = 0; j < gridcell.nbr_stands(); j++) {
  162. Stand& stand = gridcell[j];
  163. if(stand.stid == st.id)
  164. stands_frac_sum += stand.get_gridcell_fraction();
  165. }
  166. if(st_change_remain > -INPUT_RESOLUTION * 0.1 || stands_frac_sum == 0.0) {
  167. if(st_change_remain <= -INPUT_RESOLUTION && stands_frac_sum == 0.0)
  168. dprintf("\nWarning: no more stand area left of stand type %d ! Residual reduction demand %.15f ignored.\n", st.id, -st_change_remain);
  169. st_change_remain = 0.0;
  170. break;
  171. }
  172. if(young_stands_first)
  173. index = (int) gridcell.size() - 1 - i;
  174. else
  175. index = i;
  176. Stand& stand = gridcell[index];
  177. if(stand.stid == st.id) {
  178. count_st++;
  179. int first_year = max(stand.first_year, stand.clone_year);
  180. if(secondary_harvest && first_year == 0)
  181. continue;
  182. // Don't reduce stands younger than the age limit, unless this is the last stand in the loop or the initial stand has been killed
  183. if(date.year - first_year < age_limit_reduce && !reduce_all_stands && !restart)
  184. continue;
  185. // convert equal percentage of areas from all stands
  186. if(reduce_all_stands) {
  187. stand.frac_change += st_change_remain * stand.get_gridcell_fraction() / gcst.frac;
  188. stand.gross_frac_decrease -= st_change_remain * stand.get_gridcell_fraction() / gcst.frac;
  189. stand.transfer_area_st[to] = -st_change_remain * stand.get_gridcell_fraction() / gcst.frac;
  190. stand.set_gridcell_fraction(stand.get_gridcell_fraction() + st_change_remain * stand.get_gridcell_fraction() / gcst.frac);
  191. }
  192. else {
  193. if(stand.get_gridcell_fraction() > 0.0) {
  194. //all natural landcover decrease is taken from this stand
  195. if(stand.get_gridcell_fraction() >= -st_change_remain) {
  196. stand.frac_change += st_change_remain;
  197. stand.gross_frac_decrease -= st_change_remain;
  198. stand.transfer_area_st[to] -= st_change_remain;
  199. stand.set_gridcell_fraction(stand.get_gridcell_fraction() + st_change_remain);
  200. if(stand.get_gridcell_fraction() < INPUT_RESOLUTION * 0.1) {
  201. gridcell.landcover.acflux_landuse_change += stand.ccont() * stand.get_gridcell_fraction();
  202. gridcell.landcover.anflux_landuse_change += stand.ncont() * stand.get_gridcell_fraction();
  203. stand.set_gridcell_fraction(0.0);
  204. }
  205. st_change_remain = 0.0;
  206. break;
  207. }
  208. //more stands will have to be reduced
  209. else {
  210. stand.frac_change -= stand.get_gridcell_fraction();
  211. stand.gross_frac_decrease += stand.get_gridcell_fraction();
  212. stand.transfer_area_st[to] += stand.get_gridcell_fraction();
  213. st_change_remain += stand.get_gridcell_fraction();
  214. stand.set_gridcell_fraction(0.0); //will be killed below
  215. }
  216. }
  217. }
  218. }
  219. }
  220. // Restart loop and reduce oldest stand first if the rules above precluded reduction of the required area.
  221. if(st_change_remain < 0.0) {
  222. if(reduce_all_stands) {
  223. st_change_remain = 0.0;
  224. }
  225. else {
  226. if(n != SECONDARY_YOUNG)
  227. young_stands_first = false;
  228. secondary_harvest = false;
  229. restart = true;
  230. }
  231. }
  232. }
  233. }
  234. }
  235. }
  236. }
  237. else if(gcst.nstands == 1 ) {
  238. for(unsigned int i = 0; i < gridcell.size(); i++) {
  239. Stand& stand = gridcell[i];
  240. if(stand.stid == gcst.id) {
  241. bool expand_to_new_stand = gridcell.landcover.expand_to_new_stand[stand.landcover];
  242. if(gcst.gross_frac_decrease > stand.get_gridcell_fraction()) {
  243. if((gcst.gross_frac_decrease - stand.get_gridcell_fraction() > INPUT_RESOLUTION))
  244. dprintf("\nYear %d: Warning. no more stand area left of stand type %d ! Residual reduction demand %.15f ignored.\n", date.year, stand.stid, gcst.gross_frac_decrease - stand.get_gridcell_fraction());
  245. gcst.gross_frac_decrease = stand.get_gridcell_fraction();
  246. }
  247. stand.frac_change = -gcst.gross_frac_decrease;
  248. stand.gross_frac_decrease = gcst.gross_frac_decrease;
  249. stand.set_gridcell_fraction(stand.get_gridcell_fraction() + stand.frac_change);
  250. if(stand.get_gridcell_fraction() < INPUT_RESOLUTION && (!gcst.gross_frac_increase || expand_to_new_stand)) {
  251. gridcell.landcover.acflux_landuse_change += stand.ccont() * stand.get_gridcell_fraction();
  252. gridcell.landcover.anflux_landuse_change += stand.ncont() * stand.get_gridcell_fraction();
  253. stand.set_gridcell_fraction(0.0);
  254. }
  255. for(int to=0; to<nst; to++) {
  256. if(st_frac_transfer[index(st.id, to)])
  257. stand.transfer_area_st[to] = st_frac_transfer[index(st.id, to)];
  258. }
  259. }
  260. }
  261. }
  262. }
  263. }
  264. for(unsigned int j = 0; j < gridcell.nbr_stands(); j++) {
  265. Stand& stand = gridcell[j];
  266. StandType& st = stlist[stand.stid];
  267. Gridcellst& gcst = gridcell.st[stand.stid];
  268. if(!gcst.frac && stand.get_gridcell_fraction() < INPUT_RESOLUTION * 100.0) {
  269. if(stand.get_gridcell_fraction() > INPUT_RESOLUTION * 10.0)
  270. dprintf("\nYear %d: remaining stand when stand type %d fraction is 0. Residual fraction %.15f ignored. Stand killed.\n", date.year, st.id, stand.get_gridcell_fraction());
  271. gridcell.landcover.acflux_landuse_change += stand.ccont() * stand.get_gridcell_fraction();
  272. gridcell.landcover.anflux_landuse_change += stand.ncont() * stand.get_gridcell_fraction();
  273. stand.set_gridcell_fraction(0.0);
  274. }
  275. }
  276. }
  277. /// identifies which stands to expand in area
  278. /** Updates frac, frac_change, frac_old and gross_frac_increase for expanded stands
  279. * Should be preceded by a call to reduce_stands() and, optionally, transfer_to_new_stand()
  280. *
  281. * INTPUT PARAMETERS
  282. *
  283. * \param st_frac_transfer array with this year's transferred area fractions between the different stand types
  284. */
  285. void expand_stands(Gridcell& gridcell, double* st_frac_transfer) {
  286. // Updated variables are zeroed in reduce_stands()
  287. for(int i=0; i<nst; i++) {
  288. StandType& st = stlist[i];
  289. Gridcellst& gcst = gridcell.st[i];
  290. landcovertype lc = st.landcover;
  291. bool expand_to_new_stand = gridcell.landcover.expand_to_new_stand[lc];
  292. if(gcst.gross_frac_increase > 0.0) { // Not cloned stands
  293. if(!expand_to_new_stand) {
  294. for(unsigned int i = 0; i < gridcell.size(); i++) {
  295. Stand& stand = gridcell[i];
  296. if(stand.stid == st.id) {
  297. // Identify which stands to expand (not secondary stands)
  298. if(!stand.cloned) {
  299. stand.frac_change += gcst.gross_frac_increase;
  300. stand.gross_frac_increase = gcst.gross_frac_increase;
  301. stand.set_gridcell_fraction(stand.get_gridcell_fraction() + gcst.gross_frac_increase);
  302. }
  303. }
  304. }
  305. }
  306. }
  307. }
  308. }
  309. /// sets land cover transfer matrix from gross land cover change data when no input is available
  310. /** Uses rules to select preferred transfers between land covers
  311. * NB. New land cover types must be included in the preference arrays !
  312. * Also, PEATLAND, URBAN and BARREN needs to be included when using dynamic fractions for these land cover types.
  313. * INPUT PARAMETERS
  314. * \param landcoverfrac_change array with this year's difference in area fractions of the different landcovers
  315. *
  316. * OUTPUT PARAMETERS
  317. *
  318. * \param lc_frac_transfer array with this year's transferred area fractions between the different land covers
  319. */
  320. void set_lc_change_array(double landcoverfrac_change[], double lc_frac_transfer[][NLANDCOVERTYPES]) {
  321. const int NRANK = 6;
  322. int target_preference[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0};
  323. int origin_preference[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0};
  324. double receptor_remain[NLANDCOVERTYPES] = {0.0};
  325. double donor_remain[NLANDCOVERTYPES] = {0.0};
  326. int ndonor_lc = 0;
  327. int nreceptor_lc = 0;
  328. target_preference[CROPLAND][PASTURE] = 5;
  329. target_preference[CROPLAND][NATURAL] = 6;
  330. target_preference[CROPLAND][FOREST] = 4;
  331. target_preference[CROPLAND][URBAN] = 3;
  332. target_preference[CROPLAND][BARREN] = 2;
  333. target_preference[CROPLAND][PEATLAND] = 1;
  334. target_preference[PASTURE][CROPLAND] = 4;
  335. target_preference[PASTURE][NATURAL] = 6;
  336. target_preference[PASTURE][FOREST] = 5;
  337. target_preference[PASTURE][URBAN] = 3;
  338. target_preference[PASTURE][BARREN] = 1;
  339. target_preference[PASTURE][PEATLAND] = 2;
  340. target_preference[FOREST][CROPLAND] = 4;
  341. target_preference[FOREST][PASTURE] = 5;
  342. target_preference[FOREST][NATURAL] = 6;
  343. target_preference[FOREST][URBAN] = 2;
  344. target_preference[FOREST][BARREN] = 1;
  345. target_preference[FOREST][PEATLAND] = 3;
  346. target_preference[NATURAL][CROPLAND] = 4;
  347. target_preference[NATURAL][PASTURE] = 5;
  348. target_preference[NATURAL][FOREST] = 6;
  349. target_preference[NATURAL][URBAN] = 2;
  350. target_preference[NATURAL][BARREN] = 1;
  351. target_preference[NATURAL][PEATLAND] = 3;
  352. target_preference[URBAN][CROPLAND] = 5;
  353. target_preference[URBAN][PASTURE] = 4;
  354. target_preference[URBAN][FOREST] = 2;
  355. target_preference[URBAN][NATURAL] = 3;
  356. target_preference[URBAN][BARREN] = 6;
  357. target_preference[URBAN][PEATLAND] = 1;
  358. target_preference[BARREN][CROPLAND] = 4;
  359. target_preference[BARREN][PASTURE] = 5;
  360. target_preference[BARREN][FOREST] = 2;
  361. target_preference[BARREN][NATURAL] = 6;
  362. target_preference[BARREN][URBAN] = 1;
  363. target_preference[BARREN][PEATLAND] = 3;
  364. origin_preference[PASTURE][CROPLAND] = 5;
  365. origin_preference[NATURAL][CROPLAND] = 6;
  366. origin_preference[FOREST][CROPLAND] = 3;
  367. origin_preference[BARREN][CROPLAND] = 1;
  368. origin_preference[URBAN][CROPLAND] = 2;
  369. origin_preference[PEATLAND][CROPLAND] = 4;
  370. origin_preference[CROPLAND][PASTURE] = 5;
  371. origin_preference[NATURAL][PASTURE] = 6;
  372. origin_preference[FOREST][PASTURE] = 3;
  373. origin_preference[BARREN][PASTURE] = 1;
  374. origin_preference[URBAN][PASTURE] = 2;
  375. origin_preference[PEATLAND][PASTURE] = 4;
  376. origin_preference[CROPLAND][NATURAL] = 5;
  377. origin_preference[PASTURE][NATURAL] = 6;
  378. origin_preference[FOREST][NATURAL] = 4;
  379. origin_preference[BARREN][NATURAL] = 2;
  380. origin_preference[URBAN][NATURAL] = 1;
  381. origin_preference[PEATLAND][NATURAL] = 3;
  382. origin_preference[CROPLAND][FOREST] = 4;
  383. origin_preference[PASTURE][FOREST] = 5;
  384. origin_preference[NATURAL][FOREST] = 6;
  385. origin_preference[BARREN][FOREST] = 2;
  386. origin_preference[URBAN][FOREST] = 1;
  387. origin_preference[PEATLAND][FOREST] = 3;
  388. origin_preference[CROPLAND][URBAN] = 6;
  389. origin_preference[PASTURE][URBAN] = 5;
  390. origin_preference[NATURAL][URBAN] = 4;
  391. origin_preference[BARREN][URBAN] = 3;
  392. origin_preference[FOREST][URBAN] = 1;
  393. origin_preference[PEATLAND][URBAN] = 2;
  394. origin_preference[CROPLAND][PEATLAND] = 2;
  395. origin_preference[PASTURE][PEATLAND] = 4;
  396. origin_preference[NATURAL][PEATLAND] = 5;
  397. origin_preference[BARREN][PEATLAND] = 6;
  398. origin_preference[FOREST][PEATLAND] = 3;
  399. origin_preference[URBAN][PEATLAND] = 1;
  400. origin_preference[CROPLAND][BARREN] = 6;
  401. origin_preference[PASTURE][BARREN] = 5;
  402. origin_preference[NATURAL][BARREN] = 4;
  403. origin_preference[URBAN][BARREN] = 3;
  404. origin_preference[FOREST][BARREN] = 1;
  405. origin_preference[PEATLAND][BARREN] = 2;
  406. for(int i=0; i<NLANDCOVERTYPES; i++) {
  407. if(landcoverfrac_change[i] > 0.0) {
  408. receptor_remain[i] = landcoverfrac_change[i];
  409. nreceptor_lc++;
  410. }
  411. else if(landcoverfrac_change[i] < 0.0){
  412. donor_remain[i] = -landcoverfrac_change[i];
  413. ndonor_lc++;
  414. }
  415. }
  416. // Simplest cases: no ambiguities
  417. if(ndonor_lc == 1 || nreceptor_lc == 1) {
  418. for(int from=0; from<NLANDCOVERTYPES; from++) {
  419. if(landcoverfrac_change[from] < 0.0) {
  420. for(int to=0; to<NLANDCOVERTYPES; to++) {
  421. if(ndonor_lc == 1) {
  422. if(landcoverfrac_change[to] > 0.0)
  423. lc_frac_transfer[from][to] += landcoverfrac_change[to];
  424. }
  425. else if(nreceptor_lc == 1) {
  426. if(landcoverfrac_change[to] > 0.0)
  427. lc_frac_transfer[from][to] += - landcoverfrac_change[from];
  428. }
  429. }
  430. }
  431. }
  432. }
  433. else {
  434. for(int score=NRANK*2; score>0; score--) {
  435. for(int from=0; from<NLANDCOVERTYPES; from++) {
  436. if(donor_remain[from] > INPUT_RESOLUTION * 0.1) {
  437. for(int to=0; to<NLANDCOVERTYPES; to++) {
  438. // Identify receiving land covers:
  439. if(target_preference[from][to] + origin_preference[from][to] == score && receptor_remain[to] > INPUT_RESOLUTION * 0.1 && donor_remain[from] > INPUT_RESOLUTION * 0.1) {
  440. // all donor land is put into this land cover
  441. if(receptor_remain[to] >= donor_remain[from]) {
  442. lc_frac_transfer[from][to] += donor_remain[from];
  443. receptor_remain[to] -= donor_remain[from];
  444. donor_remain[from] = 0.0;
  445. break;
  446. }
  447. // transfer to more land cover types
  448. else {
  449. lc_frac_transfer[from][to] += receptor_remain[to];
  450. donor_remain[from] -= receptor_remain[to];
  451. receptor_remain[to] = 0.0;
  452. }
  453. }
  454. }
  455. }
  456. }
  457. }
  458. }
  459. }
  460. /// sets stand type transfer matrix from land cover transfer matrix when no stand type transfer input is available
  461. /** Distributes land cover transfers equally between stand types within a land cover but may use
  462. * rules to select preferred transfers between stand types/land covers (see set_lc_change_array()) if
  463. * equal_distribution is set to false.
  464. *
  465. * INPUT PARAMETERS
  466. * \param lc_frac_transfer array with this year's transferred area fractions between the different land covers
  467. *
  468. * OUTPUT PARAMETERS
  469. *
  470. * \param st_frac_transfer array with this year's transferred area fractions between the different stand types
  471. */
  472. void set_st_change_array(Gridcell& gridcell, double lc_frac_transfer[][NLANDCOVERTYPES], double* st_frac_transfer, forest_lc_frac_transfer& forest_lc_frac_transfer_s, forest_st_frac_transfer& forest_st_frac_transfer_s) {
  473. double* recip_donor_remain = new double[nst];
  474. double* recip_receptor_remain = new double[nst];
  475. double recip_lc_change[NLANDCOVERTYPES] = {0.0};
  476. double recip_lc_frac_transfer[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  477. double recip_transfer_remain[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  478. double net_transfer_remain[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  479. double net_lc_decrease[NLANDCOVERTYPES] = {0.0};
  480. for(int i=0;i<nst;i++) {
  481. recip_donor_remain[i] = 0.0;
  482. recip_receptor_remain[i] = 0.0;
  483. }
  484. // Quantify "reciprocal" lc change (gross lc change - net lc change)
  485. for(int from=0; from<NLANDCOVERTYPES; from++) {
  486. for(int to=0; to<NLANDCOVERTYPES; to++) {
  487. recip_lc_frac_transfer[from][to] = min(lc_frac_transfer[from][to], lc_frac_transfer[to][from]);
  488. recip_transfer_remain[from][to] = recip_lc_frac_transfer[from][to];
  489. recip_lc_change[from] += recip_lc_frac_transfer[from][to];
  490. net_transfer_remain[from][to] = lc_frac_transfer[from][to] - recip_lc_frac_transfer[from][to];
  491. net_lc_decrease[from] += net_transfer_remain[from][to];
  492. }
  493. }
  494. // Net land cover change
  495. int nsts[NLANDCOVERTYPES] = {0};
  496. int nsts_active[NLANDCOVERTYPES] = {0};
  497. bool multi_st = false;
  498. for(int lc=0; lc<NLANDCOVERTYPES; lc++) {
  499. for(int i=0; i<nst; i++) {
  500. if(lc == stlist[i].landcover) {
  501. nsts[lc]++;
  502. if(gridcell.st[i].frac || gridcell.st[i].frac_old)
  503. nsts_active[lc]++;
  504. }
  505. }
  506. if(nsts[lc] > 1)
  507. multi_st = true;
  508. }
  509. // Net transfers between land covers with only one stand type.
  510. for(int from=0; from<nst; from++) {
  511. StandType& st_donor = stlist[from];
  512. for(int to=0; to<nst; to++) {
  513. StandType& st_receptor = stlist[to];
  514. if(nsts[st_receptor.landcover] < 2 && nsts[st_donor.landcover] < 2 && net_transfer_remain[st_donor.landcover][st_receptor.landcover] > INPUT_RESOLUTION * 0.1) {
  515. st_frac_transfer[index(from, to)] += net_transfer_remain[st_donor.landcover][st_receptor.landcover];
  516. net_transfer_remain[st_donor.landcover][st_receptor.landcover] = 0.0;
  517. }
  518. }
  519. }
  520. // Net transfers to and from land covers that have several stand types.
  521. if(multi_st) {
  522. double net_lc_receptor_remain[NLANDCOVERTYPES] = {0.0};
  523. double net_transfer_remain_initial[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  524. for(int from=0; from<NLANDCOVERTYPES; from++) {
  525. for(int to=0; to<NLANDCOVERTYPES; to++) {
  526. net_lc_receptor_remain[to] += net_transfer_remain[from][to];
  527. net_transfer_remain_initial[from][to] = net_transfer_remain[from][to];
  528. }
  529. }
  530. double* influx_st = new double[nst];
  531. double* outflux_st = new double[nst];
  532. for(int i=0;i<nst;i++) {
  533. influx_st[i] = 0.0;
  534. outflux_st[i] = 0.0;
  535. }
  536. for(int from=0; from<nst; from++) {
  537. for(int to=0; to<nst; to++) {
  538. if(st_frac_transfer[index(from, to)] > 0.0) {
  539. influx_st[to] += st_frac_transfer[index(from, to)];
  540. outflux_st[from] += st_frac_transfer[index(from, to)];
  541. }
  542. }
  543. }
  544. // Transfer from land covers with more than one stand type:
  545. double* reduce = new double[nst];
  546. double* receptor_weight = new double[nst];
  547. for(int i=0;i<nst;i++) {
  548. reduce[i] = 0.0;
  549. receptor_weight[i] = 0.0;
  550. }
  551. double exclude_frac = 0.0;
  552. double exclude_frac_0 = 0.0;
  553. bool recalc_increase = false;
  554. for(int from=0; from<nst; from++) {
  555. StandType& st_donor = stlist[from];
  556. Gridcellst& gcst_donor = gridcell.st[from];
  557. if(gridcell.landcover.frac[st_donor.landcover])
  558. receptor_weight[from] = gcst_donor.frac / gridcell.landcover.frac[st_donor.landcover];
  559. }
  560. for(int from_lc=0; from_lc<NLANDCOVERTYPES; from_lc++) {
  561. double negative = 0.0;
  562. double reducable_frac = 0.0;
  563. bool recalc_reduce = false;
  564. if(!net_lc_decrease[from_lc])
  565. continue;
  566. for(int from=0; from<nst; from++) {
  567. StandType& st_donor = stlist[from];
  568. Gridcellst& gcst_donor = gridcell.st[from];
  569. if(st_donor.landcover == from_lc && nsts[st_donor.landcover] >=2) {
  570. double net_st_increase = net_lc_receptor_remain[st_donor.landcover] * receptor_weight[from];
  571. double net_st_decrease = net_st_increase - gcst_donor.frac_change;
  572. reduce[from] = (influx_st[from] - outflux_st[from]) + net_st_decrease;
  573. if(reduce[from] < 0.0) {
  574. negative += -reduce[from];
  575. recalc_reduce = true;
  576. }
  577. else {
  578. reducable_frac += reduce[from];
  579. }
  580. double remain = reduce[from] - gcst_donor.frac_old;
  581. if(remain > INPUT_RESOLUTION * 0.1) {
  582. exclude_frac_0 += gcst_donor.frac;
  583. exclude_frac += receptor_weight[from] * net_lc_receptor_remain[st_donor.landcover] - remain;
  584. if(net_lc_receptor_remain[st_donor.landcover])
  585. receptor_weight[from] -= remain / net_lc_receptor_remain[st_donor.landcover];
  586. reduce[from] = gcst_donor.frac_old;
  587. recalc_increase = true;
  588. }
  589. }
  590. }
  591. if(recalc_reduce) {
  592. for(int from=0; from<nst; from++) {
  593. StandType& st_donor = stlist[from];
  594. Gridcellst& gcst_donor = gridcell.st[from];
  595. if(st_donor.landcover == from_lc && reducable_frac && reduce[from] > 0.0) {
  596. reduce[from] = max(0.0, reduce[from] - negative * reduce[from] / reducable_frac);
  597. }
  598. else {
  599. reduce[from] = 0.0;
  600. }
  601. }
  602. }
  603. if(recalc_increase) {
  604. for(int from=0; from<nst; from++) {
  605. StandType& st_donor = stlist[from];
  606. Gridcellst& gcst_donor = gridcell.st[from];
  607. if(st_donor.landcover == from_lc) {
  608. if(receptor_weight[from] == gcst_donor.frac / gridcell.landcover.frac[st_donor.landcover] && gridcell.landcover.frac[st_donor.landcover] != exclude_frac_0) {
  609. receptor_weight[from] = gcst_donor.frac / (gridcell.landcover.frac[st_donor.landcover] - exclude_frac_0) * (net_lc_receptor_remain[st_donor.landcover] - exclude_frac) / net_lc_receptor_remain[st_donor.landcover];
  610. double net_st_increase = net_lc_receptor_remain[st_donor.landcover] * receptor_weight[from];
  611. double net_st_decrease = net_st_increase - gcst_donor.frac_change;
  612. reduce[from] = (influx_st[from] - outflux_st[from]) + net_st_decrease;
  613. }
  614. }
  615. }
  616. }
  617. }
  618. for(int from=0; from<NLANDCOVERTYPES; from++) {
  619. for(int to=0; to<NLANDCOVERTYPES; to++) {
  620. if(nsts[from] >=2 && net_transfer_remain[from][to] > INPUT_RESOLUTION * 0.1) {
  621. for(int i=0; i<nst; i++) {
  622. StandType& st_donor = stlist[i];
  623. if(st_donor.landcover == from && gridcell.landcover.frac_old[from]) {
  624. for(int j=0; j<nst; j++) {
  625. StandType& st_receptor = stlist[j];
  626. Gridcellst& gcst_receptor = gridcell.st[j];
  627. if(st_receptor.landcover == to && net_lc_decrease[from] && gridcell.landcover.frac[to]) {
  628. // This shouldn't happen really...
  629. if(reduce[i] < 0.0)
  630. dprintf("Year %d: st %d reduce[from] = %.15f !\n", date.get_calendar_year(), i, reduce[i]);
  631. double transfer = reduce[i] * net_transfer_remain_initial[from][to] / net_lc_decrease[from] * (gcst_receptor.frac / gridcell.landcover.frac[to]);
  632. st_frac_transfer[index(i, j)] += transfer;
  633. net_transfer_remain[from][to] -= transfer;
  634. net_lc_receptor_remain[to] -= transfer;
  635. }
  636. }
  637. }
  638. }
  639. net_transfer_remain[from][to] = 0.0;
  640. net_transfer_remain_initial[from][to] = 0.0;
  641. }
  642. }
  643. }
  644. // Transfer to land covers with more than one stand type:
  645. for(int from=0; from<nst; from++) {
  646. StandType& st_donor = stlist[from];
  647. Gridcellst& gcst_donor = gridcell.st[from];
  648. for(int to=0; to<nst; to++) {
  649. StandType& st_receptor = stlist[to];
  650. if(nsts[st_receptor.landcover] >= 2 && net_transfer_remain_initial[st_donor.landcover][st_receptor.landcover] > INPUT_RESOLUTION * 0.1
  651. && gridcell.landcover.frac_old[st_donor.landcover]) {
  652. double transfer = net_transfer_remain_initial[st_donor.landcover][st_receptor.landcover] * receptor_weight[to] * gcst_donor.frac_old / gridcell.landcover.frac_old[st_donor.landcover];
  653. st_frac_transfer[index(from, to)] += transfer;
  654. net_transfer_remain[st_donor.landcover][st_receptor.landcover] -= transfer;
  655. }
  656. }
  657. }
  658. if(influx_st)
  659. delete[] influx_st;
  660. if(outflux_st)
  661. delete[] outflux_st;
  662. if(reduce)
  663. delete[] reduce;
  664. if(receptor_weight)
  665. delete[] receptor_weight;
  666. }
  667. // Add "reciprocal" lc change (gross-net)
  668. // Also transitions from primary to secondary natural stands are done here
  669. double* donated_st = new double[nst];
  670. double* received_st = new double[nst];
  671. for(int i=0; i<nst; i++) {
  672. donated_st[i] = 0.0;
  673. received_st[i] = 0.0;
  674. }
  675. for(int from=0; from<nst; from++) {
  676. for(int to=0; to<nst; to++) {
  677. donated_st[from] += st_frac_transfer[index(from, to)];
  678. received_st[to] += st_frac_transfer[index(from, to)];
  679. }
  680. }
  681. double max_transfer_lc[NLANDCOVERTYPES] = {0.0};
  682. double max_donor_lc[NLANDCOVERTYPES] = {0.0};
  683. double max_receptor_lc[NLANDCOVERTYPES] = {0.0};
  684. for(int i=0; i<nst; i++) {
  685. StandType& st = stlist[i];
  686. Gridcellst& gcst = gridcell.st[i];
  687. max_transfer_lc[st.landcover] += min(gcst.frac - received_st[i], gcst.frac_old - donated_st[i]);
  688. max_receptor_lc[st.landcover] += gcst.frac - received_st[i];
  689. max_donor_lc[st.landcover] += gcst.frac_old - donated_st[i];
  690. }
  691. for(int lc=0; lc<NLANDCOVERTYPES; lc++) {
  692. if(!recip_lc_change[lc])
  693. continue;
  694. if(max_transfer_lc[lc] >= recip_lc_change[lc]) {
  695. for(int i=0; i<nst; i++) {
  696. StandType& st = stlist[i];
  697. Gridcellst& gcst = gridcell.st[i];
  698. if(st.landcover == lc && recip_lc_change[st.landcover] > 0.0) {
  699. double max_transfer_st = min(gcst.frac - received_st[i], gcst.frac_old - donated_st[i]);
  700. if(max_transfer_lc[st.landcover]) {
  701. recip_receptor_remain[i] = recip_lc_change[st.landcover] * max_transfer_st / max_transfer_lc[st.landcover];
  702. recip_donor_remain[i] = recip_lc_change[st.landcover] * max_transfer_st / max_transfer_lc[st.landcover];
  703. if(gcst.frac_old - donated_st[i] - recip_donor_remain[i] < -1e-13) {
  704. dprintf("Year %d: not enough room for reciprocal transfer for st %d, %g missing\n", date.get_calendar_year(), i, -(gcst.frac_old - donated_st[i] - recip_donor_remain[i]));
  705. }
  706. if(gcst.frac - received_st[i] - recip_receptor_remain[i] < -1e-13) {
  707. dprintf("Year %d: not enough room for reciprocal transfer for st %d, %g missing\n", date.get_calendar_year(), i, -(gcst.frac - received_st[i] - recip_receptor_remain[i]));
  708. }
  709. }
  710. }
  711. }
  712. }
  713. else {
  714. for(int i=0; i<nst; i++) {
  715. StandType& st = stlist[i];
  716. Gridcellst& gcst = gridcell.st[i];
  717. if(st.landcover == lc) {
  718. if(max_receptor_lc[st.landcover]) {
  719. recip_receptor_remain[i] = recip_lc_change[st.landcover] * (gcst.frac - received_st[i]) / max_receptor_lc[st.landcover];
  720. }
  721. if(max_donor_lc[st.landcover]) {
  722. recip_donor_remain[i] = recip_lc_change[st.landcover] * (gcst.frac_old - donated_st[i]) / max_donor_lc[st.landcover];
  723. }
  724. }
  725. }
  726. }
  727. }
  728. if(donated_st)
  729. delete[] donated_st;
  730. if(received_st)
  731. delete[] received_st;
  732. for(int from=0; from<nst; from++) {
  733. StandType& st_donor = stlist[from];
  734. Gridcellst& gcst_donor = gridcell.st[from];
  735. if(recip_donor_remain[from] > INPUT_RESOLUTION * 0.1) {
  736. for(int to=0; to<nst; to++) {
  737. StandType& st_receptor = stlist[to];
  738. Gridcellst& gcst_receptor = gridcell.st[to];
  739. if(recip_transfer_remain[st_donor.landcover][st_receptor.landcover] > INPUT_RESOLUTION * 0.1
  740. && recip_receptor_remain[to] > INPUT_RESOLUTION * 0.1 && recip_donor_remain[from] > INPUT_RESOLUTION * 0.1) {
  741. double donor_effective = min(recip_donor_remain[from], recip_transfer_remain[st_donor.landcover][st_receptor.landcover]);
  742. double receptor_effective = min(recip_receptor_remain[to], recip_transfer_remain[st_donor.landcover][st_receptor.landcover]);
  743. // all donor land is put into this stand type
  744. if(receptor_effective >= donor_effective) {
  745. st_frac_transfer[index(from, to)] += donor_effective;
  746. recip_receptor_remain[to] -= donor_effective;
  747. recip_donor_remain[from] -= donor_effective;
  748. recip_transfer_remain[st_donor.landcover][st_receptor.landcover] -= donor_effective;
  749. }
  750. // transfer to more stand types
  751. else {
  752. st_frac_transfer[index(from, to)] += receptor_effective;
  753. recip_donor_remain[from] -= receptor_effective;
  754. recip_receptor_remain[to] -= receptor_effective;
  755. recip_transfer_remain[st_donor.landcover][st_receptor.landcover] -= receptor_effective;
  756. }
  757. }
  758. }
  759. }
  760. }
  761. // Balance fractions for stand types within a landcover
  762. double* change_st = new double[nst];
  763. double* dif_st = new double[nst];
  764. for(int from=0; from<nst; from++) {
  765. change_st[from] = 0.0;
  766. dif_st[from] = 0.0;
  767. }
  768. for(int from=0; from<nst; from++) {
  769. for(int to=0; to<nst; to++) {
  770. change_st[from] -= st_frac_transfer[index(from, to)];
  771. change_st[to] += st_frac_transfer[index(from, to)];
  772. }
  773. }
  774. double dif_lc[NLANDCOVERTYPES] = {0.0};
  775. for(int from=0; from<nst; from++) {
  776. Gridcellst& gcst = gridcell.st[from];
  777. dif_st[from] = gcst.frac - (gcst.frac_old + change_st[from]); // fraction needed to reach target
  778. dif_lc[stlist[from].landcover] += fabs(dif_st[from]) / 2.0;
  779. }
  780. const double ROUNDING_ERROR = 1.0e-16;
  781. for(int lc=0; lc<NLANDCOVERTYPES; lc++) {
  782. // identify land cover with need for st frac adjustment
  783. if(nsts[lc] >= 2 && dif_lc[lc] >= INPUT_RESOLUTION * 0.1) {
  784. // find second land cover that has most transfers to or from the first land cover (with only one stand type with frac/frac_old > 0)
  785. bool error = true;
  786. for(int use_lc=0; use_lc<NLANDCOVERTYPES && error; use_lc++) {
  787. bool use_donor = true;
  788. if(nsts_active[use_lc] == 1 && max(lc_frac_transfer[lc][use_lc], lc_frac_transfer[use_lc][lc]) > dif_lc[lc]) {
  789. // determine whether to use transfer to or from land cover 2
  790. if(dif_lc[lc] <= lc_frac_transfer[use_lc][lc]) {
  791. bool skip = false;
  792. for(int st_help=0; st_help<nst; st_help++) {
  793. StandType& st_donor = stlist[st_help];
  794. if(st_donor.landcover == use_lc && (gridcell.st[st_help].frac || gridcell.st[st_help].frac_old)) {
  795. for(int st_dif=0; st_dif<nst; st_dif++) {
  796. StandType& st_receptor = stlist[st_dif];
  797. if(st_receptor.landcover == lc && !(dif_st[st_dif] >= 0.0 || st_frac_transfer[index(st_help, st_dif)] + dif_st[st_dif] > -ROUNDING_ERROR)) {
  798. skip = true;
  799. break;
  800. }
  801. }
  802. }
  803. }
  804. if(!skip) {
  805. use_donor = true;
  806. error = false;
  807. }
  808. }
  809. if(error && dif_lc[lc] <= lc_frac_transfer[lc][use_lc]) {
  810. bool skip = false;
  811. for(int st_dif=0; st_dif<nst; st_dif++) {
  812. StandType& st_donor = stlist[st_dif];
  813. if(st_donor.landcover == lc) {
  814. for(int st_help=0; st_help<nst; st_help++) {
  815. StandType& st_receptor = stlist[st_help];
  816. if(st_receptor.landcover == use_lc && (gridcell.st[st_help].frac || gridcell.st[st_help].frac_old) && !(dif_st[st_dif] <= 0.0 || st_frac_transfer[index(st_dif, st_help)] - dif_st[st_dif] > -ROUNDING_ERROR)) {
  817. skip = true;
  818. break;
  819. }
  820. }
  821. }
  822. }
  823. if(!skip) {
  824. use_donor = false;
  825. error = false;
  826. }
  827. }
  828. if(!error) {
  829. for(int st_dif=0; st_dif<nst; st_dif++) {
  830. StandType& st_donor = stlist[st_dif];
  831. if(st_donor.landcover == lc && fabs(dif_st[st_dif]) >= INPUT_RESOLUTION * 0.1) {
  832. for(int st_help=0; st_help<nst; st_help++) {
  833. StandType& st_receptor = stlist[st_help];
  834. if(st_receptor.landcover == use_lc && (gridcell.st[st_help].frac || gridcell.st[st_help].frac_old) && fabs(dif_st[st_dif]) >= INPUT_RESOLUTION * 0.1) {
  835. // after tests above, this should always work (overshoots < ROUNDING_ERROR are zeroed here)
  836. if(use_donor) {
  837. st_frac_transfer[index(st_help, st_dif)] += max(-st_frac_transfer[index(st_help, st_dif)], dif_st[st_dif]);
  838. dif_st[st_dif] = 0.0;
  839. }
  840. else {
  841. st_frac_transfer[index(st_dif, st_help)] -= min(st_frac_transfer[index(st_dif, st_help)], dif_st[st_dif]);
  842. dif_st[st_dif] = 0.0;
  843. }
  844. }
  845. }
  846. }
  847. }
  848. }
  849. }
  850. }
  851. if(error) {
  852. // If stand type fractions still unbalanced, try shuffling within a land cover
  853. for(int from=0; from<nst; from++) {
  854. StandType& st_donor = stlist[from];
  855. if(st_donor.landcover == lc && -dif_st[from] > INPUT_RESOLUTION * 0.1) {
  856. for(int to=0; to<nst; to++) {
  857. StandType& st_receptor = stlist[to];
  858. if(st_donor.landcover == st_receptor.landcover && -dif_st[from] > INPUT_RESOLUTION * 0.1 && dif_st[to] >= INPUT_RESOLUTION * 0.1) {
  859. double donor_effective = -dif_st[from];
  860. double receptor_effective = dif_st[to];
  861. // all donor land going to this land cover is put into this stand type
  862. if(receptor_effective >= donor_effective) {
  863. st_frac_transfer[index(from, to)] += donor_effective;
  864. lc_frac_transfer[st_donor.landcover][st_donor.landcover] += donor_effective;
  865. dif_st[from] += donor_effective;
  866. dif_st[to] -= donor_effective;
  867. }
  868. // transfer to more stand types within this land cover
  869. else {
  870. st_frac_transfer[index(from, to)] += receptor_effective;
  871. lc_frac_transfer[st_donor.landcover][st_donor.landcover] += receptor_effective;
  872. dif_st[from] += receptor_effective;
  873. dif_st[to] -= receptor_effective;
  874. }
  875. }
  876. }
  877. }
  878. }
  879. }
  880. }
  881. }
  882. if(change_st)
  883. delete[] change_st;
  884. if(dif_st)
  885. delete[] dif_st;
  886. // Remove very small transfer values derived from rounding errors in the balancing section above
  887. for(int from=0; from<nst; from++) {
  888. for(int to=0; to<nst; to++) {
  889. if(st_frac_transfer[index(from, to)] < ROUNDING_ERROR)
  890. st_frac_transfer[index(from, to)] = 0.0;
  891. }
  892. }
  893. for(int from=0; from<nst; from++) {
  894. StandType& st_donor = stlist[from];
  895. Gridcellst& gcst_donor = gridcell.st[from];
  896. for(int to=0; to<nst; to++) {
  897. StandType& st_receptor = stlist[to];
  898. Gridcellst& gcst_receptor = gridcell.st[to];
  899. double transfer_frac = st_frac_transfer[index(from, to)];
  900. if(transfer_frac > (gcst_donor.frac_old) || transfer_frac > (gcst_receptor.frac)) {
  901. if((transfer_frac - gcst_receptor.frac) > INPUT_RESOLUTION || (transfer_frac - gcst_receptor.frac) > INPUT_RESOLUTION) {
  902. dprintf("\nYear %d: st_change_array sum not compatible with st.frac_old/frac values for stand types %d and %d\n", date.year, from, to);
  903. dprintf("st_change_array=%.15f, st_donor.frac_old=%.15f, st_receptor.frac=%.15f\n", transfer_frac, gcst_donor.frac_old, gcst_receptor.frac);
  904. }
  905. if(transfer_frac > (gcst_receptor.frac) && transfer_frac < (gcst_receptor.frac + INPUT_RESOLUTION * 100.0)) {
  906. if((transfer_frac - gcst_receptor.frac) > INPUT_RESOLUTION)
  907. dprintf("Removing overshoot in array %.15f\n\n", st_frac_transfer[index(from, to)] - gcst_receptor.frac);
  908. st_frac_transfer[index(from, to)] = gcst_receptor.frac;
  909. }
  910. else if(transfer_frac > (gcst_donor.frac_old) && transfer_frac < (gcst_donor.frac_old + INPUT_RESOLUTION * 100.0)) {
  911. if((transfer_frac - gcst_donor.frac_old) > INPUT_RESOLUTION)
  912. dprintf("Removing overshoot in array %.15f\n\n", st_frac_transfer[index(from, to)] - gcst_donor.frac_old);
  913. st_frac_transfer[index(from, to)] = gcst_donor.frac_old;
  914. }
  915. }
  916. }
  917. }
  918. // Set the transfer fraction from primary stands within a stand type (simple case: equal primary/secondary ratio for all stand types in a lc->lc transfer)
  919. if(use_primary_lc_transfer) {
  920. for(int from=0; from<nst; from++) {
  921. StandType& st_donor = stlist[from];
  922. for(int to=0; to<nst; to++) {
  923. StandType& st_receptor = stlist[to];
  924. double primary_frac = 0.0;
  925. double secondary_young_frac = 0.0;
  926. if(lc_frac_transfer[st_donor.landcover][st_receptor.landcover] > 0.0) {
  927. primary_frac = forest_lc_frac_transfer_s.primary[st_donor.landcover][st_receptor.landcover] / lc_frac_transfer[st_donor.landcover][st_receptor.landcover];
  928. secondary_young_frac = forest_lc_frac_transfer_s.secondary_young[st_donor.landcover][st_receptor.landcover] / lc_frac_transfer[st_donor.landcover][st_receptor.landcover];
  929. }
  930. forest_st_frac_transfer_s.primary[index(from, to)] = primary_frac * st_frac_transfer[index(from, to)];
  931. forest_st_frac_transfer_s.secondary_young[index(from, to)] = secondary_young_frac * st_frac_transfer[index(from, to)];
  932. }
  933. }
  934. }
  935. if(recip_donor_remain)
  936. delete[] recip_donor_remain;
  937. if(recip_receptor_remain)
  938. delete[] recip_receptor_remain;
  939. }
  940. /// Handles harvest and turnover of reduced stands at landcover change.
  941. /** Sets landcover.updated to true
  942. * Stores carbon, nitrogen and water of harvested area in a temporary struct.
  943. * Should be followed by a call to stand_dynamics() to kill stands with a new area of 0
  944. *
  945. * INPUT PARAMETERS
  946. *
  947. * \param receiving_fraction sum of added area to expanding stands
  948. * \param landcover_receptor restriction of donor stands transferring to specified receptor landcover
  949. * \param landcover_donor restriction of donor stands according to landcover
  950. * \param stid_receptor restriction of donor stands transferring to specified receptor stand type
  951. * \param stid_donor restriction of stands according to stand type
  952. * \param standid restriction of donor stand
  953. *
  954. * OUTPUT PARAMETERS
  955. * \param landcover_change_transfer struct containing the following pft-specific public members:
  956. * - transfer_litter_leaf
  957. * - transfer_litter_sap
  958. * - transfer_litter_heart
  959. * - transfer_litter_root
  960. * - transfer_litter_repr
  961. * - transfer_nmass_litter_leaf
  962. * - transfer_nmass_litter_root
  963. * - transfer_nmass_litter_sap
  964. * - transfer_nmass_litter_heart
  965. * - transfer_harvested_products_slow
  966. * - transfer_harvested_products_slow_nmass
  967. * ,the following patch-level public members
  968. * - transfer_cpool_fast
  969. * - transfer_cpool_slow
  970. * - transfer_nmass_avail
  971. * - transfer_wcont_evap
  972. * - transfer_snowpack
  973. * - transfer_decomp_litter_mean
  974. * - transfer_k_soilfast_mean
  975. * - transfer_acflux_harvest
  976. * - transfer_anflux_harvest
  977. * ,the following water soil layer-specific public member:
  978. * - transfer_wcont
  979. * and the following century soil pool-specific public members:
  980. * - transfer_sompool.cmass
  981. * - transfer_sompool.fireresist
  982. * - transfer_sompool.fracremain
  983. * - transfer_sompool.ligcfrac
  984. * - transfer_sompool.nmass
  985. * - transfer_sompool.ntoc
  986. */
  987. void donor_stand_change(Gridcell& gridcell, double& receiving_fraction, landcover_change_transfer& to,
  988. int landcover_receptor = -1, int landcover_donor = -1,
  989. int stid_receptor = -1, int stid_donor = -1,
  990. int standid = -1, lc_change_harvest_params* harv_params = NULL) {
  991. int count = 0;
  992. Landcover& lc = gridcell.landcover;
  993. for(unsigned int i=0; i<gridcell.nbr_stands(); i++) {
  994. Stand& stand = gridcell[i];
  995. double donor_area = 0.0;
  996. if(stid_receptor >= 0)
  997. donor_area = stand.transfer_area_st[stid_receptor];
  998. else if(landcover_receptor >= 0)
  999. donor_area = stand.transfer_area_lc((landcovertype)landcover_receptor);
  1000. else
  1001. donor_area = stand.gross_frac_decrease;
  1002. if(donor_area > 0.0 && (stand.id == standid || standid < 0) && (stand.landcover == landcover_donor || landcover_donor < 0) && (stand.stid == stid_donor || stid_donor < 0)) {
  1003. count++;
  1004. stand.frac_temp -= donor_area;
  1005. stand.frac_temp = max(0.0, stand.frac_temp);
  1006. double scale = donor_area / receiving_fraction / (double)stand.nobj;
  1007. // Add non-living C, N and water of stand to transfer struct:
  1008. to.add_from_stand(stand, scale);
  1009. // Harvest and turnover of copies of individuals, add to transfer copy:
  1010. stand.firstobj();
  1011. while(stand.isobj) {
  1012. Patch& patch = stand.getobj();
  1013. Vegetation& vegetation = patch.vegetation;
  1014. vegetation.firstobj();
  1015. while(vegetation.isobj) {
  1016. Harvest_CN cp;
  1017. Individual& indiv = vegetation.getobj();
  1018. Patchpft& patchpft = patch.pft[indiv.pft.id];
  1019. if(indiv.has_daily_turnover())
  1020. cp.copy_from_indiv(indiv, true, false);
  1021. else
  1022. cp.copy_from_indiv(indiv, false, false);
  1023. bool lc_pool = (landcover_receptor == -1) && (stid_receptor == -1);
  1024. bool wood_harvest = false;
  1025. if(stand.landcover == NATURAL || stand.landcover == FOREST) {
  1026. if(lc_pool) {
  1027. if(stand.transfer_area_lc((landcovertype)NATURAL) || stand.transfer_area_lc((landcovertype)FOREST))
  1028. wood_harvest = true;
  1029. }
  1030. else if(landcover_receptor == NATURAL || landcover_receptor == FOREST || stid_receptor > -1 && (stlist[stid_receptor].landcover == NATURAL || stlist[stid_receptor].landcover == FOREST)) {
  1031. wood_harvest = true;
  1032. }
  1033. }
  1034. // Harvest of transferred areas:
  1035. switch (stand.landcover)
  1036. {
  1037. case CROPLAND:
  1038. harvest_crop(cp, indiv.pft, indiv.alive, indiv.cropindiv->isintercropgrass);
  1039. break;
  1040. case PASTURE:
  1041. harvest_pasture(cp, indiv.pft, indiv.alive);
  1042. break;
  1043. case NATURAL:
  1044. case FOREST:
  1045. case URBAN:
  1046. case PEATLAND:
  1047. if(harv_params) {
  1048. harvest_wood(cp, indiv.pft, indiv.alive, 1.0, harv_params->harv_eff, harv_params->res_outtake_twig, harv_params->res_outtake_coarse_root);
  1049. }
  1050. else if(wood_harvest) {
  1051. if(lc_pool) {
  1052. double lu_change_ratio = 1.0 - (stand.transfer_area_lc((landcovertype)NATURAL) + stand.transfer_area_lc((landcovertype)FOREST)) / stand.gross_frac_decrease;
  1053. // Land cover change area: frac_cut=1, harv_eff=1, res_outtake_twig=0.95, res_outtake_coarse_root=0.9
  1054. // Wood harvest area: frac_cut=1, harv_eff=pft.harv_eff, res_outtake_twig=pft.res_outtake, res_outtake_coarse_root=0
  1055. harvest_wood(cp, indiv.pft, indiv.alive, 1.0, 1.0 * lu_change_ratio + indiv.pft.harv_eff * (1.0 - lu_change_ratio), 0.95 * lu_change_ratio + indiv.pft.res_outtake * (1.0 - lu_change_ratio), 0.9 * lu_change_ratio);
  1056. }
  1057. else {
  1058. harvest_wood(cp, indiv.pft, indiv.alive, 1.0, indiv.pft.harv_eff, indiv.pft.res_outtake);
  1059. }
  1060. }
  1061. else {
  1062. harvest_wood(cp, indiv.pft, indiv.alive, 1.0, 1.0, 0.95, 0.9); // frac_cut=1, harv_eff=1, res_outtake_twig=0.95, res_outtake_coarse_root=0.9
  1063. }
  1064. break;
  1065. case BARREN: // Assuming there is nothing to harvest on barren
  1066. break;
  1067. default:
  1068. fail("Modify code to deal with landcover harvest at landcover change!\n");
  1069. }
  1070. turnover(indiv.pft.turnover_leaf, indiv.pft.turnover_root,
  1071. indiv.pft.turnover_sap, indiv.pft.lifeform, indiv.pft.landcover,
  1072. cp.cmass_leaf, cp.cmass_root, cp.cmass_sap, cp.cmass_heart,
  1073. cp.nmass_leaf, cp.nmass_root, cp.nmass_sap, cp.nmass_heart,
  1074. cp.litter_leaf,
  1075. cp.litter_root,
  1076. cp.nmass_litter_leaf,
  1077. cp.nmass_litter_root,
  1078. cp.nstore_longterm, cp.max_n_storage,
  1079. indiv.alive);
  1080. // In case any vegetation left (eg. cmass_root in pasture or grass in woodland):
  1081. kill_remaining_vegetation(cp, indiv.pft, indiv.alive, indiv.istruecrop_or_intercropgrass(), false);
  1082. //Sum added litter C & N:
  1083. to.transfer_litter_leaf[indiv.pft.id] += cp.litter_leaf * scale;
  1084. to.transfer_litter_root[indiv.pft.id] += cp.litter_root * scale;
  1085. to.transfer_litter_sap[indiv.pft.id] += cp.litter_sap * scale;
  1086. to.transfer_litter_heart[indiv.pft.id] += cp.litter_heart * scale;
  1087. to.transfer_nmass_litter_leaf[indiv.pft.id] += cp.nmass_litter_leaf * scale;
  1088. to.transfer_nmass_litter_root[indiv.pft.id] += cp.nmass_litter_root * scale;
  1089. to.transfer_nmass_litter_sap[indiv.pft.id] += cp.nmass_litter_sap * scale;
  1090. to.transfer_nmass_litter_heart[indiv.pft.id] += cp.nmass_litter_heart * scale;
  1091. // Flux to products
  1092. lc.harvest_product += cp.harvested_products_slow * donor_area / (double)stand.nobj;
  1093. lc.harvest_product_lc[stand.landcover] += cp.harvested_products_slow * donor_area / (double)stand.nobj;
  1094. lc.harvest_product_nmass += cp.harvested_products_slow_nmass * donor_area / (double)stand.nobj;
  1095. lc.harvest_product_nmass_lc[stand.landcover] += cp.harvested_products_slow_nmass * donor_area / (double)stand.nobj;
  1096. // Flux to litter
  1097. lc.acflux_harvest_wood_res += cp.dcflux_harvest_wood_res * donor_area / (double)stand.nobj;
  1098. lc.acflux_harvest_wood_res_lc[stand.landcover] += cp.dcflux_harvest_wood_res * donor_area / (double)stand.nobj;
  1099. // Flux to atm
  1100. lc.dcflux_landuse_change += cp.acflux_harvest * donor_area / (double)stand.nobj; // ecev3 - reset to 0 each day
  1101. lc.acflux_landuse_change += cp.acflux_harvest * donor_area / (double)stand.nobj;
  1102. lc.acflux_landuse_change_lc[stand.landcover] += cp.acflux_harvest * donor_area / (double)stand.nobj;
  1103. lc.anflux_landuse_change += cp.anflux_harvest * donor_area / (double)stand.nobj;
  1104. lc.anflux_landuse_change_lc[stand.landcover] += cp.anflux_harvest * donor_area / (double)stand.nobj;
  1105. // gridcell.acflux_landuse_change += -cp.debt_excess * donor_area / (double)stand.nobj;
  1106. if(ifslowharvestpool) {
  1107. to.transfer_harvested_products_slow[indiv.pft.id] += cp.harvested_products_slow * scale;
  1108. to.transfer_harvested_products_slow_nmass[indiv.pft.id] += cp.harvested_products_slow_nmass * scale;
  1109. }
  1110. vegetation.nextobj();
  1111. }
  1112. stand.nextobj();
  1113. }
  1114. }
  1115. }
  1116. }
  1117. /// Creates and kills stands at landcover change.
  1118. /** Harvest of reduced stands need to be done before with donor_stand_change().
  1119. * Should be followed by a call to receiving_stand_change() for transfer of
  1120. * carbon, nitrogen and water.
  1121. *
  1122. */
  1123. void stand_dynamics(Gridcell& gridcell) {
  1124. stlist.firstobj();
  1125. while (stlist.isobj) {
  1126. StandType& st=stlist.getobj();
  1127. Gridcellst& gcst = gridcell.st[st.id];
  1128. landcovertype lc = st.landcover;
  1129. bool expand_to_new_stand = gridcell.landcover.expand_to_new_stand[lc];
  1130. if(gcst.gross_frac_increase || gcst.gross_frac_decrease || gcst.frac == 0.0) {
  1131. double stand_sum = 0.0;
  1132. Gridcell::iterator gc_itr = gridcell.begin();
  1133. while (gc_itr != gridcell.end()) {
  1134. Stand& stand = *gc_itr;
  1135. if(stand.stid == st.id)
  1136. stand_sum += stand.get_gridcell_fraction();
  1137. ++gc_itr;
  1138. }
  1139. // first stand created
  1140. if((gcst.frac_old == 0.0 || stand_sum == 0.0) && gcst.gross_frac_increase > 0.0) {
  1141. int npatch;
  1142. if(date.year == 0)
  1143. // with npatch=0, number of patches will be chosen in the stand constructor
  1144. npatch = 0;
  1145. else
  1146. npatch = npatch_secondarystand;
  1147. Stand& stand = gridcell.create_stand_lu(st, gcst.gross_frac_increase, npatch);
  1148. }
  1149. // last stand killed
  1150. else if(gcst.frac_old > 0.0 && !gcst.frac) {
  1151. Gridcell::iterator gc_itr = gridcell.begin();
  1152. while (gc_itr != gridcell.end()) {
  1153. Stand& stand = *gc_itr;
  1154. if(stand.stid == st.id) {
  1155. gridcell.landcover.acflux_landuse_change += stand.ccont() * stand.get_gridcell_fraction();
  1156. gridcell.landcover.anflux_landuse_change += stand.ncont() * stand.get_gridcell_fraction();
  1157. gc_itr = gridcell.delete_stand(gc_itr);
  1158. }
  1159. else {
  1160. ++gc_itr;
  1161. }
  1162. }
  1163. }
  1164. else {
  1165. if(expand_to_new_stand) {
  1166. // new stand created from other landcover types (pooled option, not created in transfer_to_new_stand())
  1167. if(gcst.gross_frac_increase > 0.0) {
  1168. // Fewer patches in secondary stands if npatch_secondarystand < npatch:
  1169. Stand& stand = gridcell.create_stand_lu(st, gcst.gross_frac_increase, npatch_secondarystand);
  1170. }
  1171. }
  1172. // secondary stand killed if all of its area converted to other landcover types
  1173. if(gcst.nstands > 1 && gcst.gross_frac_decrease > 0.0) {
  1174. Gridcell::iterator gc_itr = gridcell.begin();
  1175. while (gc_itr != gridcell.end()) {
  1176. Stand& stand = *gc_itr;
  1177. if(stand.stid == st.id && stand.get_gridcell_fraction() == 0) {
  1178. gc_itr = gridcell.delete_stand(gc_itr);
  1179. if(gcst.frac < 1e-14)
  1180. gcst.frac = 0.0;
  1181. }
  1182. else {
  1183. ++gc_itr;
  1184. }
  1185. }
  1186. }
  1187. }
  1188. }
  1189. stlist.nextobj();
  1190. }
  1191. }
  1192. /// Transfers litter etc. of reduced stands to expanding stands at landcover change.
  1193. /** Transfers carbon, nitrogen and water of harvested area to expanded areas from a temporary struct.
  1194. * Is additive - can be called several times to the same landcover or standtype, if the parameter donorfrac_rel is used
  1195. *
  1196. * INPUT PARAMETERS
  1197. *
  1198. * \param landcover_change_transfer& from struct containing the following pft-specific public members:
  1199. * - transfer_litter_leaf
  1200. * - transfer_litter_sap
  1201. * - transfer_litter_heart
  1202. * - transfer_litter_root
  1203. * - transfer_litter_repr
  1204. * - transfer_nmass_litter_leaf
  1205. * - transfer_nmass_litter_root
  1206. * - transfer_nmass_litter_sap
  1207. * - transfer_nmass_litter_heart
  1208. * - transfer_harvested_products_slow
  1209. * - transfer_harvested_products_slow_nmass
  1210. * ,the following patch-level public members:
  1211. * - transfer_cpool_fast
  1212. * - transfer_cpool_slow
  1213. * - transfer_nmass_avail
  1214. * - transfer_cpool_slow
  1215. * - transfer_wcont_evap
  1216. * - transfer_snowpack
  1217. * - transfer_decomp_litter_mean
  1218. * - transfer_k_soilfast_mean
  1219. * - transfer_acflux_harvest
  1220. * - transfer_anflux_harvest
  1221. * ,the following water soil layer-specific public member:
  1222. * - transfer_wcont
  1223. * and the following century soil pool-specific public members:
  1224. * - transfer_sompool.cmass
  1225. * - transfer_sompool.fireresist
  1226. * - transfer_sompool.fracremain
  1227. * - transfer_sompool.ligcfrac
  1228. * - transfer_sompool.nmass
  1229. * - transfer_sompool.ntoc
  1230. *
  1231. * \param LCchangeCtransfer whether to transfer carbon, nitrogen and water of reduced stands to expanding stands
  1232. * \param landcover restricts receiving stands to a certain landcover, default no restriction
  1233. * \param stid restricts receiving stands to a certain standtype, default no restriction
  1234. * \param donorfrac_rel relative part of fraction (to a receiving stand) from a particular donor
  1235. * \param standid restricts receiving stands to a certain stand, default no restriction
  1236. */
  1237. void receiving_stand_change(Gridcell& gridcell, landcover_change_transfer& from,
  1238. bool LCchangeCtransfer, int landcover = -1, int stid = -1,
  1239. double donorfrac_rel = 1.0, int standid = -1) {
  1240. Gridcell::iterator gc_itr = gridcell.begin();
  1241. while (gc_itr != gridcell.end()) {
  1242. Stand& stand = *gc_itr;
  1243. if(stand.gross_frac_increase > 0.0 && (stand.landcover == landcover || landcover < 0) &&
  1244. (stand.stid == stid || stid < 0 && (stand.id == standid || standid < 0))) {
  1245. double old_frac = stand.frac_temp;
  1246. double added_frac = donorfrac_rel * stand.gross_frac_increase;
  1247. double new_frac = old_frac + added_frac;
  1248. stand.frac_temp += added_frac;
  1249. if(LCchangeCtransfer) {
  1250. stand.firstobj();
  1251. while(stand.isobj) {
  1252. Patch& patch=stand.getobj();
  1253. // add litter C & N:
  1254. for (int i=0; i<npft; i++) {
  1255. Patchpft& patchpft = patch.pft[i];
  1256. patchpft.litter_leaf = (patchpft.litter_leaf * old_frac + from.transfer_litter_leaf[i] * added_frac) / new_frac;
  1257. patchpft.litter_sap = (patchpft.litter_sap * old_frac + from.transfer_litter_sap[i] * added_frac) / new_frac;
  1258. patchpft.litter_heart = (patchpft.litter_heart * old_frac + from.transfer_litter_heart[i] * added_frac) / new_frac;
  1259. patchpft.litter_root = (patchpft.litter_root * old_frac + from.transfer_litter_root[i] * added_frac) / new_frac;
  1260. patchpft.litter_repr = (patchpft.litter_repr * old_frac + from.transfer_litter_repr[i] * added_frac) / new_frac;
  1261. patchpft.nmass_litter_leaf = (patchpft.nmass_litter_leaf * old_frac + from.transfer_nmass_litter_leaf[i] * added_frac) / new_frac;
  1262. patchpft.nmass_litter_root = (patchpft.nmass_litter_root * old_frac + from.transfer_nmass_litter_root[i] * added_frac) / new_frac;
  1263. patchpft.nmass_litter_sap = (patchpft.nmass_litter_sap * old_frac + from.transfer_nmass_litter_sap[i] * added_frac) / new_frac;
  1264. patchpft.nmass_litter_heart = (patchpft.nmass_litter_heart * old_frac + from.transfer_nmass_litter_heart[i] * added_frac) / new_frac;
  1265. if(ifslowharvestpool) {
  1266. patchpft.harvested_products_slow = (patchpft.harvested_products_slow * old_frac + from.transfer_harvested_products_slow[i] * added_frac) / new_frac;
  1267. patchpft.harvested_products_slow_nmass = (patchpft.harvested_products_slow_nmass * old_frac + from.transfer_harvested_products_slow_nmass[i] * added_frac) / new_frac;
  1268. }
  1269. }
  1270. // add soil C & N:
  1271. patch.soil.cpool_fast = (patch.soil.cpool_fast * old_frac + from.transfer_cpool_fast * added_frac) / new_frac;
  1272. patch.soil.cpool_slow = (patch.soil.cpool_slow * old_frac + from.transfer_cpool_slow * added_frac) / new_frac;
  1273. for(int i=0; i<NSOMPOOL; i++) {
  1274. patch.soil.sompool[i].cmass = (patch.soil.sompool[i].cmass * old_frac + from.transfer_sompool[i].cmass * added_frac) / new_frac;
  1275. patch.soil.sompool[i].fireresist = (patch.soil.sompool[i].fireresist * old_frac + from.transfer_sompool[i].fireresist * added_frac) / new_frac;
  1276. patch.soil.sompool[i].fracremain = (patch.soil.sompool[i].fracremain * old_frac + from.transfer_sompool[i].fracremain * added_frac) / new_frac;
  1277. patch.soil.sompool[i].ligcfrac = (patch.soil.sompool[i].ligcfrac *old_frac + from.transfer_sompool[i].ligcfrac * added_frac) / new_frac;
  1278. patch.soil.sompool[i].litterme = (patch.soil.sompool[i].litterme * old_frac + from.transfer_sompool[i].litterme * added_frac) / new_frac;
  1279. patch.soil.sompool[i].nmass = (patch.soil.sompool[i].nmass * old_frac + from.transfer_sompool[i].nmass * added_frac) / new_frac;
  1280. patch.soil.sompool[i].ntoc = (patch.soil.sompool[i].ntoc * old_frac + from.transfer_sompool[i].ntoc * added_frac) / new_frac;
  1281. }
  1282. patch.soil.nmass_avail = (patch.soil.nmass_avail * old_frac + from.transfer_nmass_avail * added_frac) / new_frac;
  1283. // ecev3 - check
  1284. const double EPS = 1.0e-16;
  1285. if (patch.soil.nmass_avail < 0.0) {
  1286. if (patch.soil.nmass_avail >= -EPS)
  1287. patch.soil.nmass_avail = 0.0; // Reset to 0 if negative but tiny
  1288. else {
  1289. dprintf("ERROR! Year %d day %d Stand %d: Negative soil.nmass_avail in receiving_stand_change: %.15f\n", date.year, date.day, patch.stand.id, patch.soil.nmass_avail);
  1290. dprintf("%.15f %.15f %.15f\n", from.transfer_nmass_avail, added_frac, new_frac);
  1291. }
  1292. }
  1293. // add other soil stuff:
  1294. for(int i=0; i<NSOILLAYER; i++) {
  1295. patch.soil.wcont[i] = (patch.soil.wcont[i] * old_frac + from.transfer_wcont[i] * added_frac) / new_frac;
  1296. }
  1297. patch.soil.wcont_evap = (patch.soil.wcont_evap * old_frac + from.transfer_wcont_evap * added_frac) / new_frac;
  1298. patch.soil.snowpack = (patch.soil.snowpack * old_frac + from.transfer_snowpack * added_frac) / new_frac;
  1299. patch.soil.snowpack_nmass = (patch.soil.snowpack_nmass * old_frac + from.transfer_snowpack_nmass * added_frac) / new_frac;
  1300. patch.soil.decomp_litter_mean = (patch.soil.decomp_litter_mean * old_frac + from.transfer_decomp_litter_mean * added_frac) / new_frac;
  1301. patch.soil.k_soilfast_mean = (patch.soil.k_soilfast_mean * old_frac + from.transfer_k_soilfast_mean * added_frac) / new_frac;
  1302. patch.soil.k_soilslow_mean = (patch.soil.k_soilslow_mean * old_frac + from.transfer_k_soilslow_mean * added_frac) / new_frac;
  1303. double aaet_5_cp[NYEARAAET] = {0.0};
  1304. patch.aaet_5.to_array(aaet_5_cp);
  1305. for(unsigned int i=0;i<from.transfer_aaet_5.size();i++)
  1306. patch.aaet_5.add((aaet_5_cp[i] * old_frac + from.transfer_aaet_5[i] * added_frac) / new_frac);
  1307. patch.soil.anfix_calc = (patch.soil.anfix_calc * old_frac + from.transfer_anfix_calc * added_frac) / new_frac;
  1308. // save individual C and N content for use in scale_indiv()
  1309. for(unsigned int i=0; i<patch.vegetation.nobj ;i++) {
  1310. Individual& indiv = patch.vegetation[i];
  1311. indiv.save_cmass_luc();
  1312. indiv.save_nmass_luc();
  1313. }
  1314. stand.nextobj();
  1315. }
  1316. // set scaling factor to be used in growth() for scaling vegetation C and N:
  1317. stand.scale_LC_change = (stand.frac_old - stand.gross_frac_decrease + min(0.0, stand.cloned_fraction)) / stand.get_gridcell_fraction();
  1318. gridcell.landcover.updated = true;
  1319. }
  1320. }
  1321. ++gc_itr;
  1322. }
  1323. }
  1324. enum {NONEWSTAND, CLONESTAND, CLONESTAND_KILLTREES, NEWSTAND_KILLALL};
  1325. /// contains rules for creation of new stands at land cover change
  1326. /** Options CLONESTAND and CLONESTAND_KILLTREES require that the receptor landcover
  1327. * allows growth of natural grass and/or tree PFTs
  1328. *
  1329. * INTPUT PARAMETERS
  1330. *
  1331. * \param landcover_donor landcover type of donor stand
  1332. * \param landcover_receptor landcover type of rexeptor stand
  1333. */
  1334. int copy_stand_type(int landcover_donor, int landcover_receptor) {
  1335. int copy_type = NONEWSTAND;
  1336. if(landcover_donor == CROPLAND) {
  1337. if(landcover_receptor == NATURAL || landcover_receptor == FOREST)
  1338. copy_type = NEWSTAND_KILLALL;
  1339. }
  1340. else if(landcover_donor == PASTURE) {
  1341. if(landcover_receptor == NATURAL || landcover_receptor == FOREST)
  1342. copy_type = CLONESTAND;
  1343. }
  1344. else if(landcover_donor == NATURAL || landcover_donor == FOREST) {
  1345. if(landcover_receptor == FOREST || landcover_receptor == NATURAL)
  1346. copy_type = CLONESTAND; // or CLONESTAND_KILLTREES/NEWSTAND_KILLALL for clearcut
  1347. }
  1348. return copy_type;
  1349. }
  1350. /// Creates unique stands from transfer events if copy_stand_type() returns >= 1 for landcovers combination
  1351. /** Either clones donor stand or creates new stand from scratch
  1352. *
  1353. * INPUT PARAMETERS
  1354. *
  1355. * \param stid_donor donor stand type
  1356. * \param stid_receptor receptor stand type
  1357. */
  1358. double transfer_to_new_stand(Gridcell& gridcell, int stid_donor, int stid_receptor) {
  1359. double cloned_area = 0.0;
  1360. for(unsigned int i=0; i<gridcell.nbr_stands(); i++) {
  1361. Stand& stand = gridcell[i];
  1362. double transfer_area = stand.transfer_area_st[stid_receptor];
  1363. if(transfer_area > 0.0 && (stand.stid == stid_donor || stid_donor < 0)) {
  1364. int copy_type = copy_stand_type(stand.landcover, stlist[stid_receptor].landcover); // NONEWSTAND, CLONESTAND, CLONESTAND_KILLTREES, NEWSTAND_KILLALL
  1365. if(copy_type) {
  1366. if(copy_type == CLONESTAND || copy_type == CLONESTAND_KILLTREES) {
  1367. Stand& new_stand = stand.clone(stlist[stid_receptor], transfer_area);
  1368. if(stand.landcover == NATURAL && (stlist[stid_receptor].naturalveg == ""))
  1369. dprintf("WARNING: cloning natural stand without allowing natural pft:s to grow in the new stand. Was this intended ?\n");
  1370. landcover_change_transfer transfer;
  1371. new_stand.cloned = true;
  1372. new_stand.gross_frac_increase = 0.0;
  1373. new_stand.frac_change = 0.0;
  1374. new_stand.cloned_fraction = transfer_area;
  1375. new_stand.origin = stand.landcover;
  1376. new_stand.firstobj();
  1377. while(new_stand.isobj) {
  1378. Patch& patch = new_stand.getobj();
  1379. Vegetation& vegetation = patch.vegetation;
  1380. vegetation.firstobj();
  1381. while(vegetation.isobj) {
  1382. Individual& indiv = vegetation.getobj();
  1383. Standpft& standpft = new_stand.pft[indiv.pft.id];
  1384. if(!standpft.active) {
  1385. // Treatment of pft individuals not allowed to grow anymore in the new stand:
  1386. // Clearcut
  1387. harvest_wood(indiv, 1.0, 1.0, 0.95, 0.9, true); // frac_cut=1, harv_eff=1, res_outtake_twig=0.95, res_outtake_coarse_root=0.9
  1388. // Grass killed, C+N goes to soil
  1389. kill_remaining_vegetation(indiv);
  1390. vegetation.killobj();
  1391. }
  1392. else if(copy_type == CLONESTAND_KILLTREES && indiv.pft.lifeform != GRASS) {
  1393. harvest_wood(indiv, 1.0, 1.0, 0.95, 0.9, true); // frac_cut=1, harv_eff=1, res_outtake_twig=0.95, res_outtake_coarse_root=0.9
  1394. vegetation.killobj();
  1395. }
  1396. else {
  1397. vegetation.nextobj();
  1398. }
  1399. }
  1400. new_stand.nextobj();
  1401. }
  1402. }
  1403. else if(copy_type == NEWSTAND_KILLALL) {
  1404. landcover_change_transfer transfer;
  1405. lc_change_harvest_params harvest_params;
  1406. // default values for natural to cropland transition
  1407. harvest_params.harv_eff = 1.0;
  1408. harvest_params.res_outtake_twig = 0.95;
  1409. harvest_params.res_outtake_coarse_root = 0.9;
  1410. donor_stand_change(gridcell, transfer_area, transfer, -1,-1, stid_receptor, -1, stand.id, &harvest_params);
  1411. Stand& new_stand = gridcell.create_stand_lu(stlist[stid_receptor], transfer_area, npatch_secondarystand);
  1412. receiving_stand_change(gridcell, transfer, true, -1, -1, 1.0, new_stand.id);
  1413. new_stand.cloned = true;
  1414. new_stand.gross_frac_increase = 0.0;
  1415. new_stand.frac_change = 0.0;
  1416. new_stand.cloned_fraction = transfer_area;
  1417. }
  1418. cloned_area += transfer_area;
  1419. stand.cloned_fraction -= transfer_area;
  1420. gridcell.st[stid_donor].gross_frac_decrease -= transfer_area;
  1421. gridcell.st[stid_donor].frac_change += transfer_area;
  1422. stand.gross_frac_decrease -= transfer_area;
  1423. stand.frac_change += transfer_area;
  1424. stand.transfer_area_st[stid_receptor] -= transfer_area;
  1425. gridcell.st[stid_receptor].gross_frac_increase -= transfer_area;
  1426. gridcell.st[stid_receptor].frac_change -= transfer_area;
  1427. if(gridcell.st[stid_donor].gross_frac_decrease < INPUT_RESOLUTION * 0.1)
  1428. gridcell.st[stid_donor].gross_frac_decrease = 0.0;
  1429. if(stand.gross_frac_decrease < INPUT_RESOLUTION * 0.1)
  1430. stand.gross_frac_decrease = 0.0;
  1431. if(stand.transfer_area_st[stid_receptor] < INPUT_RESOLUTION * 0.1)
  1432. stand.transfer_area_st[stid_receptor] = 0.0;
  1433. if(gridcell.st[stid_receptor].gross_frac_increase < INPUT_RESOLUTION * 0.1)
  1434. gridcell.st[stid_receptor].gross_frac_increase = 0.0;
  1435. }
  1436. }
  1437. }
  1438. return cloned_area;
  1439. }
  1440. /// Simulates gross landcover change by adding a specified fraction to the lc_frac_transfer array
  1441. void simulate_gross_lc_transfer(Gridcell& gridcell, double lc_frac_transfer[][NLANDCOVERTYPES]) {
  1442. // Added gross fraction transfer as percentage of the lesser of two stand types belonging to certain landcover types
  1443. double gross_lc_change_frac[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  1444. gross_lc_change_frac[CROPLAND][NATURAL] = 0.05;
  1445. gross_lc_change_frac[NATURAL][CROPLAND] = 0.05;
  1446. gross_lc_change_frac[CROPLAND][PASTURE] = 0.05;
  1447. gross_lc_change_frac[PASTURE][CROPLAND] = 0.05;
  1448. gross_lc_change_frac[PASTURE][NATURAL] = 0.05;
  1449. gross_lc_change_frac[NATURAL][PASTURE] = 0.05;
  1450. gross_lc_change_frac[FOREST][NATURAL] = 0.05;
  1451. gross_lc_change_frac[NATURAL][FOREST] = 0.05;
  1452. gross_lc_change_frac[FOREST][PASTURE] = 0.05;
  1453. gross_lc_change_frac[PASTURE][FOREST] = 0.05;
  1454. Landcover& lc = gridcell.landcover;
  1455. for(int from=0; from<NLANDCOVERTYPES; from++) {
  1456. for(int to=0; to<NLANDCOVERTYPES; to++) {
  1457. lc_frac_transfer[from][to] += gross_lc_change_frac[from][to] * min(lc.frac_old[from], lc.frac_old[to]);
  1458. }
  1459. }
  1460. }
  1461. /// Simulates gross landcover change by adding a specified fraction to the st_frac_transfer array
  1462. void simulate_gross_st_transfer(Gridcell& gridcell, double* st_frac_transfer) {
  1463. // Added gross fraction transfer as percentage of the lesser of two stand types belonging to certain landcover types
  1464. double gross_lc_change_frac[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  1465. gross_lc_change_frac[CROPLAND][NATURAL] = 0.05;
  1466. gross_lc_change_frac[NATURAL][CROPLAND] = 0.05;
  1467. gross_lc_change_frac[CROPLAND][PASTURE] = 0.05;
  1468. gross_lc_change_frac[PASTURE][CROPLAND] = 0.05;
  1469. gross_lc_change_frac[PASTURE][NATURAL] = 0.05;
  1470. gross_lc_change_frac[NATURAL][PASTURE] = 0.05;
  1471. gross_lc_change_frac[FOREST][NATURAL] = 0.05;
  1472. gross_lc_change_frac[NATURAL][FOREST] = 0.05;
  1473. gross_lc_change_frac[FOREST][PASTURE] = 0.05;
  1474. gross_lc_change_frac[PASTURE][FOREST] = 0.05;
  1475. for(int from=0; from<nst; from++) {
  1476. StandType& st_donor = stlist[from];
  1477. Gridcellst& gcst_donor = gridcell.st[from];
  1478. for(int to=0; to<nst; to++) {
  1479. StandType& st_receptor = stlist[to];
  1480. Gridcellst& gcst_receptor = gridcell.st[to];
  1481. if(gross_lc_change_frac[st_donor.landcover][st_receptor.landcover]) {
  1482. // All stand types with an area:
  1483. st_frac_transfer[index(from, to)] += gross_lc_change_frac[st_donor.landcover][st_receptor.landcover]
  1484. * min(min(gcst_donor.frac_old, gcst_donor.frac), min(gcst_receptor.frac_old, gcst_receptor.frac));
  1485. }
  1486. }
  1487. }
  1488. }
  1489. /// Called after update of st fraction and transfer values
  1490. bool check_fractions(Gridcell& gridcell, double landcoverfrac_change[], double lc_change_array[][NLANDCOVERTYPES],
  1491. double* st_change_array, bool check_lc_st_transfer = false) {
  1492. bool error = false;
  1493. // Check that landcover sum is 1:
  1494. double lc_frac_sum = 0.0;
  1495. for(int i=0; i<NLANDCOVERTYPES; i++)
  1496. lc_frac_sum += gridcell.landcover.frac[i];
  1497. if(!negligible(lc_frac_sum - 1.0, CHECK_NLIM)) {
  1498. dprintf("\nCheck 1: Year %d: landcover fraction sum: %.15f", date.get_calendar_year(), lc_frac_sum);
  1499. error = true;
  1500. }
  1501. /////////
  1502. // Check that stand type sum is 1:
  1503. double st_frac_sum = 0.0;
  1504. for(int i=0; i<nst; i++)
  1505. st_frac_sum += gridcell.st[i].frac;
  1506. if(fabs(st_frac_sum - 1.0) > INPUT_RESOLUTION * 10.0) {
  1507. dprintf("\nCheck 2: Year %d: stand type fraction sum: %.15f", date.year, st_frac_sum);
  1508. error = true;
  1509. }
  1510. /////////
  1511. // Test if the sum of gross lcc for a landcover is the same as net lcc:
  1512. double test_lc_change[NLANDCOVERTYPES] = {0.0};
  1513. for(int from=0; from<NLANDCOVERTYPES; from++) {
  1514. for(int to=0; to<NLANDCOVERTYPES; to++) {
  1515. test_lc_change[from] -= lc_change_array[from][to];
  1516. test_lc_change[to] += lc_change_array[from][to];
  1517. }
  1518. }
  1519. for(int i=0; i<NLANDCOVERTYPES; i++) {
  1520. if(fabs(test_lc_change[i] - landcoverfrac_change[i]) > INPUT_RESOLUTION * 10.0) {
  1521. dprintf("\nCheck 3: Year %d: lc_change_array sum not equal to landcoverfrac_change value for landcover %d\n", date.year, i);
  1522. dprintf("dif=%.15f", test_lc_change[i] - landcoverfrac_change[i]);
  1523. error = true;
  1524. }
  1525. }
  1526. ////////////
  1527. // Test if the sum of gross lcc for a stand type is the same as net lcc:
  1528. double *test_st_change;
  1529. test_st_change = new double[nst];
  1530. for(int i=0;i<nst;i++) {
  1531. test_st_change[i] = 0.0;
  1532. }
  1533. for(int from=0; from<nst; from++) {
  1534. for(int to=0; to<nst; to++) {
  1535. test_st_change[from] -= st_change_array[index(from, to)];
  1536. test_st_change[to] += st_change_array[index(from, to)];
  1537. }
  1538. }
  1539. for(int i=0; i<nst; i++) {
  1540. StandType& st = stlist[i];
  1541. Gridcellst& gcst = gridcell.st[i];
  1542. if(fabs(test_st_change[i] - gcst.frac_change) > INPUT_RESOLUTION * 10.0) {
  1543. dprintf("\nCheck 4: Year %d: st_change_array sum not equal to st.frac_change value for stand type %d\n", date.year, i);
  1544. dprintf("dif=%.15f", test_st_change[i] - gcst.frac_change);
  1545. error = true;
  1546. }
  1547. }
  1548. if(test_st_change)
  1549. delete[] test_st_change;
  1550. ////////////
  1551. // Test if st_change_array is compatible with st.frac_old/frac values
  1552. for(int from=0; from<nst; from++) {
  1553. StandType& st_donor = stlist[from];
  1554. Gridcellst& gcst_donor = gridcell.st[from];
  1555. for(int to=0; to<nst; to++) {
  1556. StandType& st_receptor = stlist[to];
  1557. Gridcellst& gcst_receptor = gridcell.st[to];
  1558. if(st_change_array[index(from, to)] > (gcst_donor.frac_old + INPUT_RESOLUTION) || st_change_array[index(from, to)] > (gcst_receptor.frac + INPUT_RESOLUTION)) {
  1559. dprintf("\nCheck 5: Year %d: st_change_array sum not compatible with st.frac_old/frac values for stand types %d and %d\n", date.year, from, to);
  1560. dprintf("st_change_array=%.15f, st_donor.frac_old=%.15f, st_receptor.frac=%.15f", st_change_array[index(from, to)], gcst_donor.frac_old, gcst_receptor.frac);
  1561. error = true;
  1562. }
  1563. }
  1564. }
  1565. // Test if stand type and landcover transfer matrices match each other:
  1566. if(check_lc_st_transfer) {
  1567. double lc_change[NLANDCOVERTYPES] = {0.0};
  1568. double lc_change_arr[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  1569. for(int from=0; from<nst; from++) {
  1570. StandType& st = stlist[from];
  1571. for(int to=0; to<nst; to++) {
  1572. StandType& st_dest = stlist[to];
  1573. lc_change[st.landcover] -= st_change_array[index(from, to)];
  1574. lc_change[st_dest.landcover] += st_change_array[index(from, to)];
  1575. lc_change_arr[st.landcover][st_dest.landcover] += st_change_array[index(from, to)];
  1576. }
  1577. }
  1578. for(int i=0; i<NLANDCOVERTYPES; i++) {
  1579. if(fabs(lc_change[i] - landcoverfrac_change[i]) > INPUT_RESOLUTION * 10.0) {
  1580. dprintf("\nCheck 6: Year %d: st_change_array LC sum not equal to LC value for %d\n", date.year, i);
  1581. dprintf("dif=%.15f", fabs(lc_change[i] - landcoverfrac_change[i]));
  1582. error = true;
  1583. }
  1584. }
  1585. for(int from=0; from<NLANDCOVERTYPES; from++) {
  1586. for(int to=0; to<NLANDCOVERTYPES; to++) {
  1587. if(fabs(lc_change_arr[from][to] - lc_change_array[from][to]) > INPUT_RESOLUTION * 10.0) {
  1588. dprintf("\nCheck 7: Year %d: lc_change_arr sum not equal to lc_change_array value for %d, %d\n", date.year, from, to);
  1589. dprintf("dif=%.15f", lc_change_arr[from][to] - lc_change_array[from][to]);
  1590. error = true;
  1591. }
  1592. }
  1593. }
  1594. // Test that no lc transfers are negative:
  1595. for(int from=0; from<NLANDCOVERTYPES; from++) {
  1596. for(int to=0; to<NLANDCOVERTYPES; to++) {
  1597. if(lc_change_array[from][to] < 0.0) {
  1598. dprintf("\nCheck 7b: Year %d: Negative lc_change_array[%d][%d]: %.15f\n", date.year, from, to, lc_change_array[from][to]);
  1599. error = true;
  1600. }
  1601. }
  1602. }
  1603. }
  1604. if(error)
  1605. dprintf("\n\n");
  1606. return error;
  1607. }
  1608. /// Called before updating stand.frac (reduce_stands)
  1609. bool check_fractions1(Gridcell& gridcell) {
  1610. bool error = false;
  1611. // Test if the stand type reduction is smaller than the remaining stand sum:
  1612. for(int s=0; s<nst; s++) {
  1613. StandType& st = stlist[s];
  1614. Gridcellst& gcst = gridcell.st[s];
  1615. if(gcst.frac_change < 0.0){
  1616. double stands_frac_sum = 0.0;
  1617. for(unsigned int i = 0; i < gridcell.nbr_stands(); i++) {
  1618. Stand& stand = gridcell[i];
  1619. if(stand.stid == st.id)
  1620. stands_frac_sum += stand.get_gridcell_fraction();
  1621. }
  1622. if(-gcst.frac_change - stands_frac_sum > INPUT_RESOLUTION * 10.0) {
  1623. dprintf("\nCheck 8: Year %d: stand type %d fraction reduction bigger than sum of stands\n", date.year, s);
  1624. dprintf("dif=%.15f", -gcst.frac_change - stands_frac_sum);
  1625. error = true;
  1626. }
  1627. }
  1628. }
  1629. if(error)
  1630. dprintf("\n\n");
  1631. return error;
  1632. }
  1633. /// Called after updating stand.frac and transfer values (reduce_stands)
  1634. bool check_fractions2(Gridcell& gridcell, double* st_change_array) {
  1635. bool error = false;
  1636. // Test if stand.transfer_area_st[to] sum for a stand type is equal to st_change_array[from][to)]
  1637. for(int from=0; from<nst; from++) {
  1638. StandType& st = stlist[from];
  1639. for(int to=0; to<nst; to++) {
  1640. double *test_st_change;
  1641. test_st_change = new double[nst];
  1642. for(int i=0;i<nst;i++)
  1643. test_st_change[i] = 0.0;
  1644. for(unsigned int i=0; i<gridcell.nbr_stands(); i++) {
  1645. Stand& stand = gridcell[i];
  1646. if(stand.stid == st.id)
  1647. test_st_change[to] += stand.transfer_area_st[to];
  1648. }
  1649. if(fabs(test_st_change[to] - st_change_array[index(from, to)]) > INPUT_RESOLUTION * 100.0) {
  1650. dprintf("\nCheck 9: Year %d: stand transfer area sum not equal to stand type value for stand types %d and %d\n", date.year, from, to);
  1651. dprintf("dif=%.15f", fabs(test_st_change[to] - st_change_array[index(from, to)]));
  1652. error = true;
  1653. }
  1654. if(test_st_change)
  1655. delete[] test_st_change;
  1656. }
  1657. }
  1658. ////////////
  1659. // Test if stand.frac_change for a stand is equal to stand.gross_frac_increase + stand.gross_frac_decrease:
  1660. for(unsigned int i=0; i<gridcell.nbr_stands(); i++) {
  1661. Stand& stand = gridcell[i];
  1662. if(largerthanzero(stand.frac_change - (stand.gross_frac_increase - stand.gross_frac_decrease), CHECK_NLIM+1)) {
  1663. dprintf("\nCheck 10: Year %d: frac_change is not equal to gross_frac_increase + gross_frac_decrease for stand %d\n", date.get_calendar_year(), stand.id);
  1664. dprintf("dif=%.15f\n", fabs(stand.frac_change - (stand.gross_frac_increase - stand.gross_frac_decrease)));
  1665. dprintf("frac_change=%.15f, gross_frac_increase=%.15f, gross_frac_decrease=%.15f", stand.frac_change, stand.gross_frac_increase, stand.gross_frac_decrease);
  1666. error = true;
  1667. }
  1668. }
  1669. ////////////
  1670. // Test that no stands remain when stand type fraction is 0 (after reduce_stands, so stands may remain, but with frac 0):
  1671. for(int s=0; s<nst; s++) {
  1672. StandType& st = stlist[s];
  1673. Gridcellst& gcst = gridcell.st[s];
  1674. if(gcst.frac == 0.0) {
  1675. double stands_frac_sum = 0.0;
  1676. for(unsigned int i = 0; i < gridcell.nbr_stands(); i++) {
  1677. Stand& stand = gridcell[i];
  1678. if(stand.stid == st.id)
  1679. stands_frac_sum += stand.get_gridcell_fraction();
  1680. }
  1681. if(stands_frac_sum) {
  1682. dprintf("\nCheck 11: Year %d: remaining stand when stand type %d fraction is 0\n", date.get_calendar_year(), s);
  1683. dprintf("remain=%.20f", stands_frac_sum);
  1684. error = true;
  1685. }
  1686. }
  1687. }
  1688. if(error)
  1689. dprintf("\n\n");
  1690. return error;
  1691. }
  1692. /// Called after updating stand.frac and transfer values of reduced stands (reduce_stands)
  1693. bool check_fractions3(Gridcell& gridcell) {
  1694. bool error = false;
  1695. // Test if sum of stands not equal to stand type value for stand type for reduced stands
  1696. for(int s=0; s<nst; s++) {
  1697. StandType& st = stlist[s];
  1698. Gridcellst& gcst = gridcell.st[s];
  1699. double stands_frac_sum = 0.0;
  1700. for(unsigned int i = 0; i < gridcell.nbr_stands(); i++) {
  1701. Stand& stand = gridcell[i];
  1702. if(stand.stid == st.id)
  1703. stands_frac_sum += stand.get_gridcell_fraction();
  1704. }
  1705. if(gcst.frac_change < 0.0) {
  1706. if(fabs(gcst.frac - stands_frac_sum - gcst.gross_frac_increase) > INPUT_RESOLUTION * 100.0) {
  1707. dprintf("\nCheck 12: Year %d: fraction sum of stands not equal to stand type value for stand type %d\n", date.year, s);
  1708. dprintf("dif=%.15f", fabs(gcst.frac - stands_frac_sum - gcst.gross_frac_increase));
  1709. error = true;
  1710. }
  1711. }
  1712. }
  1713. if(error)
  1714. dprintf("\n\n");
  1715. return error;
  1716. }
  1717. /// Called after creating new stands (stand_dynamics)
  1718. bool check_fractions4(Gridcell& gridcell) {
  1719. bool error = false;
  1720. // Test if sum of stands not equal to stand type value for stand type for increased or new stands
  1721. for(int s=0; s<nst; s++) {
  1722. StandType& st = stlist[s];
  1723. Gridcellst& gcst = gridcell.st[s];
  1724. double stands_frac_sum = 0.0;
  1725. for(unsigned int i = 0; i < gridcell.nbr_stands(); i++) {
  1726. Stand& stand = gridcell[i];
  1727. if(stand.stid == st.id) {
  1728. stands_frac_sum += stand.get_gridcell_fraction();
  1729. }
  1730. }
  1731. if(gcst.frac_change >= 0.0) {
  1732. if(gcst.frac && !stands_frac_sum) {
  1733. dprintf("\nCheck 13a: Year %d: no stand present when fraction is > 0 for stand type %d\n", date.year, s);
  1734. dprintf("st frac=%.15f\n\n", gcst.frac);
  1735. }
  1736. else if(fabs(gcst.frac - stands_frac_sum) > INPUT_RESOLUTION * 100.0) {
  1737. dprintf("\nCheck 13: Year %d: fraction sum of stands not equal to stand type value for stand type %d\n", date.year, s);
  1738. dprintf("dif=%.15f", fabs(gcst.frac - stands_frac_sum));
  1739. error = true;
  1740. }
  1741. }
  1742. }
  1743. if(error)
  1744. dprintf("\n\n");
  1745. return error;
  1746. }
  1747. /// Gets this year's landcover and crop area fractions, landcover transitions and checks that area changes are significant.
  1748. /** Stores changes in area fractions for the stand types
  1749. *
  1750. * OUTPUT PARAMETERS
  1751. * \param LCchangeCtransfer whether to transfer carbon, nitrogen and water of reduced stands to expanding stands
  1752. */
  1753. bool lc_changed(Gridcell& gridcell, bool& LCchangeCtransfer, InputModule* input_module) {
  1754. double stfrac_sum_old[NLANDCOVERTYPES] = {0.0};
  1755. double change_stand = 0.0;
  1756. double changeLC = 0.0;
  1757. double change_st = 0.0;
  1758. double transferred_fraction = 0.0;
  1759. double receiving_fraction = 0.0;
  1760. bool change = true;
  1761. bool force_small_change = false;
  1762. Landcover& lc = gridcell.landcover;
  1763. // Get new gridcell.landcoverfrac and/or standtype.frac from LUdata and CFTdata
  1764. input_module->getlandcover(gridcell);
  1765. // Check if any changes in land fractions have been made
  1766. for(unsigned int i = 0; i < gridcell.st.nobj; i++) {
  1767. Gridcellst& gcst = gridcell.st[i];
  1768. stfrac_sum_old[gcst.st.landcover] += gcst.frac_old;
  1769. if(fabs(gcst.frac_change) && (!gcst.frac || !gcst.frac_old))
  1770. force_small_change = true;
  1771. }
  1772. for(int i=0; i<NLANDCOVERTYPES; i++) {
  1773. changeLC += fabs(lc.frac_change[i]) / 2.0;
  1774. }
  1775. for(int lc=0; lc<NLANDCOVERTYPES; lc++) {
  1776. if(run[lc] && (!frac_fixed[lc] || !lcfrac_fixed)) {
  1777. for(unsigned int i = 0; i < gridcell.st.nobj; i++) {
  1778. Gridcellst& gcst = gridcell.st[i];
  1779. if(gcst.st.landcover == lc) {
  1780. double stfrac_change = gcst.frac_change;
  1781. if(stfrac_change < 0.0)
  1782. transferred_fraction -= stfrac_change;
  1783. else if(stfrac_change > 0.0)
  1784. receiving_fraction += stfrac_change;
  1785. if(stfrac_sum_old[lc] != 0.0) {
  1786. change_st += fabs(stfrac_change) / 2.0;
  1787. }
  1788. else {
  1789. change_st += fabs(stfrac_change);
  1790. }
  1791. change_stand += fabs(stfrac_change) / 2.0;
  1792. }
  1793. }
  1794. }
  1795. }
  1796. double change_gross_lcc = 0.0;
  1797. for(int from=0; from<NLANDCOVERTYPES; from++) {
  1798. if(run[from]) {
  1799. for(int to=0; to<NLANDCOVERTYPES; to++) {
  1800. if(run[to]) {
  1801. change_gross_lcc += lc.frac_transfer[to][from];
  1802. }
  1803. }
  1804. }
  1805. }
  1806. // if no changes, do nothing.
  1807. if(changeLC < INPUT_RESOLUTION * 0.1 && change_st < INPUT_RESOLUTION * 0.1 && change_gross_lcc < INPUT_RESOLUTION * 0.1 && !force_small_change) {
  1808. change = false;
  1809. }
  1810. // check for balance of reduced and increased stand fractions
  1811. else {
  1812. if(fabs(transferred_fraction - receiving_fraction) > 0.0001 || fabs(change_stand-receiving_fraction) > 0.0001) {
  1813. if(run[CROPLAND] && run[NATURAL] && run[PASTURE]) {
  1814. // end program if balance is expected (no landcovers inactivated)
  1815. FAIL("Transferred landcover fractions not balanced ! Abort. \n");
  1816. }
  1817. else {
  1818. // allow program to continue, but inactivate landcover change mass transfer
  1819. LCchangeCtransfer = false;
  1820. dprintf("Transferred landcover fractions not balanced !\nLandcover change fluxes not calculated.\n");
  1821. }
  1822. }
  1823. change = true;
  1824. }
  1825. return change;
  1826. }
  1827. /// Transfers harvested forest area to transition matrix
  1828. bool transfer_harvested_fractions(Gridcell& gridcell, double lc_frac_transfer[][NLANDCOVERTYPES], forest_lc_frac_transfer& forest_lc_frac_transfer_s) {
  1829. bool transfer_harvest = false;
  1830. if(gridcell.landcover.wood_harvest.prim_frac > 0.0) {
  1831. lc_frac_transfer[NATURAL][NATURAL] += gridcell.landcover.wood_harvest.prim_frac;
  1832. forest_lc_frac_transfer_s.primary[NATURAL][NATURAL] += gridcell.landcover.wood_harvest.prim_frac;
  1833. use_primary_lc_transfer = true;
  1834. transfer_harvest = true;
  1835. }
  1836. if(harvest_secondary_to_new_stand && gridcell.landcover.wood_harvest.sec_mature_frac + gridcell.landcover.wood_harvest.sec_young_frac > 0.0) {
  1837. lc_frac_transfer[NATURAL][NATURAL] += gridcell.landcover.wood_harvest.sec_mature_frac;
  1838. lc_frac_transfer[NATURAL][NATURAL] += gridcell.landcover.wood_harvest.sec_young_frac;
  1839. forest_lc_frac_transfer_s.secondary_young[NATURAL][NATURAL] += gridcell.landcover.wood_harvest.sec_young_frac;
  1840. transfer_harvest = true;
  1841. }
  1842. return transfer_harvest;
  1843. }
  1844. /// Updates all landcover, stand type and stand area fractions each year, possibly resulting in the creation and killing of stands.
  1845. /** Harvests transferred areas and transfers litter etc. of reduced stands to expanding stands and harvested matter to fluxes
  1846. * and (in the case of wood) to long-lived pools.
  1847. *
  1848. * Land cover fraction files need to be read together with gross land transition files. Net and gross land cover change is
  1849. * expected to be compatible, although rounding errors are handled by the input code.
  1850. *
  1851. * If instruction file parameter iftransfer_to_new_stand is true, new stands may be created for separate land transfers
  1852. * in transfer_to_new_stand(), according to rules in copy_stand_type().
  1853. *
  1854. * Otherwise, new transferred areas are pooled according to instruction file parameter transfer_level (0: one big pool; 1: land cover-level;
  1855. * 2: stand type-level) and rules in the gridcell arrays pool_from_all_landcovers[to] and pool_to_all_landcovers[from], set
  1856. * in the Gridcell constructor.
  1857. * New stands may then be created in stand_dynamics() from the pooled land if the gridcell array expand_to_new_stand[lc] value for the
  1858. * receiving land cover, set in the Gridcell constructor, is true (natural and forest).
  1859. * If not, soil and litter carbon and nitrogen as well as water is pooled with the receiving stands (cropland, pasture).
  1860. *
  1861. * Transferred land with living plant mass after land cover change should be put in new stands, using the stand cloning mode in
  1862. * transfer_to_new_stand(). Stands with living plant mass should not be pooled with newly transferred land, unconditionally so
  1863. * if they contain trees (pasture can be expanded without too much problems). New stands should instead be created, either in stand_dynamics()
  1864. * or transfer_to_new_stand() as described above.
  1865. */
  1866. void landcover_dynamics(Gridcell& gridcell, InputModule* input_module) {
  1867. // ecev3
  1868. int year = date.get_calendar_year();
  1869. if (false && ECEARTH && gridcell.fixedLUafter >= 0 && year >= gridcell.fixedLUafter ) {
  1870. dprintf("Return from landcover_dynamics: fixedLUfter = %d \n", year);
  1871. return;
  1872. }
  1873. double* st_frac_transfer = NULL;
  1874. forest_st_frac_transfer forest_st_frac_transfer_s(nst);
  1875. bool LCchangeCtransfer = true;
  1876. Landcover& lc = gridcell.landcover;
  1877. st_frac_transfer = new double[nst * nst];
  1878. for (int i = 0; i<nst; i++) {
  1879. gridcell.st[i].gross_frac_increase = 0.0;
  1880. gridcell.st[i].gross_frac_decrease = 0.0;
  1881. }
  1882. for(int i=0;i<NLANDCOVERTYPES;i++) {
  1883. lc.frac_change[i] = 0.0;
  1884. for(int j=0;j<NLANDCOVERTYPES;j++) {
  1885. lc.frac_transfer[i][j] = 0.0;
  1886. lc.forest_lc_frac_transfer_s.primary[i][j] = 0.0;
  1887. lc.forest_lc_frac_transfer_s.secondary_young[i][j] = 0.0;
  1888. }
  1889. }
  1890. for(int i=0;i<nst*nst;i++) {
  1891. st_frac_transfer[i] = 0.0;
  1892. forest_st_frac_transfer_s.primary[i] = 0.0;
  1893. forest_st_frac_transfer_s.secondary_young[i] = 0.0;
  1894. }
  1895. gridcell.landcover.updated = false;
  1896. for(unsigned int i=0; i<gridcell.nbr_stands(); ++i)
  1897. gridcell[i].scale_LC_change = 1.0;
  1898. bool no_changes = true;
  1899. // Transfer area to new stands during wood harvest.
  1900. if(transfer_harvested_fractions(gridcell, lc.frac_transfer, lc.forest_lc_frac_transfer_s))
  1901. no_changes = false;
  1902. // Rescale stand fractions so sum is 1
  1903. stlist.firstobj();
  1904. while (stlist.isobj) {
  1905. StandType& st = stlist.getobj();
  1906. Gridcellst& gcst = gridcell.st[st.id];
  1907. double frac_sum = 0.0;
  1908. for(unsigned int s=0;s<gridcell.nbr_stands();s++) {
  1909. Stand& stand = gridcell[s];
  1910. if(stand.stid == st.id)
  1911. frac_sum += stand.get_gridcell_fraction();
  1912. }
  1913. for(unsigned int s=0;s<gridcell.nbr_stands();s++) {
  1914. Stand& stand = gridcell[s];
  1915. if(stand.stid == st.id && frac_sum && gcst.frac)
  1916. stand.set_gridcell_fraction(stand.get_gridcell_fraction() / (frac_sum / gcst.frac));
  1917. }
  1918. stlist.nextobj();
  1919. }
  1920. // Get new landcover and stand type area fractions from input files, read transition arrays.
  1921. if(!all_fracs_const) {
  1922. // this call returns 0, causing this function to return, if no significant landcover changes this year,
  1923. // sets LCchangeCtransfer to 0 if unbalanced landcover changes (if some landcovers are inactivated), thus inactivating transfer of C and N
  1924. if(lc_changed(gridcell, LCchangeCtransfer, input_module)) {
  1925. no_changes = false;
  1926. }
  1927. }
  1928. //CLN LUCHANGE HERE
  1929. if(no_changes) {
  1930. delete[] st_frac_transfer;
  1931. return;
  1932. }
  1933. // Transfer all landcover changes to stand type transition matrix, if not already done.
  1934. if(gross_land_transfer == 3) {
  1935. // Option to read stand type transitions from file.
  1936. fail("Currently no code for option gross_land_transfer==3\n");
  1937. }
  1938. else if(gross_land_transfer == 2 && gross_input_present) {
  1939. // Option to read landcover transitions from file. Update the st_frac_transfer array.
  1940. set_st_change_array(gridcell, lc.frac_transfer, st_frac_transfer, lc.forest_lc_frac_transfer_s, forest_st_frac_transfer_s);
  1941. if(check_fractions(gridcell, lc.frac_change, lc.frac_transfer, st_frac_transfer, true))
  1942. dprintf("Fraction error after set_st_change_array()\n\n");
  1943. }
  1944. else {
  1945. // Option not to read landcover or stand type transitions or to simulate these.
  1946. // The st_frac_transfer array always needs to be updated.
  1947. const bool simulate_st = true; // gross lcc simulation at land cover level (false) or stand type level (true)
  1948. set_lc_change_array(lc.frac_change, lc.frac_transfer);
  1949. if(gross_land_transfer && !simulate_st)
  1950. simulate_gross_lc_transfer(gridcell, lc.frac_transfer);
  1951. // The lc_frac_transfer-array may contain overshoots if wood harvest with fraction transfer occurred.
  1952. double dummy;
  1953. adjust_gross_transfers(gridcell, lc.frac_change, lc.frac_transfer, lc.forest_lc_frac_transfer_s, dummy);
  1954. set_st_change_array(gridcell, lc.frac_transfer, st_frac_transfer, lc.forest_lc_frac_transfer_s, forest_st_frac_transfer_s);
  1955. if(check_fractions(gridcell, lc.frac_change, lc.frac_transfer, st_frac_transfer, true))
  1956. dprintf("Fraction error after set_st_change_array()\n\n");
  1957. if(gross_land_transfer && simulate_st)
  1958. simulate_gross_st_transfer(gridcell, st_frac_transfer);
  1959. }
  1960. check_fractions(gridcell, lc.frac_change, lc.frac_transfer, st_frac_transfer);
  1961. check_fractions1(gridcell);
  1962. for(int from=0; from<nst; from++) {
  1963. for(int to=0; to<nst; to++) {
  1964. if(st_frac_transfer[index(from, to)] > 0.0) {
  1965. gridcell.st[to].gross_frac_increase += st_frac_transfer[index(from, to)];
  1966. gridcell.st[from].gross_frac_decrease += st_frac_transfer[index(from, to)];
  1967. }
  1968. }
  1969. }
  1970. double ccont_tot_1 = gridcell.ccont();
  1971. double cflux_tot_1 = gridcell.cflux();
  1972. double ncont_tot_1 = gridcell.ncont();
  1973. double nflux_tot_1 = gridcell.nflux();
  1974. // check how many stands of each stand type exist
  1975. // identify which stands to reduce in area
  1976. // set stand variables frac, frac_old, frac_change and gross_frac_decrease for reduced stands and stand types
  1977. reduce_stands(gridcell, st_frac_transfer, forest_st_frac_transfer_s);
  1978. int error = check_fractions2(gridcell, st_frac_transfer);
  1979. error += check_fractions3(gridcell);
  1980. if(error) {
  1981. FAIL("Fraction error after reduce_stands()\n\n");
  1982. }
  1983. // Create new stands for land cover transitions when natural vegetation remains or when each transition requires a new stand:
  1984. if(iftransfer_to_new_stand) {
  1985. bool new_stand = false;
  1986. for(int from=0; from<nst; from++) {
  1987. for(int to=0; to<nst; to++) {
  1988. if(st_frac_transfer[index(from, to)] > 0.0) {
  1989. double new_stand_frac = transfer_to_new_stand(gridcell, from, to);
  1990. if(new_stand_frac) {
  1991. st_frac_transfer[index(from, to)] -= new_stand_frac;
  1992. if(st_frac_transfer[index(from, to)] < INPUT_RESOLUTION * 0.1)
  1993. st_frac_transfer[index(from, to)] = 0.0;
  1994. new_stand = true;
  1995. }
  1996. }
  1997. }
  1998. }
  1999. if(new_stand) {
  2000. int error = check_fractions(gridcell, lc.frac_change, lc.frac_transfer, st_frac_transfer);
  2001. error += check_fractions2(gridcell, st_frac_transfer);
  2002. if(error)
  2003. dprintf("Fraction error after transfer_to_new_stand()\n\n");
  2004. }
  2005. }
  2006. // set stand variables frac, frac_change and gross_frac_increase for expanding stands and stand types
  2007. expand_stands(gridcell, st_frac_transfer);
  2008. error += check_fractions(gridcell, lc.frac_change, lc.frac_transfer, st_frac_transfer);
  2009. error += check_fractions2(gridcell, st_frac_transfer);
  2010. if(error)
  2011. FAIL("Fraction error after expand_stands()\n");
  2012. // Pooling options
  2013. // transfer_level: 0: one big pool; 1: land cover-level; 2: stand type-level
  2014. if(transfer_level == 0) { // One big pool for all transfers (as in old code)
  2015. landcover_change_transfer transfer;
  2016. double receiving_fraction = 0.0;
  2017. for(int from=0; from<nst; from++) {
  2018. for(int to=0; to<nst; to++) {
  2019. if(st_frac_transfer[index(from, to)] > 0.0)
  2020. receiving_fraction += st_frac_transfer[index(from, to)];
  2021. }
  2022. }
  2023. // handle harvest and turnover of reduced stands at landcover change
  2024. donor_stand_change(gridcell, receiving_fraction, transfer);
  2025. // create and kill stands at landcover change
  2026. stand_dynamics(gridcell);
  2027. error += check_fractions4(gridcell);
  2028. if(error)
  2029. fail("Fraction error after stand_dynamics()\n");
  2030. // transfer litter etc. of reduced stands to expanding stands at landcover change
  2031. receiving_stand_change(gridcell, transfer, LCchangeCtransfer);
  2032. }
  2033. else if(transfer_level == 1) { // Landcover-level pools, larger pools available by the gridcell variables pool_from_all_landcovers and pool_to_all_landcovers (set in gridcell constructor)
  2034. landcover_change_transfer transfer_lc_2d[NLANDCOVERTYPES][NLANDCOVERTYPES];
  2035. landcover_change_transfer transfer_lc[NLANDCOVERTYPES];
  2036. landcover_change_transfer transfer_lc_from[NLANDCOVERTYPES];
  2037. double gross_landcoverfrac_increase[NLANDCOVERTYPES] = {0.0};
  2038. double gross_landcoverfrac_decrease[NLANDCOVERTYPES] = {0.0};
  2039. double gross_landcoverfrac_transfer[NLANDCOVERTYPES][NLANDCOVERTYPES] = {0.0};
  2040. for(int from=0; from<nst; from++) {
  2041. for(int to=0; to<nst; to++) {
  2042. if(st_frac_transfer[index(from, to)] > 0.0) {
  2043. gross_landcoverfrac_increase[stlist[to].landcover] += st_frac_transfer[index(from, to)];
  2044. gross_landcoverfrac_decrease[stlist[from].landcover] += st_frac_transfer[index(from, to)];
  2045. gross_landcoverfrac_transfer[stlist[from].landcover][stlist[to].landcover] += st_frac_transfer[index(from, to)];
  2046. }
  2047. }
  2048. }
  2049. // handle harvest and turnover of reduced stands at landcover change
  2050. // pool all transferred area from one landcover:
  2051. for(int from=0; from<NLANDCOVERTYPES; from++) {
  2052. if(gridcell.landcover.pool_to_all_landcovers[from]) { // alt.c
  2053. if(gross_landcoverfrac_decrease[from] > 0.0) {
  2054. donor_stand_change(gridcell, gross_landcoverfrac_decrease[from], transfer_lc_from[from], -1, from);
  2055. }
  2056. }
  2057. }
  2058. for(int to=0; to<NLANDCOVERTYPES; to++) {
  2059. double receiving_fraction = gross_landcoverfrac_increase[to];
  2060. for(int from=0; from<NLANDCOVERTYPES; from++) {
  2061. double receiving_fraction_2d = gross_landcoverfrac_transfer[from][to];
  2062. if(receiving_fraction_2d > 0.0) {
  2063. if(gridcell.landcover.pool_from_all_landcovers[to]) {
  2064. if(!gridcell.landcover.pool_to_all_landcovers[from]) { // alt.a
  2065. // Add from->to transfer to to-pool:
  2066. donor_stand_change(gridcell, receiving_fraction, transfer_lc[to], to, from);
  2067. }
  2068. else { // alt.a+c
  2069. // Add from-pool value to to-pool:
  2070. transfer_lc[to].add(transfer_lc_from[from], receiving_fraction_2d / receiving_fraction);
  2071. }
  2072. }
  2073. else {
  2074. if(!gridcell.landcover.pool_to_all_landcovers[from]) { // alt.b
  2075. // Add from->to transfer to [from][to] place in 2d-array:
  2076. donor_stand_change(gridcell, receiving_fraction_2d, transfer_lc_2d[from][to], to, from);
  2077. }
  2078. else { // alt.c
  2079. // Copy from-pool value to [from][to] place in 2d-array:
  2080. transfer_lc_2d[from][to].copy(transfer_lc_from[from]);
  2081. }
  2082. }
  2083. }
  2084. }
  2085. }
  2086. // create and kill stands at landcover change
  2087. stand_dynamics(gridcell);
  2088. error += check_fractions4(gridcell);
  2089. if(error)
  2090. fail("Fraction error after stand_dynamics()\n");
  2091. // transfer litter etc. of reduced stands to expanding stands at landcover change
  2092. for(int to=0; to<NLANDCOVERTYPES; to++) {
  2093. if(gridcell.landcover.pool_from_all_landcovers[to]) {
  2094. // Alt.a: All donor lc:s are pooled into receiving lc:s
  2095. if(gross_landcoverfrac_increase[to] > 0.0) {
  2096. receiving_stand_change(gridcell, transfer_lc[to], LCchangeCtransfer, to);
  2097. }
  2098. }
  2099. else {
  2100. for(int from=0; from<NLANDCOVERTYPES; from++) {
  2101. if(gross_landcoverfrac_transfer[from][to] > 0.0) {
  2102. if(!gridcell.landcover.pool_to_all_landcovers[from]) {
  2103. // Alt.b: Transfers from donors are independent:
  2104. receiving_stand_change(gridcell, transfer_lc_2d[from][to], LCchangeCtransfer, to, -1, gross_landcoverfrac_transfer[from][to] / gross_landcoverfrac_increase[to]);
  2105. }
  2106. else {
  2107. // Alt.c: Donor lc:s are pooled before transfer to recipients:
  2108. receiving_stand_change(gridcell, transfer_lc_from[from], LCchangeCtransfer, to, -1, gross_landcoverfrac_transfer[from][to] / gross_landcoverfrac_increase[to]);
  2109. }
  2110. }
  2111. }
  2112. }
  2113. }
  2114. }
  2115. else if(transfer_level == 2) { // Unique transfers between stand types. Stand type-level pools available by the gridcell variables pool_from_all_landcovers and pool_to_all_landcovers (set in gridcell constructor)
  2116. landcover_change_transfer* transfer_st_2d;
  2117. landcover_change_transfer* transfer_st;
  2118. landcover_change_transfer* transfer_st_from;
  2119. transfer_st_2d = new landcover_change_transfer[nst * nst];
  2120. transfer_st = new landcover_change_transfer[nst];
  2121. transfer_st_from = new landcover_change_transfer[nst];
  2122. // handle harvest and turnover of reduced stands at landcover change
  2123. // pool all transferred area from one landcover:
  2124. for(int from=0; from<nst; from++) {
  2125. StandType& st = stlist[from];
  2126. Gridcellst& gcst = gridcell.st[from];
  2127. // Pooling of donor stands within a stand type
  2128. if(gridcell.landcover.pool_to_all_landcovers[st.landcover] || !(st.landcover == NATURAL || st.landcover == FOREST)) { // alt.c
  2129. if(gcst.gross_frac_decrease > 0.0) {
  2130. donor_stand_change(gridcell, gcst.gross_frac_decrease, transfer_st_from[from], -1, -1, -1, from);
  2131. }
  2132. }
  2133. }
  2134. for(int to=0;to<nst;to++) {
  2135. StandType& st_to = stlist[to];
  2136. Gridcellst& gcst_to = gridcell.st[to];
  2137. double receiving_fraction = gcst_to.gross_frac_increase;
  2138. for(int from=0; from<nst; from++) {
  2139. StandType& st_from = stlist[from];
  2140. Gridcellst& gcst_from = gridcell.st[from];
  2141. double receiving_fraction_2d = st_frac_transfer[index(from, to)];
  2142. if(receiving_fraction_2d > 0.0) {
  2143. if(gridcell.landcover.pool_from_all_landcovers[st_to.landcover]) {
  2144. if(!(gridcell.landcover.pool_to_all_landcovers[st_from.landcover] || !(st_from.landcover == NATURAL || st_from.landcover == FOREST))) { // alt.a
  2145. // Add from->to transfer to to-pool:
  2146. donor_stand_change(gridcell, receiving_fraction, transfer_st[to], -1, -1, to, from);
  2147. }
  2148. else { // alt.a+c
  2149. // Add from-pool value to to-pool:
  2150. transfer_st[to].add(transfer_st_from[from], receiving_fraction_2d / receiving_fraction);
  2151. }
  2152. }
  2153. else {
  2154. if(!(gridcell.landcover.pool_to_all_landcovers[st_from.landcover] || !(st_from.landcover == NATURAL || st_from.landcover == FOREST))) { // alt.b
  2155. // Add from->to transfer to [from][to] place in 2d-array:
  2156. donor_stand_change(gridcell, receiving_fraction_2d, transfer_st_2d[index(from, to)], -1, -1, to, from);
  2157. }
  2158. }
  2159. }
  2160. }
  2161. }
  2162. // create and kill stands at landcover change
  2163. stand_dynamics(gridcell);
  2164. error += check_fractions4(gridcell);
  2165. if(error)
  2166. fail("Fraction error after stand_dynamics()\n");
  2167. // transfer litter etc. of reduced stands to expanding stands at landcover change
  2168. for(int to=0; to<nst; to++) {
  2169. StandType& st = stlist[to];
  2170. Gridcellst& gcst = gridcell.st[to];
  2171. if(gridcell.landcover.pool_from_all_landcovers[st.landcover]) {
  2172. // Alt.a: All donor st:s are pooled into receiving st:s
  2173. if(gcst.gross_frac_increase > 0.0) {
  2174. receiving_stand_change(gridcell, transfer_st[to], LCchangeCtransfer, -1, to);
  2175. }
  2176. }
  2177. else {
  2178. for(int from=0; from<nst; from++) {
  2179. if(st_frac_transfer[index(from, to)] > 0.0) {
  2180. StandType& st_from = stlist[from];
  2181. Gridcellst& gcst_from = gridcell.st[from];
  2182. if(!(gridcell.landcover.pool_to_all_landcovers[st_from.landcover] || !(st_from.landcover == NATURAL || st_from.landcover == FOREST))) {
  2183. // Alt.b: Transfers between donor and receptor st:s are independent:
  2184. receiving_stand_change(gridcell, transfer_st_2d[index(from, to)], LCchangeCtransfer, st.landcover, to, st_frac_transfer[index(from, to)] / gcst.gross_frac_increase);
  2185. }
  2186. else {
  2187. // Alt.c: Donor stands in a stand type are pooled before transfer to recipients:
  2188. receiving_stand_change(gridcell, transfer_st_from[from], LCchangeCtransfer, st.landcover, to, st_frac_transfer[index(from, to)] / gcst.gross_frac_increase);
  2189. }
  2190. }
  2191. }
  2192. }
  2193. }
  2194. if(transfer_st_2d)
  2195. delete[] transfer_st_2d;
  2196. if(transfer_st)
  2197. delete[] transfer_st;
  2198. if(transfer_st_from)
  2199. delete[] transfer_st_from;
  2200. }
  2201. double ccont_tot = 0.0;
  2202. double cflux_tot = gridcell.cflux();
  2203. for(unsigned int i=0; i<gridcell.size(); i++) {
  2204. Stand& stand = gridcell[i];
  2205. double ccont_stand = stand.ccont(stand.scale_LC_change);
  2206. ccont_tot += ccont_stand * stand.get_gridcell_fraction();
  2207. }
  2208. if(!negligible(cflux_tot - cflux_tot_1 + ccont_tot - ccont_tot_1, -11))
  2209. dprintf("WARNING ! C balance after lcc off by %.12f\n",cflux_tot - cflux_tot_1 + ccont_tot - ccont_tot_1);
  2210. double ncont_tot = 0.0;
  2211. double nflux_tot = 0.0;
  2212. for(unsigned int i=0; i<gridcell.size(); i++) {
  2213. Stand& stand = gridcell[i];
  2214. double ncont_stand = stand.ncont(stand.scale_LC_change);
  2215. ncont_tot += ncont_stand * stand.get_gridcell_fraction();
  2216. nflux_tot += stand.nflux() * stand.get_gridcell_fraction();
  2217. }
  2218. nflux_tot_1 = nflux_tot; // Disregard fluxes not already reset this year and the associated scaling problems
  2219. nflux_tot += gridcell.landcover.anflux_landuse_change + gridcell.landcover.anflux_harvest_slow;
  2220. // if(!negligible(nflux_tot - nflux_tot_1 + ncont_tot - ncont_tot_1, CHECK_NLIM+1))
  2221. // dprintf("WARNING ! N balance after lcc off\n");
  2222. if(st_frac_transfer)
  2223. delete[] st_frac_transfer;
  2224. }
  2225. ////////////////////////////////////////////////////////////////////////////////
  2226. // Implementation of landcover_change_transfer member functions
  2227. ////////////////////////////////////////////////////////////////////////////////
  2228. /// landcover_change_transfer constructor
  2229. landcover_change_transfer::landcover_change_transfer() {
  2230. transfer_litter_leaf = transfer_litter_sap = transfer_litter_heart = NULL;
  2231. transfer_litter_root = transfer_litter_repr = transfer_harvested_products_slow = NULL;
  2232. transfer_nmass_litter_leaf = transfer_nmass_litter_sap = transfer_nmass_litter_heart = NULL;
  2233. transfer_nmass_litter_root = transfer_harvested_products_slow_nmass = NULL;
  2234. transfer_acflux_harvest = transfer_anflux_harvest = transfer_cpool_fast = 0.0;
  2235. transfer_cpool_slow = transfer_wcont_evap = transfer_decomp_litter_mean = 0.0;
  2236. transfer_k_soilfast_mean = transfer_k_soilslow_mean = transfer_nmass_avail = 0.0;
  2237. transfer_snowpack = transfer_snowpack_nmass = transfer_anfix_calc = 0.0;
  2238. for(int i=0;i<NSOILLAYER;i++)
  2239. transfer_wcont[i] = 0.0;
  2240. for(int i=0; i<NSOMPOOL; i++)
  2241. transfer_sompool[i].ntoc = 0.0;
  2242. allocate();
  2243. }
  2244. /// landcover_change_transfer deconstructor
  2245. landcover_change_transfer::~landcover_change_transfer() {
  2246. if(transfer_litter_leaf) delete[] transfer_litter_leaf;
  2247. if(transfer_litter_sap) delete[] transfer_litter_sap;
  2248. if(transfer_litter_heart) delete[] transfer_litter_heart;
  2249. if(transfer_litter_root) delete[] transfer_litter_root;
  2250. if(transfer_litter_repr) delete[] transfer_litter_repr;
  2251. if(transfer_harvested_products_slow) delete[] transfer_harvested_products_slow;
  2252. if(transfer_nmass_litter_leaf) delete[] transfer_nmass_litter_leaf;
  2253. if(transfer_nmass_litter_sap) delete[] transfer_nmass_litter_sap;
  2254. if(transfer_nmass_litter_heart) delete[] transfer_nmass_litter_heart;
  2255. if(transfer_nmass_litter_root) delete[] transfer_nmass_litter_root;
  2256. if(transfer_harvested_products_slow_nmass) delete[] transfer_harvested_products_slow_nmass;
  2257. }
  2258. /// allocates memory for landcover_change_transfer object
  2259. void landcover_change_transfer::allocate() {
  2260. transfer_litter_leaf = new double[npft];
  2261. transfer_litter_sap = new double[npft];
  2262. transfer_litter_heart = new double[npft];
  2263. transfer_litter_root = new double[npft];
  2264. transfer_litter_repr = new double[npft];
  2265. transfer_harvested_products_slow = new double[npft];
  2266. transfer_nmass_litter_leaf = new double[npft];
  2267. transfer_nmass_litter_sap = new double[npft];
  2268. transfer_nmass_litter_heart = new double[npft];
  2269. transfer_nmass_litter_root = new double[npft];
  2270. transfer_harvested_products_slow_nmass = new double[npft];
  2271. for(int i=0;i<npft;i++) {
  2272. transfer_litter_leaf[i] = 0.0;
  2273. transfer_litter_sap[i] = 0.0;
  2274. transfer_litter_heart[i] = 0.0;
  2275. transfer_litter_root[i] = 0.0;
  2276. transfer_litter_repr[i] = 0.0;
  2277. transfer_harvested_products_slow[i] = 0.0;
  2278. transfer_nmass_litter_leaf[i] = 0.0;
  2279. transfer_nmass_litter_sap[i] = 0.0;
  2280. transfer_nmass_litter_heart[i] = 0.0;
  2281. transfer_nmass_litter_root[i] = 0.0;
  2282. transfer_harvested_products_slow_nmass[i] = 0.0;
  2283. }
  2284. }
  2285. //////////////////////////////////////////////////////////////////////////////////////////
  2286. // REFERENCES
  2287. //
  2288. // Bondeau A, Smith PC, Zaehle S, Schaphoff S, Lucht W, Cramer W, Gerten D, Lotze-Campen H,
  2289. // Müller C, Reichstein M & Smith B 2007. Modelling the role of agriculture for the
  2290. // 20th century global terrestrial carbon balance. Global Change Biology, 13:679-706.
  2291. // Lindeskog M, Arneth A, Bondeau A, Waha K, Seaquist J, Olin S, & Smith B 2013.
  2292. // Implications of accounting for land use in simulations of ecosystem carbon cycling
  2293. // in Africa. Earth Syst Dynam 4:385-407.