! =================================================================================================================================
! MODULE       : stomate_prescribe
!
! CONTACT      : orchidee-help _at_ listes.ipsl.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF         Initialize and update density, height, basal area and crown area.
!!
!!\n DESCRIPTION: None
!!
!! RECENT CHANGE(S): None
!!
!! REFERENCE(S)	:
!!
!! SVN          :
!! $HeadURL: svn://forge.ipsl.fr/orchidee/trunk/ORCHIDEE/src_stomate/stomate_prescribe.f90 $
!! $Date: 2023-11-22 20:28:54 +0100 (mer. 22 nov. 2023) $
!! $Revision: 8316 $
!! \n
!_ ================================================================================================================================

MODULE stomate_prescribe

  ! modules used:

  USE ioipsl_para
  USE xios_orchidee  
  USE time, ONLY : ts_annual_proc
  USE pft_parameters
  USE dynamic_parameters
  USE constantes
  USE function_library,    ONLY: calculate_c0_alloc, cc_to_lai, & 
       cc_to_biomass, check_vegetation_area, check_mass_balance, &
       wood_to_qmdia, Nmax, sort_circ_class_biomass, lai_to_biomass, &
       biomass_to_lai, wood_to_ba, get_printlev, weibull_class_dist, &
       calculate_rdi_boundaries

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC prescribe, prescribe_clear

  ! first call
  LOGICAL, SAVE             :: firstcall_prescribe = .TRUE.
!$OMP THREADPRIVATE(firstcall_prescribe)
  INTEGER(i_std), SAVE       :: printlev_loc       !! Local level of text output for this module
!$OMP THREADPRIVATE(printlev_loc)

CONTAINS

! =================================================================================================================================
!! SUBROUTINE   : prescribe_clear
!!
!>\BRIEF        : Set the firstcall_prescribe flag back to .TRUE. to prepare for the next simulation.
!_=================================================================================================================================

  SUBROUTINE prescribe_clear
    firstcall_prescribe=.TRUE.
  END SUBROUTINE prescribe_clear


!! ================================================================================================================================
!! SUBROUTINE   : prescribe
!!
!>\BRIEF         Works only with static vegetation and agricultural PFTs.
!!               Initialize biomass, density, presence in the first call
!!               and update them in the following.
!!
!! DESCRIPTION (functional, design, flags): \n
!! This module works only with static vegetation and agricultural PFTs.
!! In the first call, initialize density of individuals, biomass,
!! and leaf age distribution to some reasonable value. In the following calls,
!! these variables are updated following a stand replacing disturbance or
!! if the conditions are suitable for recruitment. 
!!
!! To fulfill these purposes, pipe model are used:
!! \latexonly 
!!     \input{prescribe1.tex}
!!     \input{prescribe2.tex}
!! \endlatexonly
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLES(S): ::ind, ::biomass
!!
!! REFERENCES   :
!! - Krinner, G., N. Viovy, et al. (2005). "A dynamic global vegetation model 
!!   for studies of the coupled atmosphere-biosphere system." Global 
!!   Biogeochemical Cycles 19: GB1015, doi:1010.1029/2003GB002199.
!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
!!   plant geography and terrestrial carbon cycling in the LPJ dynamic 
!!   global vegetation model, Global Change Biology, 9, 161-185.
!! - McDowell, N., Barnard, H., Bond, B.J., Hinckley, T., Hubbard, R.M., Ishii, 
!!   H., Köstner, B., Magnani, F. Marshall, J.D., Meinzer, F.C., Phillips, N., 
!!   Ryan, M.G., Whitehead D. 2002. The relationship between tree height and leaf 
!!   area: sapwood area ratio. Oecologia, 132:12–20.
!! - Novick, K., Oren, R., Stoy, P., Juang, F.-Y., Siqueira, M., Katul, G. 2009. 
!!   The relationship between reference canopy conductance and simplified hydraulic 
!!   architecture. Advances in water resources 32, 809-819.
!! - Ruger et al. 2009, J. of Ecol. doi: 10.1111/j.1365-2745.2009.01552.x  
!!
!! \n
!_ ================================================================================================================================

SUBROUTINE prescribe (&
     npts,            veget_max,                  dt,                PFTpresent, &
     everywhere,      when_growthinit,            leaf_frac,         circ_class_dist, &
     circ_class_n,    circ_class_biomass,         atm_to_bm,         forest_managed, &
     KF,              plant_status,               age,               npp_longterm, &
     lm_lastyearmax,  longevity_eff_leaf,         longevity_eff_sap, longevity_eff_root, &
     k_latosa_adapt,  light_tran_to_floor_season, lpft_replant,      species_change_map, &
     cn_leaf_init_2D, bm_sapl_2D,                 new_ind)

  !! 0. Parameters and variables declaration

  !! 0.1 Input variables

    INTEGER(i_std), INTENT(in)                             :: npts               !! Domain size (unitless)
    INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in)        :: forest_managed     !! forest management flag
    REAL(r_std), INTENT(in)                                :: dt                 !! time step (dt_days)   
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: veget_max          !! "maximal" coverage fraction of a PFT 
                                                                                 !! (LAI -> infinity) on ground. May sum to
                                                                                 !! less than unity if the pixel has
                                                                                 !! nobio area. (unitless; 0-1)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: longevity_eff_root !! Effective root turnover time that accounts
                                                                                 !! waterstress (days)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: longevity_eff_sap  !! Effective sapwood turnover time that accounts
                                                                                 !! waterstress (days)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: longevity_eff_leaf !! Effective leaf turnover time that accounts
                                                                                 !! waterstress (days)
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: lpft_replant       !! Set to true if a PFT has been clearcut
                                                                                 !! and needs to be replaced by another species
    INTEGER(i_std), DIMENSION(npts,nvm), INTENT(in)        :: species_change_map !! map which gives the PFT number that each PFT
    REAL(r_std), DIMENSION(npts,nvm),INTENT(inout)         :: light_tran_to_floor_season !! Seasonal fraction of light transmitted to canopy levels
    REAL(r_std), DIMENSION(ncirc), INTENT(in)              :: circ_class_dist    !! The probability distribution of trees 
                                                                                 !! in a circ class in case of a 
                                                                                 !! redistribution (unitless).
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)           :: cn_leaf_init_2D    !! initial leaf C/N ratio

  !! 0.2 Output variables
    REAL(r_std),DIMENSION(npts,nvm), INTENT(out)           :: new_ind            !! density of new individuals (trees m-2 day-1) 

  !! 0.3 Modified variables

    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)            :: PFTpresent         !! PFT present (0 or 1)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)        :: plant_status       !! Growth and phenological status of the plant
                                                                                 !! Different stati are listed in in constantes_var
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)        :: everywhere         !! is the PFT everywhere in the grid box or 
                                                                                 !! very localized (after its introduction) (?)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)        :: when_growthinit    !! how many days ago was the beginning of 
                                                                                 !! the growing season (days)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)        :: npp_longterm       !! "long term" net primary productivity
                                                                                 !! @tex ($gC m^{-2} year^{-1}$) @endtex
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)        :: age                !! mean age (years)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)        :: lm_lastyearmax     !! last year's maximum leaf mass for each PFT 
                                                                                 !! @tex ($gC m^{-2}$) @endtex
    REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout)  :: circ_class_n       !! Number of individuals in each circ class
                                                                                 !! @tex $(ind m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), INTENT(inout)  :: circ_class_biomass !! Biomass components of the model tree  
                                                                                 !! within a circumference class
                                                                                 !! class @tex $(g C ind^{-1})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nelements), INTENT(inout) :: atm_to_bm       !! CO2 or N taken from the atmosphere to get
                                                                                 !! the seed C and N
                                                                                 !! to create the seedlings 
                                                                                 !!  @tex (gC.m^{-2}dt^{-1})$ @endtex  
    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac       !! fraction of leaves in leaf age 
                                                                                 !! class (unitless;0-1)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: KF              !! Scaling factor to convert sapwood mass
                                                                                 !! into leaf mass (m) - this variable is 
                                                                                 !! passed to other routines
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)           :: k_latosa_adapt  !! Leaf to sapwood area adapted for long 
                                                                                 !! term water stress (m)
    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), INTENT(inout)  :: bm_sapl_2D !! Sapling biomass for the functional
                                                                                          !! allocation with a dimension for the 
                                                                                          !! circumference classes
                                                                                          !! @tex $(gC.ind^{-1})$ @endtex

  !! 0.4 Local variables
    REAL(r_std), DIMENSION(nvm)                       :: c0_alloc           !! Root to sapwood tradeoff parameter
    INTEGER(i_std)                                    :: ipts, ivm, ifm     !! index (unitless)
    INTEGER(i_std)                                    :: ipar, ilev         !! index (unitless)
    INTEGER(i_std)                                    :: icir, iele, imbc   !! index (unitless)
    INTEGER(i_std)                                    :: deb,fin, imaxt     !! index (unitless)
    REAL(r_std), DIMENSION(npts,nvm)                  :: LF                 !! Scaling factor to convert sapwood mass
                                                                            !! into root mass (unitless)
    REAL(r_std), DIMENSION(npts,nvm)                  :: lstress_fac        !! Light stress factor, based on total
                                                                            !! transmitted light (unitless, 0-1)
    REAL(r_std), DIMENSION(npts,nvm,ncirc)            :: circ_class         !! circumference of individual trees
    REAL(r_std), DIMENSION(nparts)                    :: bm_init            !! Biomass needed to initiate the next 
                                                                            !! planting @tex $(gC m^{-2})$ @endtex 
    REAL(r_std), DIMENSION(ncirc)                     :: nb_trees_i         !! Number of trees in each twentith 
                                                                            !! circumference quantile of the 
                                                                            !! distribution (ind) 
    REAL(r_std)                                       :: excedent           !! Number of trees after truncation to be 
                                                                            !! reallocated to smaller quantiles of the 
                                                                            !! distribution (ind) 
    REAL(r_std)                                       :: max_circ_init      !! Maximum initial circumferences of the 
                                                                            !! truncated exponential distribution (cm) 
    REAL(r_std), DIMENSION(ncirc)                     :: ave_tree_height    !! The height of the ideal tree in each
    REAL(r_std)                                       :: dia_init_recr      !! Initial diameter of a recruit (m). Note that 
                                                                            !! this value may differ from the input argument 
                                                                            !! dia_int which is the diameter distribution of 
                                                                            !! a newly planted PFT
    REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements):: check_intern       !! Contains the components of the internal
                                                                            !! mass balance chech for this routine
                                                                            !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nelements)        :: closure_intern     !! Check closure of internal mass balance
                                                                            !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nelements)        :: pool_start         !! Start and end pool of this routine 
                                                                            !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nelements)        :: pool_end           !! Start and end pool of this routine 
                                                                            !! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
    REAL(r_std)                                       :: sla_est            !! A first estimate of sla in case its calculation is 
                                                                            !! dynamic @tex $(m^2.gC^{-1})$ @endtex
    REAL(r_std), DIMENSION(ncirc)                     :: bypass_lm          !! leaf mass before it is set to zero
                                                                            !! for deciduous PFTs. This value is
                                                                            !! in the calculation of KF (gC tree-1)
    REAL(r_std), DIMENSION(npts,nvm)                  :: veget_max_begin    !! temporary storage of veget_max to check area conservation    
    REAL(r_std)                                       :: sapwood_density    !! Sapwood density gC m^{-3}
    REAL(r_std)                                       :: ind_init           !! Maximum number of individuals when establishing a new forest (#/m^{-2})     
    REAL(r_std)                                       :: wood_init          !! Wood mass of a sapling (g C tree^{-1})i
    REAL(r_std),DIMENSION(ncirc)                      :: ba_init            !! intitial basal area
    LOGICAL                                           :: establish          !! Logical flag to establish a new young vegetation
                                                                            !! after a stand replacing disturbance 
    LOGICAL                                           :: recruite           !! Logical flag to add saplings of a PFT under an
                                                                            !! existing canopy of the same PFT 
    REAL(r_std), DIMENSION(npts,nvm,ncirc)            :: old_circ_class_n   !! Old number of individuals in each circ class
                                                                            !! @tex $(ind m^{-2})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements) &
                                                      :: old_circ_class_biomass !! Old biomass components of the model tree  
                                                                            !! within a circumference class
                                                                            !! class @tex $(g C ind^{-1})$ @endtex
    REAL(r_std), DIMENSION(npts,nvm,nparts,nelements)  :: old_biomass       !! Old stand level biomass 
                                                                            !! tex $(gC.m^{-2})$ @endtex
    REAL(r_std), DIMENSION(ncirc)                      :: height_small      !! Height of the smallest tree in circ_class_biomass
                                                                            !! tex $(m)$ @endtex
    REAL(r_std)                                        :: dia_small         !! Diameter of the smallest tree in circ_class_biomass
                                                                            !! tex $(m)$ @endtex   
    REAL(r_std)                                        :: height_sapl       !! Height of the recruitment sapling
                                                                            !! tex $(m)$ @endtex
    REAL(r_std)                                        :: old_ind           !! density of old individuals                                        
    REAL(r_std)                                        :: cn_leaf,cn_root   !! CN ratio of leaves, and roots
                                                                            !! (gC/gN)
    REAL(r_std)                                        :: cn_wood           !! CN ratio of wood pool
                                                                            !! (gC/gN)
    REAL(r_std)                                        :: coeff 
    REAL(r_std), DIMENSION(ncirc)                      :: dummy
   
    REAL(r_std)                                        :: lambda            !! lambda of the truncated exponential
                                                                            !! distribution of initial diameters (1/cm**2)
    REAL(r_std)                                        :: tmp_ind           !! temporary number of individuals per m-2
    REAL(r_std), DIMENSION(ncirc)                      :: circ_class_n_new  !! variable used to sort circ_class_n
    REAL(r_std), DIMENSION(ncirc,nparts,nelements)     :: circ_class_biomass_new !! variable used to sort circ_class_biomass
    REAL(r_std)                                        :: k_latosa_tmp      !! Temporaray variable in the calculation of 
                                                                            !! k_latosa_adapt
    REAL(r_std), DIMENSION(nvm)                        :: n_init            !! Initial number of individuals calculated
                                                                            !! based on gradient calculation (ind m^(-2))
    REAL(r_std), DIMENSION(npts,nvm)                   :: tmp               !! Temporary variable to prepare information for xios
    REAL(r_std), DIMENSION(nvm,ncirc)                  :: dia_init          !! Initial diameter distribution of a newly planted PFT (m)
    REAL(r_std), DIMENSION(nvm)                        :: qmd_init          !! quadratic mean diameter of a newly planted PFT (m)
    REAL(r_std), DIMENSION(ncirc)                      :: circ_class_dist_loc !! local circonference class distribution (wiebull) 
