wiki:Branches/MergeOCN/Goll

Version 33 (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(:,:)

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 

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