!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*    - Calculation of annual allocation of NPP (SR PARTITION)   *!
!*    - Calculation of annual allocation of NPP of soil          *!
!*         vegetation (PARTITION_SV                              *!
!*    - Calculation of diameter at breast height (SR CALC_DBH)   *!
!*                                                               *!
!*                  Copyright (C) 1996-2018                      *!
!*     Potsdam Institute for Climate Impact Reserach (PIK)       *!
!*          Authors and contributors see AUTHOR file             *!
!*  This file is part of 4C and is licensed under BSD-2-Clause   *!
!*                   See LICENSE file or under:                  *!
!*     http://www.https://opensource.org/licenses/BSD-2-Clause   *!
!*                           Contact:                            *!
!*       https://gitlab.pik-potsdam.de/foresee/4C                *!
!*                                                               *!
!*****************************************************************!

!****************************!
!*   SUBROUTINE PARTITION   *!
!****************************!

SUBROUTINE PARTITION( p )

  !*** Declaration part ***!
  USE data_out
  USE data_par
  USE data_stand
  USE data_species
  USE data_simul

  IMPLICIT NONE

  REAL   :: lambdaf = 0.,  &      ! partitioning functions
            lambdas = 0.,  &
            lambdar = 0.,  &
            lambdac = 0.,  &
            lambdaSum = 0.,&      ! sum of the above three lambdas
            NPP = 0.,      &      ! annual NPP
            F = 0.,        &      ! state variables: foliage,
            S = 0.,        &      ! sapwood,
            H = 0.,        &      ! heartwood
            R = 0.,        &      ! fine roots,
            B = 0.,        &      ! bole height,
            Ahb = 0.,      &      ! cross sectional area heartwood at tree base
            hs = 0.,       &      ! sapwood height
            Ht = 0.,       &      ! total tree height
            Asw = 0.,      &      ! cross sectional area of sapwood in bole
            DBH = 0.,      &      ! tree diameter at breast height (DBH)
            FNew, SNew,    &      ! new states
            RNew, BNew,    &
            HtNew,  &
            HNew, Ahbnew,  &
            sigmaf = 0.,   &      ! current leaf activity rate
            sigman = 0.,   &      ! current root activity rate
            ar = 0.,       &      ! aux vars for partitioning functions
            as = 0.,       &
            ac = 0.,       &
            betar = 0.,    &
            betas = 0.,    &
            aux = 0.,      &
            Fmax,          &      ! determines whether height growth or not
            rsap,          &      ! auxiliary variable for height growth determination
            growthrate            ! height growthrate depends on relative light regime in the middle of the canopy
  REAL  ::  Sf,            &      ! senescence rates
            Ss,            &
            Sr,            &
            Gf,            &      ! growth rates
            Gs,            &
            Gr
 real ::    DBH_help
  REAL  ::  leaf_N_conc,      &   ! last years N concentration in leaves gN kgDM
            tbc_root_Ndemand, &   ! N demand for ghrowth of fine roots, branches and coarse roots g tree-1
            Nredfak               ! reduction factor for N allocation to fine roots, branches and coarse roots

  TYPE(Coh_Obj) :: p        ! pointer to cohort list

  REAL   :: term1,       &
            a1, a2, a3,  &  ! coefficients of quadratic equation
            x1 = 0.,     &
            x2 = 0.         ! solutions of quadratic equation
 
 real   ::  Fmax_old

  ! if this cohort is mistletoe infected, reduce NPP by mistletoe-specific demand
  ! demand is defined in PARTITION_MI. as mistletoe is always 1st cohort, the demand of mistletoe is calculated before the reduction here
  if (p%coh%mistletoe.eq.1) then
     p%coh%NPP = p%coh%NPP-(NPP_demand_mistletoe/p%coh%ntreea)
  endif
  ns   = p%coh%species
  F    = p%coh%x_fol
  Fmax = p%coh%Fmax
  S    = p%coh%x_sap
  R    = p%coh%x_frt
  H    = p%coh%x_hrt
  B    = p%coh%x_hbole
  NPP  = p%coh%NPP
  Ht   = p%coh%height
  Ahb  = p%coh%x_Ahb
  Sf   = p%coh%sfol
  Ss   = p%coh%ssap
  Sr   = p%coh%sfrt
  hs   = p%coh%x_hsap
  Asw  = p%coh%Asapw
  Fmax_old = Fmax
  
  DBH_help = p%coh%diam

  if (flag_end.eq.1) then
         p%coh%notViable = .TRUE.
         flag_end = 0
  end if

if(p%coh%notViable.neqv..TRUE.) then
  select case (flag_folhei)
  case (1,4)
    spar(ns)%pha = spar(ns)%pha_v1 * spar(ns)%pha_v3 *     &
        (F)**(-1-spar(ns)%pha_v3)/(spar(ns)%pha_v2+(F)**(-spar(ns)%pha_v3))**2.

  case (2) 
  
     rsap=Asw/(Asw+Ahb)
     spar(ns)%pha = 2.*spar(ns)%crown_a/(pi**0.5*(rsap*spar(ns)%pnus)**1.5*F**0.5)
  
  case (3)
    ! this version only for tests and pine trees
      spar(ns)%pha = (3500*(10.+F**0.9)-(0.9*3500.*F**0.9))/(10.+F**0.9)**2
   
  end select   ! flag_folhei 

  ! only allocate if enough NPP is available
  IF (NPP>1.0E-9) THEN
  
      select case (flag_folhei)
      case (0)
            growthrate=spar(ns)%pha*spar(ns)%pha_coeff1 + spar(ns)%pha*spar(ns)%pha_coeff2*(1./p%coh%IrelCan-1.)

      case (1,3)
            growthrate=spar(ns)%pha + spar(ns)%pha*(1./MAX(p%coh%IrelCan,0.25)-1.)
      
      case (2)
            growthrate=spar(ns)%pha + spar(ns)%pha*(1.-p%coh%IrelCan)*5.

      case (4)
            growthrate=spar(ns)%pha *0.5/MAX(p%coh%IrelCan,0.25)
      
      end select   ! flag_folhei

      sigmaf = NPP/F

         ! calculate root activity based on drought index
         ! test of a relationship which modifies fine root leaf ratio with shade tolerance:

      IF (flag_sign.eq.1 .or. flag_sign.eq.11) THEN
         term1 = spar(ns)%sigman * 10. * (((5.-spar(ns)%stol)*1.-p%coh%crown_area) / (5.-spar(ns)%stol)*1.)
         sigman = amax1(term1,spar(ns)%sigman) * p%coh%drIndAl/p%coh%nDaysGr
      ELSE
         sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
      END IF
      
      if (flag_sign .eq. 0 .or. flag_sign .eq. 1) then  
         ! auxiliary variables for fine roots
         ar = spar(ns)%pcnr * sigmaf / sigman
         betar = (Sr - R + ar*(F-Sf)) / NPP

         ! auxiliary variables for sapwood
         as = spar(ns)%prhos / spar(ns)%pnus
         aux = 2.*(B+p%coh%deltaB) + Ht
         betas = ( (as/3.)*(aux - growthrate*Sf) * (F-Sf) + Ss - S ) / NPP

         ! solve quadratic equation for lambdaf
         term1 = (1.+spar(ns)%alphac)
         a1 = term1 * as/3. * growthrate * NPP
         a2 = 1.0 + ar + term1 * as/3. * (aux + growthrate*(F-2.*Sf))
         a3 = term1*betas + betar - 1.

         x1 = (-a2 + SQRT( a2*a2 - 4.*a1*a3) ) / (2.*a1)
         x2 = (-a2 - SQRT( a2*a2 - 4.*a1*a3) ) / (2.*a1)
         lambdaf = x1
    
         if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
            lambdaf = 0.5
            lambdar = 0.5
            lambdas = 0.
            lambdac = 0.
         else

           ! calculate coefficients for sapwood and roots
             lambdar = ar * lambdaf + betar;
             lambdas = as/3. * (aux + growthrate*(F+lambdaf*NPP-2.*Sf)) * lambdaf + betas
             lambdac = spar(ns)%alphac * lambdas    

           ! check consistency of calculation, i.e. no negative values
             IF(lambdas <  0. .or. lambdas .gt. 1.) THEN
                lambdas = 0.
                lambdac = 0.
                lambdaf = (1.-betar)/(ar+1)
                lambdar = 1.-lambdaf
                
                if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
                    lambdaf = 0.5
                    lambdar = 0.5

                else if (lambdar<0) then
                   lambdar=0.
                   lambdaf=1.
                end if

             ELSE
                ! reduced allocation schemes for lamdaf<0. or lamdar<0. still to be added
                lambdaf = AMAX1( lambdaf, 0. )
                lambdar = AMAX1( lambdar, 0. )

                ! warrant that lambdaSum = 1 if balance can not be achieved this time step
                lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar
                lambdaf   = lambdaf / lambdaSum
                lambdas   = lambdas / lambdaSum
                lambdar   = lambdar / lambdaSum
                lambdac   = lambdac / lambdaSum
                lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar  ! for debugging only 

             END IF
         end if  ! lambdaf .le. 0.
  
     else   ! flag_sign = 10, 11
         ! auxiliary variables for fine roots
          ar    = spar(ns)%pcnr * sigmaf / sigman
          betar = (Sr - ar*Sf) / NPP

         ! auxiliary variables for sapwood
          as    = spar(ns)%prhos / spar(ns)%pnus
          betas = (Ss - 2.*as*hs*Sf ) / NPP

         ! auxiliary variables for coarse roots, twigs and branches
          ac    = spar(ns)%alphac

         ! linear equation system in lamda(i)
          term1   = 1. + ar + 2.*as*hs*(1+ac)
          lambdaf = 1. - (1.+ac)*betas - betar  
          lambdaf = lambdaf / term1  
          lambdar = ar * lambdaf + betar
          lambdas = 2.*as*hs * lambdaf + betas  
          lambdac = ac * lambdas    

         if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
            lambdaf = 0.5
            lambdar = 0.5
            lambdas = 0.
            lambdac = ac * lambdas    
         else

           ! calculate coefficients for sapwood and roots
             lambdar = ar * lambdaf + betar;
             lambdas = 2.*as*hs * lambdaf + betas  
             lambdac = ac * lambdas    

           ! check consistency of calculation, i.e. no negative values
             IF(lambdas <  0. .or. lambdas .gt. 1.) THEN
                lambdas = 0.
                lambdac = 0.
                lambdaf = (1.-betar)/(ar+1)
                lambdar = 1.-lambdaf
                
                if (lambdaf .le. 0. .or. lambdaf .gt. 1.) then
                    lambdaf = 0.5
                    lambdar = 0.5

                else if (lambdar<0) then
                   lambdar=0.
                   lambdaf=1.
                end if

             ELSE
                ! reduced allocation schemes for lamdaf<0. or lamdar<0. still to be added
                lambdaf = AMAX1( lambdaf, 0. )
                lambdar = AMAX1( lambdar, 0. )

                ! warrant that lambdaSum = 1 if balance can not be achieved this time step
                lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar
                lambdaf   = lambdaf / lambdaSum
                lambdas   = lambdas / lambdaSum
                lambdar   = lambdar / lambdaSum
                lambdac   = lambdac / lambdaSum
                lambdaSum = lambdaf + (1.+spar(ns)%alphac)*lambdas + lambdar  ! for debugging only 

             END IF
         end if  ! lambdaf .le. 0.
          
      endif  ! flag_sign  
  ELSE  

     lambdaf   = 0.
     lambdas   = 0.
     lambdar   = 0.

  END IF   ! IF NPP < 1.0E-09

      ! gross growth rates of compartments
      Gf = lambdaf * NPP
      Gr = lambdar * NPP
      Gs = lambdas * NPP
      p%coh%gfol = Gf
      p%coh%gfrt = Gr
      p%coh%gsap = Gs
      p%coh%x_crt  = p%coh%x_crt + Gs*spar(ns)%alphac*spar(ns)%cr_frac
      p%coh%x_tb   = p%coh%x_tb  + Gs*spar(ns)%alphac*(1.-spar(ns)%cr_frac)

      ! update of state vector
      FNew  = F + Gf - Sf
      SNew  = S + Gs - Ss
      RNew  = R + Gr - Sr
      Hnew  = H + Ss
      AhbNew= Ahb + Asw*spar(ns)%pss
     
      ! check whether height growth or not

      IF (lambdas == 0.OR.FNew<Fmax) THEN  ! treat special case where there is no height growth
        HtNew = Ht
      ELSE
      ! height growth depending on the relative light regime in the middle of the canopy
        HtNew = Ht + growthrate * (FNew-Fmax)
        Fmax=FNew
      ENDIF
        BNew = B+p%coh%deltaB

      ! copy back to original variables
      p%coh%Fmax = Fmax
      p%coh%x_fol  = FNew
      p%coh%x_sap  = SNew
      p%coh%x_frt  = RNew
      p%coh%x_hrt  = HNew
      p%coh%height = HtNew
      p%coh%x_hbole= BNew
      p%coh%x_Ahb  = AhbNew
  
  CALL CALC_DBH(BNew,Htnew,Snew,Hnew,Ahbnew,p%coh%Ahc,p%coh%ident,DBH,p%coh%dcrb,hs,Asw)
  if (flag_end.eq.1) then
      DBH = p%coh%diam
      p%coh%notViable = .TRUE.
      flag_end = 0
  end if

  ! Monitoring of current values
  if (time_out .gt. 0 .and. flag_cohout .eq. 2) then
     CALL OUT_ALL( p%coh%ident, p%coh%ntreea, NPP, DBH, growthrate,Fnew,Fmax_old,Htnew, lambdaf,lambdas,lambdar,lambdac,x1,x2)
  endif

  p%coh%x_hsap = hs
  p%coh%diam   = DBH  !  This is the new value
  p%coh%Asapw  = Asw

  p%coh%jrb = (DBH-DBH_help)*10/2
  
  if(((DBH-DBH_help)*10/2).lt.0.) p%coh%jrb = 0.
  
  ! variables required by mortality submodel
  p%coh%fol_inc  = Gf - Sf
  p%coh%bio_inc  = NPP - Sf - (1.+spar(ns)%alphac)*Ss - Sr
  p%coh%stem_inc = Gs         ! deltaH + deltaS = Ss + Gs - Ss
  p%coh%frt_inc = Gr - Sr     ! fine root increment
  p%coh%totBio = p%coh%x_fol + (1.+spar(ns)%alphac)*(p%coh%x_sap + p%coh%x_hrt) + p%coh%x_frt
  p%coh%notViable = (FNew <= 0.) .OR. (SNew <= 0.) .OR.    &
                    (RNew <= 0.) .OR. (Htnew <= Bnew)

! Nitrogen dynamics:
  leaf_N_conc = p%coh%N_fol/F

! Simple model: all (sap)wood grows with CN-ratios of branches, twigs and coarse roots.  
! When sapwood senesces N is reallocated and the new heart wood is at the level of stem CN-ratios.
! Branches, twigs and coarse roots do not senesce

! first step nitrogen related processes: N in litter, N-recallocation
  p%coh%N_pool = p%coh%N_pool + Sf/F*p%coh%N_fol*spar(ns)%reallo_fol &
                              + Sr*cpart/spar(ns)%cnr_frt*1000.* spar(ns)%reallo_frt &
                              + Ss*cpart *1000. * (1/spar(ns)%cnr_tbc - 1/spar(ns)%cnr_stem)
  p%coh%N_fol = p%coh%N_fol*(1-Sf/F)

  p%coh%litC_fol  = p%coh%litC_fol +  p%coh%ntreea * Sf * cpart
  p%coh%litC_frt  = p%coh%litC_frt +  p%coh%ntreea * Sr * cpart

  ! Species specific N content and reallocation factor (see species.par)
  ! Caution: tbc mortallity is not a litter compartment; it is assigned as heartwood
   p%coh%litN_fol  = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * (1.-spar(ns)%reallo_fol) / spar(ns)%cnr_fol
   p%coh%litN_frt  = p%coh%litN_frt +  p%coh%ntreea * Sr * cpart * (1.-spar(ns)%reallo_frt) / spar(ns)%cnr_frt

  ! second step: allocation of N to new growth
  ! before bud-break allocation to leaves is 50% of the N content of last years foliage
  tbc_root_Ndemand = Gs*cpart *kg_in_g / spar(ns)%cnr_tbc + Gr* cpart/spar(ns)%cnr_frt*kg_in_g
  IF(tbc_root_Ndemand + Gf*p%coh%med_sla*0.5 > p%coh%N_pool) THEN
     if (tbc_root_Ndemand .gt. 1E-8) then
        Nredfak = AMAX1((p%coh%N_pool-Gf*p%coh%med_sla*0.5) / tbc_root_Ndemand,0.)   ! Division by zero possible 
     else
        Nredfak = 0.
     endif
     tbc_root_Ndemand = tbc_root_Ndemand*Nredfak
  ENDIF
     p%coh%N_pool = p%coh%N_pool - tbc_root_Ndemand
     IF(p%coh%N_pool < Gf*0.5*leaf_N_conc) THEN
        p%coh%N_fol = p%coh%N_fol + p%coh%N_pool
        p%coh%N_pool = 0.
     ELSE
        p%coh%N_fol = p%coh%N_fol + Gf*0.5*leaf_N_conc
        p%coh%N_pool = p%coh%N_pool - Gf*0.5*leaf_N_conc
     ENDIF
end if

END SUBROUTINE PARTITION

!*******************************!
!*   SUBROUTINE PARTITION_SV   *!
!*******************************!

SUBROUTINE PARTITION_SV( p )

  !*** Declaration part ***!
  USE data_par
  USE data_stand
  USE data_species
  USE data_simul

  IMPLICIT NONE

  REAL   :: lambdaf = 0.,  &      ! partitioning functions
            lambdas = 0.,  &
            lambdar = 0.,  &
            NPP = 0.,      &      ! annual NPP
            F = 0.,        &      ! state variables: foliage,
            S = 0.,        &      ! sapwood,
            R = 0.,        &      ! fine roots,
            Ht = 0.,       &      ! total tree height
            FNew, SNew,    &      ! new states
            RNew,          &
            sigmaf = 0.,   &      ! current leaf activity rate
            sigman = 0.           ! current root activity rate
  REAL  ::  Sf,            &      ! senescence rates
            Ss,            &
            Sr,            &
            Gf,            &      ! growth rates
            Gs,            &
            Gr

  REAL  ::  FRsum

  REAL  ::  tbc_root_Ndemand, &   ! N demand for ghrowth of fine roots, branches and coarse roots g tree-1
            Nredfak               ! reduction factor for N allocation to fine roots, branches and coarse roots

  REAL, EXTERNAL  ::  f_lf, df_lf, ddf_lf

  INTEGER :: flag_SV_allo,  &
             rnum

  TYPE(Coh_Obj) :: p       ! pointer to cohort list

  ns   = p%coh%species
  F    = p%coh%x_fol
  S    = p%coh%x_sap
  R    = p%coh%x_frt
  NPP  = p%coh%NPP
  Ht   = p%coh%height
  Sf   = p%coh%sfol
  Ss   = p%coh%ssap
  Sr   = p%coh%sfrt

  ! choice of allocation model. 0 = constant allocation factors, 1 = allometric model
  flag_SV_allo = 1

  ! only allocate if enough NPP is available
  IF (NPP>1.0E-9) THEN

     ! calculate leaf activity based on net PS and leaf mass
     sigmaf = NPP/F

     ! calculate root activity based on drought index
! test of a relationship which modifies fine root leaf ratio with shade tolerance
      IF (flag_sign.eq.1) THEN
         sigman = amax1(spar(ns)%sigman*10*(((5.-spar(ns)%stol)*1.-p%coh%crown_area)/(5.-spar(ns)%stol)*1.),spar(ns)%sigman) * p%coh%drIndAl / p%coh%nDaysGr
      ELSE
         sigman = spar(ns)%sigman * p%coh%drIndAl / p%coh%nDaysGr
      END IF
  M_avail=(NPP+F-Sf+R-Sr+S-Ss)/kpatchsize

  IF(flag_SV_allo==0) THEN
  ! the parameters pdiam in the species.par file are used for allocation fractions
     lambdaf=spar(ns)%pdiam1
     lambdar=spar(ns)%pdiam2
     lambdas=spar(ns)%pdiam3 
  ELSE  
     FRsum=(F+R)/kpatchsize  
     CALL newt (FRsum, f_lf, df_lf, ddf_lf, 1.e-6, 100, rnum)
     IF(FRsum>M_avail .and. .not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort, rnum: ',real(rnum))
     IF(rnum==-1) THEN
        if (.not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort: ',real(p%coh%ident))
        lambdaf=0.4
        lambdar=0.4
        lambdas=0.2 
     ELSE
        lambdaf=(FRsum)/M_avail/2.
        lambdar=(FRsum)/M_avail/2.
        lambdas=1.-lambdaf-lambdar
     ENDIF
  ENDIF 

  END IF   ! IF NPP < 1.0E-09

  ! gross growth rates of compartments

  Gf = lambdaf * M_avail*kpatchsize -F +Sf
  Gr = lambdar * M_avail*kpatchsize -R +Sr
  Gs = lambdas * M_avail*kpatchsize -S +Ss

! preliminary solution for permanent seeding
  IF(lambdaf * M_avail < 1.e-4) THEN
     Gf = Gf + 1.e-4*kpatchsize
  ENDIF

  p%coh%gfol = Gf
  p%coh%gfrt = Gr
  p%coh%gsap = Gs

  ! update of state vector
  FNew  = F + Gf - Sf
  SNew  = S + Gs - Ss
  RNew  = R + Gr - Sr
  p%coh%x_fol  = FNew
  p%coh%x_sap  = SNew
  p%coh%x_frt  = RNew

  ! determine litter production from plant turnover rates
  ! first step nitrogen related processes: N in litter, N-recallocation
  p%coh%N_pool = p%coh%N_pool + Sf/F*p%coh%N_fol*spar(ns)%reallo_fol &
                              + Sr*cpart/spar(ns)%cnr_frt*1000.* spar(ns)%reallo_frt &
                              + Ss*cpart *1000. * (1/spar(ns)%cnr_tbc - 1/spar(ns)%cnr_stem)
  p%coh%N_fol = p%coh%N_fol*(1-Sf/F)

  ! Summation, due to the filling of the pool at other points as well
  p%coh%litC_fol  = p%coh%litC_fol +  p%coh%ntreea * Sf * cpart
  p%coh%litC_frt  = p%coh%litC_frt +  p%coh%ntreea * Sr * cpart

  ! New version with species specific N content and reallocation factor (see species.par) 
  ! changed to 1-reallo 
   p%coh%litN_fol  = p%coh%litN_fol + p%coh%ntreea * Sf * cpart * (1.-spar(ns)%reallo_fol) / spar(ns)%cnr_fol
   p%coh%litN_frt  = p%coh%litN_frt +  p%coh%ntreea * Sr * cpart * (1.-spar(ns)%reallo_frt) / spar(ns)%cnr_frt

  ! second step: allocation of N to new growth
  ! before bud-break allocation to leaves is 50% of the N content of last years foliage
  tbc_root_Ndemand = Gs*cpart *kg_in_g / spar(ns)%cnr_tbc + Gr* cpart/spar(ns)%cnr_frt*kg_in_g
  IF(tbc_root_Ndemand + Gf*p%coh%med_sla*0.5 > p%coh%N_pool) THEN
     if (tbc_root_Ndemand .gt. 1E-8) then
        Nredfak = AMAX1((p%coh%N_pool-Gf*p%coh%med_sla*0.5) / tbc_root_Ndemand,0.)   ! Div. by zero possible !
     else
        Nredfak = 0.
     endif
     tbc_root_Ndemand = tbc_root_Ndemand*Nredfak
  ENDIF
     p%coh%N_pool = p%coh%N_pool - tbc_root_Ndemand

    END SUBROUTINE PARTITION_SV

!*******************************!
!*   SUBROUTINE PARTITION_MI   *!
!*******************************!

SUBROUTINE PARTITION_MI( p )
  !*** Declaration part ***!
  USE data_par
  USE data_stand
  USE data_simul
  IMPLICIT NONE
  TYPE(Coh_Obj) :: p       ! pointer to cohort list
  !no partitioning, foliage mass keeps constant
  p%coh%x_fol  = p%coh%x_fol !  !FNew
  p%coh%x_sap  = 0.!SNew
  p%coh%x_frt  = 0.!RNew
END SUBROUTINE PARTITION_MI

!***************************!
! FUNCTION f_lf            *!
!***************************!

REAL FUNCTION f_lf(x)
  USE data_stand
  USE data_plant
  REAL :: x
  f_lf = ksi*x**kappa + x - M_avail 
END ! FUNCTION f_lf

!***************************!
! FUNCTION df_lf           *!
!***************************!

REAL FUNCTION df_lf(x)
  USE data_stand
  USE data_plant
  REAL :: x
  df_lf = ksi*kappa*x**(kappa-1.) + 1. 
END ! FUNCTION df_lf

!***************************!
! FUNCTION ddf_lf          *!
!***************************!

REAL FUNCTION ddf_lf(x)
  USE data_stand
  USE data_plant
  REAL :: x
  ddf_lf = ksi*kappa*(kappa-1.)*x**(kappa-2.) 
END ! FUNCTION ddf_lf

!***************************!
! SUBROUTINE CALC_DBH      *!
!***************************!
SUBROUTINE CALC_DBH(B, Ht, S, H, Ahb, Ahc, ident, dbh, dc, hs, Asw)

  !*** Declaration part ***!

USE data_par
USE data_species
USE data_simul

IMPLICIT NONE

INTEGER :: ident

REAL :: Dc         ! diameter at crown base
REAL :: B,   &     ! bole height,
        Ht,  &     ! total tree height
        S,   &     ! sapwood
        H,   &     ! heartwood
        hs,  &     ! sapwood height
        D,   &     ! stem diameter at forest floor
        DBH, &     ! tree diameter at breast height
        Ahb, &     ! cross sectional area heartwood at tree base
        Ahc, &     ! cross sectional area of heartwood at crown base
        Asw, &     ! cross sectional area of sapwood in bole
        discr, func, help, hp1, hp2,hp3, hp4
REAL :: fp, fq,  & ! coefficients of quadratic equation
        w1, w2,  & ! solutions of quadratic equation
        precision  ! criterion for acceptance of solution
real  :: sprhos    ! sapwood density [kg/cm3] 

  !*** Calculation part ***!

  precision = 1.e-5
  sprhos = spar(ns)%prhos
! calculate Diameters
      hs = (2*B +Ht)/3.
     Asw = S/(spar(ns)%prhos*hs)

! if Bole height >= height trees are dead and calculations not required
  IF(B .lt. Ht) THEN
     select case (flag_volfunc) 
     case (0)         
        D   = SQRT( (S+H)*4. / (sprhos*hs*pi) )
        IF (Ht<h_breast) THEN
           DBH = 0.0
        ELSEIF (Ht>h_breast.and.B<h_breast) then
           DBH=D-(D/(Ht-B))*(h_breast-B)
        ELSE
           DBH=D
        ENDIF

     case (1)
       D = SQRT((Ahb+Asw)*4./pi)
     ! if Bole height = 0 then there is no need to calulate Diameter at crown base and Dc = D
       IF(B.EQ.0.) THEN
          Dc = D
       ELSE
          fp  = -2. * (B/Ht) * (3.*H/(sprhos*B)-Ahb)-Ahb*(B/Ht)**2.
          fp  = -2. * B/Ht * (3.*H/(sprhos*B)-Ahb)-Ahb*(B/Ht)**2.
          fq    = ((3.*H/(sprhos) - Ahb*B) / Ht)**2.
          discr = fp**2./4.-fq
          if (abs(discr) .lt. zero) then
            discr = zero      ! avoid small values
          endif
         ! No solution
          if(discr.lt.0) then
            if (.not.flag_mult8910) then
               CALL error_mess(time,'discriminant < 0 in calc_dbh for cohort: ',real(ident))
               CALL stop_mess(time,'discriminant < 0 in calc_dbh ')
               CALL error_mess(time,'stop in calc_dbh for stand No: ',real(ip))
               CALL error_mess(time,'heart wood mass H: ',H)
               CALL error_mess(time,'bole height b: ',b)
               CALL error_mess(time,'height Ht: ',Ht)
               CALL error_mess(time,'ave. sapwood height hs: ',hs)
               CALL error_mess(time,'sapwood area Asw: ',Asw)
               CALL error_mess(time,'heartwood area at stem base Ahb: ',Ahb)
            endif
            flag_end = 1
            return
          end if

          discr = SQRT(discr)
          w1    = -fp/2. + discr
          w2    = -fp/2. - discr
1313      hp1 = SQRT(w1*Ahb)
          hp2 = (Ahb+SQRT(w1*Ahb))*B
          hp3 = (w1*Ht + (Ahb+SQRT(w1*Ahb))*B)
          help = (sprhos/3.) * (w1*Ht + (Ahb+SQRT(w1*Ahb))*B)
          func  = (sprhos/3.) * (w1*Ht + (Ahb+SQRT(w1*Ahb))*B) - H
          hp4= H* precision
          IF(abs(func) <= H * precision) THEN
            Ahc = w1
            if (.not.flag_mult8910) then
                CALL error_mess(time,' positive root is a solution in calc_dbh for cohort: ',real(ident))
                CALL error_mess(time,'stop in calc_dbh for stand No: ',real(ip))
                CALL error_mess(time,'function: ',func)
            endif
            flag_end = 1
            return
          ELSE
           func  = (sprhos/3.) * (w2*Ht + (Ahb+SQRT(w2*Ahb))*B) - H
           IF(abs(func) <= H * precision) THEN
              Ahc = w2
           ELSE
            IF(precision.LT.1e-2) THEN
               precision = precision*10.
               GOTO 1313
                if (.not.flag_mult8910) then
                   CALL error_mess(time,'no valid solution found in calc_dbh for heartwood geometry for cohort: ',real(ident))
                   CALL error_mess(time,': heart wood mass, H = ',H)
                   CALL error_mess(time,': precision requirement = ',precision)
                   CALL error_mess(time,'iteration in stand No: ',real(ip))
                endif           
            ELSE
                if (.not.flag_mult8910) then
                   CALL error_mess(time,'no valid solution found in calc_dbh for heartwood geometry for cohort: ',real(ident))
                   CALL stop_mess(time,'no valid solution found in calc_dbh for heartwood geometry')
                   CALL error_mess(time,'species No: ',real(ns))
                   CALL error_mess(time,'stop in calc_dbh for stand No: ',real(ip))
                   CALL error_mess(time,'precision requirement H*precision ',H*precision)
                   CALL error_mess(time,'heart wood mass H: ',H)
                   CALL error_mess(time,'bole height b: ',b)
                   CALL error_mess(time,'height Ht: ',Ht)
                   CALL error_mess(time,'ave. sapwood height hs: ',hs)
                   CALL error_mess(time,'sapwood area Asw: ',Asw)
                   CALL error_mess(time,'heartwood area at stem base Ahb: ',Ahb)
                endif           
               flag_end = 1
              return
            ENDIF
           ENDIF
          ENDIF
          Dc  = SQRT((Ahc+Asw)*4./pi)
       END IF
       if (Ht<=h_breast) then
           DBH = 0.0
       else if (Ht>h_breast.and.B<h_breast) then
           DBH=Dc*(Ht-h_breast)/(Ht-B)
       else
           DBH=D-(D-Dc)*h_breast/B
       end if
     end select
  ELSE
       if (.not.flag_mult8910) then
           CALL error_mess(time,'no calculation of heartwood geometry for cohort (Bole height >= height trees are dead): ',real(ident))
           CALL error_mess(time,'bole height: ',b)
           CALL error_mess(time,'height: ',Ht)
       endif
  END IF ! if B > Ht

END SUBROUTINE CALC_DBH