!_ ================================================================================================================================
 
   IF (firstcall_prescribe) THEN
      !! Initialize local printlev
      printlev_loc=get_printlev('stomate_prescribe')
      firstcall_prescribe=.FALSE.
    END IF

    IF (printlev_loc>=2) WRITE(numout,*) 'Entering stomate_prescribe'

  !! 1. Initialize biomass at first call

    !! 1.0 Initialize variables that won't be calculated for all PFTs,
    !!     to avoid a NaN if the value for bare soil is ever written.
    c0_alloc(:)=zero

    !! 1.1 Store for closing the mass balance with recruitment   
    old_biomass(:,:,:,:)= cc_to_biomass(npts,nvm,&
         circ_class_biomass(:,:,:,:,:), & 
         circ_class_n(:,:,:))
    new_ind(:,:) = zero

    !! 1.3 Initialize check for mass balance closure 
    IF (err_act.GT.1) THEN

       check_intern(:,:,:,:) = zero
       pool_start(:,:,:) = zero
       DO iele = 1,nelements

          ! atm_to_bm has as intent inout, the variable 
          ! accumulates carbon over the course of a day. 
          ! Use the difference between start and the end of 
          ! this routine
          check_intern(:,:,iatm2land,iele) = - un * &
               atm_to_bm(:,:,iele) * veget_max(:,:) * dt

          DO ipar = 1,nparts
          
             DO icir = 1,ncirc

                ! Initial biomass pool. Use circ_class_biomass
                ! because that variable is the prognostic variable
                ! biomass is syncronized to circ_class_biomass. If 
                ! many diameter classes are used, each diameter class
                ! separatly passes all criteria but when combined
                ! a mass balance problem occurs.
                pool_start(:,:,iele) = pool_start(:,:,iele) + &
                   circ_class_biomass(:,:,icir,ipar,iele) * &
                   circ_class_n(:,:,icir) * veget_max(:,:)

             ENDDO
             
          ENDDO

       ENDDO

       !! 1.3 Initialize check for area conservation
       veget_max_begin(:,:) = veget_max(:,:)

    ENDIF ! err_act.GT.1

  !! 2. Prescribe carbon and nitrogen biomass

    ! We want the vegetation to regrow after a stand replacing disturbance
    ! So loop over every point and PFT.  If there is supposed to be 
    ! vegetation here (according to veget_max) but there isn't (according 
    ! to ind), then establish a new young vegetation in that PFT.
    DO ivm = 2,nvm ! Loop over PFTs

       DO ipts = 1, npts ! Loop over pixels

          !! 2.1 Wait for the end of the year to replant
          !  If we are regrowing different species after managed forests
          !  are killed/die, we sometimes need to prevent them regrowing.
          !  This is only true if they die in the middle of the year due 
          !  to natural causes, in which case they will be replanted at 
          !  the end of the year. This loop checks to see if we regrow 
          !  this PFT now or later. We can regrow a PFT after lpft_replant 
          !  is set back to false so before that we need to process all the 
          !  information such that the correct species will be regrown. 
          !  lpft_replant is an optional variable. First test whether it is 
          !  present.
          IF(ok_change_species)THEN
             IF(lpft_replant(ipts,ivm))THEN
                CYCLE
             ENDIF
          ENDIF

          ! Debug
          IF(printlev_loc.GE.4 .AND. ipts == test_grid &
               .AND. ivm == test_pft)THEN
             WRITE(numout,*) 'PFT ',ivm
             WRITE(numout,*) 'plant_status ', plant_status(ipts,ivm)
             WRITE(numout,*) 'veget_max ', veget_max(ipts,ivm)
             WRITE(numout,*) 'recruitment_pft ',recruitment_pft(ivm)
             WRITE(numout,*) 'ok_dgvm ',  ok_dgvm
             WRITE(numout,*) 'biomass - C, ', &
                  SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1))
             WRITE(numout,*) 'biomass - N, ', &
                  SUM(SUM(circ_class_biomass(ipts,ivm,:,:,initrogen),1))
             WRITE(numout,*) 'ts_annual_proc ', ts_annual_proc
             WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(ipts,ivm,:))
             WRITE(numout,*) 'dens_target, ', dens_target(ivm)*ha_to_m2
             WRITE(numout,*) 'do_now_recruit, ',do_now_recruit
          ENDIF
          !

          ! Decide whether we have to establish new vegetation (::establish)
          ! allow recruitment (::recruite). These two local variables were
          ! introduced to allow maximal flexibility so within a single
          ! run, PFTs with and without recruitment can co-exist.
          IF( ( (plant_status(ipts,ivm) .EQ. iprescribe) .OR. & 
               (plant_status(ipts,ivm) .EQ. idead) ) .AND. &
               (SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LE. min_stomate) .AND. &
               (veget_max(ipts,ivm) .GT. min_stomate) .AND. &
               (.NOT. ok_dgvm) ) THEN

             ! The vegetation just died from replacing disturbances
             ! (two first conditions). The new vegetation could be 
             ! of the same or a different PFT than the previous 
             ! vegetation. 
             establish = .TRUE.
             recruite = .FALSE.

             ! Debug
             IF (printlev_loc>=4) THEN
                WRITE(numout,*) 'We are calling stomate_prescribe.f90 &
                     & to initialize PFT  ',ivm
             ENDIF
             !-

          ELSEIF( ( (plant_status(ipts,ivm) .NE. iprescribe) .OR. & 
               (plant_status(ipts,ivm) .NE. idead) ) .AND. &
               (SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .GT. min_stomate) .AND. &
               (veget_max(ipts,ivm) .GT. min_stomate) .AND. &
               (recruitment_pft(ivm)) .AND. &
               (.NOT. ok_dgvm) .AND. &
               (ts_annual_proc) .AND. &
               (SUM(circ_class_n(ipts,ivm,:)) .GT. dens_target(ivm)*ha_to_m2) .AND. &
               do_now_recruit) THEN

             ! There is already vegetation (first three conditions)
             ! and this specific PFT has recruitment (condition 5)
             ! Only use recruitment when the recruitment flag is TRUE
             ! (condition 7)
             ! no recruit when the stand density is below dens_target because 
             ! the model will slowly  but surely kill the forest before a new
             ! establishment
             establish = .FALSE.
             recruite = .TRUE. 

             ! Debug
             IF (printlev_loc>=4) THEN
                WRITE(numout,*) 'We are calling stomate_prescribe.f90 &
                     & for sapling recruitment in PFT  ',ivm
                WRITE(numout,*) 'Density of stem', SUM(circ_class_n(ipts,ivm,:))
             ENDIF
             !-
             
          ELSEIF( (veget_max(ipts,ivm).LT. min_stomate) .AND. &
             (SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1)) .LT. min_stomate) .AND. &
             (SUM(circ_class_n(ipts,ivm,:)) .LT. min_stomate) )THEN

             ! The PFT is not defined. Initialize the variables
             ! to ensure that all goes well when moving and merging 
             ! different age classes.
             circ_class_n(ipts,ivm,:) = zero
             circ_class_biomass(ipts,ivm,:,:,icarbon) = zero 

             ! Conditions are not fullfilled for establishing
             ! new vegetation
             establish = .FALSE.
             recruite = .FALSE.
             
          ELSEIF (veget_max(ipts,ivm) .GT. min_stomate) THEN

             ! Conditions are not fullfilled for either establishing 
             ! a new young vegetation or having recruitment under an
             ! exisisting vegetation of the same PFT.
             establish = .FALSE.
             recruite = .FALSE.

             ! Debug
             IF (printlev_loc>=4) THEN
                WRITE(numout,*) 'We are calling stomate_prescribe.f90 &
                     & but nothing should be done for PFT ',ivm
             ENDIF
             
          ELSE

             ! Unexpected condition
             CALL ipslerr_p (3,'stomate_prescribe.f90',&
                  'Unexpected condition - time to rethink the', &
                  'IF that choses between establish and recruite','')

          ENDIF

          !! 2.2 Calculate variables that will determine the size of the saplings
          IF( establish .OR. recruite )THEN

             !! 2.2.1 Calculate c0_alloc
             ! The calculation starts with c0_alloc, it only depends on PFT and 
             ! the effective longiveties.
             c0_alloc(ivm) = calculate_c0_alloc(ipts, ivm, &
                  longevity_eff_root(ipts,ivm),longevity_eff_sap(ipts,ivm))
             
             !! 2.2.2 Calculate k_latosa for newly established PFTs
             !  If a new vegetation has to be established we recalculate
             !  k_latosa from the parameter file (see note on physiological
             !  memory of k_latosa below). If we add recruits we use the KF of
             !  the existing vegetation.
             IF (establish) THEN

                ! Initialize tree level variables
                circ_class_n(ipts,ivm,:) = zero
                circ_class_biomass(ipts,ivm,:,:,:) = zero

                ! Assume that there is no physiological memory for k_latosa
                ! between different generations (this is probably true when
                ! trees are planted but seems less likely for natural
                ! regeneration. The first test that implemented memory has
                ! caused problems with the LAI from the second generation
                ! onwards. Nevertheless, we still use k_latosa_adapt. In case
                ! someone want to try again in the future, the structure of
                ! the code is ready to account for memory
                k_latosa_adapt(ipts,ivm) = k_latosa_min(ivm)

                ! Debugging
                IF (printlev_loc>=4) THEN
                   WRITE(numout,*) 'Establish a PFT (stomate_prescribe.f90):'
                   WRITE(numout,*) ' Imposing initial biomass for prescribed &
                        & trees, initial reserve mass for prescribed grasses.'
                   WRITE(numout,*) ' Declaring prescribed PFTs present.'
                   WRITE(numout,*) ' Grid point: ',ipts,' PFT Type: ',ivm
                   WRITE(numout,*) ' Establish ',establish,' Recruite ',recruite 
                ENDIF
                !-
                
             ENDIF ! establish

             !! 2.2.3 Calculate Lightstress
             !  Note that light and water stress have an opposite effect on KF. 
             !  Waterstress will decrease KF because more C should be allocated 
             !  to the roots. Light stress will increase KF because more C should 
             !  be allocated to the leaves.
             !  Lightstress varies from 0 to 1 and is calculated from the canopy 
             !  structure (veget).
             IF (establish) THEN

                ! Given that there is no vegetation at this point, the
                ! lightstress cannot be calculated and is set to 1.
                ! However, as soon as we grow a canopy it will experience
                ! light stress and so KF should be adjusted.
                lstress_fac(ipts,ivm) = un

             ELSEIF (recruite) THEN 

                ! We already have vegetation, so we should do recruitment if 
                ! we want it for this PFT. First lets compute light stress based 
                ! on the fraction of light reaching the ground (level 1),
                ! averaged over all daytime time steps (sinang>0) of the vegetative 
                ! period (gpp>0). This is based on Pgap and accounts for canopy
                ! structure. 
                lstress_fac(ipts,ivm) = light_tran_to_floor_season(ipts,ivm)

                !This value is accumulated in stomate.f90. It is accumulated
                !over an entire year and used only the last day of the year.
                !After it has been used, it can be reset. Note that 
                !age_class_dist takes place before recruitment so it has be
                !moved from one PFT to another first.


                light_tran_to_floor_season(ipts,ivm) = zero

                ! Debug
                IF(printlev_loc>=4)THEN
                   WRITE(numout,*) 'lstress_fac in stomate_prescribe is ', &
                        lstress_fac(ipts,ivm)*100, ' (%)'
                ENDIF
                !-

             ENDIF ! establish/recruite for light

          ENDIF ! establish .OR. recruite for initialization
          
          !! 2.3 Initial vegetation has to be prescribed only when the vegetation is 
          ! static
          IF (establish .OR. recruite) THEN

             !! 2.3.1 Initialize woody PFT's 
             !  Use veget_max to check whether the PFT is present. If the PFT 
             !  is present but it has no biomass, prescribe its biomass (gC m-{2}).
             IF ( is_tree(ivm) ) THEN

                ! prescribe/recruitment branching for actual process
                IF (establish) THEN

                   ! Allocation factors
                   ! Sapwood to root ratio
                   ! Following Magnani et al. 2000 "In order to decreases hydraulic 
                   ! resistance, the investment of carbon in fine roots or sapwood 
                   ! yields to the plant very different returns, both because of 
                   ! different hydraulic conductivities and because of the strong 
                   ! impact of plant height on shoot resistance. On the other hand, 
                   ! fine roots and sapwood have markedly different longevities and 
                   ! the cost of production, discounted for turnover, will differ 
                   ! accordingly. Optimal growth under hydraulic constraints requires 
                   ! that the ratio of marginal hydraulic returns to marginal annual 
                   ! cost for carbon investment in either roots or sapwood be the 
                   ! same (Bloom, Chapin & Mooney 1985; Case & Fair 1989). This
                   ! is formalized in equation (13) and further derived to obtain 
                   ! equation (17) in Magnani et al 2000. The latter is implemented 
                   ! here.Pipe_density is given in gC/m-3, convert to kg/m-3. And 
                   ! apply equation (17) in Magnani et al 2000. Note that c0_alloc 
                   ! was calculated at the start of this routine. The calculation 
                   ! itself is done in function_library

                   ! Calculate leaf area to sapwood area
                   ! To be consistent with the hydraulic limitations and pipe theory,
                   ! k_latosa is calculated from equation (18) in Magnani et al.
                   ! To do so, total hydraulic resistance and tree height need to 
                   ! known. This poses a problem as the resistance depends on the leaf 
                   ! area and the leaf area on the resistance. There is no independent 
                   ! equation and equations 12 and 18 depend on each other and 
                   ! substitution would be circular. Hence prescribed k_latosa values 
                   ! were obtained from observational records and are given in 
                   ! mtc_parameters.f90. 

                   ! The relationship between height and k_latosa as reported in 
                   ! McDowell et al 2002 and Novick et al 2009 could be implemented 
                   ! to adjust k_latosa for the height of the stand. This did NOT 
                   ! result in a realistic model behavior 
                   ! k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * &
                   !    (k_latosa_max(ivm) - (latosa_height(ivm) * &
                   !    (SUM( nb_trees_i(:) * ave_tree_height(:) ) / &
                   !       SUM( nb_trees_i(:) ))))

                   ! Also k_latosa has been reported to be a function of CO2 
                   ! concentration (Atwell et al. 2003, Tree Physiology, 23:13-21 
                   ! and Pakati et al. 2000, Global Change Biology, 6:889-897). 
                   ! This effect is not accounted for in the current code

                   ! Finally, k_latosa is also reported to be a function of 
                   ! diameter (i.e. stand thinning, Simonin et al 2006, Tree 
                   ! Physiology, 26:493-503). Here the relationship with thinning was 
                   ! interpreted as a realtionship with light stress. Note that light 
                   ! stress cannot be calculated at this time in the model because 
                   ! there is no canopy (that's why we are in prescribe!) and
                   ! so there is no lightstress. lstress was therefore set to one 
                   ! (see above). We prefered to keep this redundancy in the code 
                   ! because it makes it clear that k_latosa is calculated in the 
                   ! same way in prescribe.f90 and growth_fun_all.f90
                   ! k_latosa(ipts,ivm) = (k_latosa_adapt(ipts,ivm) + &
                   !     (lstress_fac(ipts,ivm) * &
                   !     (k_latosa_max(ivm)-k_latosa_min(ivm))))

                   ! +++CHECK+++
                   ! Ideally waterstress and lightstress should both
                   ! affect k_latosa. One of the following approaches could
                   ! be used as a starting point.
                   !  k_latosa(ipts,ivm) = k_latosa_min(ivm) + &
                   !     (wstress_fac(ipts,ivm) * lstress_fac(ipts,ivm) * &
                   !     (k_latosa_max(ivm)-k_latosa_min(ivm)))
                   !  k_latosa(ipts,ivm) = wstress_fac(ipts,ivm) * (k_latosa_min(ivm) + &
                   !      (lstress_fac(ipts,ivm) * &
                   !      (k_latosa_max(ivm)-k_latosa_min(ivm))))
                   ! ++++++++++++

                   ! Calculate conversion coefficient for sapwood area to leaf 
                   ! area 
                   ! (1) The scaling parameter between leaf and sapwood mass 
                   ! is derived from LA_ind = k_latosa * SA_ind, where 
                   ! LA_ind = leaf area of an individual, SA_ind is the sapwood 
                   ! area of an individual and k_latosa a pipe-model parameter
                   ! (2) LA_ind = Cl * sla
                   ! (3) Cs = SA_ind * height * wooddensity * tree_ff
                   ! Substitute (2) and (3) in (1)
                   ! Cl = Cs * k1 / (wooddensity * sla * tree_ff * height)
                   ! Cl = Cs*KF/height, where KF is in (m)
                   ! KF is passed to the allocation routine and it is saved in the 
                   ! restart file.
                   !+++CHECK+++
                   ! At one point it looked like a good idea to take the max of two
                   ! options but by doing so we cannot recalculate KF in phenology.
                   ! It has not been confirmed that this is really a problem. Use the
                   ! simpelest approach but leave the alternative in the code as a
                   ! suggestion of a possible solution in case something goes wrong.
!!$                   k_latosa_tmp = MAX(k_latosa_min(ivm),k_latosa_adapt(ipts,ivm) + &
!!$                        (lstress_fac(ipts,ivm) * &
!!$                        (k_latosa_max(ivm)-k_latosa_min(ivm))))
                   !+++++++++++
                   k_latosa_tmp = k_latosa_adapt(ipts,ivm) + &
                        (lstress_fac(ipts,ivm) * &
                        (k_latosa_max(ivm)-k_latosa_min(ivm)))

                   IF (sla_dyn) THEN 
                      KF(ipts,ivm) = k_latosa_tmp / &
                           (slainit(ivm) * pipe_density(ivm) * tree_ff(ivm))
                   ELSE
                      KF(ipts,ivm) = k_latosa_tmp / &
                           (sla(ivm) * pipe_density(ivm) * tree_ff(ivm))
                   END IF

                   ! Initial quadratic mean diameter of a newly planted PFT.
                   IF (ncirc .EQ. 1) THEN
                     dia_init(ivm,1) = (dia_init_min(ivm)+ dia_init_max(ivm))/2
                   ELSE
                     DO icir = 1,ncirc
                       dia_init(ivm,icir) = dia_init_min(ivm) + &
                         (icir-1) * (dia_init_max(ivm)- dia_init_min(ivm))/(ncirc-1)
                     ENDDO
                   ENDIF

                   ifm = forest_managed(ipts,ivm)
                   ! Calculate quadratic mean diameter of a newly planted stand (m)
                   ! We need to know the stand density to calculate the qmd but need
                   ! the qmd to calculate circ_class_dist. This is a circular problem.
                   ! The circularity is broken by using the middle diameter class to 
                   ! approximate the initial distribution. Note that the distribution
                   ! also depends on the forest management system.
                   icir = CEILING(ncirc/deux)

                   circ_class_dist_loc(:) = & 
                        weibull_class_dist(SQRT(dia_init(ivm,icir)**2)*m_to_cm, ivm, ifm)
                   qmd_init(ivm) = SQRT(SUM(dia_init(ivm,icir:ncirc)**2*circ_class_dist_loc(icir:ncirc)) / &
                        SUM(circ_class_dist_loc(:)))

                   ! Debug
                   IF (printlev_loc.GE.4.AND.ivm==test_pft) THEN
                      WRITE(numout,*) 'Weibull_class_dist, ipts, ivm, icir', ipts, ivm, icir
                      WRITE(numout,*) 'dia_init, forest_managed',  dia_init(ivm,icir), ifm
                      WRITE(numout,*) 'qmd_init, circ_class_dist_loc' , qmd_init(ivm), circ_class_dist_loc(:)
                   END IF
                   !-

                   ! Multiple approaches for initializing the forest have been tried.
                   ! 6) For managed forest a yield based RDI approch is used to
                   ! thin the trees at low density in order to fasten the young forest stage. 
                   ! !!WARNING!! if the RDI upper target is too low the forest cannot accumalate
                   ! enough sugar for the next budburst due to a low LAI. To avoid that we decide
                   ! to set a thresold below which a fix value of 0.4 for RDI upper target is used.
                   !+++CHECK+++
                   ! The multiplication with 0.75 is a quick fix. It just works
                   ! better. It is very arbitrary.
                   if (ifm .NE. ifm_cop) then
                     n_init(ivm) = MIN(Nmax(qmd_init(ivm) * m_to_cm, ivm) * rdi_start(ivm), &
                                       nmaxplants(ivm))
                   else
                     n_init(ivm) = MIN(Nmax(coppice_diameter(ivm) * m_to_cm, ivm) * rdi_max(ivm, ifm), & 
                                     nmaxplants(ivm))
                   end if
                   n_init(ivm) = n_init(ivm)*ha_to_m2

                   ! Debug
                   IF (printlev_loc.GT.4) THEN
                      WRITE(numout,*) 'n_init (trees/m2), ', ivm, n_init(ivm)
                   END IF
                   !-

                   DO icir = 1,ncirc
                      ! Make use of dia_init (calculated in stomate.f90 as a parameter)
                      ! do calculate the height and wood mass of the saplings of a newly
                      ! planted PFT.
                      ave_tree_height(icir) = pipe_tune2(ipts,ivm)* &
                           (dia_init(ivm,icir)**pipe_tune3(ivm))
                      wood_init = ( ave_tree_height(icir) * pi / 4. * &
                           (dia_init(ivm,icir)) ** 2. ) * pipe_density(ivm) * tree_ff(ivm)

                      CALL make_saplings(bm_sapl_2D(ipts,ivm,icir,:,:), &
                           dia_init(ivm,icir), icir, ivm, &
                           ipts, KF(ipts,ivm), c0_alloc(ivm),plant_status(ipts,ivm), &
                           establish, cn_leaf_init_2D(ipts,ivm),bypass_lm(icir))
                   ENDDO

                   ! Prescribe the newly planted population.
                   ! qmd_init is the diameter of the tree
                   ! in meters.  This is the diameter that we will use
                   ! to calculate the new distribution of trees among
                   ! the circumference classes.
                   circ_class_dist_loc(:) = weibull_class_dist(qmd_init(ivm)*m_to_cm, ivm, &
                        forest_managed(ipts,ivm))
                   circ_class_n(ipts,ivm,:) = circ_class_dist_loc(:) / SUM(circ_class_dist_loc(:)) * n_init(ivm)
                   circ_class_biomass(ipts,ivm,:,:,:) = bm_sapl_2D(ipts,ivm,:,:,:)

                   ! Debug
                   IF (printlev_loc>=4) THEN
                      DO icir=1,ncirc
                        WRITE(numout,*)'PFT, ', ivm
                         WRITE(numout,*)'icir, circ_class_n, ', icir, &
                              circ_class_n(ipts,ivm,icir)
                         WRITE(numout,*)'Initial biomass distribution - carbon'
                         WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,icarbon)
                         WRITE(numout,*) circ_class_biomass(ipts,ivm,icir,:,initrogen)
                      ENDDO
                      WRITE(numout,*)' qmd_init', qmd_init(ivm)
                      WRITE(numout,*)' n_init(ivm)', n_init(ivm)
                      WRITE(numout,*)' Maximun number of saplings', Nmax(qmd_init(ivm)*m_to_cm,ivm)
                   END IF
                   !-

                ELSEIF (recruite) THEN

                   !! Calculate the number of succesful recruits
                   ! We want to avoid adding a lot of recruits which are then killed a bit 
                   ! later in the simulation. Such an approach would also be in conflict
                   ! with the approach in which we keep a constant number of diameter
                   ! and thus have to merge the biomass of the recruits with the biomass
                   ! in the smallest diameter class. For that reason we calculate the
                   ! the number of succesful recruits that will survive and will make it
                   ! into the next diameter class. This number is expected to be rather low
                   ! and cannot be compared against observations of actual recruitment. 
                   ! It should be compared against observations of in growth.

                   ! Recruitment is based on the light stress factor [0 1] computed
                   ! above. It makes use of the functional response of recruitment
                   ! to light described by Ruger et al (2009, J. of Ecol.,
                   ! doi: 10.1111/j.1365-2745.2009.01552.x) in the 50-ha plot of Barro
                   ! Colorado Island (BCI) in Panama: 
                   ! Log10(Nrecruits)= A + B * ( Log10(lstress_fac) - MU )
                   !
                   ! Where Nrecruits is the number of recruits for an area of 25 m2.
                   ! The parameter A is the intercept (mean log10 recruits expected), 
                   ! B the slope (functional response), and MU=Log10(0.02)=-1.7 
                   ! Here the recruitment response to light can be negative (B<0),
                   ! decelerating (0<B<1), accelerating (B>1) or linear (B=1). 
                   ! MU is needed here to express the fraction of light between 0 and 1
                   ! centred on mean light 2% at BCI. The linearized (log-log) power
                   ! model is implemented here for each PFT following the original setup
                   ! of Ruger et al 2009 to ensure realistic values. IT COULD EASILY BE
                   ! SIMPLIFIED. The predicted number of recruits is for an area of
                   ! 25 m2 so we have to compute the antilogarithm and divide by 25 to
                   ! get the estimate by 1 m2 which is the unit used in ORCHIDEE.
                   IF (recruitment_alpha(ivm).LT.min_stomate .AND. &
                        recruitment_beta(ivm).LT.min_stomate) THEN
                      
                      ! Sanity check. If the users sets both values to zero
                      ! new_ind is set to zero. The simulation should be 
                      ! identical to a simulation without recruitment.
                      new_ind(ipts,ivm) = zero

                   ELSE
                      
                      ! Calculate the actual number of trees before recruitment
                      tmp_ind = SUM(circ_class_n(ipts,ivm,:))
                      
                      ! Calculate the number of recruits m2 day-1
                      ! New approach in which stand density (tmp_ind) is the main driver which
                      ! is further modulated by light. This approach simulates fewer recruits
                      ! prior to canopy closure. Als the response to a thinning depends on 
                      ! the stand structure following the thinning. This approach seems a bit
                      ! as double counting because the light should be a consequence of the stand
                      ! density.
                      new_ind(ipts,ivm) = recruitment_alpha(ivm)**(LOG10(recruitment_beta(ivm)*tmp_ind)) * &
                           SQRT(LOG10(MAX(100.*lstress_fac(ipts,ivm),un)))

                   ENDIF

                   ! Update the number of individuals
                   old_ind = tmp_ind
                   tmp_ind = tmp_ind + new_ind(ipts,ivm)

                   ! Debug
                   IF(printlev_loc>=4)THEN
                      WRITE(numout,*) 'PFT, ',ivm
                      WRITE(numout,*) 'canopy cover (LAI), ',&
                           cc_to_lai(circ_class_biomass(ipts,ivm,:,ileaf,icarbon),&
                           circ_class_n(ipts,ivm,:),ivm) 
                      WRITE(numout,*) 'used canopy_gap (%) is, ', lstress_fac(ipts,ivm)*100 
                      WRITE(numout,*) 'Original ind is, ', (tmp_ind - new_ind(ipts,ivm))*m2_to_ha
                      WRITE(numout,*) 'new_ind is, ', new_ind(ipts,ivm)*m2_to_ha
                      WRITE(numout,*) 'ind+new_ind is, ', tmp_ind*m2_to_ha
                      WRITE(numout,*) 'Dg (m) ', qmd_init(ivm)
                      ind_init = Nmax(qmd_init(ivm)*m_to_cm,ivm)
                      WRITE(numout,*) 'Nmax, ',ind_init
                   ENDIF
                   !-

                   !! Calculate the dimensions of the recruits    
                   ! The height of the saplings is prescribed and determines the reserves
                   ! which are especially important for deciduous species, which need to
                   ! survive on their reserves for the first couple of months to a whole
                   ! year because the phenology scheme requires annual mean values to get
                   ! started. Calculate sapling diameter and biomass for prescribed sapling
                   ! height (above). Note that the dimensions of a newly planted sapling and 
                   ! a recruit are not the same.
                   dia_init_recr = ( recruitment_height(ivm) / pipe_tune2(ipts,ivm) ) ** &
                        ( 1. / pipe_tune3(ivm) )
                   wood_init = (recruitment_height(ivm) * pi / 4. * (dia_init_recr) ** 2. ) * &
                        pipe_density(ivm) * tree_ff(ivm)

                   ! Use the variables from the pipe model and the allometric
                   ! relationships to calculate the different biomass components
                   ! of the saplings in each diameter class.
                   CALL make_saplings(bm_sapl_2D(ipts,ivm,1,:,:), dia_init_recr, &
                        1, ivm, ipts, KF(ipts,ivm), c0_alloc(ivm), plant_status(ipts,ivm), &
                        establish, cn_leaf_init_2D(ipts,ivm),bypass_lm(1))

                   ! Debug
                   IF (printlev_loc>=4) THEN
                      WRITE(numout,*)'Prescribed recruitment height (m): ',recruitment_height(ivm)
                      WRITE(numout,*)'Recruitment diameter (m) for prescribed height: ', dia_init
                      WRITE(numout,*)'Wood_init recruitment: ',wood_init
                      WRITE(numout,*)'ind_init ',ind_init
                      WRITE(numout,*)'bm_sapl after calling make_saplings: ',bm_sapl_2D(ipts,ivm,1,:,:)
                   ENDIF
                   !-

                   ! Update the number of individuals in the first size class. All
                   ! recruitment occurs in the circ_class with the smallest diameter
                   old_circ_class_n(ipts,ivm,:) = circ_class_n(ipts,ivm,:)
                   circ_class_n(ipts,ivm,1) = circ_class_n(ipts,ivm,1) + new_ind(ipts,ivm)                        

                   !! Determine the biomass in the first circumference class.
                   ! This is a dillution approach were an average tree is made by making use
                   ! of the biomass in the available trees and the biomass in the recruits.
                   ! This is a simple mean weighted by the number of trees. The model needs
                   ! the biomass in an individual tree (gC tree-1)
                   DO iele = 1,nelements
                      old_circ_class_biomass(ipts,ivm,:,:,iele) = &
                           circ_class_biomass(ipts,ivm,:,:,iele) 
                      circ_class_biomass(ipts,ivm,1,:,iele) = &
                           (circ_class_biomass(ipts,ivm,1,:,iele) * &
                           old_circ_class_n(ipts,ivm,1) + &
                           bm_sapl_2D(ipts,ivm,1,:,iele) * new_ind(ipts,ivm) ) / &
                           circ_class_n(ipts,ivm,1)
                   ENDDO

                   ! Debug
                   IF (printlev_loc>=4) THEN
                      WRITE(numout,*) "Canopy gap (%)", lstress_fac(ipts,ivm)*100
                      WRITE(numout,*) "Original density (N/ha): ",&
                           old_circ_class_n(ipts,ivm,:)*m2_to_ha
                      WRITE(numout,*) "New density (N/ha): ",circ_class_n(ipts,ivm,:)*m2_to_ha
                      WRITE(numout,*) "Recruited saplings per ha in first size class ", &
                           new_ind(ipts,ivm) * m2_to_ha
                      WRITE(numout,*) "Original biomass per individual", &
                           SUM(old_circ_class_biomass(ipts,ivm,:,:,icarbon),2),&
                           SUM(old_circ_class_biomass(ipts,ivm,:,:,initrogen),2)
                      WRITE(numout,*) "New biomass per individual", &
                           SUM( circ_class_biomass(ipts,ivm,:,:,icarbon),2 ),&
                           SUM( circ_class_biomass(ipts,ivm,:,:,initrogen),2 )
                      WRITE(numout,*) "Biomass of new saplings gC/sapling", &
                           bm_sapl_2D(ipts,ivm,1,:,icarbon)
                       WRITE(numout,*) "Biomass of new saplings gN/sapling", &
                           bm_sapl_2D(ipts,ivm,1,:,initrogen)
                      WRITE(numout,*) "KF ", KF(ipts,ivm)
                   ENDIF
                   !-

                ENDIF ! end prescribe/recruitment block

                IF (establish) THEN

                   ! Set leaf age classes, all leaves are current year leaves
                   leaf_frac(ipts,ivm,:) = zero
                   leaf_frac(ipts,ivm,1) = un

                   ! Set time since last beginning of growing season but only
                   ! for the first day of the whole simulation. When the model
                   ! is initialized when_growthinit is set to undef. In subsequent
                   ! time steps it should have a value. For trees without phenology
                   ! the growing season starts at the moment the PFT is prescribed
                   IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
                      when_growthinit(ipts,ivm) = 200
                   ENDIF

                ENDIF

                ! The carbon and nitrogen to build the saplings is taken from the 
                ! atmosphere, keep track of amount to calculate the C-balance 
                ! closure
                IF (establish) THEN

                   ! The whole biomass of this PFT was taken from the
                   ! atmosphere
                   atm_to_bm(ipts,ivm,:) = atm_to_bm(ipts,ivm,:) + &
                        (SUM(cc_to_biomass(npts, ivm, &
                        circ_class_biomass(ipts,ivm,:,:,:), &
                        circ_class_n(ipts,ivm,:)),1) / dt )

                ELSEIF (recruite) THEN

                   ! Only the increase in biomass of this PFT was taken
                   ! from the atmosphere
                   atm_to_bm(ipts,ivm,:) = atm_to_bm(ipts,ivm,:) + &
                        (SUM(cc_to_biomass(npts, ivm, &
                        circ_class_biomass(ipts,ivm,:,:,:), &
                        circ_class_n(ipts,ivm,:)) - &
                        old_biomass(ipts,ivm,:,:),1)  / dt )
                   
                ENDIF

                ! Check whether the order was preserved. We only want to
                ! preserve the order for the carbon biomass. Nitrogen follows
                ! carbon. We expect that the biomass of an individual tree 
                ! increases with an increasing circ class. This is probably 
                ! not important for the correct functioning of the model. So,
                ! if this turns out to be computationally expensive it is 
                ! worth trying without. We try to preserve the order as
                ! an additional mean to control the quality of the simulations.
                ! A strict order also helps to interpret the output. Hence, 
                ! if the order was not preserved we will sort the biomasses. 
                ! Note that this subroutine first checks whether the order 
                ! was preseverd.
                circ_class_biomass_new = circ_class_biomass(ipts,ivm,:,:,:)
                circ_class_n_new = circ_class_n(ipts,ivm,:)
                CALL sort_circ_class_biomass(circ_class_biomass_new, &
                     circ_class_n_new)
                circ_class_biomass(ipts,ivm,:,:,:) = circ_class_biomass_new
                circ_class_n(ipts,ivm,:) = circ_class_n_new

             ENDIF ! tree(ivm)

             !! 2.3.2 Initialize grassy PFTs
             !  Use veget_max to check whether the PFT is present. If the PFT is 
             !  present but it has no biomass, prescribe its biomass (gC m-{2}). 
             !  It is assumed that at day 1 grasses have all their biomass in the 
             !  reserve pool. The criteria exclude crops. Crops are no longer 
             !  prescribed but planted the day that begin_leaves is true.
             IF (establish) THEN

                IF ( ( .NOT. is_tree(ivm) ) .AND. &
                     ( natural(ivm) ) .AND. &
                     ( veget_max(ipts,ivm) .GT. min_stomate ) .AND. &
                     ( SUM(circ_class_biomass(ipts,ivm,1,:,icarbon) ) .LE. min_stomate ) ) THEN

                   ! Debug
                   IF(printlev_loc>=4)THEN
                      WRITE(numout,*) 'We will prescribe a new grass  '//&
                           'vegetation, the old one died'
                   ENDIF
                   ! -

                   ! For grasses we assume that an individual grass is 1 m2 of grass. 
                   ! This is set in nmaxplants in pft_parameters.f90. It could be set 
                   ! to another value but with the current code this should not have 
                   ! any meaning. The grassland does not necessarily covers the whole 
                   ! 1 m2 so adjust for the canopy cover
                   tmp_ind = nmaxplants(ivm) * ha_to_m2 * canopy_cover(ivm)

                   !+++CHECK+++
                   ! We do not deal with circumference classes for grasses
                   ! and crops, but we want to keep the arrays the same as for the trees 
                   ! so we put the young grass information into the first circumference 
                   ! class. Calculate the sapwood to leaf mass in a similar way as has 
                   ! been done for trees. For trees this approach had been justified by 
                   ! observations. For grasses such justification is not supported by 
                   ! observations but we didn't try to find it. Needs more work by 
                   ! someone interested in grasses. There might be a more elegant
                   ! solution making use of a well observed parameter. 
                   ! The mass of the structural carbon relates to the mass of the leaves 
                   ! through a prescribed parameter ::k_latosa
                    !+++CHECK+++
                   ! At one point it looked like a good idea to take the max of two
                   ! options but by doing so we cannot recalculate KF in phenology.
                   ! It has not been confirmed that this is really a problem. Use the
                   ! simpelest approach but leave the alternative in the code as a
                   ! suggestion of a possible solution in case something goes wrong.
