!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*    - Initialisation of cohorts =                              *!
!*    reads cohort information and calculates missing values     *!
!*    which are needed for stand initialisation                  *!
!*          initia                                               *!
!*          treeini                                              *!
!*          sapini                                               *!
!*          header                                               *!
!*          crown_base                                           *!
!*          crown_base_eg                                        *!
!*          fdfahc: function                                     *!
!*          ini_gener_sap                                        *!
!*          NEWTON: function numerical recipes                   *!
!*                                                               *!
!*                  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 INITIA     !
!***********************!

SUBROUTINE INITIA

! begin declaration section

      USE data_init
      USE data_par
      USE data_simul
      USE data_species
      USE data_stand
      use data_help
      IMPLICIT none

      REAL :: area  !area of database in m^2 (10000=1ha)
      INTEGER :: area_factor !factor for calculation per patch (=area/kpatchsize)
      REAL  ::   hlp_lai,share, ager
      INTEGER :: taxid, &   ! species number
                 age, &     ! tree age
                 n, &       ! number of trees
                 n_koh, &   !
                 k,     &   ! number of tree classes
                 ng_locid   ! ID stand
      INTEGER :: inunit, parunit,outunit,tmpunit,ctrlunit,listunit !units
      CHARACTER*85 zeile
      CHARACTER*80 :: infile
      CHARACTER :: source
      INTEGER :: nlines, nlines_comp, istart, fl_num, nhelp, numstand, ihelp
      INTEGER :: tax_of_BRA_id

      INTEGER,DIMENSION(:),ALLOCATABLE  :: locid_comp
      REAL rsap, cform, dummy
      REAL aux
      LOGICAL  :: select_lines
	  real standsz(10000)

      CHARACTER*40, allocatable, dimension(:)   ::helptmp
      INTEGER :: helpz

! Stand data (model initialisation)
      INTEGER baum(10),alt(10),klimid,gwa,lbanr,wgeb,lein,zei
      REAL mhoe(10),dm(10),gf(10),bon(10),en(10),psi(10)
! Parameters for missing data algorithms
      REAL p0(nspec_tree),p1(nspec_tree),p2(nspec_tree),p3(nspec_tree),p4(nspec_tree), &
       c1(nspec_tree),c2(nspec_tree),ku_a0(nspec_tree),ku_a1(nspec_tree),ku_a2(nspec_tree),&
       ku_b0(nspec_tree),ku_b1(nspec_tree),ku_b2(nspec_tree),ku_c0(nspec_tree),&
       ku_c1(nspec_tree),ku_c2(nspec_tree),wei_k1(nspec_tree),wei_k2(nspec_tree)

!     ------------------------------------------------------------------
!     INTEGER ncl  !Number of classes after classification
      integer    ncl1

      REAL dg,dmin,dmax,g,gpatch,b,c,bhd,height,hbc,hg
      REAL tot_crown_area, mixed_tot_ca, corr_la
      INTEGER pass
      REAL saquad, genDg, nbhd,x,gx,bhdmax,bhdmin,clwdth,Fint(0:100)
      REAL ku_a,ku_b,ku_c,wei_f,thdmax,p1n,p4n
	  REAL, allocatable, dimension(:)  :: nz
      REAL, allocatable, dimension(:)  :: zheigh,zbhd,zhbc
	  REAL xxr,xyr, &
           kd,  & 
           h_para,h_parb !parameter of the height function of level II sites
      INTEGER idum,anzahl, data_flag,start,baumid,dir_flag,inwahl,bz,imax
      INTEGER i,j,anzit,iz,id,icl,ios,xid,xnr,xxi,xyi, &
              bhdcl, &  !diameter classes level II
              dclmin, & !smallest diameter class level II
              ndcl, &   !amount of diameter classes of level II
              dcwdth, & !class wideness diameter classes of level II
              n_dc(30)  !stem figure of level II diameter class
      LOGICAL ehkwei, wfirst, kfirst, optimi
	   LOGICAL, allocatable, dimension(:)   ::  smaldc, bigdc
      CHARACTER*20 fnam2
      CHARACTER*5 datasets
      CHARACTER status
	  real nzsum
!     ------------------------------------------------------------------
!       ----Function----
      REAL ran0
      REAL crown_base
	  real	crown_base_eg
!     ------------------------------------------------------------------
      REAL T
      DATA T/7.0/
!     ------------------------------------------------------------------
!
! end of declaration section
!******************************************************************************
	 
	 !ncl1 = 60
     ncl1=40 
	 allocate (zheigh(ncl1), zbhd(ncl1), zhbc(ncl1), nz(ncl1))
	 allocate (smaldc(ncl1), bigdc(ncl1))
print *,' '
print *, ' *** Choice of forest stand data set: '
print *, ' 1 - Datenspeicher Waldfond'
print *, ' 2 - single tree data; classification must be performed (e.g. SILVA data)'
print *, ' 3 - Level2-data'
print *, ' 4 - already existing class file'
print *, ' 5 - FORGRA data'
print *, ' 6 - Bavarian inventory data'
WRITE(*,'(A)',advance='no') ' ***Make your choice: '
READ *, data_flag
print *,' '

clwdth=15  !set diameter class-class width 
corr_la=1. !standard value for leaf area correction in stands of high sum of crown projection areas
mixed_tot_ca=0. !sum of crown projection area for mixed stands
pass = 1 !counter for number of passes through calculation loop for mixed stands
rsap=0.3  !standard value of rsap for cases where rsap is not determined dynamically
! get unit number and open units used in all of the above cases

