!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                Stand initialisation                           *!
!*                                                               *!
!*          CONTAINS SUBROUTINES :                               *!
!*              PREPARE_STAND                                    *!
!*              internal subroutines:                            *!
!*                  SLA_INI                                      *!
!*                                                               *!
!*              CALC_INT                                         *!
!*              CALC_WEIBLA                                      *!
!*              READ_STAND (treeunit)                            *!
!*              COH_INITIAL (coh)                                *!
!*              CREATE_MISTLETOE                                 *!
!*              CREATE_SOILVEG                                   *!
!*                                                               *!
!*              CONTAiNS FUNCTIONS :                             *!
!*              SURVAGE                                          *!
!*                                                               *!
!*                  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 prepare_stand

  USE data_simul
  USE data_site
  USE data_stand
  USE data_species
  use data_climate
  use data_par
  USE data_manag

  IMPLICIT NONE

  CHARACTER      :: a, with_storage
  CHARACTER(30)  :: text
  CHARACTER(50)  :: test_stand_id
  INTEGER        :: ios,treeunit
  LOGICAL        :: exs, lstin
  INTEGER        :: help_ip, test_vf
  REAL           :: test_patchsize, xx


  REAL help_height_top  ! auxiliary var. for setting mistletoe height at uppermost crown layer
  INTEGER which_cohort
  INTEGER nr_infect_trees
  INTEGER nr_mist_per_tree
  INTEGER i
  TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

  IF(site_nr==1) THEN
      help_ip=site_nr
  ELSE
      help_ip=ip
  END IF

  pt = neu()
  anz_coh=0
  max_coh=0
  ios = -1
  nr_mist_per_tree=0

  IF(flag_stand>0) then

    exs = .false.
    stand_id = standid(help_ip)
    ! reading stand information from treefile
    inquire (File = treefile(help_ip), exist = exs)
    IF((exs .eqv. .false.) .or. (flag_stand==2)) then
      IF(exs .eqv. .false.) write(*,*) '  Stand initialization file not exists!'
      IF(flag_stand==2) write(*,*)'  Stand initialization with new file'
      write(*,'(A)',advance='no') '  Creating new file (y/n): '
      READ *, a
      IF(a.eq.'y'.or. a.eq.'Y') CALL initia

    ! planting of small trees
      if(flag_reg.eq.20) then
         call planting
         flag_reg=100
      end if
      flag_stand=1
      exs=.true.
    ENDIF
    ! read values from treefile
    IF (exs.eqv. .true.)  then
       treeunit=getunit()
       OPEN(treeunit,file=treefile(help_ip),action='read', pad='YES')  !!!! NSC implementation for ini-file
       READ(treeunit,'(A1,2X,I1,1F12.0)',iostat=ios)with_storage
       backspace treeunit
       if(with_storage .eq. 'C') then  !flag_dis .eq. 2 .and. 
          READ(treeunit,'(A1,2X,I1,1F12.0)',iostat=ios)with_storage,test_vf, test_patchsize
          call read_stand_with_nsc (treeunit)
          CLOSE(treeunit)
          kpatchsize = test_patchsize
          anz_coh       = max_coh
          coh_ident_max = anz_coh
        else
          READ(treeunit,'(I1,F12.0)',iostat=ios) test_vf, test_patchsize
  !       write(8888,*) ip, test_vf, flag_volfunc
          if(flag_multi.ne.4 .or. (flag_multi.eq.4.and.ip.eq.1) .or. (flag_multi.eq.8.and.ip.eq.1)) then
            IF(test_vf.NE.flag_volfunc) THEN
              if (.not.flag_mult8910) then
                CALL error_mess(time,"volume function in sim-file and the one used for initialisation do not match",REAL(flag_volfunc))
                CALL error_mess(time,"volume function (flag_volfunc) is set to",REAL(test_vf))
              endif
             flag_volfunc = test_vf
            ENDIF
          endif
         IF(test_patchsize .GT. 0.) THEN
           lmulti = .FALSE.
           IF(test_patchsize.NE.kpatchsize) THEN
             if (.not.flag_mult8910) then
               CALL error_mess(time,"patch size in sim-file and the one used for initialisation do not match",kpatchsize)
               CALL error_mess(time,"value in ini-file",test_patchsize)
               CALL error_mess(time,"value in sim-file",kpatchsize)
             endif
           kpatchsize = test_patchsize
           ENDIF
          ELSE
            lmulti = .TRUE.
         ENDIF
          do
           READ(treeunit,'(A)',iostat=ios) a
            IF (a .ne. '!') exit
          end do
          backspace treeunit
          ! generation of mistletoe cohort; mistletoe cohort need to be generated BEFORE tree cohorts as otherwise the light model becomes messy
          if (flag_dis.eq.1) then
            do i= 1, dis_row_nr
             if (dis_type(i) .eq. 'M') then
               if (flag_mistle.eq.0) then        !set #of mist. only once
                print *,"!! Note, implementation of mistletoe is restricted to trees of Pinus sylvestris"
                nr_mist_per_tree = dis_rel(i)
                flag_mistle=1                          ! flag indicating mistletoes
                call create_mistletoe  ! initialisation of Mistletoe
               endif
              anz_coh = max_coh
             endif
            enddo
          endif
       
          lstin = .TRUE.
          if(flag_multi.eq.4 .or. flag_multi.eq.8) stand_id = standid(help_ip)
            do while (lstin)
              IF (lmulti) THEN
               read(treeunit,*,iostat=ios) test_stand_id,  test_patchsize,text
               IF (ios .lt. 0) then
                if (.not.flag_mult8910) then
                   CALL error_mess(time,"stand identificator not found"//adjustl(stand_id)//"ip No.",real(help_ip))
                   write (*,*) '*** PREPSTAND:  program aborted'
                   write (*,*) '                stand identificator',stand_id,'  not found'
                   write (*,'(A, 2x,A)') '                 in initialisation file',treefile(help_ip)
                endif
               flag_end = 2
               return
               ENDIF
               IF (test_stand_id .ne. stand_id) THEN
                 read (treeunit,*) xx
                   do while (xx .gt. -90.0)
                    read (treeunit,*) xx
                   enddo  ! xx
                ELSE
                 lstin = .FALSE.
                 kpatchsize = test_patchsize
                 call read_stand (treeunit)
               END IF ! stand_id
              ELSE
               lstin = .FALSE.
               call read_stand (treeunit)
              END IF    ! lmulti
            end do      ! lstin
       CLOSE(treeunit)
       anz_coh       = max_coh
       coh_ident_max = anz_coh
       endif !w/o storage
    ENDIF    !exs.eqv. .true.
  END IF     !if stand >0

!if treefile not exists and not created:
IF(ios .ne. 0 .or. exs .eqv. .false.)THEN
  if (.not.flag_mult8910) PRINT *,' >>> No Stand Initialization possible '
  flag_stand=0
END IF

! Setting of height and number of mistletoe
if (flag_mistle.ne.0) then
    help_height_top=1.
    p=>pt%first
    DO WHILE (ASSOCIATED(p))
        if (p%coh%species.eq.3 .AND. p%coh%height.gt.help_height_top) then  !only on Pinus
            help_height_top=p%coh%height
            which_cohort=p%coh%ident
            nr_infect_trees=p%coh%nTreeA
        end if
        p=>p%next
    end do

    p=>pt%first
    DO WHILE (ASSOCIATED(p))
        if (p%coh%species.eq.nspec_tree+2) then
            p%coh%height  = help_height_top               !upper crown
            p%coh%x_hbole = p%coh%height-50.              !lower crown
            p%coh%nTreeA = nr_infect_trees*nr_mist_per_tree             !number of mistletoes
        end if
        if (p%coh%ident.eq.which_cohort) then             !mark uppermost tree cohort with flag mistletoe
            p%coh%mistletoe=1
        end if
        p=>p%next
    end do
end if ! end set height/number of mistletoe

! Soil Vegetation
if (flag_sveg .gt. 0) then
   call create_soilveg  ! initialisation of ground vegetation
   anz_coh = max_coh
endif

IF(flag_stand>0) CALL sla_ini
IF(flag_stand>0) CALL stand_bal_spec
CALL calc_int
CALL calc_weibla
if(flag_mg.ne.33) call overstorey

contains

SUBROUTINE sla_ini

 USE data_stand
 USE data_species

 IMPLICIT NONE
 TYPE(Coh_Obj), Pointer :: p  ! pointer to cohort list

 p => pt%first

 DO WHILE (ASSOCIATED(p))
     ns=p%coh%species
     p%coh%med_sla=spar(ns)%psla_min+spar(ns)%psla_a*0.5
     p%coh%t_leaf = p%coh%med_sla * p%coh%x_fol
     p =>p%next
 END DO
end subroutine sla_ini

end subroutine prepare_stand

!*************************************************************************

subroutine calc_int  !  calculation of intrinsic mortality rate

use data_species
implicit none
INTEGER j

do j=1,nspecies
spar(j)%intr = -log(0.01)/spar(j)%max_age
end do
end subroutine calc_int

!*************************************************************************

subroutine calc_weibla
!   calculation of parameter lamda for Weibull-distribution of sress mortality

use data_species
implicit none
INTEGER j
REAL survage

do j=1,nspecies
spar(j)%weibla = -log(0.01)/(survage(j)**weibal)
end do

end subroutine calc_weibla

!*************************************************************************

REAL function survage(ispec)
!  calculation of survival time per species depending on shade tolerance class stol

use data_species
implicit none
INTEGER :: ispec

IF(spar(ispec)%stol.eq.1) survage=20.
IF (spar(ispec)%stol.eq.2) survage=40.
IF (spar(ispec)%stol.eq.3) survage=60.
IF (spar(ispec)%stol.eq.4) survage=80.
IF (spar(ispec)%stol.eq.5) survage=100.
end function

!*************************************************************************

SUBROUTINE read_stand (treeunit)

!  Read of stand initialisation

  USE data_par
  USE data_simul
  USE data_species
  USE data_stand

  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini
  REAL           :: hdquo   ! help variable for stress initilization
  INTEGER        :: ios,treeunit, loc, i
  logical        :: treegroup_decid
  integer, dimension(5) :: decidous = (/1, 4, 5, 8, 11/)

              do
                call coh_initial (coh_ini)
                READ(treeunit,'(5f12.5,2f10.0,i7, f10.0,i7, f9.5, f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
                    coh_ini%x_Ahb, coh_ini%height, coh_ini%x_hbole, coh_ini%x_age, &
                    coh_ini%nTreeA,coh_ini%species, coh_ini%dcrb, coh_ini%diam
                IF(ios<0 .or. coh_ini%x_fol .lt. -90.0) exit

                coh_ini%nTreeD = 0.
                coh_ini%x_crt = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*spar(coh_ini%species)%cr_frac
                coh_ini%x_tb = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*(1.-spar(coh_ini%species)%cr_frac)

				!NSC Speicher initialisieren
				!hier neue version auslesen
                treegroup_decid = .False.
                do i = 1, 5
                  if (decidous(i) .eq. coh_ini%species) then
                  treegroup_decid = .True.
                  exit
                  endif
                end do
                If (treegroup_decid .eq. .True.) then
                coh_ini%x_nsc_sap = coh_ini%x_sap * decid_sap_allo * cpart !*0.5 umrechnung von kg DW zu kg C
                coh_ini%x_nsc_sap_max = coh_ini%x_nsc_sap
                coh_ini%x_nsc_tb = coh_ini%x_tb * decid_tb_allo * cpart
                coh_ini%x_nsc_tb_max = coh_ini%x_nsc_tb
                coh_ini%x_nsc_crt = coh_ini%x_crt * decid_crt_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_nsc_crt
                endif
                If (treegroup_decid .eq. .False.) then
                coh_ini%x_nsc_sap = coh_ini%x_sap * conif_sap_allo * cpart
                coh_ini%x_nsc_sap_max = coh_ini%x_nsc_sap
                coh_ini%x_nsc_tb = coh_ini%x_tb * conif_tb_allo * cpart
                coh_ini%x_nsc_tb_max = coh_ini%x_nsc_tb
                coh_ini%x_nsc_crt = coh_ini%x_crt * conif_crt_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_nsc_crt
                endif
				
                coh_ini%ident = max_coh + 1
                coh_ini%Fmax = coh_ini%x_fol
                coh_ini%x_health = 0
                coh_ini%x_hsap = 0.
                ns  = coh_ini%species
                coh_ini%N_fol=coh_ini%x_fol*spar(coh_ini%species)%ncon_fol     ! kg * mg/g --> g
                 if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.h_sapini) then
                  CALL CALC_DBH(coh_ini%x_hbole,coh_ini%height,coh_ini%x_sap,coh_ini%x_hrt,coh_ini%x_Ahb,coh_ini%Ahc,coh_ini%ident,coh_ini%diam,coh_ini%dcrb,coh_ini%x_hsap,coh_ini%asapw)
                else
                  coh_ini%x_hsap = (2*coh_ini%x_hbole + coh_ini%height)/3.
                  coh_ini%Asapw = coh_ini%x_sap/(spar(coh_ini%species)%prhos*coh_ini%x_hsap)
                end if

 ! Stress calculation
                IF (coh_ini%diam.ne. 0.) THEN
                  hdquo = coh_ini%height/ (coh_ini%diam*100)
                  IF (hdquo.gt. 1. .and. (coh_ini%x_age .gt. 10..and. coh_ini%x_age .lt.50) ) THEN
                    coh_ini%x_stress =  coh_ini%x_age/2
                  ELSE IF ( hdquo.gt. 1. .and. coh_ini%x_age .gt.50) THEN
                    coh_ini%x_stress =  coh_ini%x_age*3./7.
                  ELSE
                    coh_ini%x_stress = 0.
                  END IF
                ELSE
                  coh_ini%x_stress = 0.
                END IF  ! coh_ini

                coh_ini%x_stress = 0.
                coh_ini%nta = coh_ini%nTreeA

                IF (.not. associated(pt%first)) THEN
                  max_coh = 0
                  allocate(pt%first)
                  pt%first%coh = coh_ini
                  nullify(pt%first%next)
                ELSE
                  allocate(zeig)
                  zeig%coh = coh_ini
                  zeig%next => pt%first
                  pt%first => zeig
                END IF
                max_coh = max_coh + 1
              enddo

END SUBROUTINE read_stand

!*************************************************************************

SUBROUTINE read_stand_with_nsc (treeunit)

!  Read of stand initialisation

  USE data_par
  USE data_simul
  USE data_species
  USE data_stand


  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini
  REAL           :: hdquo   ! help variable for stress initilization
  INTEGER        :: ios,treeunit, loc, i
  logical        :: treegroup_decid
  integer, dimension(5) :: decidous = (/1, 4, 5, 8, 11/)
  character      :: a

              do
               READ(treeunit,'(A)',iostat=ios) a
               IF (a .ne. '!') exit
              end do
              backspace treeunit
              do
                call coh_initial (coh_ini)
!                READ(treeunit,'(5f12.5,2f10.0,i7, f7.0,i7, 2f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
                READ(treeunit,'(5f12.5,2f10.0,i7,f7.0,i7, 5f12.5)',iostat=ios) coh_ini%x_fol, coh_ini%x_frt, coh_ini%x_sap, coh_ini%x_hrt, &
                    coh_ini%x_Ahb, coh_ini%height, coh_ini%x_hbole, coh_ini%x_age, &
                    coh_ini%nTreeA,coh_ini%species, coh_ini%dcrb, coh_ini%diam, coh_ini%x_nsc_tb, coh_ini%x_nsc_crt, coh_ini%x_nsc_sap 

                IF(ios<0 .or. coh_ini%x_fol .lt. -90.0) exit

                coh_ini%nTreeD = 0.
                coh_ini%x_crt = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*spar(coh_ini%species)%cr_frac
                coh_ini%x_tb = (coh_ini%x_sap + coh_ini%x_hrt) * spar(coh_ini%species)%alphac*(1.-spar(coh_ini%species)%cr_frac)

				!hier nur den NSC Max-Speicher initialisieren
				!hier neue version auslesen
                treegroup_decid = .False.
                do i = 1, 5
                  if (decidous(i) .eq. coh_ini%species) then
                  treegroup_decid = .True.
                  exit
                  endif
                end do
                If (treegroup_decid .eq. .True.) then
                coh_ini%x_nsc_sap_max = coh_ini%x_sap * decid_sap_allo * cpart !*0.5 umrechnung von kg DW zu kg C
                coh_ini%x_nsc_tb_max = coh_ini%x_tb * decid_tb_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_crt * decid_crt_allo * cpart
                endif
                If (treegroup_decid .eq. .False.) then
                coh_ini%x_nsc_sap_max = coh_ini%x_sap * conif_sap_allo * cpart
                coh_ini%x_nsc_tb_max = coh_ini%x_tb * conif_tb_allo * cpart
                coh_ini%x_nsc_crt_max = coh_ini%x_crt * conif_crt_allo * cpart
                endif
!      IF(coh_ini%species==3.and.coh_ini%height.le.900) then
!      coh_ini%x_stress  =6
!      ELSE IF(coh_ini%species==1.and.coh_ini%height.le.2500) then
!      coh_ini%x_stress  =25
!      ELSE
!      coh_ini%x_stress = 0
!      ENDIF
                coh_ini%ident = max_coh + 1
                coh_ini%Fmax = coh_ini%x_fol
                coh_ini%x_health = 0
                coh_ini%x_hsap = 0.
                ns  = coh_ini%species
                coh_ini%N_fol=coh_ini%x_fol*spar(coh_ini%species)%ncon_fol     ! kg * mg/g --> g
            ! calculate diameter at breast height and pipe length by call of subroutine in partitio
                 if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.h_sapini) then
!                if (coh_ini%dcrb.eq.0..and.coh_ini%diam.eq.0..and.coh_ini%height.gt.137.) then
                  CALL CALC_DBH(coh_ini%x_hbole,coh_ini%height,coh_ini%x_sap,coh_ini%x_hrt,coh_ini%x_Ahb,coh_ini%Ahc,coh_ini%ident,coh_ini%diam,coh_ini%dcrb,coh_ini%x_hsap,coh_ini%asapw)
                else
                  coh_ini%x_hsap = (2*coh_ini%x_hbole + coh_ini%height)/3.
                  coh_ini%Asapw = coh_ini%x_sap/(spar(coh_ini%species)%prhos*coh_ini%x_hsap)
                end if

 ! Stress calculation
                IF (coh_ini%diam.ne. 0.) THEN
                  hdquo = coh_ini%height/ (coh_ini%diam*100)
                  IF (hdquo.gt. 1. .and. (coh_ini%x_age .gt. 10..and. coh_ini%x_age .lt.50) ) THEN
                    coh_ini%x_stress =  coh_ini%x_age/2
                  ELSE IF ( hdquo.gt. 1. .and. coh_ini%x_age .gt.50) THEN
                    coh_ini%x_stress =  coh_ini%x_age*3./7.
                  ELSE
                    coh_ini%x_stress = 0.
                  END IF
                ELSE
                  coh_ini%x_stress = 0.
                END IF  ! coh_ini
! provisorisch stress auf Null setzen
                coh_ini%x_stress = 0.

                coh_ini%nta = coh_ini%nTreeA

                IF (.not. associated(pt%first)) THEN
                  max_coh = 0
                  allocate(pt%first)
                  pt%first%coh = coh_ini
                  nullify(pt%first%next)
                ELSE
                  allocate(zeig)
                  zeig%coh = coh_ini
                  zeig%next => pt%first
                  pt%first => zeig
                END IF
                max_coh = max_coh + 1
              enddo

END SUBROUTINE read_stand_with_nsc

!*************************************************************************



SUBROUTINE coh_initial (coh_ini)

  USE data_simul
  USE data_soil
  USE data_stand
  USE data_species

  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini

        coh_ini%nTreeA = 0.
        coh_ini%nTreeD = 0.
        coh_ini%nTreeM = 0.
        coh_ini%nTreet = 0.
        coh_ini%nta    = 0.
        coh_ini%mistletoe = 0

        coh_ini%x_age  = 0.
        coh_ini%x_fol  = 0.
        coh_ini%x_sap  = 0.
        coh_ini%x_frt  = 0.
        coh_ini%x_hrt  = 0.
        coh_ini%x_crt  = 0.
        coh_ini%x_tb   = 0.
        coh_ini%x_hsap = 0.
        coh_ini%x_hbole= 0.
        coh_ini%x_Ahb  = 0.

        coh_ini%x_stress = 0
        coh_ini%x_health = 0

        coh_ini%bes     = 0.
        coh_ini%med_sla = 0.
        coh_ini%Fmax    = 0
        coh_ini%totBio  = 0.
        coh_ini%Dbio    = 0.
        coh_ini%height  = 0.
        coh_ini%deltaB  = 0.
        coh_ini%dcrb    = 0.
        coh_ini%diam    = 0.
        coh_ini%assi    = 0.
        coh_ini%LUE     = 0.
        coh_ini%resp    = 0.
        coh_ini%netAss  = 0.
        coh_ini%NPP     = 0.
        coh_ini%weekNPP = 0.
        coh_ini%NPPpool = 0.
        coh_ini%t_Leaf  = 0.
        coh_ini%geff    = 0.
        coh_ini%Asapw   = 0.
        coh_ini%crown_area = 0.

        coh_ini%BG        = 0.
        coh_ini%leafArea  = 0.
        coh_ini%sleafArea = 0.
        coh_ini%FPAR      = 0.
        coh_ini%antFPAR   = 0.
        coh_ini%Irel      = 0.

        coh_ini%totFPAR  = 0
        coh_ini%IrelCan  = 0
        coh_ini%botLayer = 0
        coh_ini%topLayer = 0
        coh_ini%survp    = 0.
        coh_ini%rel_fol  = 0.
        coh_ini%gfol     = 0.
        coh_ini%gfrt     = 0.
        coh_ini%gsap     = 0.
        coh_ini%sfol     = 0.
        coh_ini%sfrt     = 0.
        coh_ini%ssap     = 0.
        coh_ini%grossass = 0.
        coh_ini%maintres = 0.
        coh_ini%respsap  = 0.
        coh_ini%respfrt  = 0.
        coh_ini%respbr   = 0.

        coh_ini%height_ini = 0.
        coh_ini%ca_ini     = 0.

        coh_ini%rel_dbh_cl = 0
        coh_ini%underst    = 0

        coh_ini%fol_inc     = 0.
        coh_ini%fol_inc_old = 0.
        coh_ini%bio_inc     = 0.
        coh_ini%stem_inc    = 0.
        coh_ini%frt_inc     = 0.
        coh_ini%notViable   = .FALSE.

        coh_ini%intcap   = 0.
        coh_ini%prel     = 0.
        coh_ini%interc   = 0.
        coh_ini%prelCan  = 0.
        coh_ini%interc_st= 0.
        coh_ini%aev_i    = 0.
        coh_ini%demand   = 0.
        coh_ini%supply   = 0.
        coh_ini%watuptc  = 0.
        coh_ini%gp       = 0.
        coh_ini%drIndd   = 0.
        coh_ini%drIndPS  = 0.
        coh_ini%drIndAl  = 0.
        coh_ini%nDaysGr  = 0
        coh_ini%isGrSDay = .false.

        coh_ini%litC_fol  = 0.
        coh_ini%litC_fold = 0.
        coh_ini%litN_fol  = 0.
        coh_ini%litN_fold = 0.
        coh_ini%litC_frt  = 0.
        coh_ini%litC_frtd = 0.
        coh_ini%litN_frt  = 0.
        coh_ini%litN_frtd = 0.
        coh_ini%litC_stem = 0.
        coh_ini%litN_stem = 0.
        coh_ini%litC_tb   = 0.
        coh_ini%litC_crt  = 0.
        coh_ini%litC_tbcd = 0.
        coh_ini%litN_tb   = 0.
        coh_ini%litN_crt  = 0.
        coh_ini%litN_tbcd = 0.
        coh_ini%Nuptc_c   = 0.
        coh_ini%Nuptc_d   = 0.
        coh_ini%Ndemc_d   = 0.
        coh_ini%RedNc     = 1.
        coh_ini%N_pool    = 0.
        coh_ini%N_fol     = 0.
        coh_ini%wat_mg    = 0.   ! soley forflag_wred=9

        coh_ini%nroot   = 0
        coh_ini%shelter = 0
		coh_ini%day_bb  = 0
        coh_ini%x_nsc_sap = 0.
        coh_ini%x_nsc_tb = 0.
        coh_ini%x_nsc_crt = 0.
        coh_ini%x_nsc_sap_max = 0.
        coh_ini%x_nsc_tb_max = 0.
        coh_ini%x_nsc_crt_max = 0.
		
      if (coh_ini%species .ne. nspec_tree+2) then  ! no root allocation for mistletoe
        allocate (coh_ini%frtrel(nlay))
        allocate (coh_ini%frtrelc(nlay))
        if (flag_wred .eq. 9) then
            allocate (coh_ini%rld(nlay))
            coh_ini%rld  = 0.
		endif
		allocate (coh_ini%rooteff(nlay))
	    coh_ini%frtrel  = 0.
		coh_ini%rooteff = 0.
      end if   ! end exclude mistletoe
END SUBROUTINE coh_initial
!*************************************************************************
SUBROUTINE create_mistletoe
  USE data_plant
  USE data_simul
  USE data_species
  USE data_stand
  USE data_climate
  USE data_soil
  USE data_species
  USE data_par
  IMPLICIT NONE
  TYPE(cohort)   :: coh_ini
  real            :: help_height_top, help_height_bot
  REAL, EXTERNAL  ::  fi_lf, dfi_lf, ddfi_lf

 ! initialising of cohort of mistletoe
 call coh_initial (coh_ini)
 ! set mistletoe here to 20 m height, will be changed after, when cohorts of trees will be initialised
 help_height_top=2000 
 help_height_bot=help_height_top-50 
 ! following values are from sample calcul. of 10 year old V.austr. from Pfiz 2010
 coh_ini%ident      = max_coh + 1
 coh_ini%species    = nspec_tree+2    ! Species = species after all tree species and ground veg.
 coh_ini%nTreeA     = 1               ! #of mistletoes, to be read-in in management file
 coh_ini%nTreeD     = 0               ! dead trees
 coh_ini%nta        = coh_ini%nTreeA  ! alive trees internal calc.
 coh_ini%x_age      = 10
 coh_ini%x_fol      = mistletoe_x_fol ! fol biomass per tree [kg DW/tree], 1 Viscum (10years) see Pfiz 2010
 coh_ini%x_sap      = 0.              ! set near-zero for partitioning
 coh_ini%x_frt      = 0.              ! set near-zero for partitioning
 coh_ini%height     = help_height_top ! highest_layer   ! highest_layer of all cohorts
 coh_ini%x_hbole    = help_height_bot !
 coh_ini%med_sla    = 0.              ! average cohort specific leaf area [m2/kg] is being calculated internal
 coh_ini%Fmax       = 0               ! anual change of leaf biomass, for now: now change
 coh_ini%crown_area = 0.0189          ! max. projected crown area (m2) per individuum, calculated from Pfiz 2010
 coh_ini%t_leaf     = coh_ini%med_sla* coh_ini%x_fol      !leaf area per tree [m2]  !
 coh_ini%day_bb     = 1               ! evergreen
! no partitioning of NPP into stem/leaf etc.
! no root allocation
 allocate(zeig)
  zeig%coh = coh_ini
  zeig%next => pt%first
  pt%first => zeig
 max_coh = max_coh + 1
END SUBROUTINE create_mistletoe

    !*************************************************************************

SUBROUTINE create_soilveg

!  Read of stand initialisation

  USE data_plant
  USE data_simul
  USE data_species
  USE data_stand
  USE data_climate
  USE data_soil

  IMPLICIT NONE

  TYPE(cohort)   :: coh_ini

  real            :: lai_help, irel_help, FRsum
  integer         :: age_stand, nr, j
  integer         :: flag_SV_allo, rnum
  real			  :: troot2

  REAL, EXTERNAL  ::  fi_lf, dfi_lf, ddfi_lf
  
  age_stand = 0
  lai_help = 0.
  irel_help = 0.
  call wclas(waldtyp)

	   zeig=>pt%first
       
	   DO WHILE (ASSOCIATED(zeig))
	     ns = zeig%coh%species
	     lai_help = lai_help + zeig%coh%ntreea*zeig%coh%x_fol* spar(ns)%psla_min
		  age_stand = MAX(zeig%coh%x_age,age_stand)
		 zeig=>zeig%next
       end do

   IF((flag_stand==0 .or. age_stand .le. 5) .AND. flag_sveg ==2) THEN
      NPP_est = 10.  
   ELSE if(age_stand.le.5) then
      if(ns.eq.4) then
	     NPP_est = 5
	   else
	      NPP_est = 10.
	   end if
   ELSE if(flag_reg.ne.0) then
         NPP_est = 10
   ELSE
      lai_help = lai_help/kpatchsize
      irel_help = exp(-0.5*lai_help)
      if( svar(nspec_tree+1)%RedN .lt.0.) then
	         NPP_est = irel_help * med_rad1 * 365./100. *0.5
	   else
             NPP_est = irel_help * med_rad1 * 365./100. *0.5 * svar(nspec_tree+1)%RedN
      end if
  ENDIF

     call coh_initial (coh_ini)

       coh_ini%species = nspec_tree+1    ! numbre of species determined automatically
         ns = coh_ini%species
         flag_SV_allo=1
  IF(flag_SV_allo==0) THEN
  ! the parameters pdiam in the species.par file are used for allocation fractions
         coh_ini%x_sap = spar(ns)%pdiam3 * NPP_est/1000.*kpatchsize
		   coh_ini%x_fol = spar(ns)%pdiam1 * NPP_est/1000.*kpatchsize
		   coh_ini%x_frt = spar(ns)%pdiam2 * NPP_est/1000.*kpatchsize    
  ELSE  
     FRsum=0.8*NPP_est/1000.    ! start value as fraction of NPP in kg DM m-2   
     CALL newt (FRsum, fi_lf, dfi_lf, ddfi_lf, 0.001, 100, rnum)
     IF(rnum==-1) THEN
        if (.not.flag_mult8910) CALL error_mess(time,'no solution found for allocation for groundvegetation cohort: ',real(ns))
         coh_ini%x_sap = spar(ns)%pdiam3 * NPP_est/1000.*kpatchsize
		   coh_ini%x_fol = spar(ns)%pdiam1 * NPP_est/1000.*kpatchsize
		   coh_ini%x_frt = spar(ns)%pdiam2 * NPP_est/1000.*kpatchsize    
     ELSE
         coh_ini%x_sap = (ksi*FRsum**kappa)*kpatchsize
		   coh_ini%x_fol = (FRsum/2.)*kpatchsize
		   coh_ini%x_frt = (FRsum/2.)*kpatchsize    
     ENDIF
  ENDIF 

         coh_ini%height  = 60.
         coh_ini%x_age   = 1
         coh_ini%nTreeA  = 1
         coh_ini%ident   = max_coh + 1
         coh_ini%Fmax    = coh_ini%x_fol
         coh_ini%med_sla = spar(coh_ini%species)%psla_min + spar(coh_ini%species)%psla_a*irel_help
         coh_ini%t_leaf  = coh_ini%med_sla* coh_ini%x_fol      ! [m2]

         coh_ini%nta     = coh_ini%nTreeA
         coh_ini%ca_ini  = kpatchsize
         coh_ini%day_bb  = 100            ! assumption budding on 8.April 

! root allocation
                IF (.not. associated(pt%first)) THEN
                  max_coh = 0
                  allocate(pt%first)
                  pt%first%coh = coh_ini
                  nullify(pt%first%next)
				  call root_depth (1, pt%first%coh%species, pt%first%coh%x_age, pt%first%coh%height, pt%first%coh%x_frt, pt%first%coh%x_crt, nr, troot2, pt%first%coh%x_rdpt, pt%first%coh%nroot)
                  pt%first%coh%nroot = nr
                  do j=1,nr
                        pt%first%coh%rooteff = 1.   ! assumption for the first use
                  enddo
                  do j=nr+1, nlay
                       pt%first%coh%rooteff = 0.   ! layers with no roots
                  enddo

                ELSE
                  allocate(zeig)
                  zeig%coh = coh_ini
                  zeig%next => pt%first
                  pt%first => zeig
				  call root_depth (1, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
                  zeig%coh%nroot = nr
                  do j=1,nr
                          zeig%coh%rooteff = 1.   ! assumption for the first use
                  enddo
                   do j=nr+1, nlay
                      zeig%coh%rooteff = 0.   ! layers with no roots
                   enddo

                END IF
                max_coh = max_coh + 1

END SUBROUTINE create_soilveg

!*************************************************************************

!***************************!
! FUNCTION fi_lf           *!
!***************************!

REAL FUNCTION fi_lf(x)
  USE data_stand
  USE data_plant
  USE data_species
  REAL :: x
  fi_lf = spar(nspec_tree+1)%pss*ksi*x**kappa + (spar(nspec_tree+1)%psf+spar(nspec_tree+1)%psr)/2.*x - NPP_est/1000.
END ! FUNCTION fi_lf

!***************************!
! FUNCTION dfi_lf          *!
!***************************!

REAL FUNCTION dfi_lf(x)
  USE data_stand
  USE data_plant
  USE data_species
  REAL :: x
  dfi_lf = spar(nspec_tree+1)%pss*ksi*kappa*x**(kappa-1.) + (spar(nspec_tree+1)%psf+spar(nspec_tree+1)%psr)/2. 
END ! FUNCTION dfi_lf

!***************************!
! FUNCTION ddfi_lf         *!
!***************************!

REAL FUNCTION ddfi_lf(x)
  USE data_stand
  USE data_plant
  USE data_species
  REAL :: x
  ddfi_lf = spar(nspec_tree+1)%pss*ksi*kappa*(kappa-1.)*x**(kappa-2.) 
END ! FUNCTION ddfi_lf
