source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE/src_stomate/lpj_gap.f90 @ 4639

! 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, mortality)

    !
    ! 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
    ! mortality (fraction of trees that is dying per time step), per day in history file
    REAL(r_std), DIMENSION(npts,nvm),INTENT(out)             :: mortality

    ! 0.3 local

    ! biomass increase
    REAL(r_std), DIMENSION(npts)                             :: delta_biomass
    ! biomass increase
    REAL(r_std), DIMENSION(npts)                             :: dmortality
    ! vigour
    REAL(r_std), DIMENSION(npts)                             :: vigour
    ! natural availability, based on vigour
    REAL(r_std), DIMENSION(npts)                             :: availability
    ! indices
    INTEGER(i_std)                                           :: j,k,m
    REAL(r_std) :: ref_greff

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

    IF ( firstcall ) THEN

       firstcall = .FALSE.

    ENDIF

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

    mortality(:,:) = zero

    ref_greff =  0.035

    DO j = 2,nvm

       ! only trees

       IF ( tree(j) ) THEN

          !
          ! 1 determine availability
          !

          IF ( .NOT.  lpj_gap_const_mort  ) THEN

             !
             ! 1.1 original formulation: mortality depends on vigour
             !

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

!SZ 080806, changed to LPJ formulation according to Smith et al., 2001 

                ! 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 )

            ! note that npp_longterm is now actually longterm growth efficiency (NPP/LAI)
            ! to be fair to deciduous trees
             delta_biomass(:) = MAX( npp_longterm(:,j) - ( turnover_longterm(:,j,ileaf) + &
                  turnover_longterm(:,j,iroot) + turnover_longterm(:,j,ifruit) + & 
                  turnover_longterm(:,j,isapabove) + turnover_longterm(:,j,isapbelow) ) ,zero)

                ! scale this to the leaf surface of the tree
!!$                vigour(:) = delta_biomass(:) / (lm_lastyearmax(:,j)*sla(j)) / 70.
             vigour(:) = delta_biomass(:) / (lm_lastyearmax(:,j)*sla(j))

             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.

!SZ 080806, changed to LPJ formulation according to Smith et al., 2001 
! tuned maximal mortality to 0.05 to get realistic range of avergage age to get ~100 years at GREFF=100
! for the range of modelled annual NPP
!!$                availability(:) = min_avail / ( 1.+vigour(:)/0.17 )
                availability(:) = 0.1 / ( 1.+ref_greff*vigour(:) )

                ! 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) = MAX(min_avail,availability(:))  * dt/one_year  
!!$                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) )

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

             ENDWHERE

          ENDDO

          !
          ! 4 update number of individuals
          !

!SZ 080806, allow changing density in static case when mortality is dynamic
          IF ( control%ok_dgvm .OR. .NOT.lpj_gap_const_mort) THEN

             WHERE ( PFTpresent(:,j) )

                ind(:,j) = ind(:,j) * ( un - mortality(:,j) )

             ENDWHERE

          ENDIF
       ELSE 

          IF ( .NOT.control%ok_dgvm .AND. .NOT.lpj_gap_const_mort) THEN

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

                mortality(:,j) = 1.

             ENDWHERE
             DO k = 1, nparts

                WHERE ( PFTpresent(:,j) )

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

                ENDWHERE
             ENDDO
             
          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