ctrlunit=GETUNIT()
WRITE(*,*)site_name(ip)
OPEN (ctrlunit,FILE=TRIM(site_name(ip))//'.initctrl',STATUS='replace')
      WRITE(ctrlunit,*)'# number of trees in cohort = n trees'
      WRITE(ctrlunit,*)'# age = age'
      WRITE(ctrlunit,*)'# height = H'
      WRITE(ctrlunit,*)'# height to the base of crown = Hbc'
      WRITE(ctrlunit,*)'# breast height diameter = bhd'
      WRITE(ctrlunit,*)'# sapwood fraction of trunc cross sectional area at breast height = rsap'
      WRITE(ctrlunit,*)'# trunc diameter at tree base = D'
      WRITE(ctrlunit,*)'# trunc diameter at crown base = Dc'
      WRITE(ctrlunit,*)'# sapwood cross sectional area inside bole = Asap'
      WRITE(ctrlunit,*)'# heartwood cross sectional area at crown base = Ahc'
      WRITE(ctrlunit,*)'# heartwood cross sectional area at tree base = Ahb'
      WRITE(ctrlunit,*)'# Vol for no heartwood in crown space = Vmin'
      WRITE(ctrlunit,*)'# Vol prescribed according to empiracal volume function = Vpre'
      WRITE(ctrlunit,*)'# stem vol inherent in initialisation = Veff'
      WRITE(ctrlunit,'(A150)')'#     n trees    age         H          Hbc          bhd        rsap         D           Dc          Asap       Ahc         Ahb          Vmin        Vpre        Veff'
outunit=GETUNIT()
OPEN (outunit,FILE=TRIM(treefile(ip)),STATUS='replace')

!     ------------------------------------------------------------------
!     read in parameter for the missing-data-generator:
!     bhd-distribution from Nagel & Biging (1995),
!     crown starting height from Nagel (1995), uni-height curve according to Weimann (1980) bzw. Kuleschis (1981)
   parunit=GETUNIT()
    OPEN (parunit, FILE='input/generreg.par', STATUS='old')
    do i=1,nspec_tree
      READ (parunit,*) p0(i),p1(i),p2(i),p3(i),p4(i),c1(i),c2(i),ku_a0(i),ku_a1(i),ku_a2(i), &
           ku_b0(i),ku_b1(i),ku_b2(i),ku_c0(i),ku_c1(i),ku_c2(i),wei_k1(i),wei_k2(i)
   ENDDO
   CLOSE(parunit)
!   ---------------------------------------------------------------------
inunit=GETUNIT()

SELECT CASE(data_flag)

! ****************************************************************************
! case(1) stand generation if data source is Datenspeicher Waldfond
CASE(1)
   print *, ' Forest stand data set: Datenspeicher Waldfond'
!  preliminary: here make a choice and compile
!  datasets='singl' sets the choice of the old version which uses one single
!                   set (i.e. the first one in an input file) which contains
!                   the complete imformation for the stand in one single line
!  datasets='multi' sets the choice of a version reading a file with line by
!                   line information as in the original Datenspeicher and then
!                   writes a *.ini file for many stands with individual stand
!                   information separated by lines with stand identifiers
   print*, 'choose data set (multi/singl):'
   read(*,*) datasets
   print*, ' file name (with directory):'
   read(*,'(A)') infile
   source='D'
   standsz = 0.
   OPEN (inunit, FILE=TRIM(infile), STATUS='old')
!     ------------------------------------------------------------------
!     generating standard value out of data from the data storage unit 
!     based on estimation routine from Nagel und Biging (1995), 
!     Nagel (1995) und Gerold (1990).
!     ------------------------------------------------------------------
!
!   The following variables are read from forest inventory data: 
!   Species(baum),Age(alt),Quadratic Mean Diameter(dm),Height of tree with dm(mhoe),
!   Basal area(gf),Yield Class(bon),"Ertragsniveau"(en)
!   Additional Site variables:
!   Climate station(klimid),distance of groundwater table(gwa),soil type(lbanr),
!   forest region 'Wuchsgebiet'(wgeb),last management operation(lein), number of tree layers(zei)
!   currently not used for initialisation: xid, klimid, gwa, lbanr, wgeb, lein, bon(i), en(i)
!   lbanr (check difference to declaration!),
!   check if alt and baum can be skipped as variable names and age and species directly used
!   check idendity of hg and mhoe, dg and dm, gf and g
!     ------------------------------------------------------------------

! input of data from a dataset, first row

   IF (datasets=='singl') THEN
      READ (inunit,*)xid,klimid,lbanr,gwa,wgeb,lein, &
          zei,(baum(i), alt(i),mhoe(i),dm(i),gf(i),bon(i),en(i),i=1,zei)
      ALLOCATE(ngroups(zei))
      DO i=1,zei
           IF(baum(i).EQ.8)  ngroups(i)%taxid=1
           IF(baum(i).EQ.10) ngroups(i)%taxid=2
           IF(baum(i).EQ.11) ngroups(i)%taxid=3
           IF(baum(i).EQ.15) ngroups(i)%taxid=4
           if(baum(i).eq.12) ngroups(i)%taxid = 10
! Eucalyptus
           IF(baum(i).EQ.30)  ngroups(i)%taxid=12
           IF(baum(i).EQ.31) ngroups(i)%taxid=13
            

         IF (dm(i).eq.0) dm(i) = 0.5
         IF (mhoe(i).eq.0) mhoe(i) = 1.0
         IF (gf(i).eq.0) gf(i) = 0.25
         ngroups(i)%locid=xid
         ngroups(i)%alter=alt(i)
         ngroups(i)%mhoe=mhoe(i)
         ngroups(i)%gf=gf(i)
         ngroups(i)%dm=dm(i)
         ngroups(i)%patchsize=10000

      ENDDO
      CLOSE(inunit)
      nlines=zei
      cform=1;hlp_lai=0
    ! Initialisastion of stand data: area = 1ha
      area=10000
      area_factor=int(area/kpatchsize)
    ! read file head for description, write in ini-file
      CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
   ENDIF  !block for reading of input data DSW 'singl' = specially prepared for FORSKA

! read in stand dataEinlesen out of data storage for many stands
   IF (datasets=='multi') THEN
    select_lines=.false.
    fl_num=0


      ALLOCATE(ngroups(10000))
      numstand= 0
	  nlines=1
      ngroups%taxid=0
      ngroups%schicht=-99
      DO
        READ (inunit,*,END=3333)xid,klimid,lbanr,gwa,wgeb,lein, &
            zei,(baum(i),alt(i), psi(i), mhoe(i),dm(i),gf(i),bon(i),en(i),i=1,zei)       
	    numstand = numstand +1
		ngroups(nlines)%standsize= 0
		DO i=1,zei
           IF(baum(i).EQ.5)  ngroups(nlines)%taxid=5
           IF(baum(i).EQ.8)  ngroups(nlines)%taxid=1
           IF(baum(i).EQ.10) ngroups(nlines)%taxid=2
           IF(baum(i).EQ.11) ngroups(nlines)%taxid=3
           IF(baum(i).EQ.15) ngroups(nlines)%taxid=4
           ! the following species are preliminarily assigned
           IF(baum(i).EQ.1)  ngroups(nlines)%taxid=2   ! Abies alba
           IF(baum(i).EQ.2)  ngroups(nlines)%taxid=1   ! Acer platanoides
           IF(baum(i).EQ.3)  ngroups(nlines)%taxid=1   ! Acer pseudoplatanus
           IF(baum(i).EQ.4)  ngroups(nlines)%taxid=5   ! Alnus glutinosa
           IF(baum(i).EQ.6)  ngroups(nlines)%taxid=1   ! Carpinus betulus
           IF(baum(i).EQ.7)  ngroups(nlines)%taxid=4   ! Castanea sativa
           IF(baum(i).EQ.9)  ngroups(nlines)%taxid=4   ! Fraxinus excelsior
           IF(baum(i).EQ.12) ngroups(nlines)%taxid=5   ! Populus tremula
           IF(baum(i).EQ.13) ngroups(nlines)%taxid=4   ! Quercus petraea
           IF(baum(i).EQ.14) ngroups(nlines)%taxid=4   ! Quercus pubescencs
           IF(baum(i).EQ.16) ngroups(nlines)%taxid=1   ! Tilia cordata
           IF(baum(i).EQ.17) ngroups(nlines)%taxid=4   ! Ulmus glabra
		   iF(baum(i).EQ.21) ngroups(nlines)%taxid=10  ! Douglasie
		   iF(baum(i).EQ.22) ngroups(nlines)%taxid=6   ! Larix
		   iF(baum(i).EQ.23) ngroups(nlines)%taxid=7   ! Pinus strobus
           iF(baum(i).EQ.24) ngroups(nlines)%taxid=10   ! Douglasie
           IF (dm(i).eq.0) dm(i) = 0.5
           IF (mhoe(i).eq.0) mhoe(i) = 1.0
           IF (gf(i).eq.0) gf(i) = 0.25
           ngroups(nlines)%locid=xid
           ngroups(nlines)%alter=alt(i)
           ngroups(nlines)%mhoe=mhoe(i)
           ngroups(nlines)%gf=gf(i)
           ngroups(nlines)%dm=dm(i)
		   ngroups(nlines)%patchsize=psi(i)*10000
           ngroups(nlines)%standsize=psi(i)*10000

		 nlines=nlines+1
		   standsz(numstand) = standsz(numstand) + psi(i)*10000
        ENDDO
      ENDDO ! read loop
3333 CONTINUE
     nlines=nlines-1
     WRITE(*,*) nlines,'sets of data', numstand, 'sets of stands'
      CLOSE(inunit)
    ! read in file headder for description, write into ini-file
      cform=1;hlp_lai=0
    ! initilisation for stand data: area = stand area based on fractions of areas 
      area_factor=1
      CALL header(outunit,infile,source,cform,rsap,flag_volfunc,-99.)
      WRITE(*,*) 'number of data lines:             ', nlines
       write(*,*)'number of plots for calculations: ', fl_num
   ENDIF ! block for reading input data DSW, many lines = 'multi'
   id=1
   tmpunit=getunit()
   ihelp = 1
   istart=-99
  DO iz=1,nlines
      IF(select_lines) THEN
         DO i=1,nlines_comp
            IF(locid_comp(i)==ngroups(iz)%locid) GOTO 2233
         ENDDO  ! comparison of site id to list of sites to be selected
          CYCLE
      ENDIF  ! end of site selection
2233  CONTINUE
      WRITE(*,*) iz, nlines, ngroups(iz)%locid,ngroups(iz)%schicht
      IF(datasets=='multi'.AND.(istart.NE.ngroups(iz)%locid)) THEN
		 WRITE(outunit,*) ngroups(iz)%locid,ngroups(iz)%standsize,'stand identifier, stand area'
         ihelp = ihelp +1
		 istart=ngroups(iz)%locid
      ENDIF
      IF(datasets=='multi'.AND.ngroups(iz)%taxid==0.) THEN
  
         WRITE(*,*) 'not the right species'
         GOTO 2222
      ENDIF
      IF(datasets=='multi'.AND.ngroups(iz)%schicht==20) THEN
 ! retention trees
         age=ngroups(iz)%alter
         taxid=ngroups(iz)%taxid
         height=ngroups(iz)%mhoe
         bhd=ngroups(iz)%dm
         n_koh=ngroups(iz)%baumzahl
         hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
         CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
         GOTO 2222
      ENDIF ! end special treatment of retention trees
      IF(datasets=='multi'.AND.ngroups(iz)%dm==0.) THEN
         WRITE(4444,*)'data insufficient for: ',ngroups(iz)%locid,' line: ',iz
         GOTO 2222
      ENDIF
 
      IF(datasets=='multi'.AND.ngroups(iz)%mhoe<h_sapini*0.01 .or. ngroups(iz)%gf.eq.0.) THEN
          aux = ngroups(iz)%standsize/10000.
         height=ngroups(iz)%mhoe
         n_koh=ngroups(iz)%baumzahl* aux
          age=ngroups(iz)%alter        
		  taxid = ngroups(iz)%taxid
		  WRITE(4444,*)'sapling init needed for: ',ng_locid,' line: ',iz
		call ini_gener_sap(outunit, taxid,age,height,n_koh)
          GOTO 2222
       ENDIF
	  optimi=.false.
      anzahl= 0;start=1
  
	  allocate(helptmp(10000000))
	  helptmp = ' '
 ! generation of single trees out of population mean values
    DO
	  helptmp = ' '
        IF((start==1).or.(.not.optimi))THEN
           T =7
           anzahl=0
           start=0
           wfirst=.true.
           kfirst=.true.
           WRITE(*,*)ngroups(iz)%locid,ngroups(iz)%patchsize
           age=ngroups(iz)%alter
           dg=ngroups(iz)%dm                    !quadratic mean diameter
           hg=ngroups(iz)%mhoe                  !corresponding height to dg
           taxid=ngroups(iz)%taxid              !species
           g=ngroups(iz)%gf                     !basal area/ha
           gpatch=g/area_factor                 !basal area/patch
           IF (datasets=='multi')  gpatch=g*ngroups(iz)%standsize/10000.
   ! selection of uni-height curve: Beech, Spruce, Oak calculated according to Weimann,
   ! other species of tree according to Kuleschis (vergl. Gerold 1990)
        IF (taxid==3.OR.taxid==5) THEN
            ehkwei=.false.
        ELSE
            ehkwei=.true.
        ENDIF
        IF ((dg-T).lt. 3.0)  THEN
           T=dg-4.0
        IF (T.lt.0.3) T=0.3
           ENDIF
   ! Estimation of Dmax out of dg (Gerold 1990)
        Dmax=8.2+1.8*dg-0.01*dg**2
           IF (dg.le.2) Dmax=dg+2
! calculation for the Weibull-distribution function
!     in case b or c are calcuted too small, p1 and p4 respectively have to be modified
          p1n=p1(taxid)
          IF (p1n.lt.((1.0001-p0(taxid))/Dg))  p1n=(1.0001-p0(taxid))/Dg
          b=p0(taxid)+p1n*Dg
          p4n=p4(taxid)
          IF (p4n.lt.((1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax)) p4n=(1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax
          c=p2(taxid)+p3(taxid)*Dg+p4n*Dmax

          anzit=0
          thdmax=5.0
      ENDIF ! end of introductory calculation and repetitions without optimisation

      genDg=0
      nbhd=0
      saquad=0
      bhdmax=0
      bhdmin=100
      clwdth=0
      gx=0
      idum=1
      x=0

!----------------------------
! generation of single trees
     DO
         IF (gx.gt.gpatch) exit
		 x = ran0(idum)
	     bhd=b*((T/b)**c-log(1.-x))**(1./c)
      if ( bhd.ge. 0.5*Dg) then
         IF (bhd.gt.bhdmax) bhdmax=bhd
         IF (bhd.lt.bhdmin) bhdmin=bhd
         IF ((.not. optimi) .and. (bhd.gt.(1.5*dmax))) bhd=1.5*dmax

!***height calculation according to uni-height curve
         IF (ehkwei) THEN
!           uni-height curve of Weimann (1980)
            IF (wfirst) THEN
               wei_f=wei_k1(taxid)+wei_k2(taxid)*hg
               wfirst=.false.
            ENDIF
            IF (bhd.ge.(dg-hg/2.)) THEN
               height=hg+wei_f*(log(hg-dg+bhd)-log(hg))
            ELSE
               height=(hg+wei_f*(log(hg/2.)-log(hg))-1.3)*(bhd/(dg-hg/2.))**0.5+1.3
            ENDIF
         ELSE
!           uni-height curve of Kuleschis (1981)
            IF (kfirst) THEN
               ku_a=1-(ku_a0(taxid)+ku_a1(taxid)*dg+ku_a2(taxid)*dg**2)
               ku_b=ku_b0(taxid)+ku_b1(taxid)*dg+ku_b2(taxid)*dg**2
               ku_c=ku_c0(taxid)+ku_c1(taxid)*dg+ku_c2(taxid)*dg**2
               kfirst=.false.
            ENDIF
               height=hg*(ku_a+(ku_b/(bhd+dg/2.))*dg+(ku_c/(bhd+dg/2.)**2)*dg**2)
         ENDIF
		 if(taxid.eq.10) then
! height curve after Bwinpro Douglas fir
		    height = 1.3 +(hg-1.3)*exp(-(0.199651*dg+4.63277655)*((1/bhd) - (1/dg)))
		 end if
		 if(taxid.eq.12.or. taxid.eq.13) then
! Medhurst et al. 1999
		    height = 3.665629*bhd**0.541
		  end if

! solution for small stands; tree dimensions below 3 m = rubbish
               IF (height.gt.(bhd*3.)) height=bhd*3.
               IF (height.lt.1.35) height=1.35+bhd
		 if(taxid.eq.12.or. taxid.eq.13) then
! Eucalyptus
		    hbc = crown_base_eg(height, bhd)
		 else	      
            hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
		 end if

         IF ((height-hbc).lt. 0.5) hbc= height - 0.5
         write(helptmp(nbhd+1), '(3f7.1,2i7)') bhd,height,hbc,age,taxid 
         gx=gx+1E-4*pi*(bhd/2.)**2
         nbhd=nbhd+1
         anzahl=anzahl+1
         saquad=saquad+bhd**2
       end if  ! BHD test
      ENDDO ! single tree calculation
!---calculates the generated Dg and test deviations of Dg and Dmax of the population value.
!   if deviation greater 20% a fittinf of the parameters acording to the Weibull-distribution is done
!   the standard generation is repeated in several iterations.
!---the optimisation can be shut off with optimi=.false.

      genDg=SQRT(saquad/nbhd)
      IF((.not. optimi) .or. (Dg .lt. 7)) exit
        IF(ABS(genDg-Dg).gt.(Dg/10.).or.(bhdmax-Dmax).gt. (Dmax/thdmax)) THEN
         IF (ABS(genDg-Dg).gt.(Dg/10.))THEN
            p1n=p1n*Dg/genDg
            IF (p1n.lt.((1.0001-p0(taxid))/Dg)) p1n=(1.0001-p0(taxid))/Dg
            b=p0(taxid)+p1n*Dg
         ELSE
            p4n=p4n*Dmax/bhdmax
            IF (p4n.lt.((1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax)) &
                    p4n=(1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax
            c=p2(taxid)+p3(taxid)*Dg+p4n*Dmax
         ENDIF
         anzahl=anzahl-Int(nbhd)
         anzit=anzit+1
         IF (anzit.ge.50) THEN
            IF (thdmax.eq.2) THEN
               print *,'id/zei: ',id,iz,' Optimization not  successful. Biased STAND.INI will be generated'
               optimi=.false.
            ELSE
               anzit=0
               thdmax=2.0
               b=p0(taxid)+p1(taxid)*Dg
               c=p2(taxid)+p3(taxid)*Dg+p4(taxid)*Dmax
            ENDIF
         ENDIF
      ELSE
         exit
      ENDIF

    ENDDO
! end of generation of single trees

 ! classification of single values in diameter cohorts
    clwdth=1+AINT((bhdmax-bhdmin)/ncl1)  !calculation of class widths
 !   write(4444,*) 'clwdth', clwdth, bhdmax, bhdmin, ncl1
    DO i=1,ncl1
       nz(i)=0
       zbhd(i)=0
       zheigh(i)=0
       zhbc(i)=0
    ENDDO
    DO j=1,nbhd
      read(helptmp(j), *) bhd,height,hbc,age,taxid         
        IF(height<1.3)   WRITE(4444,*)'bhd ',bhd,'height ',height,'art ',taxid
       icl=INT(bhd/clwdth)+1
           IF(icl.gt.ncl1) icl=ncl1
           nz(icl)=nz(icl)+1         !addition stem numbre of diameter classes
           zbhd(icl)=zbhd(icl)+bhd   !sum of diametes of diameter calsses
           zheigh(icl)=zheigh(icl)+height   !sum of height value of classes
           zhbc(icl)=zhbc(icl)+hbc   !sum of crown starting height of classes

    ENDDO

    deallocate(helptmp)
    tot_crown_area=0.
    DO i=1,ncl1
       IF (nz(i).ne.0) THEN
          bhd=zbhd(i)/nz(i)
          height=zheigh(i)/nz(i)
          hbc=zhbc(i)/nz(i)
          n_koh=NINT(nz(i)/area_factor)
          tot_crown_area=tot_crown_area+n_koh*PI*(MIN(spar(taxid)%crown_a*bhd+spar(taxid)%crown_b,spar(taxid)%crown_c))**2
       ENDIF
    ENDDO

    IF(tot_crown_area>1.1*kpatchsize) THEN
       corr_la=kpatchsize/tot_crown_area
    ELSE
       corr_la=1.
    ENDIF

    DO i=1,ncl1
       IF (nz(i).ne.0) THEN
          bhd=zbhd(i)/nz(i)
          height=zheigh(i)/nz(i)
          hbc=zhbc(i)/nz(i)
          n_koh=NINT(nz(i)/area_factor)
    !     --- 4C-specific calculations:
          IF(height<1.3)   WRITE(4444,*)ngroups(iz)%locid,'bhd ',bhd,'height ',height,'art ',taxid
          IF(height*100<h_sapini) THEN
             CALL sapini(outunit,taxid, height,hbc, n_koh,age)
             WRITE(4444,*)ngroups(iz)%locid,bhd,taxid
          ELSE
             CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
          ENDIF
        ENDIF
    ENDDO     !cohort loop
2222 CONTINUE
     if(datasets=='multi') then
	   IF (iz.ne.nlines.AND.datasets=='multi'.AND.(istart.NE.ngroups(iz+1)%locid)) WRITE(outunit,*) '-99.9'
	 end if
2244 CONTINUE
  ENDDO        !line loop
  CLOSE(outunit)
  CLOSE(ctrlunit)
  RETURN

! ****************************************************************************
! case(6) stand generation if data source is from Bavarian inventories
CASE(6)
   print *, ' Forest stand data set: Bavarian inventories'
   infile='/data/safe/4C/4C_input/stand/Bayernw.dat'
   source='B'
   OPEN (inunit, FILE=TRIM(infile), STATUS='old')

      listunit=GETUNIT()
      OPEN (listunit, FILE='/home/lasch/4c/v0.99e1/input/koord.txt', STATUS='old')

!     ------------------------------------------------------------------
!     generated standard values of data from data storage based on
!     estimation routines of Nagel and Biging (1995), Nagel (1995) and
!     Gerold (1990).
!     ------------------------------------------------------------------
!
!   The following variables are read from forest inventory data:
!   Species(baum),Age(alt),Quadratic Mean Diameter(dm),Height of tree with dm(mhoe),
!   Basal area(gf),Yield Class(bon),"Ertragsniveau"(en)
!      
!     ------------------------------------------------------------------

! read in stad data of multiple stands out of records
    select_lines=.true.
	datasets='multi'
    fl_num=0
      IF(select_lines) THEN
         READ(listunit,*)nlines_comp
         ALLOCATE(locid_comp(nlines_comp))
         DO i=1,nlines_comp    ! reading list of sites to be initialised
            READ(listunit,*) locid_comp(i)
         ENDDO  ! end reading list of sites to be initialised
      ENDIF  ! end of reading file with sites to be selected
      IF(select_lines) CLOSE(listunit)
      CALL assign_BAY
      CALL init_plenter_param
      READ (inunit,*)
      READ (inunit,*)nlines
      ALLOCATE(ngroups(nlines))
         istart=1
         READ(inunit,*) dummy, dummy, dummy,  ngroups(1)%locid, dummy, &
           ngroups(1)%schicht, ngroups(1)%BRAid, dummy, dummy, ngroups(1)%alter, &
           dummy, dummy, ngroups(1)%dm, ngroups(1)%mhoe, ngroups(1)%baumzahl, &
           ngroups(1)%gf, ngroups(1)%volume, dummy
         ngroups(1)%taxid=tax_of_BRA_id(ngroups(1)%BRAid)
		 ngroups(1)%standsize=40000
         IF(ngroups(1)%alter==0.OR.ngroups(1)%mhoe==0.OR.ngroups(1)%dm==0.OR.ngroups(1)%volume==0.OR.ngroups(1)%gf==0) CALL data_gap_fill_DSW(1)
      DO i=2,nlines
         READ(inunit,*) dummy, dummy, dummy,  ngroups(i)%locid, dummy, &
           ngroups(i)%schicht, ngroups(i)%BRAid, dummy, dummy, ngroups(i)%alter, &
           dummy, dummy, ngroups(i)%dm, ngroups(i)%mhoe, ngroups(i)%baumzahl, &
           ngroups(i)%gf, ngroups(i)%volume, dummy
         WRITE(*,*) 'set no', i, 'read'
         ngroups(i)%taxid=tax_of_BRA_id(ngroups(i)%BRAid)
		 ngroups(i)%standsize=40000
         ! preliminary solution: larches mapped to pine
         IF(ngroups(i)%taxid==6) ngroups(i)%taxid=3
         IF(ngroups(i)%taxid==0) THEN
 
         ELSE
            IF(ngroups(i)%alter==0.OR.ngroups(i)%mhoe==0.OR.ngroups(i)%dm==0.OR.ngroups(i)%gf==0) THEN
            WRITE(7333,*)'set ',i,'not enough data or below 1.3 m height'
!			CALL data_gap_fill_DSW(i)
            ENDIF
         ENDIF
         IF(ngroups(i)%locid.NE.ngroups(istart)%locid) THEN
		    istart=i
		    fl_num=fl_num+1
         ENDIF
      ENDDO ! readin loop for multi data-set
      CLOSE(inunit)
    ! read file headder for description, write in ini-file
      cform=1;hlp_lai=0
    ! initialisation of stand records: area =
    !                     stand area calculated according to partial areas.
      area_factor=1
      CALL header(outunit,infile,source,cform,rsap,flag_volfunc,-99.)
   id=1
   WRITE (fnam2,'(a,i1,a)') 'schicht',id,'.tmp'
   tmpunit=getunit()
   istart=-99
   DO iz=1,nlines

    ng_locid = ngroups(iz)%locid
    taxid=ngroups(iz)%taxid

      IF(select_lines) THEN
         DO i=1,nlines_comp
            IF(locid_comp(i)==ng_locid) GOTO 2255
         ENDDO  ! comparison of site id to list of sites to be selected
         CYCLE
      ENDIF  ! end of site selection
2255  CONTINUE

      IF(datasets=='multi'.AND.(istart.NE.ng_locid)) THEN

         WRITE(outunit,*) ng_locid,ngroups(iz)%standsize,'stand identifier, stand area'
         istart=ng_locid
         aux = ngroups(iz)%standsize/10000.
      ENDIF
      IF(datasets=='multi'.AND.taxid==0.) THEN
         ! solution for bushes must be found 
         WRITE(*,*) 'not the right species'
         GOTO 2277
      ENDIF
		 IF(ngroups(iz)%baumzahl<30.AND.ngroups(iz)%baumzahl>0) ngroups(iz)%schicht=5
      IF(datasets=='multi'.AND.ngroups(iz)%schicht==5) THEN
         ! retention trees can be directly initialized since they are not distributed onto different height cohorts 
         WRITE(4444,*) 'single type ',ngroups(iz)%schicht
         age=ngroups(iz)%alter
         height=ngroups(iz)%mhoe
         bhd=ngroups(iz)%dm
         n_koh=ngroups(iz)%baumzahl*aux
         hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
         CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)

         GOTO 2277
      ENDIF ! end special treatment of retention trees
      IF(datasets=='multi'.AND.ngroups(iz)%dm==0.and.ngroups(iz)%mhoe>h_sapini*0.01) THEN
         WRITE(4444,*)'data insufficient for: ',ng_locid,' line: ',iz
         GOTO 2277
      ENDIF
      IF(datasets=='multi'.AND.ngroups(iz)%mhoe<h_sapini*0.01) THEN
          height=ngroups(iz)%mhoe
          n_koh=ngroups(iz)%baumzahl* aux
          age=ngroups(iz)%alter
		call ini_gener_sap(outunit, taxid,age,height,n_koh)
          GOTO 2277
      ENDIF

        T=7
           age=ngroups(iz)%alter
           dg=ngroups(iz)%dm                    !quadratic mean diameter
           hg=ngroups(iz)%mhoe                  !corresponding height to dg
           g=ngroups(iz)%gf                     !basal area/ha
           gpatch=g*4.                          !basal area/patch
		   bz=ngroups(iz)%baumzahl*4.           !tree numbre/patch
!		   clwdth=dg/20.
           clwdth=dg/5
           
   ! selection of uni-height curve: beech, spruce, oak calculation according to Weimann,
   ! other species of trees after Kuleschis (vergl. Gerold 1990)
        IF (taxid==3.OR.taxid==5) THEN
            ehkwei=.false.
        ELSE
            ehkwei=.true.
        ENDIF
        ! zuweisen der PArameterwerte für Einheitshöhenkurve 
         IF (ehkwei) THEN
         ! uni-height curve from Weimann (1980)
               wei_f=wei_k1(taxid)+wei_k2(taxid)*hg
         ELSE
         ! uni-height curve from Kuleschis (1981)
               ku_a=1-(ku_a0(taxid)+ku_a1(taxid)*dg+ku_a2(taxid)*dg**2)
               ku_b=ku_b0(taxid)+ku_b1(taxid)*dg+ku_b2(taxid)*dg**2
               ku_c=ku_c0(taxid)+ku_c1(taxid)*dg+ku_c2(taxid)*dg**2
         ENDIF
           IF ((dg-T).lt. 3.0)  THEN
           T=dg-4.0
        IF (T.lt.0.3) T=0.3
           ENDIF
   ! Estimation of Dmax from dg (Gerold 1990)
        Dmax=8.2+1.8*dg-0.01*dg**2
           IF (dg.le.2) Dmax=dg+2
! Calculation of parameter for Weibull-distribution 
!     in case b or c is calculated too small, 
!     p1 and p4 respectively have to be modified
          p1n=p1(taxid)
          IF (p1n.lt.((1.0001-p0(taxid))/Dg))  p1n=(1.0001-p0(taxid))/Dg
          b=p0(taxid)+p1n*Dg

        Dmin = 0.1*Dg
		IF(Dg>70) Dmin = 2.*Dg - Dmax
          p4n=p4(taxid)
          IF (p4n.lt.((1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax)) p4n=(1.0005-p2(taxid)-p3(taxid)*Dg)/Dmax
          c=p2(taxid)+p3(taxid)*Dg+p4n*Dmax
          anzit=0
          thdmax=5.0

      helpz=0
      DO
  	     imax=INT((Dmax-Dmin)/clwdth)
          if(imax.gt.30) then
                 imax= 30
                 clwdth= (Dmax-Dmin)/30.
          end  if
         if(helpz.gt.50) goto 2277
         helpz= helpz + 1

         Fint(0)=0.
	     gx=0.
	     bhd=Dmin+0.5*clwdth
	     DO i = 1,imax
	        Fint(i)=1-exp(-((bhd-Dmin)/b)**c)
		    gx=gx+(Fint(i)-Fint(i-1))*bhd**2
		    bhd=bhd+clwdth
         END DO
	     gx=gx*PI/4*1e-4*bz
	     IF(ABS(gx-gpatch)>0.02*gpatch) THEN
            IF(gx>gpatch) THEN
		       c=c*gpatch/gx
            ELSE
			   IF(Dmin<0.8*Dg) THEN
			      Dmin=Dmin*1.05
               ELSE
     	          c=c*gx/gpatch
               ENDIF
		    ENDIF
	     ELSE
	        EXIT
	     ENDIF
      END DO
      bhd=Dmin+0.5*clwdth
      DO i = 1,imax
		 n_koh=NINT((Fint(i)-Fint(i-1))*bz)
         !***calculate height according to uni-height curve
         IF (ehkwei) THEN
         ! uni-height curve from Weimann (1980)
            IF (bhd.ge.(dg-hg/2.)) THEN
               height=hg+wei_f*(log(hg-dg+bhd)-log(hg))
            ELSE
               height=(hg+wei_f*(log(hg/2.)-log(hg))-1.3)*(bhd/(dg-hg/2.))**0.5+1.3
            ENDIF
         ELSE
         ! uni-height curve from Kuleschis (1981)
               height=hg*(ku_a+(ku_b/(bhd+dg/2.))*dg+(ku_c/(bhd+dg/2.)**2)*dg**2)
         ENDIF
         ! solution for small stands; tree dimensions below 3 m = rubbish
               IF (height.gt.(bhd*3.)) height=bhd*3.
               IF (height.lt.1.35) height=1.35+bhd
         hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
         IF ((height-hbc).lt. 0.5) hbc= height - 0.5
         CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
          if(fail.eq.1)  write(4444,*) 'negative root in newton', ng_locid,iz
		 bhd=bhd+clwdth
      END DO

2277 CONTINUE
     IF (iz.ne.nlines.AND. datasets=='multi'.AND.(istart.NE.ngroups(iz+1)%locid)) WRITE(outunit,*) '-99.9'
2266 CONTINUE
  ENDDO        !sign loop
  CLOSE(outunit)
  CLOSE(ctrlunit)
  RETURN

CASE(2)
334 CONTINUE
   CALL assign_DSW
   inwahl=0
   source='S'
   PRINT *, 'If you want to use SILVA data, type: 1'
   PRINT *, 'If you want to use levelII data from Sachsen, type: 2'
   PRINT *, 'If you want to use single tree data with tree class information, type: 3'
   PRINT *, ' if you want to use data like level II single tree data and generate one tree cohorts, type: 4'
   READ(*,*) inwahl
   IF (inwahl<1.OR.inwahl>4) THEN
      WRITE(*,*) 'You should use integer 1, 2,3 or 4 for the choice of data source'
      GOTO 334
   ENDIF
333 CONTINUE
   IF (inwahl==1) PRINT *, ' Forest stand data set: SILVA (classification must be performed)'
   IF (inwahl==2) PRINT *, ' Forest stand data set: levelII Sachsen (classification must be performed)'
   IF (inwahl==3) PRINT *, ' Forest stand data set: single tree data with tree type information (classification must be performed)'
   IF (inwahl==4) PRINT *, ' Forest stand data set: single tree data without clissification'
   WRITE(*,'(A)')
   WRITE(*,'(A)')' Do you want to read the input file from '
   WRITE(*,'(A)')' 1 - the Standard 4C stand directory on data/safe/4C/4C_input/stand'
   WRITE(*,'(A)')' 2 - or do you want to specify another directory?'
   WRITE(*,'(A)',advance='no') ' ***Make your choice: '
   READ(*,*) dir_flag
   IF(dir_flag.EQ.1) THEN
      WRITE(*,'(A)',advance='no')' Input file: '
      READ (*,'(A)') infile
   ELSEIF(dir_flag.EQ.2) THEN
      WRITE(*,'(A)',advance='no')' Input directory and file: '
      READ (*,'(A)') infile
   ELSE
      WRITE(*,*) 'You should use integer 1 or 2 for the choice of the input mode. Please try again!'
      GOTO 333
   ENDIF
337 CONTINUE
   cform=1;hlp_lai=0
   IF(dir_flag.EQ.1) OPEN (inunit,FILE='/data/safe/4C/4C_input/stand/'//trim(infile),STATUS='old')
   IF(dir_flag.EQ.2) OPEN (inunit,FILE=trim(infile),STATUS='old')
! initialising for stand records: area = 1ha
   area=10000
   IF(inwahl==2.OR.inwahl==3.OR.inwahl==4) THEN
! class width
      clwdth=1  !set diameter of classes width
      READ(inunit,'(a85)')zeile
      READ(inunit,*) area
      READ(inunit,'(a85)')zeile
   ENDIF
  area_factor = 1.
  kpatchsize = area

! read in file headder for descriptions, write in ini-file
   CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)

! classification of single values into diameter cohorts
   IF(inwahl==1) THEN
      READ(inunit,'(a85)')zeile
      READ(inunit,'(a85)')zeile
   ENDIF

335 CONTINUE
   DO i=1,ncl1
      nz(i)=0
      zbhd(i)=0
      zheigh(i)=0
      zhbc(i)=0
   ENDDO

nhelp=0
   DO
      IF(inwahl==1) READ(inunit,*,IOSTAT=ios)xnr,baumid,bhd,height,hbc,kd,xxr,xyr,xxi,xyi
      IF(inwahl==2.or.inwahl.eq.4) THEN
         READ(inunit,*,IOSTAT=ios)xnr,taxid,bhd,height,hbc,age
		 nhelp = nhelp+1
		 if(bhd.le.10) bhd=11.
         bhd=bhd/10.
         IF(hbc>-99.99.AND.hbc<-99.8) THEN
            hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
            IF(height-hbc<0.5) CALL error_mess(time,"crown to shallow in tree",REAL(xnr))
         ENDIF
      ENDIF
      IF(inwahl==3) THEN
         READ(inunit,*,IOSTAT=ios)xnr,taxid,bhd,height,hbc,ager,status
         IF(taxid>=100) taxid=tax_of_BRA_id(taxid)
         age = INT(ager)
         bhd=bhd/10.
         IF(hbc>-99.99.AND.hbc<-99.8) THEN
            hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
            IF(height-hbc<0.5) CALL error_mess(time,"crown to shallow in tree",REAL(xnr))
            IF((height-hbc)/height<0.5) hbc=0.5*height
            IF(bhd<=3.) hbc=0.
         ENDIF
      ENDIF
      IF (ios<0) exit
      IF (xnr==-9999) exit
      IF (inwahl==4) exit
      icl=INT(bhd/clwdth)+1
      IF(inwahl.eq.4.or.(inwahl==3.AND.status.NE.'F'.AND.status.NE.'Z'.AND.status.NE.'V'.and.status.NE.'H'.and.status.NE.'U'.and. status.NE.'B'))THEN
      ELSE
         IF(icl.gt.ncl1) icl=ncl1
         nz(icl)=nz(icl)+1       !sum stem numbre of diameter class 
         zbhd(icl)=zbhd(icl)+bhd  !sum up the diameters of a class
         zheigh(icl)=zheigh(icl)+height    !sum up height value of a class 
         zhbc(icl)=zhbc(icl)+hbc  !sum up crown startin height of a class
      ENDIF
   ENDDO
  nzsum=sum(nz)
  IF(inwahl.ne.4) THEN
   tot_crown_area=0.
   DO i=1,ncl1
      IF (nz(i).ne.0) THEN
         bhd=zbhd(i)/nz(i)
         height=zheigh(i)/nz(i)
         hbc=zhbc(i)/nz(i)
         if(hbc<0.025) hbc = 0.
         if(hbc>=0.025.and.hbc<0.05) hbc =0.05
         n_koh=NINT(nz(i)/area_factor)
         IF(inwahl==1) THEN
            SELECT CASE(baumid)
               CASE(5)
                 taxid=1
               CASE(1)
                 taxid=2
               CASE(3)
                 taxid=3
               CASE default
                 taxid=99
            END select
         ENDIF
         tot_crown_area=tot_crown_area+n_koh*PI*(MIN(spar(taxid)%crown_a*bhd+spar(taxid)%crown_b,spar(taxid)%crown_c))**2
      ENDIF
   ENDDO

    IF(tot_crown_area>1.1*kpatchsize) THEN
       corr_la=kpatchsize/tot_crown_area
    ELSE
       corr_la=1.
    ENDIF
    IF(pass==1) THEN
       mixed_tot_ca = mixed_tot_ca + tot_crown_area
    ELSE
       corr_la=kpatchsize/mixed_tot_ca
    ENDIF        

   DO i=1,ncl1
      IF (nz(i).ne.0) THEN
         bhd=zbhd(i)/nz(i)
         height=zheigh(i)/nz(i)
         hbc=zhbc(i)/nz(i)
         if(hbc<0.025) hbc = 0.
         if(hbc>=0.025.and.hbc<0.05) hbc =0.05
         n_koh=NINT(nz(i)/area_factor)
         IF(inwahl==1) THEN
            SELECT CASE(baumid)
               CASE(5)
                 taxid=1
               CASE(1)
                 taxid=2
               CASE(3)
                 taxid=3
               CASE default
                 taxid=99
            END select
         ENDIF
! --- 4C-specific calculation:
         WRITE(*,*) 'call :', taxid,bhd,height,hbc,nz(i),n_koh
         IF( height<(h_sapini/100.)) then
            call sapini(outunit,taxid, height, hbc, n_koh,age)
         ELSE
            CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
         ENDIF
       ENDIF
   ENDDO
   else if(xnr.ne.-999) then
      n_koh = 1
      print*, 'xnr:', xnr
      IF( height<(h_sapini/100.)) then
            call sapini(outunit,taxid, height, hbc, n_koh,age)
         ELSE
            CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
         ENDIF
   end if

   IF (xnr==-9999) GOTO 335
   if(inwahl==4.and.xnr==-999) then
      CLOSE(inunit)
      CLOSE(outunit)
     CLOSE(ctrlunit)
     RETURN
   end if
   if(inwahl==4) goto 335
   CLOSE(inunit)
   CLOSE(outunit)
   IF(mixed_tot_ca>1.1*kpatchsize .AND. pass == 1) THEN
      OPEN (outunit,FILE=TRIM(treefile(ip)),STATUS='replace')
      pass = 2
      GOTO 337
   ENDIF
   CLOSE(ctrlunit)
   RETURN

CASE(3)
444  print *, ' Forest stand data set: Level2-Daten'
   source='L'
   WRITE(*,'(A)')
   WRITE(*,'(A)')' Do you want to read the input file from '
   WRITE(*,'(A)')' 1 - the Standard 4C stand directory on data/safe/4C/4C_input/stand'
   WRITE(*,'(A)')' 2 - or do you want to specify another directory?'
   WRITE(*,'(A)',advance='no') ' ***Make your choice: '
   READ(*,*) dir_flag
   IF(dir_flag.EQ.1) THEN
      WRITE(*,'(A)',advance='no')' Input file: '
      READ (*,'(A)') infile
   ELSEIF(dir_flag.EQ.2) THEN
      WRITE(*,'(A)',advance='no')' Input directory and file: '
      READ (*,'(A)') infile
   ELSE
      WRITE(*,*) 'You should use integer 1 or 2 for the choice of the input mode. Please try again!'
      GOTO 444
   ENDIF
   cform=1;hlp_lai=0
   IF(dir_flag.EQ.1) OPEN (inunit,FILE='/data/safe/4C/4C_input/stand/'//trim(infile),STATUS='old')
   IF(dir_flag.EQ.2) OPEN (inunit,FILE=trim(infile),STATUS='old')
!------------------------------------------------------------------
! Read in level II data according to diamter classes
   READ(inunit,'(a85)')zeile
   READ(inunit,'(a85)')zeile
   READ(inunit,'(a85)')zeile
   READ(inunit,*)age,taxid,area, rsap, &
     dclmin, & !smallest diameter of experimentation patches
     ndcl,  &  !amount diameter class
     dcwdth    !class width
   READ(inunit,*)h_para,h_parb, &  !parameter of height function after Lockow
     (n_dc(i),i=1,ndcl) !stem numbre per diameter class
   close(inunit)
   clwdth=dcwdth

!   ---------------------------------------------------------------------

! current patch size = value specified by kpatchsize
   area_factor=int(area/kpatchsize)

! read in file headder for desciption, write into ini-file
  CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
   DO i=1,ncl1
     nz(i)=0
     zbhd(i)=0
     zheigh(i)=0
     zhbc(i)=0
   ENDDO

   bhdcl=dclmin
   DO i=1,ndcl
      bhd=bhdcl
      height=h_para*(0.01*bhd)**h_parb  !height function after regression from Lockow
      hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
      IF ((height-hbc).lt. 0.5) hbc= height - 0.5
      icl=INT(bhd/clwdth)+1
      IF(icl.gt.ncl1) icl=ncl1
      nz(icl)=nz(icl)+n_dc(i)      !sum stem numbre of diameter class 
      zbhd(icl)=zbhd(icl)+bhd*n_dc(i) !sum up diameters of a class 
      zheigh(icl)=zheigh(icl)+height*n_dc(i) !sum up height values of a class
      zhbc(icl)=zhbc(icl)+hbc*n_dc(i) !sum up crown starting height of a class
      bhdcl=bhdcl+dcwdth
   ENDDO

  smaldc(1)=.false.
   DO i=1,ncl1
      IF (smaldc(i)) THEN
         IF (i<ncl1) smaldc(i+1)=.true.
      ELSE
         IF (i<ncl1) smaldc(i+1)=.false.
         n_koh=NINT(nz(i)/area_factor)
         IF (n_koh>0) THEN
           IF (i<ncl1) smaldc(i+1)=.true.
         ENDIF
      ENDIF
   ENDDO

   bigdc(ncl1)=.false.
   DO i=ncl1,1,-1
      IF (bigdc(i)) THEN
         IF (i>1) bigdc(i-1)=.true.
      ELSE
         IF (i>1) bigdc(i-1)=.false.
         n_koh=NINT(nz(i)/area_factor)
         IF (n_koh>0) THEN
            IF (i>1) bigdc(i-1)=.true.
         ENDIF
      ENDIF
   ENDDO

   DO i=1,ncl1
      IF (nz(i).ne.0) THEN
         n_koh=NINT(nz(i)/area_factor)
         IF (n_koh==0) THEN   !if no trees in cohorte, shift trees to next class
            zbhd(i+1)=zbhd(i+1)+zbhd(i) !add diameter to sum of next class
            zheigh(i+1)=zheigh(i+1)+zheigh(i) !add height to sum of next class
            zhbc(i+1)=zhbc(i+1)+zhbc(i) !add hbc to sum of next class
            nz(i+1)=nz(i+1)+nz(i)  !add trees to next class
            nz(i)=0   !empty class
         ELSE
            bhd=zbhd(i)/nz(i)
            height=zheigh(i)/nz(i)
            hbc=zhbc(i)/nz(i)
!   --- 4C-specific calculations:
           CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
         ENDIF
      ENDIF
      IF (.not.bigdc(i+1)) exit
   ENDDO

   DO j=ncl1,(i+1),-1
      IF (nz(j).ne.0) THEN
         n_koh=NINT(nz(j)/area_factor)
         IF (n_koh==0) THEN    !if no trees in cohorte, shift trees to next class
            zbhd(j-1)=zbhd(j-1)+zbhd(j) !add diameter to sum of next class
            zheigh(j-1)=zheigh(j-1)+zheigh(j) !add height to sum of next class
            zhbc(j-1)=zhbc(j-1)+zhbc(j)  !add hbc to sum of next class
            nz(j-1)=nz(j-1)+nz(j)  !add trees to next class
            nz(j)=0 !empty class
         ELSE
            bhd=zbhd(j)/nz(j)
            height=zheigh(j)/nz(j)
            hbc=zhbc(j)/nz(j)
!   --- 4C-specific calculation:
           CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
         ENDIF
      ENDIF
      IF (.not. smaldc(i)) exit
   ENDDO
   CLOSE(outunit)
   CLOSE(ctrlunit)
   RETURN

CASE(4)

   WRITE(*,*) 'Do you want to use the standard procedure - type: S'
   WRITE(*,*) 'or Manfred Lexers input format - type: L'
   READ(*,*) source
   WRITE(*,'(A)',advance='no')' Input file: '
   READ(*,'(A)') infile
   cform=1;hlp_lai=0
   IF(flag_volfunc.EQ.0) THEN
      WRITE(*,'(A)',advance='no')' Input form factor (Default in 4C = 1): '
      READ *, cform
   ENDIF
   OPEN (inunit,FILE=TRIM(infile),STATUS='old')

   ! read in data from input-file
   IF (source=='S') THEN
      READ(inunit,*)source, taxid, rsap
      READ(inunit,*) area
      READ(inunit,*,END=10)n,k,age
      area_factor = 1.
      CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)

      !read in data
      DO i=1,k
         READ(inunit,*,END=10)bhd,height,share,hbc
            IF(hbc>-99.99.AND.hbc<-99.8) THEN
               hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
            END IF
          n_koh = NINT(share*n)
         CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
      ENDDO
   ELSE
      READ(inunit,*) area
       kpatchsize= area
      CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
      !read in data
      DO
         READ(inunit,*,iostat=ios)bhd,taxid,height,n_koh,age
         if(ios < 0) exit
         IF(height.ne.0 .AND. n_koh.ne.0) then
           IF(height<h_sapini*0.01) then
              CALL ini_gener_sap(outunit,taxid,age,height,n_koh)
           else
             hbc=crown_base(height,c1(taxid),c2(taxid),bhd)
             CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
           end if
         ENDIF

      ENDDO
   ENDIF
10 continue

PRINT*, 'Bestandesblattfläche (pro ha): ', hlp_lai*area_factor
   CLOSE(inunit)
   CLOSE(outunit)
   CLOSE(ctrlunit)

! FORGRA data input

CASE(5)
WRITE(*,'(A)',advance='no')' Input file: '
   READ(*,'(A)') infile
   cform=1;hlp_lai=0
   IF(flag_volfunc.EQ.0) THEN
      WRITE(*,'(A)',advance='no')' Input form factor (Default in 4C = 1): '
      READ *, cform
   ENDIF
   OPEN (inunit,FILE=TRIM(infile),STATUS='old')

   ! read in data from input file
   READ(inunit,*)source, rsap
   READ(inunit,*) area
   READ(inunit,*,END=20)n,k
   area_factor=int(area/kpatchsize)

   CALL header(outunit,infile,source,cform,rsap,flag_volfunc,kpatchsize)
   !read in data
   DO i=1,k
      READ(inunit,*,END=20)bhd,height,share,hbc,age,taxid
      n_koh=NINT(share*n/area_factor)
      IF(height<h_sapini) THEN
          CALL sapini(outunit,taxid, height,hbc, n_koh,age)
      ELSE
          CALL treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
      ENDIF
   ENDDO
20 CONTINUE

   CLOSE(outunit)
   CLOSE(ctrlunit)

CASE default

    PRINT *,' False number'
    RETURN

END select
WRITE(*,*) 'initialisation terminated'
	 deallocate (zheigh, zbhd, zhbc, nz)
	 deallocate (smaldc, bigdc)
     if (allocated(locid_comp))deallocate(locid_comp)

END subroutine initia

!****************************!
!*   SUBROUTINE TREEINI     *!
!****************************!

SUBROUTINE treeini(outunit,ctrlunit,taxid,source,bhd,height,hbc,n_koh,cform,rsap,age,hlp_lai,corr_la)
!  Species (taxid) must be handed over (Beech 1, Spruce 2, Pine 3, Oak 4)
!  Source is specifying data source
!  height and hbc are read in meter and is converted later to cm
!  n_koh numbre of trees in a cohort
!  -------------------------------------------------------------------------
  USE data_init
  USE data_par
  USE data_simul
  USE data_species
  USE data_stand
  USE data_help
  IMPLICIT none

!   ----VARIABLEN---
  REAL    :: bhd,height,hbc,hlp_lai,hfd,vd,VS,Vg,k1,k2,k3,hm,Ahc,Veff,dbc,corr_la
  REAL     ::  swheight,stembio,afol,asap,dbase, dcbase,volratio,d1,d2,h1,h2,a1,b0, x_ges
  INTEGER :: taxid, &   ! species number
             age, &     ! tree age
             n_koh
  INTEGER :: outunit,ctrlunit !units
  CHARACTER*85 zeile
  CHARACTER(75):: infile
  CHARACTER :: source
  REAL rsap, cform, sicrsap, lifrac, rsapfit
  INTEGER taumax, ring

!  function
      REAL newton

   sicrsap=rsap
! since the fraction of wood which is sapwood generally is not measured at the
! plots for which the model is initialized, it needs to be approximated
! the following rsap initialisation has been fitted to a pine run at Kienhorst
   rsapfit=1.-1.544e-8*age**4+4.343e-6*age**3-3.359e-4*age**2-4.557e-4*age
! estimation of rsap from average diameter increase
! attention: age of tree when first ring has been grown at 1.3 m must be estimated
! for the time being this is set to 5
! If hbc < h_breast, rsap and Asap (below) have to be calculated at lower height
   hm=height
   height=height*100
   hbc=hbc*100
   lifrac=1.-spar(taxid)%pss
   IF(age>6) THEN
      IF(hbc<h_breast) THEN
         taumax=age-INT(hbc/h_breast*5.)
      ELSE
         taumax=age-5
      ENDIF
      rsap=0.
      DO ring = 0,taumax-1
         rsap=rsap+exp(ring*log(lifrac))*(2.*(taumax-ring)-1.)
      END DO
      rsap=rsap/taumax**2
   ELSE
      rsap=1.
   ENDIF
   rsap=rsap*corr_la
!   --- calculate height of Sapwood-Pipes and stem-mass
   swheight=2.*hbc/3.+height/3.

if(taxid.ne.12. .and. taxid.ne.13) then
      if(taxid.eq.10) then
  ! after BWINpro , Bergel 1974
      hfd = (-200.31914/(height*bhd*bhd))+(0.8734/bhd) - 0.0052*log(bhd*bhd) + 7.3594/(height*bhd) + 0.46155    
   else
      k1=par_S(taxid,1)+par_S(taxid,2)*log(bhd)+par_S(taxid,3)*log(bhd)**2
      k2=par_S(taxid,4)+par_S(taxid,5)*log(bhd)+par_S(taxid,6)*log(bhd)**2
      k3=par_S(taxid,7)+par_S(taxid,8)*log(bhd)+par_S(taxid,9)*log(bhd)**2
      hfd=exp(k1+k2*log(hm)+k3*log(hm)**2)
   end if
! vd volume with SILVA equations
   vd=(hfd*pi*bhd**2)/40000
else
! Eucalyptus, Binkley et al 2002
   vd = 0.00005447*bhd**1.921157*(height/100)**0.950581
! Stape et. al 2010 Fkt. VER
   vd = (0.027*bhd**2.221*(height/100)**0.625)/500
! Stape et al 2010 Fkt ARA
   vd = (0.004*bhd**1.959*(height/100)**1.512)/500
end if

! vs volume with Eberswalde equations
  if(taxid.eq.3)   vs = exp(parEBW(10,1)+parEBW(10,2)*log(bhd)+parEBW(10,3)*log(hm))
   IF(taxid==3) vd = vs
   IF(flag_volfunc.EQ.0) THEN
      IF(source.ne.'S') stembio= swheight*spar(taxid)%prhos*cform*pi*(bhd/2.)**2
      IF(source.eq.'S') THEN
         stembio=vd*spar(taxid)%prhos*1000000
         bhd= SQRT(stembio*4/(swheight*spar(taxid)%prhos*cform*pi))
      ENDIF
!   --- seperation of sap wood and heartwood and sap wood cross section
      x_Ahb= 0.
      x_sap=rsap*stembio
      x_hrt=(1-rsap)*stembio
      asap=rsap*pi*(bhd/2.)**2
      !   --- estimation of leafe matter and leave area
      x_fol=asap*spar(taxid)%pnus
      afol=x_fol*(spar(taxid)%psla_min+0.5*spar(taxid)%psla_a)
      hlp_lai=hlp_lai+afol*n_koh
      !   --- fine root matter roughly estimated
      x_frt=x_fol
      IF(n_koh>0) WRITE(outunit,'(5f12.5,2f10.0,3i7)')x_fol,x_frt,x_sap,x_hrt,x_Ahb,height,hbc,age,n_koh,taxid
   ELSEIF(flag_volfunc.EQ.1) THEN
      IF (hbc>h_breast.AND.hbc<h_breast+h_bo_br_diff) hbc=h_breast
      IF (hbc==h_breast) dbc=bhd
      IF (hbc<h_breast) THEN
         dbc=bhd/height*(h_breast-hbc)+bhd   ! dbc = diameter at base of the crown
         asap=PI/4.*dbc**2.*rsap
      ELSE
         asap=PI/4.*bhd**2.*rsap !change Martin bhd>>dbc as written ins description and rsap weg
      ENDIF

      rsap = asap/((pi*bhd*bhd)/4)
      x_sap=spar(taxid)%prhos*asap*swheight
      ! first guess for start values of Ahc
       IF (hbc<=h_breast) THEN
          Ahc=PI/4.*dbc**2.-asap
          x_Ahb=PI/4.*(dbc*age/taumax)**2.-asap
       ELSE
          Ahc=PI/4.*bhd**2.*(1.-rsap)*0.04
          Ahc=Newton(Ahc,asap,bhd,hbc,height,Vd)

          if(fail.eq.1) return
          x_Ahb=PI/4.*((bhd-(4./PI*(asap+Ahc))**0.5*h_breast/hbc)/(1.-h_breast/hbc))**2-asap
      ENDIF
      ! Vg for test purposes = volume if no heartwood in crown space 
      Vg=1./3.*height*asap+2./3.*hbc*asap+1./3.*hbc*x_Ahb
      !   --- seperation of sap wood and heartwood and splitting of sap wood cross section
      stembio=spar(taxid)%prhos*(1./3.*height*(asap+Ahc)+1./3.*hbc*(2.*asap+x_ahb+(x_ahb*ahc)**0.5))
      volratio=1.0

   if(infile=='input/bwi2_blmwert1.prn') then 
      !Spruce
      if(taxid.eq.2)then
      !after Wirth et al. 2002 Tree physiology
       b0=-2.83958
       d1=2.55203
       d2=-0.14991
       h1=-0.19172
       h2=0.25739
       a1=-0.08278
       volratio=(exp(b0+d1*log(bhd)+d2*(log(bhd))**2+h1*log(height/100)+h2*(log(height/100))**2+a1*log(age+0.01)))/stembio
      endif
      !Pine
      if(taxid.eq.3)then
      !after Zianis et al. 2005 Silva Fennica EFI BEFs Europe
      volratio=exp(-2.6768+7.5939*(bhd/(bhd+13))+0.0151*height/100+0.8799*log(height/100))/stembio
      endif
      !for douglas fir (correction after bartelink 1996, forest ecol. manag.)
      if(taxid.eq.10)then
      volratio=exp(-3.229+1.901*log(bhd)+0.807*log(height/100))/stembio
      endif
   end if       
      x_sap=x_sap*volratio
      x_hrt=stembio*volratio-x_sap
      x_ges=x_hrt+x_sap
      x_Ahb=x_Ahb*volratio
      asap=asap*volratio

      if (x_hrt/x_ges .gt. 0.5 .and. taxid .eq. 2 .and. age .gt. 100) then !query too heigh heart wood percentage
         x_hrt=0.5*stembio*volratio
         x_sap=0.5*stembio*volratio
      endif
      if (x_hrt/x_ges .gt. 0.35 .and. taxid .eq. 3 .or. taxid .eq. 10) then !query too heigh heart wood percentage
         x_hrt=0.35*stembio*volratio
         x_sap=0.65*stembio*volratio
      endif
      Veff=(1./3.*height*(asap+Ahc)+1./3.*hbc*(2.*asap+x_ahb+(x_ahb*ahc)**0.5))*0.000001
      dbase = ((x_Ahb+asap)*4./PI)**0.5
      dcbase = ((Ahc+asap)*4./PI)**0.5
       WRITE(ctrlunit,'(2I5, 12F12.5)') n_koh,age,height,hbc,bhd,rsap,dbase,dcbase,asap,ahc,x_ahb,Vg/1000000,Vd,Veff

    !   --- estimation leaf matter and leaf area

      x_fol=asap*spar(taxid)%pnus*volratio
      afol=x_fol*(spar(taxid)%psla_min+0.5*spar(taxid)%psla_a)
      hlp_lai=hlp_lai+afol*n_koh
    !   --- fine root matter rough estimate
      x_frt=x_fol
      IF(n_koh>0) WRITE(outunit,'(5f12.5,2f10.0,3i7, 2f12.5)')x_fol,x_frt,x_sap,x_hrt,x_Ahb,height,hbc,age,n_koh,taxid, dcbase,bhd
   ENDIF
END subroutine treeini

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    SUBROUTINE SAPINI                !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  initilization of seedling cohorts with given height according to relations used in growth_seed

SUBROUTINE sapini(outunit,taxid, height, hbc, n_koh,iage)
   USE data_species
   USE data_stand
   use data_help
   IMPLICIT none
   REAL :: height,hbc,hhelp 
   INTEGER  :: outunit,n_koh ,taxid,iage
   REAL     :: x1,x2,xacc,shelp
   real     :: rtflsp, sapwood 

external sapwood
external rtflsp
! Shootbiomass kg  from height (cm), originally x_sap [mg]
hhelp = height * 100.

IF (taxid.ne.2) x_sap = exp(( LOG(hhelp)-LOG(spar(taxid)%pheight1))/spar(taxid)%pheight2)/1000000.
    IF (taxid.eq.2) THEN
        x1 = 1.
        x2 = 2.
        xacc=(1.0e-10)*(x1+x2)/2
! solve  equation for calculation of sapwood from height; determine root
        heihelp = hhelp
        hnspec = taxid
        shelp=rtflsp(sapwood,x1,x2,xacc)
        x_sap = (10**shelp)/1000000         !  transformation mg ---> kg
    ENDIF

! leaf matter
x_fol = (spar(taxid)%seeda*(x_sap** spar(taxid)%seedb))   ![kg]

! fine root matter rough estimate
x_frt = x_fol

! cross sectional area of heartwood
x_ahb = 0.
x_hrt = 0.

 IF(n_koh>0) WRITE(outunit,'(5f12.5,2f10.0,3i7)')x_fol,x_frt,x_sap,x_hrt,x_Ahb,hhelp,hbc,iage,n_koh,taxid
END subroutine sapini

FUNCTION ran0(idum)
   INTEGER idum,IA,IM,IQ,IR,MASK
   REAL ran0,AM
   PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,MASK=123459876)
   INTEGER kran
   idum=ieor(idum,MASK)
   kran=idum/IQ
   idum=IA*(idum-kran*IQ)-IR*kran
   IF (idum.lt.0) idum=idum+IM
   ran0=AM*idum
   idum=ieor(idum,MASK)
   RETURN
END
!  (C) Copr. 1986-92 Numerical Recipes Software 0)+0143$!-.

SUBROUTINE header(outunit,infile,source,cform,rsap,flag_volfunc,patchsize)
! write file headder into ini-file
  INTEGER    ::   outunit, flag_volfunc
  REAL       ::   rsap, cform, patchsize
  CHARACTER(75)  :: infile
  CHARACTER  ::   source

  WRITE(outunit,'(I1,1F12.0,A32)')flag_volfunc,patchsize,' ! = volume function, patch size'
  WRITE(outunit,'(A15,A1,A13,A80)') '! data source= ',source,' source file= ',infile
  WRITE(outunit,'(A57)') '! sapwood fraction and form factor now dynamic per cohort '
  WRITE(outunit,'(a37)')'! 4C Tree Initialization File (Stand)'
  WRITE(outunit,'(a1)')'!'
  WRITE(outunit,'(a51)')'! contains the following data (single tree values):'
  WRITE(outunit,'(a1)')'!'
  WRITE(outunit,'(a31)')'! x_fol:   foliage biomass (kg)'
  WRITE(outunit,'(a33)')'! x_frt:   fine root biomass (kg)'
  WRITE(outunit,'(a31)')'! x_sap:   sapwood biomass (kg)'
  WRITE(outunit,'(a33)')'! x_hrt:   heartwood biomass (kg)'
  WRITE(outunit,'(a65)')'! x_Ahb:   cross sectional area of heartwood at stem base (cm**2)'
  WRITE(outunit,'(a27)')'! height:  tree height (cm)'
  WRITE(outunit,'(a27)')'! x_hbole: bole height (cm)'
  WRITE(outunit,'(a27)')'! x_age:   tree age (years)'
  WRITE(outunit,'(a26)')'! n:       number of trees'
  WRITE(outunit,'(a35)')'! sp:      species (integer number)'
  WRITE(outunit,'(a33)')'! DC:      diameter at crown base'
  WRITE(outunit,'(a37)')'! DBH:      diameter at breast height'
  WRITE(outunit,'(a1)')'!'
  WRITE(outunit,'(a120)')'!     x_fol       x_frt       x_sap       x_hrt       x_Ahb     height   x_hbole  x_age      n     sp     DC         DBH'
END subroutine header

  FUNCTION crown_base(height,c1,c2,bhd)
IMPLICIT NONE
REAL crown_base
REAL height,bhd,c1,c2
     !---  estimate crown starting height according to Nagel (1995)
     crown_base=height*(1.-exp(-1.*(c1+c2*height/bhd)**2))
END function crown_base

 Function crown_base_eg(height,bhd)
 IMPLICIT NONE
 real crown_base_eg
 real height, bhd

! after Nutto etal.  2006
 crown_base_eg= -5.12 -0.407*bhd + 1.193*height
 if ( crown_base_eg.lt. 0.) crown_base_eg = 0.
END function crown_base_eg

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
SUBROUTINE fdfahc(X,F,DF,asap,bhd,hbc,height,Vd,J)
USE data_par
USE data_simul
use data_help
IMPLICIT none
REAL X,F,DF,asap,bhd,hbc,height,Vd,C,dCdX
INTEGER J
  fail=0
   IF (asap+X.LE.0) THEN
       WRITE(*,*) 'negative root at calculation C in fdfahc, program will stop'
       STOP
   ENDIF
   C=(bhd-(4./PI*(asap+X))**0.5*h_breast/hbc)/(1.-h_breast/hbc)
   dCdX=(-h_breast)/hbc/(1.-h_breast/hbc)/(4./PI*(asap+X))**0.5*2./PI
   IF (X*(PI/4.*C**2.-asap).LE.0) THEN
       fail=1
       return
   ENDIF
   F=1./3.*height*(asap+X)+1./3.*hbc*(asap+PI/4.*C**2.+(X*(PI/4.*C**2.-asap))**0.5)-Vd*1000000.
   DF=1./3.*(height+hbc*PI/2.*C*dCdX+hbc*0.5/(X*(PI/4.*C**2.-asap))**0.5*(PI/4.*C**2+X*PI/2.*C*dCdX-asap))
END subroutine fdfahc

FUNCTION NEWTON(X,asap,bhd,hbc,height,Vd)
use data_help
IMPLICIT NONE
REAL newton
REAL F,DF,X,DX,asap,bhd,hbc,height,Vd
INTEGER J,stepmax
! Newton is to be called with a start value for X
! a subroutine NEWFDF is to be included in the main program which
! calculates the value of the function and its derivative at X and
! returns them in the variables F and DF
      PARAMETER (stepmax=5000)
      DO 7 J=1,stepmax
        CALL fdfAhc(X,F,DF,asap,bhd,hbc,height,Vd,J)
         if(fail.eq.1) return
        IF(DF.EQ.0.0) THEN
          DX=0.01*X
        ELSE
          DX=F/DF
        ENDIF
        Newton=X
        IF(DX.GT.X) DX=X/2.
        X=X-DX
        IF(ABS(DX).LT.0.0005) RETURN
7     END DO
END

SUBROUTINE ini_gener_sap(outunit,taxid,age,pl_height, nplant)
 USE data_stand
 USE data_par
 USE data_species
 USE data_soil
 USE data_help
 USE data_plant
 USE data_manag
 IMPLICIT NONE
 integer    :: nplant,       &
               taxid,        &
               nclass,       &
               i,nr,         &
               age,          &
               outunit
 real       :: pl_height,    &
               height,       &
               hhelp,        &
               hbc,          &
               sdev,         &
               help,         &
               nstot
real         :: rtflsp, sapwood
real         :: hmin_est       ! empirical estimated minimum height

real, dimension(:), allocatable   :: hei,     &
                                     nschelp
integer,dimension(:),allocatable  :: nsc

external sapwood
external rtflsp

sdev = hsdev(taxid)
if (nplant.eq.0) nplant= numplant(taxid)
height = pl_height*100
if(height .lt. 100) then
  hmin_est = height - height*0.2
else
   hmin_est = height - height*0.1
end if
if(nplant.eq.1) hmin_est = height
  nclass= nint((height+2*sdev) - hmin_est) + 1
if(nplant.eq.1) nclass =1
  if(nplant.lt.200) nclass=1
  allocate(hei(nclass))
  allocate(nschelp(nclass))
  allocate(nsc(nclass))
   nstot = 0
   help =  (1/(sqrt(2*pi)*sdev))
   do i = 1, nclass
! height per class
          hei(i) = hmin_est + (i-1)
          nschelp(i) = help*exp(-((hei(i)-height)**2)/(2*(sdev)**2))
          nstot = nstot + nschelp (i)
   end do

! scaling of plant number per cohort
   do i = 1,nclass
         nsc(i) = nint((nschelp(i)*nplant/nstot) + 0.5)
   end do
 if(nplant.eq.1) nsc(1) = nplant
   do i = 1,nclass
       hhelp = hei(i)*0.01
        hbc=0
      call sapini(outunit,taxid, hhelp, hbc,nsc(i),age)
   end do
END SUBROUTINE ini_gener_sap
