! recalculate vegetation cover and LAI
!
! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/lpj_cover.f90,v 1.9 2010/04/06 15:44:01 ssipsl Exp $
! IPSL (2006)
!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
MODULE lpj_cover

  ! modules used:

  USE ioipsl
  USE stomate_constants
  USE constantes_veg

  IMPLICIT NONE

  ! private & public routines

  PRIVATE
  PUBLIC cover

CONTAINS

  SUBROUTINE cover (npts, cn_ind, ind, biomass, &
       veget_max, veget_max_old, veget, lai, litter, carbon)

    !
    ! 0 declarations
    !

    ! 0.1 input

    ! Domain size
    INTEGER(i_std), INTENT(in)                                   :: npts
    ! crown area (m**2) per ind.
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)          :: cn_ind
    ! density of individuals (1/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)          :: ind
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)          :: veget_max_old

    ! 0.2 modified fields
    ! biomass (gC/(m**2 of ground))
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)   :: biomass
    ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)       :: veget_max

    ! 0.3 output

    ! fractional coverage on ground, taking into account LAI (=grid-scale fpc)
    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)       :: veget
    ! leaf area index OF AN INDIVIDUAL PLANT
    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)         :: lai

    ! metabolic and structural litter, above and below ground (gC/(m**2 of ground))
    REAL(r_std),DIMENSION(npts,nlitt,nvm,nlevs), INTENT(inout)         :: litter
    !  carbon pool: active, slow, or passive,(gC/(m**2 of ground))
    REAL(r_std),DIMENSION(npts,ncarb,nvm), INTENT(inout)               :: carbon

    ! 0.4 local

    ! index
    INTEGER(i_std)                                         :: i,j

    ! Litter dilution (gC/m²)
    REAL(r_std),DIMENSION(npts,nlitt,nlevs)                            :: dilu_lit
    ! Soil Carbondilution (gC/m²)
    REAL(r_std),DIMENSION(npts,ncarb)                                  :: dilu_soil_carbon

    ! conversion vectors
    REAL(r_std),DIMENSION(nvm)                                         :: delta_veg
    ! vecteur de conversion
    REAL(r_std)                                                        :: delta_veg_sum

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

    !
    ! 1 If the vegetation is dynamic, calculate new maximum vegetation cover for
    !   natural plants
    !

    IF ( control%ok_dgvm ) THEN

       veget_max(:,ibare_sechiba) = 1.

       DO j = 2,nvm

          IF ( natural(j) ) THEN

             veget_max(:,j) = ind(:,j) * cn_ind(:,j)

          ENDIF

          veget_max(:,ibare_sechiba) = veget_max(:,ibare_sechiba) - veget_max(:,j)

       ENDDO

       veget_max(:,ibare_sechiba) = MAX( veget_max(:,ibare_sechiba), zero )

    ENDIF

    DO i = 1, npts         
       ! Generation of the conversion vector

       delta_veg(:) = veget_max(i,:)-veget_max_old(i,:)
       delta_veg_sum = SUM(delta_veg,MASK=delta_veg.LT.zero)

       dilu_lit(i,:,:) = zero
       dilu_soil_carbon(i,:) = zero
       DO j=1, nvm
          IF ( delta_veg(j) < -min_stomate ) THEN 
             dilu_lit(i,:,:)=  dilu_lit(i,:,:) - delta_veg(j)*litter(i,:,j,:) / delta_veg_sum
             dilu_soil_carbon(i,:)=  dilu_soil_carbon(i,:) - delta_veg(j) * carbon(i,:,j) / delta_veg_sum
          ENDIF
       ENDDO

       DO j=1, nvm
          IF ( delta_veg(j) > min_stomate) THEN

             ! Dilution of reservoirs

             ! Litter
             litter(i,:,j,:)=(litter(i,:,j,:) * veget_max_old(i,j) + dilu_lit(i,:,:) * delta_veg(j)) / veget_max(i,j)

             ! Soil carbon
             carbon(i,:,j)=(carbon(i,:,j) * veget_max_old(i,j) + dilu_soil_carbon(i,:) * delta_veg(j)) / veget_max(i,j)

          ENDIF
          !SZ correct biomass to conserve mass since it's defined on veget_max
          IF(j.GE.2.AND.veget_max_old(i,j).GT.min_stomate.AND.veget_max(i,j).GT.min_stomate) THEN
             biomass(i,j,:)=biomass(i,j,:)*veget_max_old(i,j)/veget_max(i,j)
          ENDIF

       ENDDO
    ENDDO

    !
    ! 2 Calculate LAI
    !   The LAI is defined on the space covered by the crown of the plant.
    !    ( biomass / veget_max ) is in gC/(m**2 covered by the crown)
    !
    !MM in Soenke code but not in merge version ; must keep that ??
!!$    DO j = 2,nvm
!!$       lai(:,j) = biomass(:,j,ileaf,icarbon)*sla(j)
!!$    ENDDO

    !
    ! 3 calculate grid-scale fpc (foliage protected cover)
    !

    DO j = 2,nvm
       DO i = 1, npts
          IF (lai(i,j) == val_exp) THEN               
             veget(i,j) = veget_max(i,j)
          ELSE
             veget(i,j) = veget_max(i,j) * ( un - exp( - lai(i,j) * ext_coeff(j) ) )
          ENDIF
       ENDDO
    ENDDO
    !
    veget(:,ibare_sechiba) = un
    DO j = 2,nvm
       veget(:,ibare_sechiba) = veget(:,ibare_sechiba) - veget(:,j)
    ENDDO

  END SUBROUTINE cover

END MODULE lpj_cover