| Version 58 (modified by dgoll, 9 years ago) (diff) |
|---|
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:
- 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 - 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
