! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/lpj_crown.f90,v 1.12 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_crown
  !---------------------------------------------------------------------
  !- calculate individual crown area from stem mass.
  !---------------------------------------------------------------------
  USE ioipsl
  USE stomate_data
  USE constantes
  USE pft_parameters
  !-
  IMPLICIT NONE
  !-
  ! private & public routines
  !-
  PRIVATE
  PUBLIC crown
  !-
CONTAINS
  !-
  !===
  !-
  SUBROUTINE crown &
       &  (npts, PFTpresent, ind, biomass, veget_max, cn_ind, height)
    !---------------------------------------------------------------------
    ! 0 declarations
    !-
    ! 0.1 input
    !-
    ! Domain size
    INTEGER(i_std),INTENT(in) :: npts
    ! Is pft there
    LOGICAL,DIMENSION(npts,nvm),INTENT(in) :: PFTpresent
    ! density of individuals (1/(m**2 of ground))
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: ind
    ! biomass (gC/(m**2 of ground))
    REAL(r_std),DIMENSION(npts,nvm,nparts),INTENT(in) :: biomass
    !-
    ! 0.2 modified fields
    !-
    ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
    !-
    REAL(r_std),DIMENSION(npts,nvm),INTENT(inout) :: veget_max
    !-
    ! 0.3 output
    !-
    ! crown area (m**2) per ind.
    !-
    REAL(r_std),DIMENSION(npts,nvm),INTENT(out) :: cn_ind
    !-
    ! height of vegetation (m)
    !-
    REAL(r_std),DIMENSION(npts,nvm),INTENT(inout) :: height
    !-
    ! 0.4 local
    !-
    ! wood mass of an individual
    !-
    REAL(r_std),DIMENSION(npts) :: woodmass
    !-
    ! index
    !-
    INTEGER(i_std) :: j
    !-
    ! stem diameter
    !-
    REAL(r_std),DIMENSION(npts) :: dia
    REAL(r_std),DIMENSION(nvm) :: height_presc_12
    !---------------------------------------------------------------------
    !-
    ! 1 initializations
    !-
    ! 1.1 check if DGVM activated
    !-
    IF (.NOT.control%ok_dgvm) THEN
       STOP 'crown: not to be called with static vegetation.'
    ENDIF
    !-
    ! 1.2 initialize output to zero
    !-
    cn_ind(:,:) = 0.0
    ! no convertion, just cop
    height_presc_12(1:nvm) = height_presc(1:nvm)
    !-
    ! 2 calculate (or prescribe) crown area
    !-
    DO j = 2,nvm
       IF (tree(j)) THEN
          !-----
          !---- 2.1 trees
          !-----
          IF (natural(j)) THEN
             !------ 2.1.1 natural
             WHERE (PFTpresent(:,j) .AND.ind(:,j).GT.min_stomate)
                !-------- 2.1.1.1 calculate individual wood mass
                woodmass(:) = &
                     &         (biomass(:,j,isapabove)+biomass(:,j,isapbelow) &
                     &         +biomass(:,j,iheartabove)+biomass(:,j,iheartbelow))/ind(:,j)
                !-------- 2.1.1.2 stem diameter (pipe model)
                dia(:) = (woodmass(:)/(pipe_density*pi/4.*pipe_tune2)) &
                     &                **(1./(2.+pipe_tune3))
                !-------- 2.1.1.3 height
                height(:,j) = pipe_tune2*(dia(:)**pipe_tune3)
                WHERE (height(:,j) > height_presc_12(j))
                   dia(:) = (height_presc_12(j)/pipe_tune2)**(1./pipe_tune3)
                   height(:,j) = height_presc_12(j)
                ENDWHERE
                !-------- 2.1.1.4 crown area: for large truncs, crown area cannot
                !--------         exceed a certain value, prescribed through maxdia.
                cn_ind(:,j) = pipe_tune1*MIN(dia(:),maxdia(j))**pipe_tune_exp_coeff
             ENDWHERE
          ELSE
             !------ 2.1.2 tree is agricultural - stop
             STOP 'crown: cannot treat agricultural trees.'
          ENDIF
       ELSE
          !-----
          !---- 2.2 grasses
          !-----
          WHERE (PFTpresent(:,j))
             !------ 2.2.1 an "individual" is 1 m**2 of grass
             cn_ind(:,j) = 1.
          ENDWHERE
       ENDIF
       !---
       !-- 2.3 recalculate vegetation cover if natural
       !       ind and cn_ind are 0 if not present
       !---
       !SZ isn't this physically inconsistent with the assumptions of sechiba??
       ! the actual LPJ formulation calculates fpc from maximum LAI, and fpar from actual LAI=veget
       IF (natural(j).AND.control%ok_dgvm) THEN
          veget_max(:,j) = ind(:,j) * cn_ind(:,j)
       ENDIF
    ENDDO
    !-------------------
  END SUBROUTINE crown
  !-
  !===
  !-
END MODULE lpj_crown