source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/lpj_cover.f90 @ 257

! 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_data
  USE pft_parameters

  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, turnover_daily, bm_to_litter)

    !
    ! 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
    ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground at beginning of time step
    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
    ! Turnover rates (gC/(m**2 of ground)/day)
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)          :: turnover_daily
    ! conversion of biomass to litter (g/m**2 / day
    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)          :: bm_to_litter

    ! 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(inout)         :: 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,k,m

    ! 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,reduct
    ! vecteur de conversion
    REAL(r_std)                                                        :: delta_veg_sum,diff,sr
    REAL(r_std), DIMENSION(npts)                                       :: frac_nat,sum_vegettree,sum_vegetgrass
    REAL(r_std), DIMENSION(npts)                                       :: sum_veget_natveg

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

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

    IF ( control%ok_dgvm ) THEN

       ! some initialisations
       frac_nat(:) = un
       sum_veget_natveg(:) = zero
       sum_vegettree(:) = zero
       sum_vegetgrass(:) = zero

       veget_max(:,ibare_sechiba) = un

       DO j = 2,nvm

          IF ( natural(j) ) THEN

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

          ELSE
             !fraction occupied by agriculture needs to be substracted for the DGVM
             !this is used below to constrain veget for natural vegetation, see below
             frac_nat(:) = frac_nat(:) - veget_max(:,j)

          ENDIF

       ENDDO

       DO i = 1, npts 

          IF (sum_veget_natveg(i) .GT. frac_nat(i) .AND. frac_nat(i) .GT. min_stomate) THEN

             DO j = 2,nvm
                IF( natural(j) ) THEN
                   veget_max(i,j) =  veget_max(i,j) * frac_nat(i) / sum_veget_natveg(i)
                ENDIF
             ENDDO

          ENDIF
       ENDDO

       DO j = 2,nvm
          veget_max(:,ibare_sechiba) = veget_max(:,ibare_sechiba) - veget_max(:,j)
       ENDDO
       veget_max(:,ibare_sechiba) = MAX( veget_max(:,ibare_sechiba), zero )

       ! 1.3 calculate carbon fluxes between PFTs to maintain mass balance
       !

       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

             IF(j.GE.2.AND.veget_max_old(i,j).GT.min_stomate.AND.veget_max(i,j).GT.min_stomate) THEN

                ! Correct biomass densities (i.e. also litter fall) to conserve mass 
                ! since it's defined on veget_max

                biomass(i,j,:)=biomass(i,j,:)*veget_max_old(i,j)/veget_max(i,j)
                turnover_daily(i,j,:)=turnover_daily(i,j,:)*veget_max_old(i,j)/veget_max(i,j)
                bm_to_litter(i,j,:)=bm_to_litter(i,j,:)*veget_max_old(i,j)/veget_max(i,j)

             ENDIF

          ENDDO
       ENDDO
    ENDIF

    !
    ! 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 ??
!NV, MM : we keep those comments for compatibility with CMIP5 computations.
!! They have to be uncommented avec CMIP5 versions in the trunk !
!!$    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
             IF ( control%ok_dgvm ) THEN
!!$SZneed to check this - this formulation will cause 100% veget, otherwise there will always 
!!$ be some percent bare ground
                veget(i,j) = ind(i,j) * cn_ind(i,j)  * ( un - EXP( - lai(i,j) * ext_coeff(j) ) )
             ELSE
                veget(i,j) = veget_max(i,j) * ( un - EXP( - lai(i,j) * ext_coeff(j) ) )
             ENDIF
          ENDIF

          ! check sums of fpc for natural vegetation (see correction below!) in dynamic mode
          IF ( control%ok_dgvm ) THEN

             IF(natural(j))THEN
                IF(tree(j)) THEN
                   sum_vegettree(i)=sum_vegettree(i)+veget(i,j)
                ELSE 
                   sum_vegetgrass(i)=sum_vegetgrass(i)+veget(i,j)
                ENDIF
             ENDIF

          ENDIF
       ENDDO
    ENDDO


    ! 3.1 correct gridscale fpc for dynamic vegetation
!!$SZ, this part should be obsolete now that veget_max is forced to 1.0
!!$ nevertheless maintained just for savety. Whoever wants to test
!!$ whether this works without is invited to do so.

    ! in the DGVM mode, we can arrive at a sum of veget slighly exceeding 1.0,
    ! because mainly of grass dynamics...
    ! In this case, we devide the fpar over natural vegetation first such that 
    ! grasses are shadowed by trees, and in the theoretically impossible case that
    ! this is not sufficient, reduce proportionally all veget's.
    !
    IF ( control%ok_dgvm ) THEN

       DO i = 1,npts

          diff=sum_vegettree(i)+sum_vegetgrass(i)-frac_nat(i)
          reduct(:) = 0.
          ! ordinary case, the reason too much grasses 
          ! reduce grass veget to match the maximum
          IF (diff .GT. 0. ) THEN

             IF (sum_vegetgrass(i).GT.min_stomate) THEN
                sr=0.
                DO j=2,nvm
                   IF(natural(j).AND..NOT.tree(j)) THEN
                      reduct(j)=-MIN(diff,sum_vegetgrass(i))*veget(i,j)/sum_vegetgrass(i)
                      sr=sr+reduct(j)
                   ENDIF
                ENDDO
                diff=diff+sr
             ENDIF

          ENDIF

          ! this is theoretically impossible, since trees can only occupy 95%,
          ! but better be save than sorry
          IF (diff .GT. min_stomate ) THEN

             IF (sum_vegettree(i).GT.min_stomate) THEN
                sr=0.
                DO j=2,nvm
                   IF(natural(j).AND.tree(j)) THEN
                      reduct(j)=-MIN(diff,sum_vegettree(i))*veget(i,j)/sum_vegettree(i)
                      sr=sr+reduct(j)
                   ENDIF
                ENDDO
                diff=diff+sr 
             ENDIF

          ENDIF

!!$          ! tell user if the problem could not be resolved
!!$          ! in theory the model should stop here!
!!$          IF (diff .GT. min_stomate ) THEN
!!$
!!$             write(numout,*) 'ATT, DGVM!: veget exceeds bareground without vegetation left'
!!$             write(numout,*) 'ATT, DGVM!: is this a bug? cell: ',i
!!$             write(numout,*) 'ATT, DGVM!: veget ',veget(i,:)
!!$
!!$          ENDIF

          ! finally, implement the reduction. (reduc is negative!)
          veget(i,:)=veget(i,:)+reduct(:)

       ENDDO

    ENDIF

    veget(:,ibare_sechiba) = un
    DO j = 2,nvm
       veget(:,ibare_sechiba) = veget(:,ibare_sechiba) - veget(:,j)
    ENDDO

  END SUBROUTINE cover

END MODULE lpj_cover