!!$                   k_latosa_tmp = MAX(k_latosa_min(ivm),k_latosa_adapt(ipts,ivm) + &
!!$                        (lstress_fac(ipts,ivm) * &
!!$                        (k_latosa_max(ivm)-k_latosa_min(ivm))))
                   !+++++++++++
                   k_latosa_tmp = k_latosa_adapt(ipts,ivm) + lstress_fac(ipts,ivm) * &
                        (k_latosa_max(ivm)-k_latosa_min(ivm))

                   IF (sla_dyn) THEN 
                      KF(ipts,ivm) = k_latosa_tmp / &
                           (slainit(ivm) * tree_ff(ivm) * pipe_density(ivm))
                   ELSE
                      KF(ipts,ivm) = k_latosa_tmp / &
                           (sla(ivm) * tree_ff(ivm) * pipe_density(ivm))
                   END IF

                   ! Calculate leaf to root area    
                   LF(ipts,ivm) = c0_alloc(ivm) * KF(ipts,ivm)

                   ! Debug
                   IF(printlev_loc>=4.AND. test_pft == ivm)THEN
                      WRITE(numout,*) 'KF, ', KF(ipts,ivm)
                      WRITE(numout,*) 'LF, c0_alloc, ', c0_alloc(ivm) * KF(ipts,ivm), &
                           c0_alloc(ivm)
                   ENDIF
                   !-
                   qmd_init(ivm) = zero
                   ! initialize everything to make sure there are not random values 
                   ! floating around for ncirc = 1
                   bm_sapl_2D(ipts,ivm,1,:,:) = val_exp

                   ! Similar as for trees, the initial height of the vegetation was 
                   ! defined. By using lai_to_biomass dynamic sla is accounted for. 
                   bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) = lai_to_biomass(&
                        height_init(ivm)/lai_to_height(ivm),ivm)

                   ! Use allometric relationships to define the root mass based on 
                   ! leaf mass. Some 'sapwood mass' is needed to store the reserves. 
                   ! An arbitrairy fraction of 5% was used 
                   bm_sapl_2D(ipts,ivm,1,iroot,icarbon) = bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) / &
                        LF(ipts,ivm)
                   bm_sapl_2D(ipts,ivm,1,isapabove,icarbon) = bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) / &
                        KF(ipts,ivm)

                   ! Some of the biomass components that exist for trees are undefined 
                   ! for grasses
                   bm_sapl_2D(ipts,ivm,1,isapbelow,icarbon) = zero
                   bm_sapl_2D(ipts,ivm,1,iheartabove,icarbon) = zero
                   bm_sapl_2D(ipts,ivm,1,iheartbelow,icarbon) = zero
                   bm_sapl_2D(ipts,ivm,1,ifruit,icarbon) = zero
                   bm_sapl_2D(ipts,ivm,1,icarbres,icarbon) = zero

                   ! Pools that are defined in the same way for trees and grasses
                   bm_sapl_2D(ipts,ivm,1,ilabile,icarbon) = labile_to_total * &
                        (bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) + &
                        fcn_root(ivm) * bm_sapl_2D(ipts,ivm,1,iroot,icarbon) + fcn_wood(ivm) * &
                        (bm_sapl_2D(ipts,ivm,1,isapabove,icarbon) + &
                        bm_sapl_2D(ipts,ivm,1,isapbelow,icarbon) + &
                        bm_sapl_2D(ipts,ivm,1,icarbres,icarbon)))

                   ! Allocate the nitrogen
                   cn_leaf=cn_leaf_init_2D(ipts,ivm) 
                   cn_wood=cn_leaf/fcn_wood(ivm) 
                   cn_root=cn_leaf/fcn_root(ivm) 

                   !+++CHECK+++
                   ! See comments at a similmar block of code for the
                   ! tree saplings
                   ! Use the C/N ratio to calculate the N content for
                   ! all the biomass components.
                   bm_sapl_2D(ipts,ivm,1,:,initrogen) = &
                        bm_sapl_2D(ipts,ivm,1,:,icarbon) / cn_root

                   ! Overwrite the N content for the leaves and wood
                   ! components. Note that only sapabove is defined the
                   ! other wood components are zero
                   bm_sapl_2D(ipts,ivm,1,ileaf,initrogen) = &
                        bm_sapl_2D(ipts,ivm,1,ileaf,icarbon) / cn_leaf
                   bm_sapl_2D(ipts,ivm,1,isapabove,initrogen) = &
                        bm_sapl_2D(ipts,ivm,1,isapabove,icarbon) / cn_wood
                   !++++++++++++

                   ! Avoid deciduous PFTs to have leaves out at establishment
                   ! Whether the saplings have leaves or don't have leaves the first 
                   ! year doesn't really matter. Either the approach is correct in the 
                   ! northern hemisphere or in the southern hemisphere. Anyhow, a 
                   ! spin-up is recommended to avoid issues with the initial conditions
                   IF ( pheno_type(ivm) .NE. 1 ) THEN

                      ! Not evergreen. Deciduous PFTs now need to survive half a year 
                      ! before bud burst. To ensure survival there are several options: 
                      ! (a) either the height of the initial vegetation is increased 
                      ! (this results in more reserves) or (b) the reserves could be
                      ! increased. The second option may result in result in numerical 
                      ! issues further down the code as the optimal reserve level is 
                      ! calculated from the other biomass pools

                      ! Notice that the sapwood is not put into reserves.  Sapwood is
                      ! important in phenology when the amount of leaves and roots
                      ! are determined, so it needs to be already present.
                      bm_sapl_2D(ipts,ivm,1,icarbres,:) = bm_sapl_2D(ipts,ivm,1,icarbres,:) + &
                           bm_sapl_2D(ipts,ivm,1,ileaf,:) + &
                           bm_sapl_2D(ipts,ivm,1,iroot,:)
                      bm_sapl_2D(ipts,ivm,1,ileaf,:) = zero
                      bm_sapl_2D(ipts,ivm,1,iroot,:) = zero

                      ! When there are no leaves, the crop/grass is dormant
                      plant_status(ipts,ivm) = idormant

                   ELSE

                      ! Initilize carbohydrate reserves for evergreen PFTs
                      bm_sapl_2D(ipts,ivm,1,icarbres,:) = zero 

                      ! Evergreen plants always have a canopy
                      plant_status(ipts,ivm) = icanopy

                   ENDIF

                   ! Debug
                   IF (printlev_loc>=4) THEN
                      DO iele=1,nelements
                         WRITE(numout,*) ' young grass biomass (gC):',ivm,ipts
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,ileaf,:,) ',&
                              bm_sapl_2D(ipts,ivm,1,ileaf,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,ispabove,:) ',&
                              bm_sapl_2D(ipts,ivm,1,isapabove,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,isapbelow,:) ',&
                              bm_sapl_2D(ipts,ivm,1,isapbelow,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,iheartabove,:) ',&
                              bm_sapl_2D(ipts,ivm,1,iheartabove,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,iheartbelow,:) ',&
                              bm_sapl_2D(ipts,ivm,1,iheartbelow,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,iroot,:) ',&
                              bm_sapl_2D(ipts,ivm,1,iroot,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,ifruit,:)', &
                              bm_sapl_2D(ipts,ivm,1,ifruit,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,icarbres,:)', &
                              bm_sapl_2D(ipts,ivm,1,icarbres,iele)
                         WRITE(numout,*) '  bm_sapl_2D(ipts,:,1,ilabile,:)', &
                              bm_sapl_2D(ipts,ivm,1,ilabile,iele)
                      ENDDO
                   ENDIF
                   ! -

                   ! circ_class_biomass (gC tree-1). This 
                   ! allows to more easily check the mass balance closure
                   circ_class_biomass(ipts,ivm,1,:,:) = bm_sapl_2D(ipts,ivm,1,:,:)
                   circ_class_n(ipts,ivm,1) = tmp_ind

                   ! Set leaf age classes -> all leaves will be current year leaves
                   leaf_frac(ipts,ivm,:) = zero
                   leaf_frac(ipts,ivm,1) = un

                   ! Set time since last beginning of growing season but only
                   ! for the first day of the whole simulation. When the model
                   ! is initialized when_growthinit is set to undef. In subsequent
                   ! time steps it should have a value.
                   IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
                      when_growthinit(ipts,ivm) = 200
                   ENDIF

                   ! The carbon and nitrogen biomass to build the saplings is taken 
                   ! from the atmosphere, keep track of amount to calculate the mass 
                   ! balance closure
                   DO iele = 1,nelements
                      atm_to_bm(ipts,ivm,iele) = atm_to_bm(ipts,ivm,iele) + &
                           ((SUM(circ_class_biomass(ipts,ivm,1,:,iele))*&
                           circ_class_n(ipts,ivm,1)) / dt)
                   ENDDO

                   ! Debug
                   IF (printlev_loc>=4) THEN
                      DO iele=1,nelements
                         WRITE(numout,*) ' root to sapwood tradeoff (LF) : ', c0_alloc(ivm)
                         WRITE(numout,*) ' grass sapling biomass: ',bm_sapl_2D(ipts,ivm,1,:,iele)
                      ENDDO
                   ENDIF
                   ! -

                ELSEIF (.NOT. natural(ivm)) THEN

                   ! Initialize croplands - else leaves_begin will never become true
                   ! Set time since last beginning of growing season but only
                   ! for the first day of the whole simulation. When the model
                   ! is initialized when_growthinit is set to undef. In subsequent
                   ! time steps it should have a value.
                   IF (when_growthinit(ipts,ivm) .EQ. undef) THEN
                      when_growthinit(ipts,ivm) = 200
                   ENDIF
                   qmd_init(ivm) = zero
                   ! Debug
                   IF (printlev_loc>=4) WRITE(numout,*) 'Initialize crops',ivm
                   !-

                ENDIF ! .NOT. tree(ivm)

             ENDIF ! establish (for .NOT. is_tree(ivm)) 

             !! 2.3.3 Declare PFT present 
             !  Now that the PFT has biomass it should be declared 'present'
             !  everywhere in that grid box. Assign some additional properties
             IF (establish) THEN
                PFTpresent(ipts,ivm) = .TRUE.
                everywhere(ipts,ivm) = un
                age(ipts,ivm) = zero
                npp_longterm(ipts,ivm) = npp_longterm_init
                lm_lastyearmax(ipts,ivm) = SUM(bm_sapl_2D(ipts,ivm,:,ileaf,icarbon) * &
                     circ_class_n(ipts,ivm,:))
             ENDIF

          ENDIF ! establish or recruite

       ENDDO ! loop over pixels

    ENDDO ! loop over PFTs


 !! 3. Check numerical consistency of this routine

    IF (err_act.GT.1) THEN

       ! 3.1 Check surface area
       CALL check_vegetation_area("stomate_prescribe", npts, veget_max_begin, &
            veget_max,'pft')

       ! 3.2 Mass balance closure 
       ! 3.2.1 Calculate final biomass
       pool_end(:,:,:) = zero 
       DO ipar = 1,nparts
          DO iele = 1,nelements
             DO icir = 1,ncirc
                pool_end(:,:,iele) = pool_end(:,:,iele) + &
                     (circ_class_biomass(:,:,icir,ipar,iele) * &
                     circ_class_n(:,:,icir) * veget_max(:,:))
             ENDDO
          ENDDO
       ENDDO

       ! 3.2.2 Calculate mass balance
       ! Common processes
       DO iele=1,nelements
          check_intern(:,:,iatm2land,iele) = check_intern(:,:,iatm2land,iele) + &
               atm_to_bm(:,:,iele) * veget_max(:,:) * dt
          check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
               pool_start(:,:,iele)) 
       ENDDO

       closure_intern(:,:,:) = zero
       DO imbc = 1,nmbcomp
          DO iele=1,nelements
             ! Debug
             IF (printlev_loc>=3) THEN              
                IF(iele .EQ. 1) WRITE(numout,"(9999(G12.5,:,','))") &
                  'imbc, test_pft, check_intern', test_grid , imbc, &
                  test_pft, check_intern(test_grid,test_pft,imbc,iele)
                ENDIF
             !-
             closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
                  check_intern(:,:,imbc,iele)
          ENDDO
       ENDDO

       ! 3.3 Check mass balance closure
       CALL check_mass_balance("stomate_prescribe", closure_intern, npts, &
            pool_end, pool_start, veget_max, 'pft')
      
    ENDIF ! err_act .GT. 1

    IF(printlev_loc>=3) WRITE(numout,*) 'Leaving stomate_prescribe.f90'

  END SUBROUTINE prescribe



