! Phenology
!
! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate_phenology.f90,v 1.11 2010/04/06 16:06:34 ssipsl Exp $
! IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
MODULE stomate_phenology

  ! modules used:

  USE ioipsl
  USE stomate_data
  USE constantes
  USE pft_parameters

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC phenology,phenology_clear

  ! first call
  LOGICAL, SAVE                                              :: firstcall = .TRUE.
  LOGICAL, SAVE                                              :: firstcall_hum = .TRUE.
  LOGICAL, SAVE                                              :: firstcall_moi = .TRUE.
  LOGICAL, SAVE                                              :: firstcall_humgdd = .TRUE.
  LOGICAL, SAVE                                              :: firstcall_moigdd = .TRUE.

CONTAINS

  SUBROUTINE phenology_clear
    firstcall=.TRUE.
    firstcall_hum=.TRUE.
    firstcall_moi = .TRUE.
    firstcall_humgdd = .TRUE.
    firstcall_moigdd = .TRUE.
  END SUBROUTINE phenology_clear

  SUBROUTINE phenology (npts, dt, PFTpresent, &
       veget_max, &
       tlong_ref, t2m_month, t2m_week, gpp, &
       maxmoiavail_lastyear, minmoiavail_lastyear, &
       moiavail_month, moiavail_week, &
       gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
       senescence, time_lowgpp, time_hum_min, &
       biomass, leaf_frac, leaf_age, &
       when_growthinit, co2_to_bm, lai)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                        :: npts
    ! time step in days
    REAL(r_std), INTENT(in)                                     :: dt
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)                  :: PFTpresent
    ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: veget_max
    ! "long term" 2 meter reference temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: tlong_ref
    ! "monthly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: t2m_month
    ! "weekly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: t2m_week
    ! daily gross primary productivity (gC/(m**2 of ground)/day)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gpp
    ! last year's maximum moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: maxmoiavail_lastyear
    ! last year's minimum moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: minmoiavail_lastyear
    ! "monthly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_month
    ! "weekly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: moiavail_week
    ! growing degree days, threshold -5 deg C
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gdd_m5_dormance
    ! growing degree days, since midwinter
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: gdd_midwinter
    ! number of chilling days since leaves were lost
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ncd_dormance
    ! number of growing days, threshold -5 deg C
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: ngd_minus5
    ! is the plant senescent? (only for deciduous trees - carbohydrate reserve)
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)                  :: senescence
    ! duration of dormance (d)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: time_lowgpp
    ! time elapsed since strongest moisture availability (d)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: time_hum_min

    ! 0.2 modified fields

    ! biomass (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)     :: biomass
    ! fraction of leaves in leaf age class
    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_frac
    ! leaf age (days)
    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_age
    ! how many days ago was the beginning of the growing season
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: when_growthinit
    ! co2 taken up (gC/(m**2 of total ground)/day)
    !NV passge 2D
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                 :: co2_to_bm

    ! 0.3 output

    ! leaf area index
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: lai

    ! 0.4 local

    ! are we allowed to decalre the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm)                              :: allow_initpheno
    ! biomass we would like to have
    REAL(r_std), DIMENSION(npts)                                :: bm_wanted
    ! biomass we use (from carbohydrate reserve or from atmosphere)
    REAL(r_std), DIMENSION(npts)                                :: bm_use
    ! minimum leaf mass (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts)                                :: lm_min
    ! does the leaf age distribution have to be reset?
    LOGICAL(r_std), DIMENSION(npts)                             :: age_reset
    ! indices
    INTEGER(i_std)                                              :: i,j,m
    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm)                              :: begin_leaves

    REAL(r_std), DIMENSION(npts,nvm)                          :: histvar

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering phenology'

    !
    ! 1 first call
    !

    IF ( firstcall ) THEN

       WRITE(numout,*) 'phenology:'

       WRITE(numout,*) '   > take carbon from atmosphere if carbohydrate' // &
            ' reserve too small: ', always_init

       WRITE(numout,*) '   > minimum time since last beginning of a growing' // &
            ' season (d): ', min_growthinit_time

       firstcall = .FALSE.

    ENDIF

    !
    ! 2 various things
    !

    !
    ! 2.1 allow detection of the beginning of the growing season if dormance was
    !     long enough and last beginning of growing season was a sufficiently 
    !     long time ago
    !

    allow_initpheno(:,ibare_sechiba) = .FALSE. 
    DO j = 2,nvm

       WHERE ( ( time_lowgpp(:,j) .GE. lowgpp_time(j) ) .AND. &
            ( when_growthinit(:,j) .GT. min_growthinit_time )          )
          allow_initpheno(:,j) = .TRUE.
       ELSEWHERE
          allow_initpheno(:,j) = .FALSE.
       ENDWHERE

    ENDDO

    WHERE(allow_initpheno)
       histvar=un
    ELSEWHERE
       histvar=zero
    ENDWHERE
    CALL histwrite (hist_id_stomate, 'ALLOW_INITPHENO', itime, histvar, npts*nvm, horipft_index)

    !
    ! 2.2 increase counter: how many days ago was the beginning of the growing season
    !     Needed for allocation
    !

    when_growthinit(:,:) = when_growthinit(:,:) + dt

    !
    ! 3 Check biometeorological conditions
    !

    ! default: phenology does not start
    begin_leaves(:,:) = .FALSE.

    ! different kinds of phenology
    
    ! used in all the differents models of phenology DS 17112010
    t_always = ZeroCelsius + t_always_add

    DO j = 2,nvm

       SELECT CASE ( pheno_model(j) )

       CASE ( 'hum' )

          CALL pheno_hum (npts, j, PFTpresent, allow_initpheno, &
               moiavail_month, moiavail_week, &
               maxmoiavail_lastyear, minmoiavail_lastyear, &
               begin_leaves)

       CASE ( 'moi' )

          CALL pheno_moi (npts, j, PFTpresent, allow_initpheno, &
               time_hum_min, &
               moiavail_month, moiavail_week, &
               begin_leaves)


       CASE ( 'ncdgdd' )

          CALL pheno_ncdgdd (npts, j, PFTpresent, allow_initpheno, &
               ncd_dormance, gdd_midwinter, &
               t2m_month, t2m_week, begin_leaves)

       CASE ( 'ngd' )

          CALL pheno_ngd (npts, j, PFTpresent, allow_initpheno, ngd_minus5, &
               t2m_month, t2m_week, begin_leaves)

       CASE ( 'humgdd' )

          CALL pheno_humgdd (npts, j, PFTpresent, allow_initpheno, gdd_m5_dormance, &
               maxmoiavail_lastyear, minmoiavail_lastyear, &
               tlong_ref, t2m_month, t2m_week, &
               moiavail_week, moiavail_month, &
               begin_leaves)

       CASE ( 'moigdd' )

          CALL pheno_moigdd (npts, j, PFTpresent, allow_initpheno, gdd_m5_dormance, &
               time_hum_min, &
               tlong_ref, t2m_month, t2m_week, &
               moiavail_week, moiavail_month, &
               begin_leaves)

       CASE ( 'none' )

          ! no action

       CASE default

          WRITE(numout,*) 'phenology: don''t know how to treat this PFT.'
          WRITE(numout,*) '  number:',j
          WRITE(numout,*) '  phenology model: ',pheno_model(j)
          STOP

       END SELECT

    ENDDO

    WHERE(begin_leaves)
       histvar=un
    ELSEWHERE
       histvar=zero
    ENDWHERE
    CALL histwrite (hist_id_stomate, 'BEGIN_LEAVES', itime, histvar, npts*nvm, horipft_index)

    !
    ! 4 leaves start to grow if meteorological conditions are favourable and if
    !   leaf regrowth is allowed (cf also turnover)
    !

    DO j = 2,nvm

       age_reset(:) = .FALSE.

       DO i = 1, npts

          IF ( begin_leaves(i,j) ) THEN

             lm_min(i) = lai_initmin(j) / sla(j)

             ! do we have to put a minimum biomass into the leaves?

             IF ( biomass(i,j,ileaf) .LT. lm_min(i) ) THEN

                !
                ! 4.1 determine how much biomass we can use
                !

                bm_wanted(i) = 2. * lm_min(i)

                ! eventually take the missing carbon from the atmosphere and
                ! put it into carbohydrate reserve

                IF ( always_init .AND. ( biomass(i,j,icarbres) .LT. bm_wanted(i) ) ) THEN
                   !NV passage 2D
                   co2_to_bm(i,j) = co2_to_bm(i,j) + ( bm_wanted(i) - biomass(i,j,icarbres) ) / dt

                   biomass(i,j,icarbres) = bm_wanted(i)

                ENDIF

                bm_use(i) = MIN( biomass(i,j,icarbres), bm_wanted(i) )

                !
                ! 4.2 dispatch that biomass on leaves and roots
                !

                biomass(i,j,ileaf) = biomass(i,j,ileaf) + bm_use(i) / 2.

                biomass(i,j,iroot) = biomass(i,j,iroot) + bm_use(i) / 2.

                !
                ! 4.3 decrease reservoir biomass
                !

                biomass(i,j,icarbres) = biomass(i,j,icarbres) - bm_use(i)

                !
                ! 4.4 decide whether we have to reset then leaf age distribution
                !     (done later for better vectorization)
                !

                age_reset(i) = .TRUE.

             ENDIF  ! leaf mass is very low

             !
             ! 4.5 reset counter: start of the growing season
             !

             when_growthinit(i,j) = 0.0

          ENDIF    ! start of the growing season

       ENDDO      ! loop over grid points

       !
       ! 4.6 reset leaf age distribution where necessary
       !     simply say that everything is in the youngest age class
       !

       ! 4.6.1 fractions

       WHERE ( age_reset(:) )
          leaf_frac(:,j,1) = un
       ENDWHERE
       DO m = 2, nleafages
          WHERE ( age_reset(:) )
             leaf_frac(:,j,m) = zero
          ENDWHERE
       ENDDO

       ! 4.6.2 ages

       DO m = 1, nleafages
          WHERE ( age_reset(:) )
             leaf_age(:,j,m) = zero
          ENDWHERE
       ENDDO

    ENDDO        ! loop over PFTs


    IF (bavard.GE.4) WRITE(numout,*) 'Leaving phenology'

  END SUBROUTINE phenology

  !
  ! ==============================================================================
  ! Phenology: begins if "weekly" soil humidity starts to exceed a certain threshold
  !            value. This value depends on last year's max and min humidity ...
  !            Always initiate growing season if soil moisture exceeds a certain threshold.
  !

  SUBROUTINE pheno_hum (npts, j, PFTpresent, allow_initpheno, &
       moiavail_month, moiavail_week, &
       maxmoiavail_lastyear, minmoiavail_lastyear, &
       begin_leaves)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! PFT index
    INTEGER(i_std), INTENT(in)                               :: j
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: PFTpresent
    ! are we allowed to decalre the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: allow_initpheno
    ! "monthly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_month
    ! "weekly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_week
    ! last year's maximum moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: maxmoiavail_lastyear
    ! last year's minimum moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: minmoiavail_lastyear

    ! 0.2 output

    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)              :: begin_leaves

    ! 0.3 local

    REAL(r_std)                                              :: moiavail_always
    ! first call
    REAL(r_std), DIMENSION(npts)                             :: availability_crit
    ! index
    INTEGER(i_std)                                           :: i

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering hum'

    !
    ! Initializations
    !

    !
    ! 1.1 messages
    !

    IF ( firstcall_hum ) THEN

       WRITE(numout,*) 'pheno_hum:'
       WRITE(numout,*) '   > moisture availability above which moisture tendency doesn''t matter: '
       WRITE(numout,*) '         trees:', moiavail_always_tree
       WRITE(numout,*) '         grasses:', moiavail_always_grass

       firstcall_hum = .FALSE.

    ENDIF

    !
    ! 1.2 initialize output
    !

    begin_leaves(:,j) = .FALSE.

    !
    ! 1.3 check the prescribed critical value
    !

    IF ( hum_frac(j) .EQ. undef ) THEN

       WRITE(numout,*) 'hum: hum_frac is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    !
    ! 1.4 critical moisture availability above which we always detect the beginning of the
    !     growing season.
    !

    IF ( tree(j) ) THEN
       moiavail_always = moiavail_always_tree
    ELSE
       moiavail_always = moiavail_always_grass
    ENDIF

    !
    ! 2 PFT has to be there and start of growing season must be allowed
    !

    DO i = 1, npts

       IF ( PFTpresent(i,j) .AND. allow_initpheno(i,j) ) THEN

          ! critical availability: depends on last year's max and min.

          availability_crit(i) = minmoiavail_lastyear(i,j) + hum_frac(j) * &
               ( maxmoiavail_lastyear(i,j) - minmoiavail_lastyear(i,j) )

          ! the favorable season starts if the "monthly" moisture availability is still quite
          ! low, but the "weekly" availability is already higher (as it reacts faster).
          ! If monthly moisture availability is high enough, also initiate growing season if
          ! this has not happened yet.

          IF ( ( ( moiavail_week(i,j)  .GE. availability_crit(i) ) .AND. &
               ( moiavail_month(i,j) .LT. moiavail_week(i,j) )   ) .OR. &
               ( moiavail_month(i,j) .GE. moiavail_always )                ) THEN
             begin_leaves(i,j) = .TRUE.
          ENDIF

       ENDIF        ! PFT there and start of growing season allowed

    ENDDO

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving hum'

  END SUBROUTINE pheno_hum

  !
  ! ==============================================================================
  ! Phenology: begins if moisture minium was a sufficiently long time ago.
  !            Additionally, "weekly" soil humidity must be higher that "monthly" soil
  !            humidity.
  !

  SUBROUTINE pheno_moi (npts, j, PFTpresent, allow_initpheno, &
       time_hum_min, &
       moiavail_month, moiavail_week, &
       begin_leaves)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! PFT index
    INTEGER(i_std), INTENT(in)                               :: j
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: PFTpresent
    ! are we allowed to decalre the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: allow_initpheno
    ! time elapsed since strongest moisture availability (d)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: time_hum_min
    ! "monthly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_month
    ! "weekly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_week

    ! 0.2 output

    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)              :: begin_leaves

    ! 0.3 local

    ! moisture availability above which moisture tendency doesn't matter
    REAL(r_std)                                              :: moiavail_always
    ! index
    INTEGER(i_std)                                           :: i

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering moi'

    !
    ! Initializations
    !

    !
    ! 1.1 messages
    !

    IF ( firstcall_moi ) THEN

       WRITE(numout,*) 'pheno_moi:'
       WRITE(numout,*) '   > moisture availability above which moisture tendency doesn''t matter: '
       WRITE(numout,*) '         trees:', moiavail_always_tree
       WRITE(numout,*) '         grasses:', moiavail_always_grass

       firstcall_moi = .FALSE.

    ENDIF

    !
    ! 1.2 initialize output
    !

    begin_leaves(:,j) = .FALSE.

    !
    ! 1.3 check the prescribed critical value
    !

    IF ( hum_min_time(j) .EQ. undef ) THEN

       WRITE(numout,*) 'moi: hum_min_time is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    !
    ! 1.4 critical moisture availability above which we always detect the beginning of the
    !     growing season.
    !

    IF ( tree(j) ) THEN
       moiavail_always = moiavail_always_tree
    ELSE
       moiavail_always = moiavail_always_grass
    ENDIF

    !
    ! 2 PFT has to be there and start of growing season must be allowed
    !

    DO i = 1, npts

       IF ( PFTpresent(i,j) .AND. allow_initpheno(i,j) ) THEN

          ! the favorable season starts if the moisture minimum was a sufficiently long 
          ! time ago and if the "monthly" moisture availability is lower than the "weekly"
          ! availability (this means that soil moisture is increasing).
          ! If monthly moisture availability is high enough, also initiate growing season if
          ! this has not happened yet.

          IF  ( ( ( moiavail_week(i,j) .GT. moiavail_month(i,j) ) .AND. &
               ( time_hum_min(i,j) .GT. hum_min_time(j) )    ) .OR. &
               ( moiavail_month(i,j) .GE. moiavail_always )                     ) THEN
             begin_leaves(i,j) = .TRUE.
          ENDIF

       ENDIF        ! PFT there and start of growing season allowed

    ENDDO

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving moi'

  END SUBROUTINE pheno_moi

  !
  ! ==============================================================================
  ! Phenology: leaves are put on if gdd exceeds a critical value.
  !            Additionally, there has to be at least some moisture.
  !            Set gdd to undef if beginning of the growing season detected.
  !

  SUBROUTINE pheno_humgdd (npts, j, PFTpresent, allow_initpheno, gdd, &
       maxmoiavail_lastyear, minmoiavail_lastyear, &
       tlong_ref, t2m_month, t2m_week, &
       moiavail_week, moiavail_month, &
       begin_leaves)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! PFT index
    INTEGER(i_std), INTENT(in)                               :: j
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: PFTpresent
    ! are we allowed to decalre the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: allow_initpheno
    ! growing degree days, calculated since leaves have fallen
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: gdd
    ! last year's maximum moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: maxmoiavail_lastyear
    ! last year's minimum moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: minmoiavail_lastyear
    ! "long term" 2 meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: tlong_ref
    ! "monthly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_month
    ! "weekly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_week
    ! "weekly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_week
    ! "monthly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_month

    ! 0.2 output

    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)              :: begin_leaves

    ! 0.3 local

    ! moisture availability above which moisture tendency doesn't matter
    REAL(r_std)                                              :: moiavail_always
    ! critical moisture availability
    REAL(r_std), DIMENSION(npts)                             :: moiavail_crit
    ! long term temperature, C
    REAL(r_std), DIMENSION(npts)                             :: tl
    ! critical GDD
    REAL(r_std), DIMENSION(npts)                             :: gdd_crit
    ! index
    INTEGER(i_std)                                           :: i

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering humgdd'

    !
    ! 1 Initializations
    !

    !
    ! 1.1 messages
    !

    IF ( firstcall_humgdd ) THEN

       WRITE(numout,*) 'pheno_humgdd:'
       WRITE(numout,*) '   > moisture availability above which moisture tendency doesn''t matter: '
       WRITE(numout,*) '         trees:', moiavail_always_tree
       WRITE(numout,*) '         grasses:', moiavail_always_grass
       WRITE(numout,*) '   > monthly temp. above which temp. tendency doesn''t matter: ', &
            t_always

       firstcall_humgdd = .FALSE.

    ENDIF

    !
    ! 1.2 initialize output
    !

    begin_leaves(:,j) = .FALSE.

    !
    ! 1.3 check the prescribed critical values
    !

    IF ( ANY(pheno_gdd_crit(j,:) .EQ. undef) ) THEN

       WRITE(numout,*) 'humgdd:pheno_gdd_crit is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    IF ( hum_frac(j) .EQ. undef ) THEN

       WRITE(numout,*) 'humgdd: hum_frac is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    !
    ! 1.4 critical moisture availability above which we always detect the beginning of the
    !     growing season.
    !

    IF ( tree(j) ) THEN
       moiavail_always = moiavail_always_tree
    ELSE
       moiavail_always = moiavail_always_grass
    ENDIF

    !
    ! 2 PFT has to be there, start of growing season must be allowed, 
    !   and gdd has to be defined
    !

    DO i = 1, npts

       IF ( PFTpresent(i,j) .AND. allow_initpheno(i,j) .AND. &
            ( gdd(i,j) .NE. undef )                           ) THEN

          ! is critical gdd reached and is temperature increasing?
          ! be sure that at least some humidity

          moiavail_crit(i) = minmoiavail_lastyear(i,j) + hum_frac(j) * &
               ( maxmoiavail_lastyear(i,j) - minmoiavail_lastyear(i,j) )

          tl(i) = tlong_ref(i) - ZeroCelsius
          gdd_crit(i) = pheno_gdd_crit(j,1) + tl(i)*pheno_gdd_crit(j,2) + &
               tl(i)*tl(i)*pheno_gdd_crit(j,3)

          IF ( ( gdd(i,j) .GE. gdd_crit(i) ) .AND. &
               ( ( t2m_week(i) .GT. t2m_month(i) ) .OR. &
               ( t2m_month(i) .GT. t_always )          ) .AND. &
               ( ( ( moiavail_week(i,j)  .GE. moiavail_crit(i) ) .AND. &
               ( moiavail_month(i,j) .LT. moiavail_crit(i) )        ) .OR. &
               ( moiavail_month(i,j) .GE. moiavail_always )                   ) )  THEN
             begin_leaves(i,j) = .TRUE.
          ENDIF

       ENDIF        ! PFT there and start of growing season allowed

    ENDDO

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving humgdd'

  END SUBROUTINE pheno_humgdd

  !
  ! ==============================================================================
  ! Phenology: leaves are put on if gdd exceeds a critical value.
  !            Additionally, a certain time must have elapsed since the moisture minimum.
  !            Set gdd to undef if beginning of the growing season detected.
  !

  SUBROUTINE pheno_moigdd (npts, j, PFTpresent, allow_initpheno, gdd, &
       time_hum_min, &
       tlong_ref, t2m_month, t2m_week, &
       moiavail_week, moiavail_month, &
       begin_leaves)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! PFT index
    INTEGER(i_std), INTENT(in)                               :: j
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: PFTpresent
    ! are we allowed to decalre the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: allow_initpheno
    ! growing degree days, calculated since leaves have fallen
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: gdd
    ! time elapsed since strongest moisture availability (d)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: time_hum_min
    ! "long term" 2 meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: tlong_ref
    ! "monthly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_month
    ! "weekly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_week
    ! "weekly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_week
    ! "monthly" moisture availability
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: moiavail_month

    ! 0.2 output

    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)              :: begin_leaves

    ! 0.3 local

    ! moisture availability above which moisture tendency doesn't matter
    REAL(r_std)                                              :: moiavail_always
    ! long term temperature, C
    REAL(r_std), DIMENSION(npts)                             :: tl
    ! critical GDD
    REAL(r_std), DIMENSION(npts)                             :: gdd_crit
    ! index
    INTEGER(i_std)                                           :: i

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering moigdd'

    !
    ! 1 Initializations
    !

    !
    ! 1.1 messages
    !

    IF ( firstcall_moigdd ) THEN

       WRITE(numout,*) 'pheno_moigdd:'
       WRITE(numout,*) '   > moisture availability above which moisture tendency doesn''t matter: '
       WRITE(numout,*) '         trees:', moiavail_always_tree
       WRITE(numout,*) '         grasses:', moiavail_always_grass
       WRITE(numout,*) '   > monthly temp. above which temp. tendency doesn''t matter: ', &
            t_always

       firstcall_moigdd = .FALSE.

    ENDIF

    !
    ! 1.2 initialize output
    !

    begin_leaves(:,j) = .FALSE.

    !
    ! 1.3 check the prescribed critical values
    !

    IF ( ANY(pheno_gdd_crit(j,:) .EQ. undef) ) THEN

       WRITE(numout,*) 'moigdd: pheno_gdd_crit is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    IF ( hum_min_time(j) .EQ. undef ) THEN

       WRITE(numout,*) 'moigdd: hum_min_time is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    !
    ! 1.4 critical moisture availability above which we always detect the beginning of the
    !     growing season.
    !

    IF ( tree(j) ) THEN
       moiavail_always = moiavail_always_tree
    ELSE
       moiavail_always = moiavail_always_grass
    ENDIF

    !
    ! 2 PFT has to be there, start of growing season must be allowed, 
    !   and gdd has to be defined
    !

    DO i = 1, npts

       IF ( PFTpresent(i,j) .AND. allow_initpheno(i,j) .AND. &
            ( gdd(i,j) .NE. undef )                           ) THEN

          ! is critical gdd reached and is temperature increasing?
          ! has enough time gone by since moisture minimum and is moisture increasing?

          tl(i) = tlong_ref(i) - ZeroCelsius
          gdd_crit(i) = pheno_gdd_crit(j,1) + tl(i)*pheno_gdd_crit(j,2) + &
               tl(i)*tl(i)*pheno_gdd_crit(j,3)

          IF ( ( gdd(i,j) .GE. gdd_crit(i) ) .AND. &
               ( ( t2m_week(i) .GT. t2m_month(i) ) .OR. &
               ( t2m_month(i) .GT. t_always )          ) .AND. &
               ( ( ( time_hum_min(i,j) .GT. hum_min_time(j) ) .AND. &
               ( moiavail_week(i,j) .GT. moiavail_month(i,j) )            ) .OR. &
               ( moiavail_month(i,j) .GE. moiavail_always )                         ) )  THEN
             begin_leaves(i,j) = .TRUE.
          ENDIF

       ENDIF        ! PFT there and start of growing season allowed

    ENDDO

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving moigdd'

  END SUBROUTINE pheno_moigdd


  !
  ! ==============================================================================
  ! Phenology: leaves are put on if a certain relationship between ncd since leaves were 
  !            lost (number of chilling days) and gdd since midwinter (growing degree 
  !            days) is fulfilled
  !

  SUBROUTINE pheno_ncdgdd (npts, j, PFTpresent, allow_initpheno, &
       ncd_dormance, gdd_midwinter, &
       t2m_month, t2m_week, begin_leaves)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! PFT index
    INTEGER(i_std), INTENT(in)                               :: j
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: PFTpresent
    ! are we allowed to declare the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: allow_initpheno
    ! number of chilling days since leaves were lost
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: ncd_dormance
    ! growing degree days since midwinter
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: gdd_midwinter
    ! "monthly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_month
    ! "weekly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_week

    ! 0.2 output

    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)              :: begin_leaves

    ! 0.3 local

    ! index
    INTEGER(i_std)                                           :: i
    ! critical gdd
    REAL(r_std)                                              :: gdd_min

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering ncdgdd'

    !
    ! 1 Initializations
    !

    !
    ! 1.1 initialize output
    !

    begin_leaves(:,j) = .FALSE.

    !
    ! 1.2 check the prescribed critical values
    !

    IF ( ncdgdd_temp(j) .EQ. undef ) THEN

       WRITE(numout,*) 'ncdgdd: ncdgdd_temp is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    !
    ! 2 PFT has to be there and start of growing season must be allowed
    !

    DO i = 1, npts

       IF ( PFTpresent(i,j) .AND. allow_initpheno(i,j) .AND. &
            ( gdd_midwinter(i,j) .NE. undef ) .AND. &
            ( ncd_dormance(i,j) .NE. undef )                  ) THEN

          ! critical gdd

          gdd_min = ( gddncd_ref / exp(gddncd_curve*ncd_dormance(i,j)) - gddncd_offset )

          ! has the critical gdd been reached and are temperatures increasing?

          IF ( ( gdd_midwinter(i,j) .GE. gdd_min ) .AND. &
               ( t2m_week(i) .GT. t2m_month(i) ) ) THEN
             begin_leaves(i,j) = .TRUE.
             gdd_midwinter(i,j)=undef
          ENDIF

       ENDIF        ! PFT there and start of growing season allowed

    ENDDO

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving ncdgdd'

  END SUBROUTINE pheno_ncdgdd

  !
  ! ==============================================================================
  ! Phenology: leaves are put on if ngd (number of growing days, defined as
  !            days with t>-5 deg C) exceeds a critical value.
  !

  SUBROUTINE pheno_ngd (npts, j, PFTpresent, allow_initpheno, ngd, &
       t2m_month, t2m_week, begin_leaves)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! PFT index
    INTEGER(i_std), INTENT(in)                               :: j
    ! PFT exists
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: PFTpresent
    ! are we allowed to declare the beginning of the growing season?
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)               :: allow_initpheno
    ! growing degree days
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: ngd
    ! "monthly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_month
    ! "weekly" 2-meter temperatures (K)
    REAL(r_std), DIMENSION(npts), INTENT(in)                 :: t2m_week

    ! 0.2 output

    ! signal to start putting leaves on
    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)              :: begin_leaves

    ! 0.3 local

    ! index
    INTEGER(i_std)                                           :: i

    ! =========================================================================

    IF (bavard.GE.3) WRITE(numout,*) 'Entering ngd'

    !
    ! Initializations
    !

    !
    ! 1.1 initialize output
    !

    begin_leaves(:,j) = .FALSE.

    !
    ! 1.2 check the prescribed critical value
    !

    IF ( ngd_crit(j) .EQ. undef ) THEN

       WRITE(numout,*) 'ngd: ngd_crit is undefined for PFT',j
       WRITE(numout,*) 'We stop.'
       STOP

    ENDIF

    !
    ! 2 PFT has to be there and start of growing season must be allowed
    !

    DO i = 1, npts

       IF ( PFTpresent(i,j) .AND. allow_initpheno(i,j) ) THEN

          ! is critical ngd reached and are temperatures increasing?

          IF ( ( ngd(i,j) .GE. ngd_crit(j) ) .AND. &
               ( t2m_week(i) .GT. t2m_month(i) )        ) THEN
             begin_leaves(i,j) = .TRUE.
          ENDIF

       ENDIF        ! PFT there and start of growing season allowed

    ENDDO

    IF (bavard.GE.4) WRITE(numout,*) 'Leaving ngd'

  END SUBROUTINE pheno_ngd

END MODULE stomate_phenology