= Daniel's page = CNP-Dev version based on MERGE-OCN revision 2698 = Bugs fixed = == '''1. Mass conservation issue: stomate_growth_fun_all.f90''' == The carbon being allocated to biomass pools must no be substracted before nutrient limitation of allocation is computed. This can be fixed by following: {{{ !DSGdebug_01 ! do NOT this now, do it after nutrient limitation on allocation is considered in bm_alloc_tot(ipts,j) ! ! Update the labile carbon pool ! biomass(ipts,j,ilabile,icarbon) = biomass(ipts,j,ilabile,icarbon) - & ! bm_alloc_tot(ipts,j) !DSGdebug_01 !! 3.10 Maintenance respiration }}} {{{ ! The calculation of ::resp_maint is solely based on the demand i.e. ! given the biomass and the condition of the plant, how much should be ! respired. It is not sure that this demand can be satisfied i.e. the ! calculated maintenance respiration may exceed the available carbon !DSGdebug_01 ! DEFAULT CASE: There is no deficit which must be subtracted from labile deficit = zero !DSGdebug_01 IF ( bm_alloc_tot(ipts,j) - resp_maint(ipts,j) .LT. zero ) THEN }}} [...] {{{ ! Not enough carbon to pay the deficit, the individual ! is going to die at the end of this day bm_alloc_tot(ipts,j) = bm_alloc_tot(ipts,j) + & biomass(ipts,j,icarbres,icarbon) biomass(ipts,j,icarbres,icarbon) = zero ! Truncate the maintenance respiration to the available carbon resp_maint(ipts,j) = bm_alloc_tot(ipts,j) !DSGdebug_01 ! There is no deficit which must be subtracted from labile deficit = zero !DSGdebug_01 ENDIF }}} [...] {{{ ! Final ::resp_maint is know bm_alloc_tot(ipts,j) = bm_alloc_tot(ipts,j) - resp_maint(ipts,j) !DSGdebug_01 ! Subtracted the deficit from labile pool biomass(ipts,j,ilabile,icarbon) = biomass(ipts,j,ilabile,icarbon) - & (resp_maint(ipts,j) + deficit) !DSGdebug_01 !! 3.11 Growth respiration ! Calculate total growth respiration and update allocatable carbon ! Growth respiration is a tax on productivity, not actual allocation ! Total growth respiration has be calculated before the allocation ! takes place because the allocation itself is not linear. After ! the allocation has been calculated, growth respiration can be ! calculated for each biomass component separatly. The unit of ! resp_growth is gC m-2 dt-1 resp_growth(ipts,j) = frac_growthresp(j) * MAX(zero, bm_alloc_tot(ipts,j)) bm_alloc_tot(ipts,j) = bm_alloc_tot(ipts,j) - resp_growth(ipts,j) !DSGdebug_01 biomass(ipts,j,ilabile,icarbon) = biomass(ipts,j,ilabile,icarbon) - & resp_growth(ipts,j) !DSGdebug_01 }}} [...] {{{ !======================================================= ! Block from OCN but I did not find it in DOFOCO ??? ! ! 5.1 retrieve allocated biomass from labile pool (nitrogen, or new allocation) ! !DSGdebug_01 ! This is the right spot to remove bm_alloc_tot: biomass(:,:,ilabile,icarbon) = biomass(:,:,ilabile,icarbon) - bm_alloc_tot(:,:) !DSGdebug_01 biomass(:,:,ilabile,initrogen) = biomass(:,:,ilabile,initrogen) - n_alloc_tot(:,:) }}} === UPDATE [01/06/2015]=== {{{ !! 5.2.7 Don't grow wood, use C to fill labile pool ELSEIF ( (.NOT. grow_wood) .AND. (b_inc_tot .GT. min_stomate) ) THEN ! Calculate the C that needs to be distributed to the ! labile pool. The fraction is proportional to the ratio ! between the total allocatable biomass and the unallocated ! biomass per tree (b_inc now contains the unallocated ! biomass). At the end of the allocation scheme bm_alloc_tot ! is substracted from the labile biomass pool to update the ! biomass pool (biomass(:,:,ilabile) = biomass(:,:,ilabile) - ! bm_alloc_tot(:,:)). At that point, the scheme puts the ! unallocated b_inc into the labile pool. What we ! want is that the unallocated fraction is removed from ! ::bm_alloc_tot such that only the allocated C is removed ! from the labile pool. b_inc_tot will be moved back into ! the labile pool in 5.2.11 bm_alloc_tot(ipts,j) = bm_alloc_tot(ipts,j) - b_inc_tot !DSGdebug_01a We didn't remove bm_alloc_tot yet from labile so no !need to correct labile pool for changes in bm_alloc_tot (which is residual) !DSGdebug_01a biomass(ipts,j,ilabile,icarbon) = & !DSGdebug_01a biomass(ipts,j,ilabile,icarbon) + b_inc_tot }}} [...] {{{ !! 5.3.7 Don't grow wood, use C to fill labile pool ELSEIF ( (.NOT. grow_wood) .AND. (b_inc_tot .GT. min_stomate) ) THEN ! Calculate the C that needs to be distributed to the ! labile pool. The fraction is proportional to the ratio ! between the total allocatable biomass and the unallocated ! biomass per tree (b_inc now contains the unallocated ! biomass). At the end of the allocation scheme bm_alloc_tot ! is substracted from the labile biomass pool to update the ! biomass pool (biomass(:,:,ilabile) = biomass(:,:,ilabile) - ! bm_alloc_tot(:,:)). At that point, the scheme puts the ! unallocated b_inc into the labile pool. What we ! want is that the unallocated fraction is removed from ! ::bm_alloc_tot such that only the allocated C is removed ! from the labile pool. b_inc_tot will be moved back into ! the labile pool in 5.2.11 bm_alloc_tot(ipts,j) = bm_alloc_tot(ipts,j) - b_inc_tot !DSGdebug_01a We didn't remove bm_alloc_tot yet from labile so no !need to correct labile pool for changes in bm_alloc_tot (which is residual) !DSGdebug_01a biomass(ipts,j,ilabile,icarbon) = biomass(ipts,j,ilabile,icarbon) + & !DSGdebug_01a b_inc_tot }}} [...] {{{ ELSE ! Not enough carbon to pay the deficit ! There is likely a bigger problem somewhere in ! this routine WRITE(numout,*) 'WARNING 23: PFT, ipts: ',j,ipts CALL ipslerr_p (3,'growth_fun_all',& 'WARNING 23: numerical problem overspending ',& 'when trying to account for unallocatable C ','') ENDIF ELSE !DSGdebug_01a We didn't remove bm_alloc_tot yet from labile so no !need to correct labile pool for changes in bm_alloc_tot (which is residual) !DSGdebug_01a ! Move the unallocated carbon back into the labile pool !DSGdebug_01a biomass(ipts,j,ilabile,icarbon) = & !DSGdebug_01a biomass(ipts,j,ilabile,icarbon) + residual(ipts,j) ENDIF }}} === UPDATE [03/06/2015]=== {{{ ! Not enough carbon to pay the deficit ! There is likely a bigger problem somewhere in ! this routine WRITE(numout,*) 'WARNING 11: PFT, ipts: ',j,ipts CALL ipslerr_p (3,'growth_fun_all',& 'WARNING 11: numerical problem overspending ',& 'when trying to account for unallocatable C ','') ENDIF ELSE ! Move the unallocated carbon back into the labile pool !DSGdebug_01b: no need to as we did not yet subtracted bm_alloc_tot !DSGdebug_01b biomass(ipts,j,ilabile,icarbon) = & !DSGdebug_01b biomass(ipts,j,ilabile,icarbon) - residual(ipts,j) ENDIF }}} == '''2. Mass conservation issue: stomate_phenology.f90''' == The nutrient demand must be calculated AFTER the final Cl_init and Cr_init are known. This can be fixed by: {{{ ! The biomass available to use is set to be the minimum of the biomass of ! the labile pool (if carbon not taken from the atmosphere), and ! the wanted biomass. bm_use(i) = MIN( biomass(i,j,ilabile,icarbon) + biomass(i,j,icarbres,icarbon), & bm_wanted(i) ) !DSGdebug_02 ! the nutrients need to support the biomass: !DSGdebug_02 bm_wanted_n(i) = (Cl_init + Cr_init*fcn_root(j))/cn_leaf_prescribed(j) }}} [...] {{{ ! In case nitrogen or phosphorus is not sufficiently available ! downregulate new leaf biomass to respect leaf stoichiometry; ! DSG: this violates the ratio used to calculate the ! leave-root-sapwood relationships: is this OK? !DSGdebug_02: moved after Cl_init and Cr_init are updated ! the nutrients need to support the biomass: bm_wanted_n(i) = (Cl_init + Cr_init*fcn_root(j))/cn_leaf_prescribed(j) !DSGdebug_02 }}} == '''3. Parameter value issue: constantes_mtc.f90''' == The parameter k_latosa_max and k_latosa_min were initially designed for tree PFT only. However these variables are also used for grass PFT (stomate_growth_fun_all.f90). Therefore these parameter cannot be set to undef. Parameter set to value of tree PFT. {{{ !DSGdebug_03 REAL(r_std), PARAMETER, DIMENSION(nvmc) :: k_latosa_max_mtc = & !! Maximum leaf-to-sapwood area ratio as defined in McDowell et al & (/ undef, 5000., 5000., 5000., 3000., 5000., 5000., & !! 2002, Oecologia and compiled in Hickler et al 2006, Appendix S2 ! & 5000., 5000., undef, undef, undef, undef /) !! (unitless) & 5000., 5000., 5000., 5000., 5000., 5000. /) !! (unitless) REAL(r_std), PARAMETER, DIMENSION(nvmc) :: k_latosa_min_mtc = & !! Minimum leaf-to-sapwood area ratio as defined in McDowell et al & (/ undef, 1500., 1500., 1500., 1000., 1500., 1500., & !! 2002, Oecologia and compiled in Hickler et al 2006, Appendix S2 ! & 1500., 1500., undef, undef, undef, undef /) !! (unitless) & 1500., 1500., 1500., 1500., 1500., 1500. /) !! (unitless) !DSGdebug_03 }}} == '''4. Stoichiometry issue: stomate_turnover.f90''' == There must be Nitrogen losses from abovesap wood when there are Carbon losses to ensure CN ratio. Here I assume that sapwood nitrogen can be recycled by plants according to leaves (this can be discussed). {{{ dturnover(:) = biomass(:,ivm,iroot,initrogen) * leaf_frac(:,ivm,ilage) * turnover_rate(:) biomass(:,ivm,ilabile,initrogen) = biomass(:,ivm,ilabile,initrogen) + recycle_root(ivm) * dturnover(:) biomass(:,ivm,iroot, initrogen) = biomass(:,ivm,iroot, initrogen) - dturnover(:) turnover(:,ivm,iroot, initrogen) = turnover(:,ivm,iroot,initrogen) + ( un - recycle_root(ivm) ) * dturnover(:) dturnover(:) = biomass(:,ivm,ifruit,initrogen) * leaf_frac(:,ivm,ilage) * turnover_rate(:) biomass(:,ivm,ifruit, initrogen) = biomass(:,ivm,ifruit,initrogen) - dturnover(:) turnover(:,ivm,ifruit,initrogen) = turnover(:,ivm,ifruit,initrogen) + dturnover(:) !DSGdebug_04 dturnover(:) = biomass(:,ivm,isapabove,initrogen) * leaf_frac(:,ivm,ilage) * turnover_rate(:) biomass(:,ivm,ilabile,initrogen) = biomass(:,ivm,ilabile,initrogen) + recycle_leaf(ivm) * dturnover(:) biomass(:,ivm,isapabove,initrogen) = biomass(:,ivm,isapabove,initrogen) - dturnover(:) turnover(:,ivm,isapabove,initrogen) = turnover(:,ivm,isapabove,initrogen) + ( un - recycle_leaf(ivm) ) * dturnover(:) !DSGdebug_04 }}} == '''5. Negative biomass: stomate_growth_fun_all.f90''' == In case the variable "residual" is not zero and nutrient availability reduces bm_alloc_tot to a value smaller than the value of residual, the allocated biomass gets negative. This can be fixed by correcting the biomass to be allocated by the residual before calculating nutrient limitation. {{{ ! Move the unallocated carbon back into the labile pool biomass(ipts,j,ilabile,icarbon) = & biomass(ipts,j,ilabile,icarbon) + residual(ipts,j) ENDIF !DSGdebug_05 ! correct the biomass to be allocated by the residual. ! this has to be done HERE as we need to known the actual biomass ! being allocated for the calculation of nutrient limitation bm_alloc_tot(ipts,j) = bm_alloc_tot(ipts,j) - residual(ipts,j) !DSGdebug_05 }}} [...] {{{ !DSGdebug_05 ! bm_alloc(:,j,k,icarbon) = f_alloc(:,j,k) * (bm_alloc_tot(:,j) - residual(:,j)) bm_alloc(:,j,k,icarbon) = f_alloc(:,j,k) * bm_alloc_tot(:,j) ! residual was already removed from bm_alloc_tot !DSGdebug_05 }}} Further changes needed regarding residual to ensure mass conservation: 1. The residual has to be initialized with zero: {{{ ! If npp is not initialized, bare soil value is never set. npp(:,:) = zero !DSGdebug_05a ! Not having this results in an unitilized error !DSGdebug_05a ! with valgrid, but I can't figure out why. It always !DSGdebug_05a ! seems to be set before being used. ! DSG: Of course you get errors, when a variable was never set but is used ! (for example in case first_call=true) ! I set it to zero, as it should be zero in the aformentioned case !DSGdebug_05a residual(:,:) = val_exp !DSGdebug_05a residual(:,:) = zero !DSGdebug_05a }}} 2. In case of nutrient down regulating bm_alloc the labile pool has to be corrected for not allocated residual {{{ !DSGdebug_01 biomass(:,:,ilabile,:) = biomass(:,:,ilabile,:) - alloc_tot(:,:,:) !DSGdebug_01 !DSGdebug_05a WHERE ((residual(:,:).GT. zero).AND.(residual(:,:) .GT. bm_alloc_tot(:,:))) biomass(:,:,ilabile,icarbon) = biomass(:,:,ilabile,icarbon) - (residual(:,:) - bm_alloc_tot(:,:)) END WHERE !DSGdebug_05a }}} == '''6. Mass conservation issue: stomate_phenology.f90''' == When calculating the share of new biomass to leaves (Cl_init) and roots (Cr_init) of each respective circ_class, the biomass has to be divided by the number of tree per class (circ_class_n); which wasn't done. {{{ ! Distribute the biomass over the leaves and roots (gC tree-1) ! Use Cl_init to estimate the share for each circumference class ! leaf biomass = FK * Cs / height (allometric relationship) ! root biomass = KF / LF * Cs / height ! Convert from gC m-2 to gC tree-1 ! +++CHECK+++ ! Cl_init + Cr_init can exceed bm_use. bm_use should be used in these equations ! DSG: I cannot confirm the statement in the line before ! but the equation were wrong, the biomass has to be divided by numbers of tree (circ_class_n) circ_class_biomass(i,j,l,ileaf,icarbon) = circ_class_biomass(i,j,l,ileaf,icarbon) + & Cl_init * ( KF(i,j) * Cs_tree(l) / height(i,j,l) * circ_class_n(i,j,l) ) / & !DSGdebug_06 SUM( KF(i,j) * Cs_tree(:) / height(i,j,:) * circ_class_n(i,j,:) ) SUM( KF(i,j) * Cs_tree(:) / height(i,j,:) * circ_class_n(i,j,:) ) / circ_class_n(i,j,l) !DSGdebug_06 circ_class_biomass(i,j,l,iroot,icarbon) = circ_class_biomass(i,j,l,iroot,icarbon) + & Cr_init * ( KF(i,j) * Cs_tree(l) / height(i,j,l) * circ_class_n(i,j,l) ) / & !DSGdebug_06 SUM( KF(i,j) * Cs_tree(:) / height(i,j,:) * circ_class_n(i,j,:) ) SUM( KF(i,j) * Cs_tree(:) / height(i,j,:) * circ_class_n(i,j,:) ) / circ_class_n(i,j,l) !DSGdebug_06 !++++++++++++ }}} == '''7. Mass conservation issue: stomate_phenology.f90''' == see in code comment for explanation {{{ ! Calculate the available biomass in roots, sapwood and leaves (gC ind-1) IF ( biomass(i,j,ileaf,icarbon) .LT. min_stomate .AND. biomass(i,j,iroot,icarbon) .LT. min_stomate) THEN Cs_grass = biomass(i,j,isapabove,icarbon) ELSE WRITE(6,*) 'There is leaf or root carbon that should not be here, something could have gone wrong in senescence' !DSGdebug_07: that doesn't justify to violate mass conservation Cs_grass = biomass(i,j,isapabove,icarbon) ENDIF }}} {{{ ! Cl_init and Cr_init were recalculated to properly account for bm_use (the available C) IF (j==test_pft) WRITE (6,*) 'Cr_init (end)=',Cr_init !DSGdebug_07 biomass(i,j,ileaf,icarbon) = Cl_init !DSGdebug_07 biomass(i,j,iroot,icarbon) = Cr_init biomass(i,j,ileaf,icarbon) = biomass(i,j,ileaf,icarbon) + Cl_init biomass(i,j,iroot,icarbon) = biomass(i,j,iroot,icarbon) + Cr_init !DSGdebug_07 }}} == 8. Parameter value issue: constantes_mtc.f90 == The value of k_sap for crops and grasses was set to "undef". New value of 3.e-4 (personal communication Sebastiaan Luyssaert). {{{ REAL(r_std), PARAMETER, DIMENSION(nvmc) :: k_sap_mtc = & !! Maximal sapwood specific conductivity. Values compiled in T. Hickler & (/ undef, 50., 10., 8., 5., 30., 8., & !! et al. 2006. @tex $(m^{2} s^{-1} MPa^{-1})$ @endtex !DSGdebug_08 & 20., 8., undef, undef, undef, undef /)*1.e-4 & 20., 8., 3., 3., 3., 3. /)*1.e-4 !DSGdebug_08 }}} == 9. Mass conservation issue: stomate_turnover.f90 [CORRECTED 09/06/2015] === The fraction of nitrogen recycled is calculated wrong when updating circ_class_biomass (*,ilabile,introgen). This can be fixed by: CORRECTION: "(ivm)" was missing for the variable "recycle_leaf" and "recycle_root" {{{ !nitrogen circ_class_biomass(:,ivm,icirc,ileaf,initrogen) = circ_class_biomass(:,ivm,icirc,ileaf,initrogen) * & (un - leaf_frac(:,ivm,ilage) * turnover_rate(:)) circ_class_biomass(:,ivm,icirc,iroot,initrogen) = circ_class_biomass(:,ivm,icirc,iroot,initrogen) * & (un - leaf_frac(:,ivm,ilage) * turnover_rate(:)) circ_class_biomass(:,ivm,icirc,ifruit,initrogen) = circ_class_biomass(:,ivm,icirc,ifruit,initrogen) * & (un - leaf_frac(:,ivm,ilage) * turnover_rate(:)) ! add the recycled nutrients to labile pool: circ_class_biomass(:,ivm,icirc,ilabile,initrogen) = & circ_class_biomass(:,ivm,icirc,ilabile,initrogen) & + circ_class_biomass(:,ivm,icirc,ileaf,initrogen) * & ! add from leaves !DSGdebug_12 (un - leaf_frac(:,ivm,ilage) * turnover_rate(:) * recycle_leaf(ivm)) & ( leaf_frac(:,ivm,ilage) * turnover_rate(:) * recycle_leaf) & + circ_class_biomass(:,ivm,icirc,iroot,initrogen) * & ! add from roots !DSGdebug_12 (un - leaf_frac(:,ivm,ilage) * turnover_rate(:) * recycle_root(ivm)) ( leaf_frac(:,ivm,ilage) * turnover_rate(:) * recycle_root) }}} [...] {{{ circ_class_biomass(:,ivm,icirc,ileaf,iphosphorus) = circ_class_biomass(:,ivm,icirc,ileaf,iphosphorus) * & (un - dt / leaffall(ivm)) circ_class_biomass(:,ivm,icirc,iroot,iphosphorus) = circ_class_biomass(:,ivm,icirc,iroot,iphosphorus) * & (un - dt / leaffall(ivm)) ! add the recycled nutrients to labile pool: circ_class_biomass(:,ivm,icirc,ilabile,initrogen) = & circ_class_biomass(:,ivm,icirc,ilabile,initrogen) & + circ_class_biomass(:,ivm,icirc,ileaf,initrogen) * & ! add from leaves !DSGdebug_12 (un - dt / leaffall(ivm) * recycle_leaf) & ( dt / leaffall(ivm)) * recycle_leaf(ivm) & + circ_class_biomass(:,ivm,icirc,iroot,initrogen) * & ! add from roots !DSGdebug_12 (un - dt / leaffall(ivm) * recycle_root) ( dt / leaffall(ivm)) * recycle_root(ivm) }}} === UPDATE 05/06/2015 === stomate_growth_fun_all.f90 Fruits were forgotten in these calculations {{{ ! Allocation of icarbres at the tree level (gC tree-1) circ_class_biomass(ipts,j,l,icarbres,m) = temp_share * & biomass(ipts,j,icarbres,m) / circ_class_n(ipts,j,l) IF((m .EQ. initrogen).OR.(m .EQ. iphosphorus)) THEN circ_class_biomass(ipts,j,l,ileaf,m) = temp_share * & biomass(ipts,j,ileaf,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,isapabove,m) = temp_share * & biomass(ipts,j,isapabove,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,isapbelow,m) = temp_share * & biomass(ipts,j,isapbelow,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,iheartabove,m) = temp_share * & biomass(ipts,j,iheartabove,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,iheartbelow,m) = temp_share * & biomass(ipts,j,iheartbelow,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,iroot,m) = temp_share * & biomass(ipts,j,iroot,m) / circ_class_n(ipts,j,l) !DSGdebug_12a circ_class_biomass(ipts,j,l,ifruit,m) = temp_share * & biomass(ipts,j,ifruit,m) / circ_class_n(ipts,j,l) !DSGdebug_12a ENDIF ELSE circ_class_biomass(ipts,j,l,ilabile,m) = zero circ_class_biomass(ipts,j,l,icarbres,m) = zero IF((m .EQ. initrogen).OR.(m .EQ. iphosphorus)) THEN circ_class_biomass(ipts,j,l,ileaf,m) = zero circ_class_biomass(ipts,j,l,isapabove,m) = zero circ_class_biomass(ipts,j,l,isapbelow,m) = zero circ_class_biomass(ipts,j,l,iheartabove,m) = zero circ_class_biomass(ipts,j,l,iheartbelow,m) = zero circ_class_biomass(ipts,j,l,iroot,m) = zero !DSGdebug_12a circ_class_biomass(ipts,j,l,ifruit,m) = zero !DSGdebug_12a ENDIF ENDIF }}} === UPDATE: 09/06/2015 ==== This might not be needed, but it decreases the risk the variable turnover is not set or set to a wrong value. {{{ WHERE ( senescence(:,ivm) ) ! Turnover at the stand level (gC m-2) turnover(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) * dt / leaffall(ivm) turnover(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) * dt / leaffall(ivm) turnover(:,ivm,ileaf,initrogen) = biomass(:,ivm,ileaf,initrogen) * dt / leaffall(ivm) turnover(:,ivm,iroot,initrogen) = biomass(:,ivm,iroot,initrogen) * dt / leaffall(ivm) !DSGdebug_12b: just to be sure... ELSEWHERE turnover(:,ivm,ileaf,icarbon) = zero turnover(:,ivm,ileaf,initrogen) = zero !DSGdebug_12b ENDWHERE }}} == 10. Parameter value issue: pft_parameters.f90 == The value for leaffall is set to undef for a couple of PFT, but it is used in stomate_turnover.f90 for every PFT. I fixed this by setting all undef to 10. like all the other PFTs: {{{ !Config Key = LEAFFALL !Config Desc = length of death of leaves, tabulated !Config if = OK_STOMATE !Config Def = undef, 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10., 10. !Config Help = !Config Units = [days] CALL getin_p('LEAFFALL',leaffall) }}} == 11. Mass conservation issue: stomate_turnover.f90 == The mass of biomass and circ_class_biomass starts to diverge after leaf fall in stomate_turnover.f90. I rewrote the code in a more fail-safe format and thereby fixed the bug. I don't know where the error is in the original piece of code. {{{ !! 4.3 Loop over leaf age classes DO ilage = 1, nleafages turnover_rate(:) = zero WHERE ( leaf_age(:,ivm,ilage) .GT. leaf_age_crit(:,ivm) / deux) turnover_rate(:) = MIN( 0.99_r_std, dt / ( leaf_age_crit(:,ivm) * & ( leaf_age_crit(:,ivm) / leaf_age(:,ivm,ilage) )**quatre ) ) ENDWHERE IF (is_tree(ivm)) THEN !DSGdebug_13 ! Stand level turnover (gC m-2) !DSGdebug_13 ! Leaves !DSGdebug_13 dturnover(:) = biomass(:,ivm,ileaf,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:) !DSGdebug_13 biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - dturnover(:) !DSGdebug_13 turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + dturnover(:) !DSGdebug_13 !DSGdebug_13 ! save leaf mass change !DSGdebug_13 delta_lm(:,ilage) = - dturnover(:) !DSGdebug_13 !DSGdebug_13 ! Roots !DSGdebug_13 dturnover(:) = biomass(:,ivm,iroot,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:) !DSGdebug_13 biomass(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) - dturnover(:) !DSGdebug_13 turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + dturnover(:) !DSGdebug_13 !DSGdebug_13 ! Fruit !DSGdebug_13 dturnover(:) = biomass(:,ivm,ifruit,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:) !DSGdebug_13 biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - dturnover(:) !DSGdebug_13 turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + dturnover(:) !DSGdebug_13 !DSGdebug_13 !DSGdebug_13 ! Stand level turnover (gC m-2) !DSGdebug_13 ! Leaves !DSGdebug_13 dturnover(:) = biomass(:,ivm,ileaf,initrogen) * leaf_frac(:,ivm,ilage) * turnover_rate(:) !DSGdebug_13 biomass(:,ivm,ilabile,initrogen) = biomass(:,ivm,ilabile,initrogen) + & !DSGdebug_13 recycle_leaf(ivm) * dturnover(:) !DSGdebug_13 !DSGdebug_13 biomass(:,ivm,ileaf,initrogen) = biomass(:,ivm,ileaf,initrogen) - dturnover(:) !DSGdebug_13 turnover(:,ivm,ileaf,initrogen) = turnover(:,ivm,ileaf,initrogen) + dturnover(:) * (1. - recycle_leaf(ivm)) !DSGdebug_13 !DSGdebug_13 !DSGdebug_13 ! Roots !DSGdebug_13 dturnover(:) = biomass(:,ivm,iroot,initrogen) * leaf_frac(:,ivm,ilage) * turnover_rate(:) !DSGdebug_13 biomass(:,ivm,ilabile,initrogen) = biomass(:,ivm,ilabile,initrogen) + & !DSGdebug_13 recycle_root(ivm) * dturnover(:) !DSGdebug_13 !DSGdebug_13 biomass(:,ivm,iroot,initrogen) = biomass(:,ivm,iroot,initrogen) - dturnover(:) !DSGdebug_13 turnover(:,ivm,iroot,initrogen) = turnover(:,ivm,iroot,initrogen) + dturnover(:) * (1. - recycle_root(ivm)) !DSGdebug_13 !DSGdebug_13 !DSGdebug_13 !DSGdebug_13 ! Fruit !DSGdebug_13 dturnover(:) = biomass(:,ivm,ifruit,initrogen) * leaf_frac(:,ivm,ilage) * turnover_rate(:) !DSGdebug_13 biomass(:,ivm,ifruit,initrogen) = biomass(:,ivm,ifruit,initrogen) - dturnover(:) !DSGdebug_13 turnover(:,ivm,ifruit,initrogen) = turnover(:,ivm,ifruit,initrogen) + dturnover(:) !DSGdebug_13 !DSGsimplify DO ielem = 1,nelements ! Stand level turnover (gC m-2) for ! - leaves, ... dturnover(:) = biomass(:,ivm,ileaf,ielem) * leaf_frac(:,ivm,ilage) * turnover_rate(:) biomass(:,ivm,ileaf,ielem) = biomass(:,ivm,ileaf,ielem) - dturnover(:) IF (ielem == icarbon) THEN ! no recycling turnover(:,ivm,ileaf,ielem) = turnover(:,ivm,ileaf,ielem) + dturnover(:) ! save leaf mass change; carbon only delta_lm(:,ilage) = - dturnover(:) ELSEIF (ielem == initrogen) THEN ! recycle biomass(:,ivm,ilabile,ielem) = biomass(:,ivm,ilabile,ielem) + recycle_leaf(ivm) * dturnover(:) turnover(:,ivm,ileaf, ielem) = turnover(:,ivm,ileaf, ielem) + ( un - recycle_leaf(ivm) ) * dturnover(:) ELSEIF (ielem == iphosphorus) THEN ! recycle biomass(:,ivm,ilabile,ielem) = biomass(:,ivm,ilabile,ielem) + p_recycle_leaf(ivm) * dturnover(:) turnover(:,ivm,ileaf, ielem) = turnover(:,ivm,ileaf, ielem) + ( un - p_recycle_leaf(ivm) ) * dturnover(:) ELSE STOP END IF ! - roots, ... dturnover(:) = biomass(:,ivm,iroot,ielem) * leaf_frac(:,ivm,ilage) * turnover_rate(:) biomass(:,ivm,iroot,ielem) = biomass(:,ivm,iroot,ielem) - dturnover(:) IF (ielem == icarbon) THEN ! no recycling turnover(:,ivm,iroot,ielem) = turnover(:,ivm,iroot,ielem) + dturnover(:) ! save root mass change; carbon only delta_lm(:,ilage) = - dturnover(:) ELSEIF (ielem == initrogen) THEN ! recycle biomass(:,ivm,ilabile,ielem) = biomass(:,ivm,ilabile,ielem) + recycle_root(ivm) * dturnover(:) turnover(:,ivm,iroot, ielem) = turnover(:,ivm,iroot, ielem) + ( un - recycle_root(ivm) ) * dturnover(:) ELSEIF (ielem == iphosphorus) THEN ! recycle biomass(:,ivm,ilabile,ielem) = biomass(:,ivm,ilabile,ielem) + p_recycle_root(ivm) * dturnover(:) turnover(:,ivm,iroot, ielem) = turnover(:,ivm,iroot, ielem) + ( un - p_recycle_root(ivm) ) * dturnover(:) ELSE STOP END IF ! ... and fruits. dturnover(:) = biomass(:,ivm,ifruit,ielem) * leaf_frac(:,ivm,ilage) * turnover_rate(:) biomass(:,ivm,ifruit,ielem) = biomass(:,ivm,ifruit,ielem) - dturnover(:) turnover(:,ivm,ifruit,ielem) = turnover(:,ivm,ifruit,ielem) + dturnover(:) ENDDO !DSGsimplify !DSGdebug_13 }}} == 12. Mass conservation issue: stomate_turnover.f90 == When updating the circ_class_biomass(*,ilabile,initrogen) with nitrogen from recyling the exisiting nitrogen in was forgotten. {{{ ! Tree level (gC tree-1) circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = zero circ_class_biomass(:,ivm,icirc,iroot,icarbon) = zero circ_class_biomass(:,ivm,icirc,ifruit,icarbon) = zero circ_class_biomass(:,ivm,icirc,ileaf ,initrogen) = zero circ_class_biomass(:,ivm,icirc,iroot ,initrogen) = zero circ_class_biomass(:,ivm,icirc,ifruit,initrogen) = zero circ_class_biomass(:,ivm,icirc,ileaf ,iphosphorus) = zero circ_class_biomass(:,ivm,icirc,iroot ,iphosphorus) = zero circ_class_biomass(:,ivm,icirc,ifruit,iphosphorus) = zero circ_class_biomass(:,ivm,icirc,ilabile ,initrogen) = & !DSGdebug_14 circ_class_biomass(:,ivm,icirc,ilabile ,initrogen) & !DSGdebug_14 + circ_class_biomass(:,ivm,icirc,ileaf ,initrogen) * recycle_leaf(ivm) & + circ_class_biomass(:,ivm,icirc,iroot ,initrogen) * recycle_root(ivm) }}} == 13. Stoichiometry issue: stomate_growth_fun_all.f90 == In case of limited growth (alloc_tot(ipts,j,icarbon).LT.min_stomate) all *_alloc must be set to zero: {{{ ELSE ! no biomass growth; thus set allocated nutrients to zero alloc_tot(ipts,j,:) = zero !DSGdebug_16 bm_alloc(ipts,j,:,:) = zero !DSGdebug_16 }}} == 14. Mass conservation issue: stomate_growth_fun_all.f90 == All the circ_class_biomass pools have to updated after the actual growth (accounting for nutrient limitation) is known. {{{ !DSGdebug_16 IF((m .EQ. initrogen).OR.(m .EQ. iphosphorus)) THEN ! we need to synchronize all the circ_class pools if it is ! carbon or nutrients: circ_class_biomass(ipts,j,l,ileaf,m) = temp_share * & biomass(ipts,j,ileaf,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,isapabove,m) = temp_share * & biomass(ipts,j,isapabove,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,isapbelow,m) = temp_share * & biomass(ipts,j,isapbelow,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,iheartabove,m) = temp_share * & biomass(ipts,j,iheartabove,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,iheartbelow,m) = temp_share * & biomass(ipts,j,iheartbelow,m) / circ_class_n(ipts,j,l) circ_class_biomass(ipts,j,l,iroot,m) = temp_share * & biomass(ipts,j,iroot,m) / circ_class_n(ipts,j,l) !DSGdebug_12a circ_class_biomass(ipts,j,l,ifruit,m) = temp_share * & biomass(ipts,j,ifruit,m) / circ_class_n(ipts,j,l) !DSGdebug_12a !DSGdebug_16 ENDIF }}} == 15. Mass conservation issue: stomate_turnover.f90 == The calculation of the resorption from ileaf and iroot to ilabile uses the stocks of ileaf and iroot. Therefore the calculation have to be done before the stocks of ileaf and iroot are calculated. {{{ ! Turnover at tree level (gC tree-1) DO icirc = 1,ncirc WHERE ( senescence(:,ivm) ) !DSGdebug_16: has to be updated before circ_class_biomass is updated: ! add the recycled nutrients to labile pool: circ_class_biomass(:,ivm,icirc,ilabile,initrogen) = & circ_class_biomass(:,ivm,icirc,ilabile,initrogen) & + circ_class_biomass(:,ivm,icirc,ileaf,initrogen) * & ! add from leaves !DSGdebug_12 (un - dt / leaffall(ivm) * recycle_leaf) & dt / leaffall(ivm) * recycle_leaf(ivm) & + circ_class_biomass(:,ivm,icirc,iroot,initrogen) * & ! add from roots !DSGdebug_12 (un - dt / leaffall(ivm) * recycle_root) dt / leaffall(ivm) * recycle_root(ivm) circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = circ_class_biomass(:,ivm,icirc,ileaf,icarbon) * & (un - dt / leaffall(ivm) ) circ_class_biomass(:,ivm,icirc,iroot,icarbon) = circ_class_biomass(:,ivm,icirc,iroot,icarbon) * & (un - dt / leaffall(ivm)) circ_class_biomass(:,ivm,icirc,ileaf,initrogen) = circ_class_biomass(:,ivm,icirc,ileaf,initrogen) * & (un - dt / leaffall(ivm)) circ_class_biomass(:,ivm,icirc,iroot,initrogen) = circ_class_biomass(:,ivm,icirc,iroot,initrogen) * & (un - dt / leaffall(ivm)) !DSGdebug_16 ! add the recycled nutrients to labile pool: !DSGdebug_16 circ_class_biomass(:,ivm,icirc,ilabile,initrogen) = & !DSGdebug_16 circ_class_biomass(:,ivm,icirc,ilabile,initrogen) & !DSGdebug_16 + circ_class_biomass(:,ivm,icirc,ileaf,initrogen) * & ! add from leaves !DSGdebug_16 !DSGdebug_12 (un - dt / leaffall(ivm) * recycle_leaf) & !DSGdebug_16 dt / leaffall(ivm) * recycle_leaf(ivm) & !DSGdebug_16 + circ_class_biomass(:,ivm,icirc,iroot,initrogen) * & ! add from roots !DSGdebug_16 !DSGdebug_12 (un - dt / leaffall(ivm) * recycle_root) !DSGdebug_16 dt / leaffall(ivm) * recycle_root(ivm) }}} [...] {{{ !! 6.1 For deciduous trees: next to leaves, also fruits and fine roots are dropped ! For deciduous trees: next to leaves, also fruits and fine roots are dropped: fruit ::biomass(:,ivm,ifruit) ! and fine root ::biomass(:,ivm,iroot) carbon pools are set to zero. IF ( is_tree(ivm) .AND. ( senescence_type(ivm) .NE. 'none' ) ) THEN ! Check whether we shed the remaining leaves. The condition depends on biomass(ileaf) ! so first calculate the sheding at the tree level. because biomass(ileaf) is not ! changed at the tree level the same statement can then be used at the stand level. DO icirc = 1,ncirc ! check whether we shed the remaining leaves WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. senescence(:,ivm) .AND. & ( biomass(:,ivm,ileaf,icarbon) .LT. (lai_initmin(ivm) / 2.)/sla(ivm) ) ) ! Tree level (g tree-1) !DSGdebug_16: the recycling has to be done before the source !pools are set to zero, as these pool are part of the equation: ! relocate nutrient circ_class_biomass(:,ivm,icirc,ilabile ,initrogen) = & !DSGdebug_14 circ_class_biomass(:,ivm,icirc,ilabile ,initrogen) & !DSGdebug_14 + circ_class_biomass(:,ivm,icirc,ileaf ,initrogen) * recycle_leaf(ivm) & + circ_class_biomass(:,ivm,icirc,iroot ,initrogen) * recycle_root(ivm) circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = zero circ_class_biomass(:,ivm,icirc,iroot,icarbon) = zero circ_class_biomass(:,ivm,icirc,ifruit,icarbon) = zero circ_class_biomass(:,ivm,icirc,ileaf ,initrogen) = zero circ_class_biomass(:,ivm,icirc,iroot ,initrogen) = zero circ_class_biomass(:,ivm,icirc,ifruit,initrogen) = zero !DSGdebug_16 circ_class_biomass(:,ivm,icirc,ilabile ,initrogen) = & !DSGdebug_16 !DSGdebug_14 !DSGdebug_16 circ_class_biomass(:,ivm,icirc,ilabile ,initrogen) & !DSGdebug_16 !DSGdebug_14 !DSGdebug_16 + circ_class_biomass(:,ivm,icirc,ileaf ,initrogen) * recycle_leaf(ivm) & !DSGdebug_16 + circ_class_biomass(:,ivm,icirc,iroot ,initrogen) * recycle_root(ivm) }}} = '''Mass conservation checks''' = Mass closure given by: mass_before + mass_change = mass_after == '''stomate_lpj.f90''' == {{{ !DSG mass conservation ======================================== mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3) !! 5. Grow new biomass - respiration, npp and allocation ! Call the allometry based allocation (based on Sitch et al 2003 and Zaehle et al 2010) CALL growth_fun_all (npts, dt_days, veget_max, PFTpresent, & senescence, when_growthinit, moiavail_growingseason, t2m_week, & gpp_daily, gpp_week, resp_maint_part, resp_maint, & resp_growth, npp_daily, biomass, age, & leaf_age, leaf_frac, use_reserve, t_photo_stress, & lab_fac, lai_target, ind, rue_longterm, & circ_class_n, circ_class_biomass, c0_alloc, cn_leaf_season, np_leaf_season, & KF, n_uptake_daily, p_uptake_daily) !DSG mass conservation ============================================ mass_change(:,:,icarbon) = npp_daily(:,:)*dt_days mass_change(:,:,initrogen) = SUM(n_uptake_daily(:,:,:),DIM=3) mass_change(:,:,iphosphorus) = p_uptake_daily(:,:) }}} {{{ !DSG mass conservation ======================================== mass_before(:,:,:) = SUM(biomass(:,:,:,:),DIM=3) CALL phenology_prognostic (npts, dt_days, PFTpresent, veget_max, & tlong_ref, t2m_month, t2m_week, gpp_daily, & maxmoiavail_lastyear, minmoiavail_lastyear, moiavail_month, moiavail_week, & gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, & senescence, time_hum_min, biomass, leaf_frac, & leaf_age, when_growthinit, co2_to_bm, circ_class_n, & circ_class_biomass, ind, c0_alloc, KF) !DSG mass conservation ============================================ mass_change(:,:,icarbon) = zero mass_change(:,:,initrogen) = zero mass_change(:,:,iphosphorus) = zero }}}