!! ================================================================================================================================
!! SUBROUTINE   : make_sapling
!!
!>\BRIEF Given the dimensions and the allocation factors, calculate the biomass of a sapling        
!!
!! DESCRIPTION (functional, design, flags):
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLES(S): ::bm_sapl, ::KF
!!
!! REFERENCES   :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

   SUBROUTINE make_saplings(bm_sapl_2D_loc, tree_dia, &
                           icir, ivm, ipts, KF_loc, c0_alloc_loc, plant_status_loc, &
                           establish, cn_leaf_init_2D_loc,bypass_lm_loc)

    !! 0. Parameters and variables declaration

    !! 0.1 Input variables
    REAL(r_std), INTENT(in)                            :: tree_dia           !! Tree diameter (m)
    INTEGER(i_std), INTENT(in)                         :: ivm, icir, ipts    !! index (unitless)
    REAL(r_std), INTENT(in)                            :: c0_alloc_loc       !! Root to sapwood tradeoff parameter
    LOGICAL, INTENT(in)                                :: establish          !! yes: establish a new young vegetation. 
                                                                             !! No: add saplings under an existing vegetation
    REAL(r_std), INTENT(in)                            :: cn_leaf_init_2D_loc!! initial leaf C/N ratio 
    REAL(r_std), INTENT(in)                            :: KF_loc             !! Scaling factor to convert sapwood mass
                                                                             !! into leaf mass (m) - this variable is 
                                                                             !! passed to other routines

    !! 0.2 Output variables
    REAL(r_std), INTENT(out)                           :: bypass_lm_loc      !! leaf mass before it is set to zero
                                                                             !! for deciduous PFTs. This value is
                                                                             !! in the calculation of KF (gC tree-1)

    !! 0.3 Modified variables
    REAL(r_std), DIMENSION(:,:),INTENT(inout)          :: bm_sapl_2D_loc     !! Sapling biomass for the functional 
    REAL, INTENT(inout)                                :: plant_status_loc   !! Growth and phenological status of the plant 
                                                                             !! Different stati defined in constants
    !! 0.4  Local variables
    INTEGER(i_std)                                    :: iele                !! index (unitless) 
    REAL(r_std)                                       :: cn_leaf,cn_root     !! CN ratio of leaves and roots (gC/gN)
    REAL(r_std)                                       :: cn_wood             !! CN ratio of wood pool (gC/gN)        
    REAL(r_std)                                       :: tree_height         !! Tree height (m)

