! $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.
  !- SZ, I've put the woodmass calculation out of this routine
  !      because after the very first establishment, woodmass
  !      could not be calculated here as veget_max = zero and 
  !      d_ind not known...
  !---------------------------------------------------------------------
  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, woodmass_ind, 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
    ! woodmass of the individual, needed to calculate crownarea in lpj_crownarea
    REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: woodmass_ind
    !-
    ! 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 .AND. lpj_gap_const_mort) THEN
       STOP 'crown: not to be called with static vegetation.'
    ENDIF
    !-
    ! 1.2 initialize output to zero
    !-
    cn_ind(:,:) = zero
    ! 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)
             WHERE (PFTpresent(:,j) .AND.woodmass_ind(:,j).GT.min_stomate)
!!$SZ note that woodmass_ind needs to be defined on the individual, hence
!!$ biomass*veget_max/ind, not as stated here, correction MERGE
!!$!-------- 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)) &
                dia(:) = (woodmass_ind(:,j)/(pipe_density*pi/4.*pipe_tune2)) &
                     &                **(1./(2.+pipe_tune3))
                !-------- 2.1.1.3 height
                height(:,j) = pipe_tune2*(dia(:)**pipe_tune3)
!!$SZ: The constraint on height has nothing to do with LPJ (for that purpose there's dia_max
!!$ cannot see why this is necessary - it also blurrs the output, hence I leave it commented
!!$                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) = un
          ENDWHERE
       ENDIF
       !---
       !-- 2.3 recalculate vegetation cover if natural
       !       ind and cn_ind are 0 if not present
       !---
!!$SZ: since now all state variables are defined on veget_max it is very
!!$ dangerous to change this several times in stomate_lpj, as then GPP, turnover and allocated 
!!$ biomass are not defined on the same space! Hence, veget_max is now kept constant
!!$ and updated at the end of stomate_lpj in lpj_cover.f90
!!$ Eventually, this routine should only be called once at the beginning and the end of stomate_lpj
!!$ or prefereably cn_ind made a saved state variable!
!!$    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