! gap routine - place for new plants
!
! Death rate of trees is estimated by evaluating their vigour (based on npp).
! For large availabilities, lifetime is 50 years (!?).
! Age of stands is not considered, although availability death rate should probably
! depend on age.
!
! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/lpj_gap.f90,v 1.10 2009/01/06 15:01:25 ssipsl Exp $
! IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
MODULE lpj_gap

  ! modules used:

  USE ioipsl
  USE stomate_constants
  USE constantes_veg
  USE parallel

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC gap,gap_clear

  ! first call
  LOGICAL, SAVE                                           :: firstcall = .TRUE.

CONTAINS


  SUBROUTINE gap_clear
    firstcall = .TRUE.
  END SUBROUTINE gap_clear

  SUBROUTINE gap (npts, dt, &
       npp_longterm, turnover_longterm, lm_lastyearmax, &
       PFTpresent, biomass, ind, bm_to_litter)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                     :: npts
    ! Time step (days)
    REAL(r_std), INTENT(in)                                  :: dt
    ! "long term" net primary productivity (gC/(m**2 of ground)/year)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: npp_longterm
    ! "long term" turnover rate (gC/(m**2 of ground)/year)
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)     :: turnover_longterm
    ! last year's maximum leaf mass, for each PFT (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)            :: lm_lastyearmax

    ! 0.2 modified fields

    ! Is pft there
    LOGICAL, DIMENSION(npts,nvm), INTENT(in)            :: PFTpresent
    ! biomass (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)  :: biomass
    ! Number of individuals / (m**2 of ground)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)         :: ind
    ! biomass taken away (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)  :: bm_to_litter

    ! 0.3 local

    ! which kind of mortality
    LOGICAL, SAVE                                           :: constant_mortality
    ! biomass increase
    REAL(r_std), DIMENSION(npts)                             :: delta_biomass
    ! vigour
    REAL(r_std), DIMENSION(npts)                             :: vigour
    ! natural availability, based on vigour
    REAL(r_std), DIMENSION(npts)                             :: availability
    ! mortality (fraction of trees that is dying per time step), per day in history file
    REAL(r_std), DIMENSION(npts,nvm)                        :: mortality
    ! indices
    INTEGER(i_std)                                           :: j,k

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

    IF ( firstcall ) THEN

       firstcall = .FALSE.

       !Config  Key  = LPJ_GAP_CONST_MORT
       !Config  Desc = constant tree mortality
       !Config  Def  = y
       !Config  Help = If yes, then a constant mortality is applied to trees. 
       !Config         Otherwise, mortality is a function of the trees' 
       !Config         vigour (as in LPJ).

       constant_mortality = .TRUE.
       CALL getin_p('LPJ_GAP_CONST_MORT', constant_mortality)     
       WRITE(numout,*) 'gap: constant mortality:', constant_mortality

    ENDIF

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

    mortality(:,:) = zero

    DO j = 2,nvm

       ! only trees

       IF ( tree(j) ) THEN

          !
          ! 1 determine availability
          !

          IF ( .NOT. constant_mortality ) THEN

             !
             ! 1.1 original formulation: mortality depends on vigour
             !

             WHERE ( PFTpresent(:,j) .AND. ( lm_lastyearmax(:,j) .GT. min_stomate ) )

                ! how much did the tree grow per year?

                delta_biomass(:) = &
                     MAX( npp_longterm(:,j) - ( turnover_longterm(:,j,ileaf) + &
                     turnover_longterm(:,j,iroot) + turnover_longterm(:,j,ifruit) ), &
                     0._r_std )

                ! scale this to the leaf surface of the tree

                vigour(:) = delta_biomass(:) / (lm_lastyearmax(:,j)*sla(j)) / 70.

             ELSEWHERE

                vigour(:) = 0.0

             ENDWHERE

             WHERE ( PFTpresent(:,j) )

                ! note that availability is never above 0.02, i.e. lifetime of 50 years when very
                ! low vigour.

                availability(:) = 0.02 / ( 1.+vigour(:)/0.17 )

                ! Mortality (fraction per time step).
                ! In the original DGVM, mortality was set to zero if there was strong fire
                ! perturbation.
                ! This has been de-activated since the npp is not influenced by fire,
                ! as opposed to the original DGVM. Instead, mortality is simply
                ! equal to the availability, modulated by the time step.
                ! Exact formulation: mor = 1. - ( 1. - availability ) ** (dt/one_year)
                ! approximation ok as availability < 0.02 << 1

                mortality(:,j) = availability(:) * dt/one_year

             ENDWHERE

          ELSE

             !
             ! 1.2 Alternative version: Constant mortality
             !

             WHERE ( PFTpresent(:,j) )

                mortality(:,j) = dt/(residence_time(j)*one_year)

             ENDWHERE

          ENDIF

          !
          ! 2 Special for the DGVM:
          !   mortality is one if npp is zero or negative.
          !

          IF ( control%ok_dgvm ) THEN

             WHERE ( PFTpresent(:,j) .AND. ( npp_longterm(:,j) .LE. min_stomate ) )

                mortality(:,j) = 1.

             ENDWHERE

          ENDIF

          !
          ! 3 update biomass, create litter
          !

          DO k = 1, nparts

             WHERE ( PFTpresent(:,j) )

                bm_to_litter(:,j,k) = bm_to_litter(:,j,k) + mortality(:,j) * biomass(:,j,k)

                biomass(:,j,k) = biomass(:,j,k) * ( 1. - mortality(:,j) )

             ENDWHERE

          ENDDO

          !
          ! 4 update number of individuals
          !

          IF ( control%ok_dgvm ) THEN

             WHERE ( PFTpresent(:,j) )

                ind(:,j) = ind(:,j) * ( 1. - mortality(:,j) )

             ENDWHERE

          ENDIF

       ENDIF       ! only trees

    ENDDO         ! loop over pfts

    !
    ! 5 history
    !

    ! output in fraction of trees that dies/day.
    ! exact formulation: 1. - ( 1. - mortality ) ** (1./dt)
    mortality = mortality / dt

    CALL histwrite (hist_id_stomate, 'MORTALITY', itime, &
         mortality, npts*nvm, horipft_index)

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

  END SUBROUTINE gap

END MODULE lpj_gap