!_ ================================================================================================================================

      ! The woody biomass is contained in four components. Thus, 
      ! wood_init = isapabove + isapbelow + iheartabove + iheartbelow. 
      ! Given that isapbelow = isapbelow and iheartabove = 
      ! iheartbelow =  bm_sapl_heartabove*isapabove. If 
      ! bm_sapl_heartbelow = bm_sapl_heartabove = 0.2, then 
      ! isapabove = wood_init/2.4
          
      IF ((dia_init_max(ivm)+dia_init_min(ivm))/2 .GT. 0.05) THEN
        CALL ipslerr_p(3,'Initial diameter of sapling was set to more than 5cm',&
        'The current estimate of a sapling biomass does not account &
        for heartwood','If you want to use larger diameters add a &
        heartwood estimate in the code','')
        ! bm_sapl_2D_loc(isapabove,icarbon) =
        ! tree_dia**(2+pipe_tune3(ivm))*pi/4* &
        !   2*tree_ff(ivm)*pipe_density(ivm)*pipe_tune2(ipts,ivm) / &
        !  (2. + bm_sapl_heartabove + bm_sapl_heartbelow)
        !bm_sapl_2D_loc(isapbelow,icarbon) = bm_sapl_2D_loc(isapabove,icarbon)
        !bm_sapl_2D_loc(iheartabove,icarbon) =  bm_sapl_heartabove * &
        !   bm_sapl_2D_loc(isapabove,icarbon)
        !bm_sapl_2D_loc(iheartbelow,icarbon) =  bm_sapl_heartbelow * &
        !   bm_sapl_2D_loc(isapbelow,icarbon)

      ELSE

         bm_sapl_2D_loc(iheartabove,icarbon) =  0.0001
         bm_sapl_2D_loc(iheartbelow,icarbon) =  0.0001
         bm_sapl_2D_loc(isapabove,icarbon) = tree_dia**(2+pipe_tune3(ivm))*pi/4* &
            tree_ff(ivm)*pipe_density(ivm)*pipe_tune2(ipts,ivm) - &
            bm_sapl_2D_loc(iheartabove,icarbon)
         bm_sapl_2D_loc(isapbelow,icarbon) = bm_sapl_2D_loc(isapabove,icarbon)

      ENDIF  

      ! Use the allometric relationships to calculate initial leaf 
      ! and root mass. Note that wstress should be accounted for
      ! through c0_alloc. There was an inconsistency in the
      ! original calculation (OCN) - most pools were in gN but 
      ! leaves were in gC. A correction is proposed. Given that 
      ! icarbes was just set to zero the equation for ilabile is 
      ! a bit overkill and icarbes could be omitted.
      tree_height =  ((((bm_sapl_2D_loc(isapabove,icarbon)+ &
           bm_sapl_2D_loc(isapbelow,icarbon)+bm_sapl_2D_loc(iheartabove,icarbon)+ &
           bm_sapl_2D_loc(iheartbelow,icarbon))/(tree_ff(ivm)* &
           pipe_density(ivm))*4/pi)**(pipe_tune3(ivm)/2)) &
               *pipe_tune2(ipts,ivm))**(1/(pipe_tune3(ivm)/2+1))

      bm_sapl_2D_loc(ileaf,icarbon) = ( bm_sapl_2D_loc(isapabove,icarbon) + &
           bm_sapl_2D_loc(isapbelow,icarbon) ) * KF_loc /  tree_height 
      bm_sapl_2D_loc(iroot,icarbon) = bm_sapl_2D_loc(ileaf,icarbon) / &
           ( KF_loc * c0_alloc_loc )
      bm_sapl_2D_loc(icarbres,icarbon)=zero
      bm_sapl_2D_loc(ifruit,icarbon) = zero
      bm_sapl_2D_loc(ilabile,icarbon) =  labile_to_total * &
           (bm_sapl_2D_loc(ileaf,icarbon)  + &
           fcn_root(ivm) * bm_sapl_2D_loc(iroot,icarbon) + &
           fcn_wood(ivm) * (bm_sapl_2D_loc(isapabove,icarbon) + &
           bm_sapl_2D_loc(isapbelow,icarbon) + &
           bm_sapl_2D_loc(icarbres,icarbon)))      

      ! Debug
      IF (printlev_loc>=4) THEN 
         WRITE(numout,*) ' +++++ CHECK INSIDE make_sapling ++++ '
         WRITE(numout,*) 'Circumference class: ',icir,' PFT type: ',ivm
         WRITE(numout,*) ' root to sapwood tradeoff p :', c0_alloc_loc
         WRITE(numout,*) ' sapling height, sapling diameter ', &
              tree_height , ( tree_height/pipe_tune2(ipts,ivm) ) ** &
              ( 1. / pipe_tune3(ivm) )
         WRITE(numout,*) 'pipe_density, ',pipe_density(ivm)
      ENDIF
      ! -

      ! Allocate the nitrogen
      cn_leaf=cn_leaf_init_2D_loc
      cn_wood=cn_leaf/fcn_wood(ivm)
      cn_root=cn_leaf/fcn_root(ivm)

      ! Use the C/N ratio to calculate the N content for
      ! all the biomass components.
      bm_sapl_2D_loc(:,initrogen) = &
           bm_sapl_2D_loc(:,icarbon) / cn_root

      ! Overwrite the N content for the leaves and wood
      ! components.
      bm_sapl_2D_loc(ileaf,initrogen) = &
           bm_sapl_2D_loc(ileaf,icarbon) / cn_leaf
      bm_sapl_2D_loc(isapabove,initrogen) = &
           bm_sapl_2D_loc(isapabove,icarbon) / cn_wood
      bm_sapl_2D_loc(isapbelow,initrogen) = &
           bm_sapl_2D_loc(isapbelow,icarbon) / cn_wood
      bm_sapl_2D_loc(iheartabove,initrogen) = &
           bm_sapl_2D_loc(iheartabove,icarbon) / cn_wood
      bm_sapl_2D_loc(iheartbelow,initrogen) = &
           bm_sapl_2D_loc(iheartbelow,icarbon) / cn_wood

      ! Store the leaf mass so it can be used in the calculation of KF
      ! when using a dynamic sla.
      bypass_lm_loc = bm_sapl_2D_loc(ileaf,icarbon)

      ! Avoid deciduous PFTs to have leaves out at establishment
      ! Whether the saplings have leaves or don't have leaves the
      ! first year doesn't really matter. Either the approach is
      ! correct in the northern hemisphere or in the southern
      ! hemisphere. Anyhow, a spin-up is needed to avoid issues
      ! with the initial conditions
      IF ( pheno_type(ivm) .NE. 1 ) THEN
         
         ! Not evergreen. Deciduous PFTs now need to survive an extra
         ! year before bud burst. To ensure survival there are several
         ! options: (a) either the height of the initial vegetation is
         ! increased (this results in more reserves) or (b) the reserves
         ! could be increased. The second option may result in numerical
         ! issues further down the code as the optimal reserve level is
         ! calculated from the other biomass pools. Also some initial
         ! tests showed that higher reserves simply resulted in more
         ! respiration.
         ! Tip: use taller trees to start with if they die between
         ! prescribe and leaf onset
         bm_sapl_2D_loc(icarbres,:) = (bm_sapl_2D_loc(icarbres,:) + &
              bm_sapl_2D_loc(ileaf,:) + bm_sapl_2D_loc(iroot,:))
         bm_sapl_2D_loc(ileaf,:) = zero
         bm_sapl_2D_loc(iroot,:) = zero
      
         ! When deciduous trees have no leaves they are dormant
         ! If the code is used for recruitment use the actual 
         ! plant_status of the overstorey trees for the saplings
         IF (establish) THEN
            plant_status_loc = idormant
         ENDIF
         
      ELSE
         
         ! Initilize carbohydrate reserves for evergreen PFTs
         bm_sapl_2D_loc(icarbres,:) = zero

         ! Evergreen trees always have a canopy. If the code
         ! is used for recruitment use the actual plant_status 
         ! of the overstorey trees for the saplings
         IF (establish) THEN
            plant_status_loc = icanopy
         ENDIF
    
      ENDIF

      ! Debug
      IF (printlev_loc>=4) THEN
         DO iele=1,nelements
            ! This is the only use of icir, ivm, and ipts in this routine.
            ! This could be placed outside the routine, but it's a bit
            ! cleaner here.
            WRITE(numout,*) ' sapling biomass (gC):',icir,ivm,ipts
            WRITE(numout,*) '  bm_sapl_2D_loc(ileaf,:,) ',&
                 bm_sapl_2D_loc(ileaf,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(ispabove,:) ',&
                 bm_sapl_2D_loc(isapabove,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(isapbelow,:) ',&
                 bm_sapl_2D_loc(isapbelow,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(iheartabove,:) ',&
                 bm_sapl_2D_loc(iheartabove,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(iheartbelow,:) ',&
                 bm_sapl_2D_loc(iheartbelow,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(iroot,:) ',&
                 bm_sapl_2D_loc(iroot,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(ifruit,:)', &
                 bm_sapl_2D_loc(ifruit,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(icarbres,:)', &
                 bm_sapl_2D_loc(icarbres,iele)
            WRITE(numout,*) '  bm_sapl_2D_loc(ilabile,:)', &
                 bm_sapl_2D_loc(ilabile,iele)
         ENDDO
      ENDIF
      
   END SUBROUTINE make_saplings
  
END MODULE stomate_prescribe
