| Version 37 (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
!! 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
[...]
! 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
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
