!*****************************************************************!
!*                                                               *!
!*                     4C (FORESEE)                              *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines:                               *!
!*             PREPARE_SITE and PREPARE_CLIMATE                  *!
!*                                                               *!
!* Contains subroutines:                                         *!
!*                                                               *!
!* PREPARE_SITE:                                                 *!
!*            preparation of site specific simulation parameters *!
!*                                                               *!
!*            contains internal subroutines:                     *!
!*            SITEMENU: choice of inputs                         *!
!*            EDITFILE: edit filenames                           *!
!*            READSOIL: input of soil parameter                  *!
!*            READCN: input of C-N-parameter                     *!
!*            READVALUE: input of start values for               *!
!*                       soil water and C-N-modeling             *!
!*            ALLOC_SOIL: allocate soil variables                *!
!*            STAND_BAL_INI: allocate and init stand variables   *!
!*            CONTROL_FILE: saving all parameters                *!
!*                          and start conditions for each site   *!
!*                                                               *!
!* READDEPO:  reading deposition data                            *!
!* READREDN:  reading values of redN                             *!
!* READLIT:   reading initialisation data of litter fractions    *!
!*                                                               *!
!* PREPARE_CLIMATE: reading of site specific climate input data  *!
!*                  from file                                    *!
!*                  contains internal subroutines:               *!
!*                  READ_DWD                                     *!
!*                  READ_CLI                                     *!
!*                  CLIMFILL                                     *!
!*                                                               *!
!* STORE_PARA:      multi run - restore of changed parameter     *!
!*                                                               *!
!*                  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_site

! input of site specific data

  use data_climate
  use data_inter
  use data_manag
  use data_mess
  use data_out
  use data_par
  use data_simul
  use data_site
  use data_soil
  use data_soil_cn
  use data_species
  use data_stand
  use data_tsort
  use data_frost


implicit none
integer i,ios,help, help_ip
character a
character :: text
character(10)  :: helpsim, text2
logical:: ex=.TRUE.
real parerr
real, external :: avg_sun_incl
character(100)  :: helpx

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' prepare_site'

WRITE(helpsim,'(I4)') anz_sim
read(helpsim,*) anh


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

! Initialization of climate data
IF (flag_clim==1 .or. ip==1 .or. flag_multi .eq.5) THEN
    call prepare_climate
END IF

if (flag_end .gt. 0) return
ios=0; help=0
do
  if (ip==1 .and. flag_mult9) then
   if (flag_trace) write (unit_trace, '(I4,I10,A,I3,A5,L5)') iday, time_cur, ' prepare_site ip=',ip,'  ex=',ex  
    call readspec
    call readsoil  ! reading soil parameter
    IF (flag_end .gt.0) return
    if (flag_soilin .eq. 0) call readvalue  ! Initialization of simulation start values for soil layers
    
    ! biochar
    if (flag_bc .gt. 0) call bc_appl
    
    ! Deposition data
    call readdepo

    ! Input redN
    if (flag_multi .ne. 4 .or. flag_multi .ne. 8 ) call readredN

    flag_mult9 = .FALSE.
  else
   if (flag_trace) write (unit_trace, '(I4,I10,A,I3,A5,L5)') iday, time_cur, ' prepare_site ip=',ip,'  ex=',ex   

    ! Deposition data
    call readdepo

    select case (flag_multi)
    case (1,6)
       call readspec
       if (flag_soilin .eq. 0) call readvalue  ! Initialization of simulation start values for soil layers
       call readredN           ! Input redN
       call readsoil  ! reading soil parameter   
       
       do
          jpar = jpar + 1
             if (vpar(jpar) .gt. -90.0) then
                helpx = simpar(jpar)
                call store_para(vpar(jpar), helpx, parerr)
                IF (parerr .ne. 1.) then
                   CALL  error_mess(time,'parameter variation: '//trim(simpar(jpar))//' not found',vpar(jpar))
                   write (*,*) '***  parameter variation: ', trim(simpar(jpar)), ' not found,   see error log'
                endif
             else
                exit
             endif
       enddo
    
    ! biochar
      if (flag_bc .gt. 0 .or. flag_decomp .gt. 100) call bc_appl

    case (2,4)


       call readsoil  ! reading soil parameter
       if (flag_soilin .eq. 0) call readvalue  ! Initialization of simulation start values for soil layers

    case (5)
       call readspec   
       call readsoil  
       if (flag_soilin .eq. 0) call readvalue  ! Initialization of simulation start values for soil layers
       call readredN                           ! Input redN

    case (7)
       call assign_co2par      
       call readsoil                           ! reading soil parameter
       if (flag_soilin .eq. 0) call readvalue  ! Initialization of simulation start values for soil layers
       call readredN                           ! Input redN

    case (8, 9, 10)

       call readsoil         ! reading soil parameter
       IF (flag_end .gt.0) return
       call readredN         ! Input redN or test resp.
    end select
  endif
exit
enddo

! Setting flag_inth and prec_stad_red from flag_int
if (flag_int .lt. 1000) then
    flag_inth = flag_int
else
  ! Conversion character ==> number and vice versa  
    write (helpsim,'(I4)') flag_int
    text2 = helpsim(2:2)
    read (text2,*) flag_inth
    text2 = helpsim(3:4)
    read (text2,*) prec_stand_red
endif

if (.not.flag_mult8910) then
    unit_soil = getunit()
    open (unit_soil,file=trim(dirout)//trim(site_name(help_ip))//'_soil.ini'//anh,status='replace')
    WRITE (unit_soil,'(2A)') '! Soil initialisation, site name: ',site_name(help_ip)
endif

call stand_bal_ini !allocation of stand summation variables

! Initialization of CO2
call assign_co2par
! Initialisation litter compartments
call readlit
! Initialization of soil model with profile data
call soil_ini    ! Aufruf ohne s_cn_ini
! Initialization disturbances
IF (flag_dis .eq. 1 .or. flag_dis .eq. 2) CALL dist_ini
! Initialization of stand
call prepare_stand
IF (flag_end .gt.0) return
! calculation of latitude in radians
xlat = lat/90.*pi*0.5
! calculation of average sun inclination
avg_incl = AVG_SUN_INCL(lat)   ! degrees
beta=avg_incl*PI/180           ! radians
! read externally prescribed bud burst days
CALL readbudb
! Initialization management
IF(flag_mg.ne.0.and. flag_mg.ne.5) call manag_ini
IF(flag_mg.eq. 5) then
    thin_dead = 1
    allocate(thin_flag1(nspec_tree))
	thin_flag1  = 0
end if	  
! Initialization of output file per site
call prep_out
call stand_balance
call CROWN_PROJ
call standup
call root_ini   ! initialisation of root distribution
call s_cn_ini

! Initialization of soil temperature model with stand data
call s_t_ini

! control file for saving simulation environment
! output of first Litter-Input at start
if(flag_mult8910 .and. (anz_sim .gt. 1)) then
    continue
else 
    IF ((ip .eq. 1 .or. flag_multi .eq. 1 .or. flag_multi .eq. 6) .and. (time_out .ne. -2) ) call control_file
endif

! hand over of the litter-initialising
call litter
if ((flag_decomp .eq. 20) .or. (flag_decomp .eq. 21)) then
    call testfile(valfile(ip),ex)
    if (ex .eqv. .true.) then
       ios = 0
       unit_litter = getunit()
       open(unit_litter,file=valfile(ip),status='old',action='read')
       if (flag_multi .ne. 9) print *,' *** Open file of litter input data ',valfile(ip),'...'
           do
              read(unit_litter,*) text
              IF(text .ne. '!')then
                 backspace(unit_litter);exit
              endif
           enddo
    endif
endif
call cn_inp

! read flux data
if (flag_eva .gt.10) call evapo_ini
      ! yearly output
      IF (time_out .gt. 0) THEN
         IF (mod(time,time_out) .eq. 0) CALL outyear (1)
         IF (mod(time,time_out) .eq. 0) CALL outyear (2)
      ENDIF

    contains

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

subroutine readsoil ! Input of soil parameter

use data_par
use data_soil_t
use data_site

implicit none

integer :: inunit, helpnl, helpnr, ihelp
real helpgrw, hlong, hlat
character :: text
character(30) :: hor, boart, helpid

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' readsoil'

! Setting of flag_surf from flag_cond
select case (flag_cond)

    case (0,1,2,3)
       flag_surf = 0

    case (10,11,12,13)
       flag_surf = 1

    case (30,31,32,33)
       flag_surf = 3
       
end select

! Setting of flag_bc from flag_decomp
if (flag_decomp .ge. 100) then
    flag_decomp = flag_decomp - 100
    flag_bc = 1
else
    flag_bc = 0
endif

call testfile(sitefile(ip),ex)
IF (ex .eqv. .true.) then
  inunit = getunit()
  ios=0
  open(inunit,file=sitefile(ip),iostat=ios,status='old',action='read')
  if (.not.flag_mult8910) then
      print *,'***** Reading soil parameter from file ',sitefile(ip),'...'
      write (unit_err, *) 'Soil parameter from file   ',trim(sitefile(ip))
  endif

  do
     read(inunit,*) text
     IF(text .ne. '!')then
        backspace(inunit)
        exit
     endif
  enddo
     
  if (flag_multi .eq. 8.or. flag_multi.eq.5.or. flag_mult910) then  
      read(inunit,*) text
     IF((text .eq. 'N') .or. (text .eq. 'n'))then
         flag_soilin = 3
     else
         flag_soilin = 2
         backspace(inunit)
     endif
  else
     read(inunit,*) text
     IF((text .eq. 'N') .or. (text .eq. 'n'))then
         flag_soilin = 1
     else
         flag_soilin = 0
         backspace(inunit)
     endif
     soilid(ip) = valfile(ip)
  endif
  if ((text .eq. 'S') .or. (text .eq. 's'))then
      flag_soilin = 4
      read(inunit,*) text
  endif
  if (.not.flag_mult8910) then
      write (unit_err, *) 'Soil identity number       ', trim(soilid(ip))
      write (unit_err, *) 'Climate ID                 ', trim(clim_id(ip))
  endif
  
  if (flag_soilin .eq. 1 .or. flag_soilin .ge. 3) then
    flag_hum = 1
  endif  
  
  if (flag_cond .ge. 40) then
    flag_hum = 0
  endif
  
  select case (flag_soilin)

  case (0,1)   ! single files f. j. site

     read (inunit,*,iostat=ios) long
     read (inunit,*,iostat=ios) lat
     read (inunit,*,iostat=ios) nlay
     read (inunit,*,iostat=ios) nroot_max
     read (inunit,*,iostat=ios) helpgrw
     
     if (helpgrw .gt. 1) then
        grwlev = helpgrw
     else
        fakt   = helpgrw
        grwlev = 1000.
     endif
     
     read (inunit,*,iostat=ios) w_ev_d
     read(inunit,*,iostat=ios) k_hum     ! mineralization constants of humus
     read(inunit,*,iostat=ios) k_hum_r
     read(inunit,*,iostat=ios) k_nit     ! nitrification constant

     IF(help==0) call alloc_soil
         read (inunit,*,iostat=ios) text
         select case (flag_soilin)
         case (0)      ! old input structure
             do i = 1, nlay
               read (inunit,*,iostat=ios) text
               read (inunit,*,iostat=ios) thick(i),pv_v(i),dens(i),f_cap_v(i), &
                     wilt_p_v(i),spheat(i),phv(i),wlam(i)
             end do
             skelv = 0.
     
         case(1)      ! new input structure
             do i = 1, nlay
               read (inunit,*,iostat=ios) helpnr, thick(i),pv_v(i),f_cap_v(i),wilt_p_v(i), &
                     dens(i),spheat(i),phv(i),wlam(i),skelv(i), sandv(i),clayv(i),humusv(i),&
                     C_hum(i), N_hum(i),NH4(i),NO3(i)
                     if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
                       if (phv(i) .le. 0.01) phv(i)=6.0         ! if flag_wurz 4 or 6 is used for calculation a pH-value is assumed
                     endif
             end do
         end select  ! flag_soilin (0,1)
     
    if (.not.flag_mult8910) print *, ' '
     IF (ios .ne.0) then
          print *,' >>>FORESEE message: Error during reading soil data!'
          WRITE(*,'(A)',advance='no') '  Stop program (y/n)? '
          read *, a
          IF ( a .eq. 'y' .or. a .eq. 'Y') then
              print *, '  STOP program!'
              stop
          endif
          IF (help==1)  call dealloc_soil
          print *,'  Check your input choice!!!'
      endif   ! ios

  case (2)   ! all sites are read from one file; old structure
     
     ios = 0
     do while (ios .eq. 0)
         read (inunit,*,iostat=ios) helpid, helpnl, helpnr
         if (trim(soilid(ip)) .ne. trim(helpid)) then
            do i = 1, helpnl
               read (inunit,*,iostat=ios) helpid
            enddo
         else
            nlay      = helpnl
            nroot_max = helpnr
            if (help==0) call alloc_soil
            do i = 1, nlay
               read (inunit,*,iostat=ios) helpnl, hor, boart, depth(i), thick(i),pv_v(i),dens(i), &
                    f_cap_v(i), wilt_p_v(i), spheat(i),phv(i),wlam(i), &
                    C_hum(i), N_hum(i), NH4(i), NO3(i), temps(i)
            enddo
            lat = latitude(ip)
			grwlev = gwtable(ip)
            exit
         endif
     enddo

     IF (ios .lt. 0) then
       if (.not.flag_mult8910) print *,' >>>FORESEE message: soil_id ', soilid(ip), ' not found'
       if (.not.flag_mult8910) print *,'  Check your input choice!!!'
       if (help==1)  call dealloc_soil
       CALL error_mess(time,"soil identificator not found "//adjustl(soilid(ip))//" ip No. ",real(help_ip))
       flag_end = 5
       return
     ENDIF   ! ios

     skelv = 0.

  case (3)   ! all sites are read from one file; new structure
     
     ios = 0
     do while (ios .eq. 0)
         read (inunit,*,iostat=ios) helpid, helpnl, helpnr
         if (trim(soilid(ip)) .ne. trim(helpid)) then
            do i = 1, helpnl
               read (inunit,*,iostat=ios) helpid
            enddo
         else
            nlay      = helpnl
            nroot_max = helpnr
            if (help==0) call alloc_soil
             do i = 1, nlay
               read (inunit,*,iostat=ios) helpnr, hor, boart, depth(i), thick(i),pv_v(i),f_cap_v(i), &
                     wilt_p_v(i),dens(i),spheat(i),phv(i),wlam(i),skelv(i), sandv(i), &
                     clayv(i),humusv(i),C_hum(i), N_hum(i),NH4(i),NO3(i)
                     if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
                       if (phv(i) .le. 0.01) phv(i)=6.0         ! if flag_wurz 4 or 6 is used for calculation a pH-value is assumed
                     endif
             end do
            lat = latitude(ip)
			grwlev = gwtable(ip)
            exit
         endif
     enddo
     IF (ios .lt. 0) then
       if (.not.flag_mult8910) print *,' >>>FORESEE message: soil_id ', soilid(ip), ' not found'
       if (.not.flag_mult8910) print *,'  Check your input choice!!!'
       if (help==1)  call dealloc_soil
       CALL error_mess(time,"soil identificator not found "//adjustl(soilid(ip))//"ip No.",real(help_ip))
       flag_end = 5
       return
     ENDIF   ! ios

  case (4)   ! one file several sites

     if (.not.flag_mult8910) print *,'   Reading soil model parameter from soil type file... ', soilid(ip)

     ios = 0
     do while (ios .eq. 0)
         read (inunit,*,iostat=ios) helpid
         if (trim(soilid(ip)) .ne. trim(helpid)) then
            read (inunit,*,iostat=ios) text
            read (inunit,*,iostat=ios) text
            read (inunit,*,iostat=ios) helpnl
            do i = 1, helpnl+6
               read (inunit,*,iostat=ios) boart
            enddo
            read (inunit,*,iostat=ios) boart
         else
             read (inunit,*,iostat=ios) hlong
             read (inunit,*,iostat=ios) hlat
             read (inunit,*,iostat=ios) nlay
             read (inunit,*,iostat=ios) nroot_max
             read (inunit,*,iostat=ios) helpgrw
             if (flag_multi .eq. 8.or. flag_multi.eq.5.or. flag_mult910) then
                if (abs(latitude(ip)) .gt. 90.) lat = latitude(ip)
			    grwlev = gwtable(ip)
             else
                 if (helpgrw .gt. 1) then
                    grwlev = helpgrw
                 else
                    fakt   = helpgrw
                    grwlev = 1000.
                 endif
                long = hlong
                lat = hlat
             endif  
             read (inunit,*,iostat=ios) w_ev_d
             read(inunit,*,iostat=ios) k_hum     ! mineralization constants of humus
             read(inunit,*,iostat=ios) k_hum_r
             read(inunit,*,iostat=ios) k_nit     ! nitrification constant

             IF(help==0) call alloc_soil
 
             read (inunit,*,iostat=ios) text
              do i = 1, nlay
               read (inunit,*,iostat=ios) helpnr, thick(i),pv_v(i),f_cap_v(i),wilt_p_v(i), &
                     dens(i),spheat(i),phv(i),wlam(i),skelv(i), sandv(i),clayv(i),humusv(i),&
                     C_hum(i), N_hum(i),NH4(i),NO3(i)
                     if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
                       if (phv(i) .le. 0.01) phv(i)=6.0         ! if flag_wurz 4 or 6 is used for calculation a pH-value is assumed
                     endif
             end do
             IF (ios .ne.0) then
                  print *,' >>>FORESEE message: Error during reading soil data!'    
                  print *, ' Program stopped!'
                  IF (help==1)  call dealloc_soil
                  flag_end = 7
                  return
              endif   ! ios
           exit
         endif
     enddo

     if (.not.flag_mult8910) print *, ' '
     IF (ios .lt. 0) then
       if (.not.flag_mult8910) then
           print *,' >>>FORESEE message: soil_id ', soilid(ip), ' not found'
           print *,'  Check your input choice!!!'
       endif
       if (help==1)  call dealloc_soil
       CALL error_mess(time,"soil identificator not found "//adjustl(soilid(ip))//"ip No.",real(help_ip))
       flag_end = 5
       return
     ENDIF   ! ios

  end select  ! flag_soilin
  close(inunit)
endif    ! ex

if (nroot_max .lt. 0) then
    do i=1, nlay
        if (C_hum(i) .gt. zero) nroot_max = i
    enddo
endif
if (.not.flag_mult8910) then
    write (unit_err, *) 'Latitude                ',lat
    write (unit_err,*)
endif

end subroutine readsoil

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

subroutine readvalue ! Input of cn-parameters and start values for soil model

integer :: inunit   
character :: text

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' readvalue'

call testfile(valfile(ip),ex)
IF (ex .eqv. .true.) then
   ios = 0
   inunit = getunit()
   open(inunit,file=valfile(ip),status='old',action='read')
   if (.not.flag_mult8910) print *,' *** Reading initial soil values from file ',valfile(ip),'...'
   do
      read(inunit,*) text
      IF(text .ne. '!')then
         backspace(inunit);exit
      endif
   enddo
 ! Soil temperature
   read(inunit,*,iostat=ios) text
   read(inunit,*,iostat=ios) (temps(i),i=1,nlay)
 ! C-content of humus
    read(inunit,*,iostat=ios) text
    read(inunit,*,iostat=ios) (C_hum(i),i=1,nlay)
 ! N-content of humus
    read(inunit,*,iostat=ios) text
    read(inunit,*,iostat=ios) (N_hum(i),i=1,nlay)
 ! NH4-content
    read(inunit,*,iostat=ios) text
    read(inunit,*,iostat=ios) (NH4(i),i=1,nlay)
 ! NO3-content
    read(inunit,*,iostat=ios) text
    read(inunit,*,iostat=ios) (NO3(i),i=1,nlay)
endif

IF (ios .ne. 0) then
print *,' >>>FORESEE message: Error during reading start values!'
WRITE(*,'(A)',advance='no') '  Stop program (y/n)? '
read *, a
  IF ( a .eq. 'y' .or. a .eq. 'Y') then
  print *, '  STOP program!'
  stop
  ELSE
  call dealloc_soil
  print *,'  Check your input choice!!!'
  end if
endif
close(inunit)

end subroutine readvalue

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

subroutine alloc_soil
use data_soil_t
use data_soil

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' alloc_soil'

help=0
allocate(thick(nlay))
allocate(mid(nlay))
allocate(depth(nlay))
allocate(pv(nlay))
allocate(pv_v(nlay))
allocate(dens(nlay))
allocate(f_cap_v(nlay))
allocate(field_cap(nlay))
allocate(wilt_p(nlay))
allocate(wilt_p_v(nlay))
allocate(vol(nlay))
allocate(quarzv(nlay))
allocate(sandv(nlay))
allocate(BDopt(nlay))
allocate(clayv(nlay))
allocate(siltv(nlay))
allocate(humusv(nlay))
allocate(fcaph(nlay))
allocate(wiltph(nlay))
allocate(pvh(nlay))
allocate(dmass(nlay))
allocate(skelv(nlay))
allocate(skelfact(nlay))
allocate(spheat(nlay))
allocate(phv(nlay))
allocate(wlam(nlay))
allocate(wats(nlay))
allocate(watvol(nlay))
allocate(wat_res(nlay))
wat_res = 0.

allocate(perc(nlay))
allocate(wupt_r(nlay))
allocate(wupt_ev(nlay))
allocate(s_drought(nlay))
allocate(root_fr(nlay))
!allocate(dp_rfr(nlay))
allocate(temps(nlay))
allocate (C_opm(nlay))
allocate (C_hum(nlay))
allocate (C_opmfrt(nlay))
allocate (C_opmcrt(nlay))
allocate (N_opm(nlay))
allocate (N_hum(nlay))
allocate (N_opmfrt(nlay))
allocate (N_opmcrt(nlay))
allocate (NH4(nlay))
allocate (NO3(nlay))
allocate (Nupt(nlay))
allocate (Nmin(nlay))
allocate (rmin_phv(nlay))
allocate (rnit_phv(nlay))
allocate (cnv_opm(nlay))
allocate (cnv_hum(nlay))
allocate(slit(nspecies))
allocate(slit_1(nspecies))

if (flag_bc .gt. 0) then
    allocate (C_bc(nlay))
    allocate (N_bc(nlay))
    C_bc = 0.
    N_bc = 0.
endif

do i=1,nspecies
slit(i)%C_opm_frt = 0.
slit(i)%N_opm_frt = 0.
slit(i)%C_opm_crt = 0.
slit(i)%N_opm_crt = 0.
slit(i)%C_opm_tb = 0.
slit(i)%N_opm_tb = 0.
slit(i)%C_opm_stem = 0.
slit(i)%N_opm_stem = 0.
enddo

nlay2 = nlay+2
mfirst = 1

allocate (sh(mfirst:nlay2))
allocate (sv(mfirst:nlay2))
allocate (sb(mfirst:nlay2))
allocate (sbt(mfirst:nlay2))
allocate (t_cb(mfirst:nlay2))
allocate (t_cond(mfirst:nlay2))
allocate (h_cap(mfirst:nlay2))
allocate (sxx(mfirst:nlay2))
allocate (svv(mfirst:nlay2))
allocate (svva(mfirst:nlay2))
allocate (soh(mfirst:nlay2))
allocate (son(mfirst:nlay2+1))
help=1
C_opm = 0
allocate(fr_loss(nlay))
allocate(redis(nlay))

end subroutine alloc_soil

!------------------------------------------------------------------
subroutine stand_bal_ini

use data_stand

implicit none

integer i   

allocate(diam_class(num_class, nspecies)); diam_class=0
allocate(diam_class_t(num_class, nspecies)); diam_class_t=0
allocate(diam_class_h(num_class,nspecies)); diam_class_h=0
allocate(diam_class_age(num_class,nspecies)); diam_class_age=0
allocate(diam_class_mvol(num_class,nspecies)); diam_class_mvol=0
allocate(diam_classm(num_class,nspecies)); diam_classm=0
allocate(diam_classm_h(num_class,nspecies)); diam_classm_h=0
allocate(height_class(num_class)); height_class =0

! array of potential litter (dead stems and twigs/branches for the next years 
allocate(dead_wood(nspec_tree))
do i = 1,nspec_tree
    allocate(dead_wood(i)%C_tb(lit_year))
    allocate(dead_wood(i)%N_tb(lit_year))
    allocate(dead_wood(i)%C_stem(lit_year))
    allocate(dead_wood(i)%N_stem(lit_year))
    dead_wood(i)%C_tb = 0.
    dead_wood(i)%N_tb = 0.
    dead_wood(i)%C_stem = 0.
    dead_wood(i)%N_stem = 0.
enddo

end subroutine stand_bal_ini

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

subroutine control_file ! saving simulation parameter and start conditions for each site
real buckdepth  
character(8) actdate
character(10) acttime
character(150) site_help
integer help_ip,  j
TYPE(Coh_Obj), Pointer :: help_coh  ! pointer to cohort list

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

! Write soil initialisation file
if (flag_mult8910) then
    site_help = site_name1
else
    site_help = site_name(help_ip)
endif

if (.not.flag_mult8910 .or. (flag_mult8910 .and. anh .eq. "1") .or. (flag_mult8910 .and. time_out .gt. 0.)) then
    if (.not.flag_mult8910) then
        WRITE (unit_soil,'(26A)') 'Layer',' Depth(cm)',' F-cap(mm)',' F-cap(Vol%)','   Wiltp(mm)', &
          ' Wiltp(Vol%)',' Pore vol.',' Skel.(Vol%)',' Density','  Spheat','      pH','    Wlam',    &
          ' Water(mm)',' Water(Vol%)',' Soil-temp.',' C_opm g/m2', &
          ' C_hum g/m2',' N_opm g/m2',' N_hum g/m2',' NH4 g/m2',' NO3 g/m2','  humus part',' d_mass g/m2', '  Clay','  Silt','  Sand'
        do i = 1,nlay
        WRITE (unit_soil,'(I5,2F10.2,3F12.2,F10.2,F12.2,4F8.2,F10.2,F12.2, 5F11.2,2F9.4,2E12.4, 3F6.1)') i,depth(i),field_cap(i),f_cap_v(i),wilt_p(i), &
              wilt_p_v(i),pv_v(i), skelv(i)*100., dens(i),spheat(i),phv(i),wlam(i),   &
              wats(i),watvol(i),temps(i),c_opm(i),c_hum(i),n_opm(i), n_hum(i),nh4(i),no3(i),humusv(i),dmass(i), clayv(i)*100., siltv(i)*100., sandv(i)*100.

        end do
    endif

    ! Write control file
    call date_and_time(actdate, acttime)
    unit_ctr = getunit()
    open(unit_ctr,file=trim(dirout)//trim(site_help)//'.ctr'//anh,status='replace')
    WRITE(unit_ctr,'(2A)')     '*** Site name: ',site_name(help_ip)
    WRITE(unit_ctr,'(2A)')     '    Appendix ' ,anh
    WRITE(unit_ctr,'(A,F7.2)') '  Longitude: ', long
    WRITE(unit_ctr,'(A,F7.2)') '  Latitude:  ', lat
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(10A)') '               ----   Version: v2.2   ---- '  
    WRITE(unit_ctr,'(10A)') '  Date: ',actdate(7:8),'.',actdate(5:6),'.',actdate(1:4), &
                            '  Time: ',acttime(1:2),':',acttime(3:4)
    WRITE(unit_ctr,'(A,A)') '   Simulation control file:      ',trim(simfile)
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Data files:'
    IF(flag_clim==1)then
    WRITE(unit_ctr,'(A,A)') ' Climfile:                   ',trim(climfile(ip))
    ELSE
    WRITE(unit_ctr,'(A,A)') ' Climfile:                   ',trim(climfile(1))
    endif
    WRITE(unit_ctr,'(A,A)') ' Sitefile:                   ',trim(sitefile(help_ip))
    WRITE(unit_ctr,'(A,A)') ' Start value file:           ',trim(valfile(help_ip))

    ! Initialization of stand
    IF( flag_multi==3 .OR. (site_nr>1 .AND. flag_stand>0) ) THEN
    WRITE(unit_ctr,'(A,A)') ' Stand initialization:       ',trim(treefile(ip))
    ELSE IF( ip==1 .AND. flag_stand>0) THEN
    WRITE(unit_ctr,'(A,A)') ' Stand initialization:       ',trim(treefile(ip))
    ELSE IF (flag_stand==0) THEN
    WRITE(unit_ctr,'(A,A)') ' Stand initialization:       none'
    endif
    IF (lmulti) WRITE(unit_ctr,'(A,A)') ' Stand identificator: ', adjustl(standid(ip))
    WRITE(unit_ctr,*) ' '
    IF(flag_mg.ne.0 .and. flag_mg.ne.5) then
    WRITE(unit_ctr,'(A,A)') ' Management control file:   ',trim(manfile(ip))
    ELSE
    WRITE(unit_ctr,'(A)') ' Management:                 none'
    endif
    WRITE(unit_ctr,'(A,A)') ' Deposition file:            ',trim(depofile(ip))
    WRITE(unit_ctr,'(A,A)') ' N reduction file:           ',trim(redfile(ip))
    WRITE(unit_ctr,'(A,A)') ' Litter initialisation file: ',trim(litfile(ip))
    if (flag_stat .gt. 0) WRITE(unit_ctr,'(A,A)') ' File with measurements:     ',trim(mesfile(1))
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Soil description '
    WRITE(unit_ctr,'(A,I3)') ' Number of soil layers:    ',nlay
    WRITE(unit_ctr,'(A,I3)') ' Number of rooting layers: ',nroot_max
    WRITE(unit_ctr,'(A,I3)') ' Ground water from layer:  ',nlgrw
    WRITE(unit_ctr,'(A,F5.1)') ' Evaporation depth (cm):   ',w_ev_d
    call bucket(bucks_100, bucks_root, buckdepth)
    buckdepth = buckdepth/100
    WRITE(unit_ctr,'(A,F5.2,A,F7.2)') ' Bucket size (mm), ', buckdepth,' m depth: ',bucks_100
    WRITE(unit_ctr,'(A,F7.2)') ' Bucket size (mm) of rooting zone: ',bucks_root
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Soil water conditions'
    WRITE(unit_ctr,'(12A)') 'Layer ','Depth(cm) ','F-cap(mm) ','F-cap(Vol%) ','Wiltp(mm) ', &
      'Wiltp(Vol%) ','Pore vol. ','Density   ','Spheat  ','pH-value  ','   Wlam','   skel. '
    do i = 1,nlay
    WRITE(unit_ctr,'(I5,12F10.2)') i,depth(i),field_cap(i),f_cap_v(i),wilt_p(i), &
       wilt_p_v(i),pv_v(i),dens(i),spheat(i),phv(i),wlam(i),skelv(i)
    end do
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Soil initial values'
    WRITE(unit_ctr,'(9A)') 'Layer ','Water-cont. ','Soil-temp. ','C_opm     ', &
      'C_hum     ','N_opm     ','N_hum     ','NH4-cont. ','NO3-cont '
    do i=1,nlay
    WRITE(unit_ctr,'(I5, 2F10.2, 6F10.4)') i,wats(i),temps(i),c_opm(i),c_hum(i),n_opm(i), &
      n_hum(i),nh4(i),no3(i)
    end do
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '       N_tot       C_tot     N_antot    N_humtot    C_humtot   C_opm_fol   C_opm_tb   C_opm_frt   C_opm_crt   C_opm_stem '
    WRITE(unit_ctr,'(10F12.4)') N_tot, C_tot, N_an_tot, N_hum_tot, C_hum_tot, C_opm_fol, C_opm_tb, C_opm_frt, C_opm_crt, C_opm_stem
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)',advance='no') 'Mineralization constant of humus - humus layer (k_hum): '
    WRITE(unit_ctr,'(F10.5)') k_hum
    WRITE(unit_ctr,'(A)',advance='no') 'Mineralization constant of humus - mineral soil (k_hum_r): '
    WRITE(unit_ctr,'(F10.5)') k_hum_r
    WRITE(unit_ctr,'(A)',advance='no') 'Nitrification constant (k_nit): '
    WRITE(unit_ctr,'(F10.5)') k_nit
    WRITE(unit_ctr,*) ' '
    if (flag_bc .gt.0) then
    WRITE(unit_ctr,'(A)') '*** Biochar application '
    WRITE(unit_ctr,'(A)') '    year  C-content(%)  C/N-ratio  depth  mass(kg/ha dry mass)'
      do j = 1, n_appl_bc
      WRITE(unit_ctr,'(I7,F14.1, F11.1, I7, F18.1)') &
            y_bc(j), cpart_bc(j), cnv_bc(j), bc_appl_lay(j), C_bc_appl(j)      
      enddo
    WRITE(unit_ctr,'(F10.5)')  
    endif
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Stand initialisation'
    WRITE(unit_ctr,'(A)')'  Coh      x_fol       x_frt       x_sap       x_hrt       x_Ahb     height   x_hbole  x_age      n     sp       DC          DBH'
    help_coh => pt%first
    DO WHILE (ASSOCIATED(help_coh))
    WRITE(unit_ctr,'(I5,5f12.5,2f10.0,i7,f7.0,i7, 2f12.5)') help_coh%coh%ident, help_coh%coh%x_fol, help_coh%coh%x_frt, help_coh%coh%x_sap, help_coh%coh%x_hrt, &
              help_coh%coh%x_Ahb, help_coh%coh%height, help_coh%coh%x_hbole, help_coh%coh%x_age, &
              help_coh%coh%nTreeA,help_coh%coh%species, help_coh%coh%dcrb, help_coh%coh%diam
    help_coh => help_coh%next
    END DO
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Simulation control'
    WRITE(unit_ctr,'(A66,I4)') 'Run option:                                                       ',flag_multi
    WRITE(unit_ctr,'(A66,I4)') 'Start year:                                                       ',time_b
    WRITE(unit_ctr,'(A66,I4)') 'Number of simulation years - year:                                ', year
    WRITE(unit_ctr,'(A60,F12.1)') 'Patch size [m²] - kpatchsize:                                     ',kpatchsize
    WRITE(unit_ctr,'(A60,F12.1)') 'Thickness of leaf layers - dz:                                    ',dz
    WRITE(unit_ctr,'(A66,I4)') 'Time step for photosynthesis calculations (days) -    ns_pro:     ',ns_pro
    WRITE(unit_ctr,'(A66,I4)') 'Mortality (0-OFF,1-ON stress, 2- ON stress+intr) -  flag_mort:    ',flag_mort
    WRITE(unit_ctr,'(A66,I4)') 'Regeneration (0-OFF,1-ON, 2-weekly growth of seedl.) - flag_reg:  ',flag_reg
    WRITE(unit_ctr,'(A66,I4)') 'use FORSKA for regeneration (0-OFF,1-ON) - flag_forska:           ',flag_lambda
    WRITE(unit_ctr,'(A66,I4)') 'Stand initialization (0-no,1-from *.ini,2-generate) - flag_stand: ',flag_stand
    WRITE(unit_ctr,'(A66,I4)') 'Ground vegetation initialization (0-no,1-generate) - flag_sveg:     ',flag_sveg
    WRITE(unit_ctr,'(A66,I4)') 'Stand management (0-no,1-yes, 2 - seed once) - flag_mg:           ',flag_mg
    WRITE(unit_ctr,'(A66,I4)') 'Disturbance (0-OFF, 1-ON ) - flag_dis:                            ',flag_dis
    WRITE(unit_ctr,'(A66,I4)') 'Light absoption algorithm (1,2,3,4) - :                           ',flag_light
    WRITE(unit_ctr,'(A66,I4)') 'Foliage-height relationship (0,1) - flag_folhei:                  ',flag_folhei
    WRITE(unit_ctr,'(A66,I4)') 'Volume function trunc (0,1) - flag_volfunc:                       ',flag_volfunc
    WRITE(unit_ctr,'(A66,I4)') 'Respiration model (0-0.5*NPP,1-organ specific) - flag_resp:       ',flag_resp
    WRITE(unit_ctr,'(A66,I4)') 'Limitation (0-NO,1-water, 2-N, 3-water+N) - flag_limi:            ',flag_limi
    WRITE(unit_ctr,'(A66,I4)') 'Flag for decomposition model - flag_decomp:                       ',flag_decomp
    WRITE(unit_ctr,'(A66,I4)') 'Root spec. activity (0-const,1-varying) - flag_sign:              ',flag_sign
    WRITE(unit_ctr,'(A66,I4)') 'Water uptake function soil (1,2,3,4) - flag_wred:                 ',flag_wred
    WRITE(unit_ctr,'(A66,I4)') 'Root distribution - flag_wurz:                                    ',flag_wurz
    WRITE(unit_ctr,'(A66,I4)') 'Heat conductance - flag_cond:                                     ',flag_cond
    WRITE(unit_ctr,'(A66,I4)') 'Interception - flag_int:                                          ',flag_int
    WRITE(unit_ctr,'(A66,I4)') 'Evapotranspiration - flag_eva:                                    ',flag_eva
    WRITE(unit_ctr,'(A66,I4)') 'CO2 (0-constant,1-historic increase,2-step change)- flag_co2:     ',flag_co2
    WRITE(unit_ctr,'(A66,I4)') 'Sort flag - flag_sort:                                            ',flag_sort
    WRITE(unit_ctr,'(A66,I4)') 'wpm flag - flag_wpm:                                              ',flag_wpm
    WRITE(unit_ctr,'(A66,I4)') 'Analysis of measurements - flag_stat:                             ',flag_stat
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A66,A)')  'Species parameter file:                                           ',trim(specfile(help_ip))
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '*** Species parameter description'
    WRITE(unit_ctr,'(A51,I4)')   ' Species number: ', nspecies
    WRITE(unit_ctr,'(A51,I4)')   ' Number of tree species: ', nspec_tree
    WRITE(unit_ctr,*) '              ********** '
    WRITE(unit_ctr,'(A25,A9,2X,A30)') 'Short Name', '  Spec-Nr', 'Latin Name                   '
    WRITE(unit_ctr,*) ' '
    do i=1,nspecies
    WRITE(unit_ctr,'(A25,I9,2X,A30)') trim(spar(i)%species_short_name), i, spar(i)%species_name
    enddo
    WRITE(unit_ctr,*) '              ********** '
    WRITE(unit_ctr,'(A51,15A16)')   ' Species name: ', (trim(spar(i)%species_short_name),i=1,nspecies)
    WRITE(unit_ctr,1010) ' Maximal age - max_age: ', (spar(i)%max_age,i=1,nspecies)
    WRITE(unit_ctr,1010) ' Stress rec. time - yrec: ', (spar(i)%yrec,i=1,nspecies)
    WRITE(unit_ctr,1010) ' Shade tolerance - stol: ', (spar(i)%stol,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Extinction coeff - pfext: ', (spar(i)%pfext,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Root activity rate - sigman: ', (spar(i)%sigman,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Respiration coeff - respcoeff: ', (spar(i)%respcoeff,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Growth resp. par. - prg: ', (spar(i)%prg,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Maint.resp.par./sapwood - prms: ', (spar(i)%prms,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Maint.resp.par./fineroot - prmr: ', (spar(i)%prmr,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Senesc.par. foliage - psf: ', (spar(i)%psf,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Senesc.par. sapwood - pss: ', (spar(i)%pss,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Senesc.par. fineroot - psr: ', (spar(i)%psr,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N/C ratio of biomass - pcnr: ', (spar(i)%pcnr,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N concentration of foliage - ncon_fol: ', (spar(i)%ncon_fol,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N concentration of fine roots - ncon_frt: ', (spar(i)%ncon_frt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N concentration of coarse roots - ncon_crt: ', (spar(i)%ncon_crt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N concentration of twigs and branches - ncon_tbc: ', (spar(i)%ncon_tbc,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N concentration of stemwood - ncon_stem: ', (spar(i)%ncon_stem,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Reallocation parameter of foliage - reallo_fol: ', (spar(i)%reallo_fol,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Reallocation parameter of fine root - reallo_frt: ', (spar(i)%reallo_frt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Ratio of coarse wood - alphac: ', (spar(i)%alphac,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Coarse root fraction of coarse wood - cr_frac: ', (spar(i)%cr_frac,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Sapwood density - prhos: ', (spar(i)%prhos,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Proport.const.(pipe mod.) - pnus: ', (spar(i)%pnus,i=1,nspecies)
    IF(flag_folhei==0) THEN
    WRITE(unit_ctr,1000) ' Height growth parameter - pha: ', (spar(i)%pha,i=1,nspecies)
    ELSEIF(flag_folhei==1) THEN
    WRITE(unit_ctr,1000) ' Height growth par. 1 - pha_v1: ', (spar(i)%pha_v1,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Height growth par. 2 - pha_v2: ', (spar(i)%pha_v2,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Height growth par. 3 - pha_v3: ', (spar(i)%pha_v3,i=1,nspecies)
    ELSE
    WRITE(unit_ctr,'(A51,I3)')     ' non valid flag value - flag_folhei : ',flag_folhei
    ENDIF
    WRITE(unit_ctr,1000) ' Height growth parameter coeff 1 - pha_coeff1: ', (spar(i)%pha_coeff1,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Height growth parameter coeff 2 - pha_coeff2: ', (spar(i)%pha_coeff2,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Crown radius - DBH ratio parameter a - crown_a: ', (spar(i)%crown_a,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Crown radius - DBH ratio parameter b - crown_b: ', (spar(i)%crown_b,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Crown radius - DBH ratio parameter c - crown_c: ', (spar(i)%crown_c,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Minimum specific leaf area - psla_min: ', (spar(i)%psla_min,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Light dep. specific leaf area - psla_a: ', (spar(i)%psla_a,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Efficiency parameter - phic: ', (spar(i)%phic,i=1,nspecies)
    WRITE(unit_ctr,1000) ' N content - pnc: ', (spar(i)%pnc,i=1,nspecies)
    WRITE(unit_ctr,1000) ' kco2_25: ', (spar(i)%kCO2_25,i=1,nspecies)
    WRITE(unit_ctr,1000) ' ko2_25: ', (spar(i)%kO2_25,i=1,nspecies)
    WRITE(unit_ctr,1000) ' CO2/O2 specif. value - pc_25: ', (spar(i)%pc_25,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Q10_kco2: ', (spar(i)%q10_kCO2,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Q10_ko2: ', (spar(i)%q10_kO2,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Q10_pc: ', (spar(i)%q10_pc,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Rd to Vm ratio - pb: ', (spar(i)%pb,i=1,nspecies)

    WRITE(unit_ctr,1000) ' PIM: Inhibitor min temp. - PItmin: ', (spar(i)%PItmin,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Inhibitor opt temp. - PItopt: ', (spar(i)%PItopt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Inhibitor max temp. - PItmax: ', (spar(i)%PItmax,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Inhibitor scaling factor - PIa: ', (spar(i)%PIa,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Promotor min temp. - PPtmin: ', (spar(i)%PPtmin,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Promotor opt temp. - PPtopt: ', (spar(i)%PPtopt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Promotor max temp. - PPtmax: ', (spar(i)%PPtmax,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Promotor scaling factor - PPa: ', (spar(i)%PPa,i=1,nspecies)
    WRITE(unit_ctr,1000) ' PIM: Promotor scaling factor - PPb: ', (spar(i)%PPb,i=1,nspecies)
    WRITE(unit_ctr,1000) ' CSM: chilling base temp. - CSTbC: ', (spar(i)%CSTbC,i=1,nspecies)
    WRITE(unit_ctr,1000) ' CSM: base temp. - CSTbT: ', (spar(i)%CSTbT,i=1,nspecies)
    WRITE(unit_ctr,1000) ' CSM: scaling factor - CSa: ', (spar(i)%CSa,i=1,nspecies)
    WRITE(unit_ctr,1000) ' CSM: scaling factor - CSb: ', (spar(i)%CSb,i=1,nspecies)
    WRITE(unit_ctr,1000) ' TSM: base temp. - LTbT: ', (spar(i)%LTbT,i=1,nspecies)
    WRITE(unit_ctr,1000) ' TSM: critical temperature sum - LTcrit: ', (spar(i)%LTcrit,i=1,nspecies)
    WRITE(unit_ctr,1010) ' TSM: start day after 1.11. - Lstart: ', (spar(i)%Lstart,i=1,nspecies)
    WRITE(unit_ctr,1000) ' usefd pheno model - Phmodel: ', (spar(i)%Phmodel,i=1,nspecies)
    WRITE(unit_ctr,1000) ' End day for phenology - end_bb: ', (spar(i)%end_bb,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Fpar_mod - fpar_mod: ', (spar(i)%fpar_mod,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Intercep.cap. - ceppot_spec: ', (spar(i)%ceppot_spec,i=1,nspecies)
    WRITE(unit_ctr,1000) ' photosynthesis response to N-limitation - Nresp: ',  (spar(i)%Nresp,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Regeneration flag - regflag: ',  (spar(i)%regflag,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Seedrate: ',  (spar(i)%seedrate,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Seedmass: ',  (spar(i)%seedmass,i=1,nspecies)
    WRITE(unit_ctr,1000) ' Standard dev. of seedrate - seedsd: ',  (spar(i)%seedsd,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - seeda: ',  (spar(i)%seeda,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - seedb: ',  (spar(i)%seedb,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - pheight1: ',  (spar(i)%pheight1,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - pheight2: ',  (spar(i)%pheight2,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - pheight3: ',  (spar(i)%pheight3,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - pdiam1: ',  (spar(i)%pdiam1,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - pdiam2: ',  (spar(i)%pdiam2,i=1,nspecies)
    WRITE(unit_ctr,1000) ' all. parameter - pdiam3: ',  (spar(i)%pdiam3,i=1,nspecies)
    WRITE(unit_ctr,1000) ' decomp. parameter foliage - k_opm_fol: ',  (spar(i)%k_opm_fol,i=1,nspecies)
    WRITE(unit_ctr,1000) ' synth. parameter foliage - k_syn_fol: ',  (spar(i)%k_syn_fol,i=1,nspecies)
    WRITE(unit_ctr,1000) ' decomp. parameter fine roots - k_opm_frt: ',  (spar(i)%k_opm_frt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' synth. parameter fine roots - k_syn_frt: ',  (spar(i)%k_syn_frt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' decomp. parameter coarse roots - k_opm_crt: ',  (spar(i)%k_opm_crt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' synth. parameter coarse roots - k_syn_crt: ',  (spar(i)%k_syn_crt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' decomp. parameter twigs/branches - k_opm_tb: ',  (spar(i)%k_opm_tb,i=1,nspecies)
    WRITE(unit_ctr,1000) ' synth. parameter twigs/branches - k_syn_tb: ',  (spar(i)%k_syn_tb,i=1,nspecies)
    WRITE(unit_ctr,1000) ' decomp. parameter stem - k_opm_stem: ',  (spar(i)%k_opm_stem,i=1,nspecies)
    WRITE(unit_ctr,1000) ' synth. parameter dtem - k_syn_stem: ',  (spar(i)%k_syn_stem,i=1,nspecies)

    WRITE(unit_ctr,1000)
    WRITE(unit_ctr,1000) ' spec_rl: ',  (spar(i)%spec_rl,i=1,nspecies)
    WRITE(unit_ctr,1000) ' tbase: ',  (spar(i)%tbase,i=1,nspecies)
    WRITE(unit_ctr,1000) ' topt: ',  (spar(i)%topt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' bdmax_coef: ',  (spar(i)%bdmax_coef,i=1,nspecies)
    WRITE(unit_ctr,1000) ' porcrit_coef: ',  (spar(i)%porcrit_coef,i=1,nspecies)
    WRITE(unit_ctr,1000) ' ph_opt_max: ',  (spar(i)%ph_opt_max,i=1,nspecies)
    WRITE(unit_ctr,1000) ' ph_opt_min: ',  (spar(i)%ph_opt_min,i=1,nspecies)
    WRITE(unit_ctr,1000) ' ph_max: ',  (spar(i)%ph_max,i=1,nspecies)
    WRITE(unit_ctr,1000) ' ph_min : ',  (spar(i)%ph_min ,i=1,nspecies)
    WRITE(unit_ctr,1000) ' v_growth: ',  (spar(i)%v_growth,i=1,nspecies)

    WRITE(unit_ctr,1000)
    WRITE(unit_ctr,1000) ' C/N ratio of foliage - cnr_fol: ', (spar(i)%cnr_fol,i=1,nspecies)
    WRITE(unit_ctr,1000) ' C/N ratio of fine roots - cnr_frt: ', (spar(i)%cnr_frt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' C/N ratio of coarse roots - cnr_crt: ', (spar(i)%cnr_crt,i=1,nspecies)
    WRITE(unit_ctr,1000) ' C/N ratio of twigs and branches - cnr_tbc: ', (spar(i)%cnr_tbc,i=1,nspecies)
    WRITE(unit_ctr,1000) ' C/N ratio of stemwood - cnr_stem: ', (spar(i)%cnr_stem,i=1,nspecies)

    WRITE(unit_ctr,1000)
    WRITE(unit_ctr,1000) ' Reduction factor - RedN: ', (svar(i)%RedN, i=1,nspecies)

    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,'(A)') '******   Model parameter   ******'
    WRITE(unit_ctr,1020) 'Optimum ratio of ci to ca [-] - Lambda: ',lambda
    WRITE(unit_ctr,1020) 'Molar mass of carbon [g/mol] - Cmass: ',Cmass
    WRITE(unit_ctr,1020) 'Minimum conductance [mol/(m2*d)] - gmin: ',gmin
    WRITE(unit_ctr,1020) 'Shape of PS response curve - ps: ',ps
    WRITE(unit_ctr,1020) 'Slope of N function at 20 °C [g(N) (mymol s-1)-1] - pn: ',pn
    WRITE(unit_ctr,1020) 'Minimum N content [g/g] - nc0: ',nc0
    WRITE(unit_ctr,1020) 'C3 quantum efficiency - qco2: ',qco2
    WRITE(unit_ctr,1020) 'Scaling parameter - qco2a: ',qco2a
    WRITE(unit_ctr,1020) 'Partial pressure of oxygen (kPa) - o2: ',o2
    WRITE(unit_ctr,1020) 'Atmospheric CO2 content (mol/mol) - co2: ',co2_st
    WRITE(unit_ctr,1020) 'Albedo of the canopy - pfref: ',pfref
    WRITE(unit_ctr,1020) 'Part of C in biomass [-] - cpart: ',cpart
    WRITE(unit_ctr,1020) 'Ratio of molecular weights of water and air - rmolw: ',rmolw
    WRITE(unit_ctr,1020) 'Universal gas constant [J/mol/K] = [Pa/m3/K] - R_gas: ',R_gas
    WRITE(unit_ctr,1020) 'von Karman''s constant [-] - c_karman: ',c_karman
    WRITE(unit_ctr,1020) 'Specific heat of air at const. pressure [J/g/K] - c_air: ',c_air
    WRITE(unit_ctr,1020) 'Psychrometer constant [hPa/K] - psycro: ',psycro
    WRITE(unit_ctr,1020) 'Breast height for inventory measurements [cm] - h_breast: ',h_breast
    WRITE(unit_ctr,1020) 'Height for sapling allometry - h_sapini: ',h_sapini
    WRITE(unit_ctr,1020) 'Min. diff. b. height of crown base and breast height- h_bo_br_diff: ',h_bo_br_diff
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p1_co2: ',p1_co2
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p2_co2: ',p2_co2
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p3_co2: ',p3_co2
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p4_co2: ',p4_co2
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of CO2 scenario - p5_co2: ',p5_co2
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p1_co2h: ',p1_co2h
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p2_co2h: ',p2_co2h
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p3_co2h: ',p3_co2h
    WRITE(unit_ctr,1020) 'Parameter variable for calculation of historical CO2 scenario - p4_co2h: ',p4_co2h
    WRITE(unit_ctr,1020) 'Threshold of air temperature for snow accumulation [°C] - temp_snow:    ',temp_snow
    WRITE(unit_ctr,1020) 'Parameter for calculation of transpiration demand  - alfm:              ',alfm
    WRITE(unit_ctr,1020) 'Parameter for calculation of transpiration demand [mol/(m2*d)] - gpmax: ',gpmax 
    WRITE(unit_ctr,1020) 'Parameter for growing degree day calculation - thr_gdd: ',thr_gdd
    
    IF (flag_multi==2) THEN
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,*) 'runs with climate scenarios produced by adding summands to every daily temperature'
    WRITE(unit_ctr,*) 'and modifying every single precipitation value by a multiplicand'
    WRITE(unit_ctr,*) 'run ident      deltaT     delta P factor'
    ENDIF

    !   mangament parameter adaptation management

    IF (flag_mg.eq.2. .and. flag_reg .eq. 0) then
    WRITE(unit_ctr,*) ' '
    WRITE(unit_ctr,*) '***Managment parameter case flag_mg = 2 (user specified)  ***'
    WRITE(unit_ctr,'(A35,4F15.5)')  'height for management control(cm)', ho1,ho2,ho3,ho4
    WRITE(unit_ctr,'(A35,5I15)')  'management flags thr1-thr5' , thr1,thr2, thr3,thr4,thr5
    WRITE(unit_ctr,'(A35,F15.5)')   'height for directional felling', thr6
    WRITE(unit_ctr,'(A35,I15)')   'measure at rotation', thr7
    WRITE(unit_ctr,'(A35,I15)')   'regeneration measure', mgreg
    WRITE(unit_ctr,'(A35,F15.5)')  'lower/upper limit of height(cm)', limit
    WRITE(unit_ctr,'(A35,I15)')  'number of years between thinning',thinstep
    WRITE(unit_ctr,'(A35,F15.5)') 'rel. value for directional felling', direcfel
    WRITE(unit_ctr,'(A35,5F15.5)')'number of Zielbaeume(spec.)', (zbnr(i),i=1,nspec_tree)
    WRITE(unit_ctr,'(A35,5F15.5)')'rel. value for tending of pl.',(tend(i), i =1,nspec_tree)
    WRITE(unit_ctr,'(A35,5I15)')'rotation ',(rot(i), i =1,nspec_tree)
    WRITE(unit_ctr,'(A35,5I15)')'age of nat./pl. regeneration',(regage(i), i =1,nspec_tree)
    end IF

    IF (flag_multi .ne. 2.and. flag_mg.ne.2 .and. flag_reg .eq.0) close(unit_ctr)
endif  !  flag_mult8910

1000 FORMAT (A51,15 F16.5)
1010 FORMAT (A51,15 I16)
1020 FORMAT(A70,F15.5)

end subroutine control_file

end subroutine prepare_site

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

SUBROUTINE readbudb

use data_simul
use data_species
use data_stand

implicit none

   DO ns=1,nspecies
      IF(spar(ns)%phmodel==4) THEN
         WRITE(*,*) 'Please type the day of budburst for 4C species number ',ns,':'
         READ(*,*) svar(ns)%ext_daybb
      ENDIF
   ENDDO

END subroutine readbudb

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

SUBROUTINE readdepo

use data_climate
use data_depo
use data_out
use data_simul
use data_site

implicit none

character  text
integer  hx, unit_dep, i,j,ios, ii
!integer  realrec
integer  id,im,iy,itz1, itz2, hyear1, hyear2, hyear3, hy
logical  ex
real     hNO, hNH

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' readdepo'

if (.not.allocated(NOd)) allocate (NOd (1:366,1:year))
if (.not.allocated(NHd)) allocate (NHd (1:366, 1:year))

! for areal usage standard/constant deposition is set as concentration:
if (flag_multi .eq. 8 .or. flag_mult910) then
    flag_depo = 2
    NOd = NOdep(ip)       ! concentration mg/l
    NHd = NHdep(ip)       ! concentration mg/l
    return
endif

NOd = 0.
NHd = 0.

if (.not.flag_mult8910)  print *
inquire (File = depofile(ip), exist = ex)  ! test whether file exists
  IF(ex .eqv. .false.) then
    if (.not.flag_mult8910) then
        hx = 0
        print *,' >>>FORESEE message: Cannot find deposition data - all data set to zero!'
        CALL error_mess(hx,'Cannot find deposition data - all data set to ',REAL(hx))
    endif
  else
    if (.not.flag_mult8910) print *, ' >>>FORESEE message: Now reading DEPOSITION data from file, please wait...'
!   now read data from file
    unit_dep = getunit()
    OPEN (unit_dep,FILE=depofile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')

    flag_depo = 1
    read(unit_dep,*) text
    select case (text)
    case ('C', 'c')     ! concentrations mg/l
        flag_depo = 2
        read(unit_dep,*) text

    case ('Y', 'y')   ! Yearly constant deposition mg/m2
        flag_depo = 3
        read(unit_dep,*) text

    case ('A', 'a')   ! Annual sum of deposition g/m2
        flag_depo = 4
        read(unit_dep,*) text
    
    end select

    do
      IF (text .ne. '!') then
         backspace(unit_dep)
     exit
      endif
      read(unit_dep,*) text
    enddo

!    assignment of dates
!    fill in missing values with current values until current date
!    fill in missing values at the end
        hyear1 = 0
        hyear2 = 0
        hyear3 = 1
        itz1 = 1
        itz2 = 1
        select case (flag_depo)
        case(4)
            do while ((ios .eq. 0) .and. (hyear1 .lt. year))
                read(unit_dep,*,iostat=ios) iy, hNO, hNH
                if (ios .eq.0) then
                    if (iy .gt. time_b+year) then
                        hyear1 = year
                    else
                        hyear1 = iy - time_b + 1
                    endif
                    if ((hyear1 .le. year) .and. (hyear1 .gt. 0)) then        ! save from simulation start year onwards
                        do i = 1,366
                            NOd(i,hyear1) = hNO * 1000./366.                  ! report of year [g/m2] in daily values [mg/m2] 
                            NHd(i,hyear1) = hNH * 1000./366.
                        enddo
                        hy = hyear1-1
                        do while ((hy .gt. hyear2) .and. (hy .gt. 0))  
                            do i = 366, 1, -1
                                NOd(i,hy) = hNO * 1000./366.
                                NHd(i,hy) = hNH * 1000./366.
                            enddo
                            hy = hy - 1
                        enddo
                        hyear2 = hyear1
                     endif   ! 0 < hyear1 < year
               else  ! ios .ne. 0
                    if (hyear1 .le. 0) then
                        hyear1 = 1
                        hyear2 = 1
                    endif
                    continue
                endif  ! ios = 0
            enddo
                  
        case default
            do while ((ios .eq. 0) .and. (hyear1 .lt. year))
                read(unit_dep,*,iostat=ios) id,im,iy, hNO, hNH
                if (ios .eq.0) then
                    call daintz(id,im,iy,itz1) 
                    if (iy .gt. time_b+year) then
                        hyear1 = year
                    else
                        hyear1 = iy - time_b + 1
                    endif
                    if ((hyear1 .le. year) .and. (hyear1 .gt. 0)) then        ! save from simulation start year onwards 
                        NOd(itz1,hyear1) = hNO 
                        NHd(itz1,hyear1) = hNH 

                        select case (flag_depo)
                        case(1,2)                     
                            if (hyear1 .eq. hyear3) then
                                if (itz1 .gt. 1) then
                                    do i = itz1-1, itz2, -1
                                        NOd(i,hyear1) = hNO
                                        NHd(i,hyear1) = hNH
                                    enddo
                                endif
                            else
                                if (itz2 .lt. recs(hyear3)) then
                                    do i = itz2+1, recs(hyear3)
                                        NOd(i,hyear3) = hNO
                                        NHd(i,hyear3) = hNH
                                    enddo
                                endif
                                itz2 = 1
                                if (itz1 .gt. 1) then
                                    do i = itz1-1, itz2, -1
                                        NOd(i,hyear1) = hNO
                                        NHd(i,hyear1) = hNH
                                    enddo
                                endif
                                hy = hyear1-1
                                do while ((hy .gt. hyear3) .and. (hy .gt. 0))
                                    do i = 366, 1, -1
                                        NOd(i,hy) = hNO
                                        NHd(i,hy) = hNH
                                    enddo
                                    hy = hy - 1
                                enddo
                            endif   ! hyear1 .eq. hyear3
                            hyear3 = hyear1
                            itz2   = itz1
                            hyear2 = hyear3
         
                        case(3)                     ! fill in of constant year values
                            do i = 1,366
                                NOd(i,hyear1) = hNO
                                NHd(i,hyear1) = hNH
                            enddo
                            hy = hyear1-1
                            do while ((hy .gt. hyear2) .and. (hy .gt. 0))  
                                do i = 366, 1, -1
                                    NOd(i,hy) = hNO
                                    NHd(i,hy) = hNH
                                enddo
                                hy = hy - 1
                            enddo
                            hyear2 = hyear1
                            itz2 = 366
                       end select   ! flag_depo 1-3
                        
                    endif   ! 0 < hyear1 < year
                else  ! ios .ne. 0
                    if (hyear1 .le. 0) then
                        hyear1 = 1
                        hyear2 = 1
                    endif
                    continue
                endif  ! ios = 0
            enddo
        end select  ! flag_depo    

!   fill in of the missing data at the end
        select case (flag_depo)
        case (3)
            if (hyear1 .lt. year) then
                hy = hyear1+1
                do while (hy .le. year)  
                    do i = 366, 1, -1
                        NOd(i,hy) = hNO
                        NHd(i,hy) = hNH
                    enddo
                    hy = hy + 1
                enddo
            else    ! if date is outside the simulation period, it will be completly filled in
              do j = 1, year
                do i = 1, 366   
                    NOd(i,j) = hNO
                    NHd(i,j) = hNH
                enddo
              enddo          
            endif
        
        case default
            if (hyear2 .le. year) then
              if (itz2 .lt. recs(hyear2)) then
                if (.not.flag_mult8910) then
                    hx = iy
                    CALL error_mess(hx,' Not enough data records in deposition file, iostat = ',REAL(ios))
                    WRITE (unit_err,*) ' >>>FORESEE message: Fill next values with same data '
                    WRITE (unit_err,'(A,2I4,A,2I4)')'                      from internal simulation time', itz2, hyear2, '   to', recs(hyear2), year
                endif
                do j = hyear2, year
                    ii = 1
                    if (j .eq. hyear2) ii = itz2
                    do i = ii, 366
                        NOd(i,j) = hNO
                        NHd(i,j) = hNH
                    enddo
                enddo 
              else
                hy = hyear2+1
                do while (hy .le. year)  
                    do i = 366, 1, -1
                        NOd(i,hy) = hNO
                        NHd(i,hy) = hNH
                    enddo
                    hy = hy + 1
                enddo
              endif
            else    ! if date is outside the simulation period, it will be completly filled in
              do j = 1, year
                do i = 1, 366   
                    NOd(i,j) = hNO
                    NHd(i,j) = hNH
                enddo
              enddo          
            endif
        end select
    close (unit_dep)
  endif

  write (*,*)

END subroutine readdepo


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

SUBROUTINE readredN

use data_out
use data_site
use data_species
use data_stand
use data_simul

implicit none

character  text
integer  hx, unit_red, i,ios
logical  ex

if (.not.flag_mult8910) print *
if (flag_multi .lt. 8) then
  inquire (File = "./input/.", exist = ex)  ! test whether file exists
  inquire (File = redfile(ip), exist = ex)  ! test whether file exists
  IF(ex .eqv. .false.) then
    print *,' >>>FORESEE message: Cannot find data of RedN - internal calculation'
    hx = 0
    CALL error_mess(hx,'Cannot find data of RedN - internal calculation ',REAL(hx))
  else
    print *, ' >>>FORESEE message: Now reading RedN data from file, please wait...'
    unit_red = getunit()
    OPEN (unit_red,FILE=redfile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')

    DO
      READ(unit_red,*) text
      IF (text .ne. '!') THEN
         backspace(unit_red)
         EXIT
      ENDIF
    ENDDO

    read (unit_red,*,iostat=ios) (svar(i)%RedN, i=1,nspecies)
    close (unit_red)
  endif   ! ex
else
  do i = 1, nspecies
    svar(i)%RedN = RedN_list(i, ip)
  enddo
endif   ! flag_multi

IF(flag_limi==0 .OR. flag_limi==1) THEN
   DO i=1,nspecies
      svar(i)%RedN = 1.
   ENDDO
ENDIF

do i = 1,nspecies
    if (svar(i)%RedN .lt. 0) then       ! no values; internal calculation
        if (flag_multi .lt. 8) then
            print *,' >>>FORESEE message: Cannot find data of RedN - internal calculation for', spar(i)%species_short_name
            write (unit_err, '(A,I3,1X,A)') 'Cannot find data of RedN - internal calculation for species ',i, spar(i)%species_short_name
        endif
        flag_redn = .TRUE.
    endif
enddo

   if (.not.flag_mult8910) write (*,*)

END subroutine readredN

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

SUBROUTINE readlit

!use data_climate
use data_out
use data_soil_cn
use data_species
use data_stand
use data_simul

implicit none

character  text
integer  unit_lit, i,ios
integer  nspec_lit
logical  ex
real     help, hx
real, dimension(22) :: helpin

flag_lit = 0
if (flag_mult8910) then
    inquire (File = litfile(1), exist = ex)  ! test whether file exists
else
    print *
    inquire (File = litfile(ip), exist = ex)  ! test whether file exists
endif
  IF(ex .eqv. .false.) then
    if (.not.flag_mult8910) then
        print *,' >>>FORESEE message: Cannot find data of litter initialisation - internal calculation'
        hx = 0.
        write (unit_err,*)
        write (unit_err,*) 'Cannot find data of litter initialisation - internal calculation '
    endif
  else
    if (.not.flag_mult8910) print *, ' >>>FORESEE message: Now reading litter initialisation data from file, please wait...'
!   now read data from file
    unit_lit = getunit()
    OPEN (unit_lit,FILE=litfile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')

    do
      read(unit_lit,*) text
      IF (text .ne. '!') then
         backspace(unit_lit)
         exit
      endif
    enddo

    helpin = 0.
    slit%C_opm_fol = 0.
    read (unit_lit,*) nspec_lit
    read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_fol, i=1,nspec_lit)
    read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_tb , i=1,nspec_lit)
    read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_frt(1), i=1,nspec_lit)
    read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_crt(1), i=1,nspec_lit)
    read (unit_lit,*,iostat=ios) text, (slit(i)%C_opm_stem,i=1,nspec_lit)
    flag_lit = 1

    help = 0.
    hx   = 1.
    do i=1,nspecies
       if (slit(i)%C_opm_fol .gt. 0) then       
         totfol_lit  = totfol_lit  + slit(i)%C_opm_fol
         totfrt_lit  = totfrt_lit  + slit(i)%C_opm_frt(1)
         tottb_lit   = tottb_lit   + slit(i)%C_opm_tb
         totcrt_lit  = totcrt_lit  + slit(i)%C_opm_crt(1)
         totstem_lit = totstem_lit + slit(i)%C_opm_stem
       else
          hx = -1.
       endif
    enddo
    help = totfol_lit

    if ((help .gt. 0.) .or. (hx .gt. 0) .and. .not.flag_mult8910) then
       CALL error_mess(0,'Using data of litter initialisation from file '//trim(litfile(ip)), hx)    
    else
      ! no values; internal calculation of litter initialisation
        if (.not.flag_mult8910) then
           print *,' >>>FORESEE message: No data of litter initialisation - internal calculation'
           hx = 0.
           CALL error_mess(0,'No data of litter initialisation - internal calculation ', hx)
        endif
       flag_lit = 0
    endif
    close (unit_lit)
  endif   ! ex

  if (.not.flag_mult8910) write (*,*)

END subroutine readlit

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

subroutine prepare_climate         
! read climate file

use data_climate
use data_out
use data_simul
use data_stand

implicit none

type clifile     ! new data type for all climate parameters
  integer :: day,mon,ye
  real :: m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11
end type clifile
type (clifile), allocatable,dimension(:,:) :: climall !variable for data type climfile
character(1) c
character :: text
integer :: i,j,ios, unit_cli
integer :: realrec = 0
integer :: repflag = 0
logical :: ex

if (.not.flag_mult8910) then
    print *, ' '
    print *, '  Input of climate data:    '
endif

call testfile(climfile(ip),ex) !input filename, test whether file exists
 IF(ex .eqv. .false.) then
 print *,' >>>FORESEE message: Cannot find climate data - program STOP!'
 stop
 endif
 if (.not.flag_mult8910) print *, ' >>>FORESEE message: Now reading CLIMATE data from file, please wait...'
!now read data from file
unit_cli = getunit()
OPEN (unit_cli,FILE=climfile(ip),IOSTAT=ios,STATUS='OLD',ACTION='READ')
allocate (recs (1:year))
allocate (dd (1:366,1:year));allocate (mm (1:366, 1:year))
allocate (yy (1:year));allocate (tp (-2:366,1:year))
allocate (hm (0:366,1:year));allocate (prc (0:366,1:year))
allocate (prs (0:366,1:year));allocate (rd (0:366,1:year))
allocate (tn (0:366,1:year))
allocate (tx (0:366,1:year))
allocate (vp (0:366,1:year))
allocate (sdu (0:366,1:year))
allocate (wd (0:366,1:year))
allocate (sde (0:366,1:year))
allocate (bw (0:366,1:year))

dd = -99.9
mm = -99.9
yy = -99.9
tn = -99.9
tx = -99.9
wd = -99.9    ! wind initialisation

IF (index(climfile(ip),'.cli') .ne. 0) then
    flag_climtyp = 1
    do
        read(unit_cli,*) text
        IF (text .ne. '!') then
           IF (text .eq. 'N') then
              flag_climtyp = 2
           else IF(text.eq.'T') then
              flag_climtyp = 3
           else
              backspace(unit_cli)
              exit
           endif
        endif
    enddo
else if (index(climfile(ip),'.txt') .ne. 0) then
    flag_climtyp = 4
else
   flag_climtyp = 5  
end IF
call read_cli
close(unit_cli)
if (flag_end .gt. 0) return


IF (realrec < year .and. repflag == 0) then
   year = realrec
else
   IF (repflag == 1) then
      call climfill
   end IF
end IF
med_rad1 = 0.
do j = 1, year-1
  tp(0,j+1) = tp(recs(j),j)
  tp(-1,j+1)= tp(recs(j)-1,j)
  tp(-2,j+1)= tp(recs(j)-2,j)
  hm(0,j+1) = hm(recs(j),j);prc(0,j+1) = prc(recs(j),j);prs(0,j+1) = prs(recs(j),j)
  rd(0,j+1) = rd(recs(j),j)
  wd(0,j+1) = wd(recs(j),j)
  bw(0,j+1) = bw(recs(j),j)
  vp(0,j+1) = vp(recs(j),j)
  sdu(0,j+1) = sdu(recs(j),j)
  sde(0,j+1) = sde(recs(j),j)
  tx(0,j+1) = tx(recs(j),j)
  tn(0,j+1) = tn(recs(j),j)

  if( yy(j) .eq.time_b) then
    do i=1, recs(j)

       med_rad1 = med_rad1 + rd(i, j)
    end do
    med_rad1 = med_rad1/recs(1)

  end if
end do
tp(-2,1) = tp(1,1); tp(-1,1) = tp(1,1); tp(0,1) = tp(1,1)
hm(0,1) = hm(1,1);prc(0,1) = prc(1,1);prs(0,1) = prs(1,1)
rd(0,1) = rd(1,1)
wd(0,1)=wd(1,1)
vp(0,1) = vp(1,1)
bw(0,1) = bw(1,1)
tn(0,1) = tn(1,1)
tx(0,1) = tx(1,1)
sdu(0,1) =sdu(1,1)
sde(0,1) = sde(1,1)

contains

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

subroutine read_dwd

character(3) text
integer help, help1, help2, help3
allocate (climall (0:366,1:year))

j=1
c = 'n'
do
   IF (j > year) then
      realrec = year
      exit
   end IF
   if (.not.flag_mult8910) print *, 'Year ',j
   read(unit_cli,*) text
   if(text.ne.'ta ') then
      backspace(unit_cli)
   end if

   do i = 1, 366
      read (unit_cli,*,IOSTAT=ios) climall(i,j)

     help2 = climall(i,j)%day
     help3 = climall(i,j)%mon
     help =  climall(i,j)%ye
     help1 = climall(i-1,j)%ye
     if (help.eq.2099 .and.help1.eq.2100.and. i.eq.366) then
      end if
   end do
  IF (climall(365,j)%ye == climall(366,j)%ye) then
      recs(j) = 366
   else
      backspace unit_cli
      climall(366,j)%day = 0
      climall(366,j)%mon = 0
      climall(366,j)%ye = 0
      recs(j) = 365
	  help = help-1
   end IF
   IF (j < year .and. ios < 0 .and. c .eq. 'n') then
      realrec = j
      if (.not.flag_mult8910) then 
          print *, ' >>>FORESEE message: Not enough climate data records in file!'
          call error_mess(0,'read_cli: Not enough data records in climate file; number of complete years: ',real(realrec))
          write(unit_err,'(A,I5)')' read_cli: Fill next values with same from first year, day: ',i_exit
          write(unit_err,'(A,I5)')' read_cli: Fill next values with same data up to years:     ',year
          repflag = 1
          exit
      endif
    else if(j.eq.year.and.ios < 0) then
	    realrec = year
		exit
   end IF

   j=j+1
    if(help.lt.time_b) j = j-1

end do

do j = 1, realrec
  yy(j) = climall(1,j)%ye
  do i = 1, recs(j)
    dd(i,j) = climall(i,j)%day
    mm(i,j) = climall(i,j)%mon
    tx(i,j) = climall(i,j)%m1
    tp(i,j) = climall(i,j)%m2
    tn(i,j) = climall(i,j)%m3

    prc(i,j) = climall(i,j)%m4
    hm(i,j) = climall(i,j)%m5
    prs(i,j) = climall(i,j)%m6
    vp(i,j) = climall(i,j)%m7
    sdu(i,j) = climall(i,j)%m8
    bw(i,j) = climall(i,j)%m9
    rd(i,j) = climall(i,j)%m10

    wd(i,j) = climall(i,j)%m11
  end do
end do

close(9)
deallocate (climall) 
end subroutine read_dwd

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

subroutine read_cli

implicit none

integer  :: testtext, hp
character(11) :: text2
character(4) :: text

testtext=0
c      = 'n'
j      = 1
hp = 0
read(unit_cli,'(A)') text2
hp = index(text2,'.')  
backspace(unit_cli)

do
    IF(j > year) exit

    select case(flag_climtyp)
    case (1)
      do i=1,366
        if (hp .gt. 0) then
            read(unit_cli,*,iostat=ios) text2,tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j)
            text = text2(1:2)
            write (text,'(A)') text2(1:2)
            read (text,*) dd(i,j)
            write (text,'(A)') text2(4:5)
            read (text,*) mm(i,j)
            write (text,'(A)') text2(7:10)
            read (text,*) yy(j)
        else
            read(unit_cli,*,iostat=ios) dd(i,j),mm(i,j),yy(j),tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j)
        endif  ! hp
        i_exit = i
        if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
            recs(j) = i
            write (*,*) 'Year ',j, yy(j)
            realrec = j
            if (j .eq. year) ios = -10 
            exit
        endif
        if (ios .ne. 0) exit
      end do

    case (2)
      do i=1,366
        read(unit_cli,*) dd(i,j),mm(i,j),yy(j),&
                                    tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j),wd(i,j)
        i_exit = i
        if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
            recs(j) = i
            write (*,*) 'Year ',j, yy(j)
            realrec = j
            if (j .eq. year) ios = -10 
            exit
        endif
        if (ios .ne. 0) exit
      end do

    case (3)
      do i=1,366
        if (hp .gt. 0) then
            read(unit_cli,*,iostat=ios) text2, &
                                    tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j),wd(i,j), tx(i,j),tn(i,j)
            text = text2(1:2)
            write (text,'(A)') text2(1:2)
            read (text,*) dd(i,j)
            write (text,'(A)') text2(4:5)
            read (text,*) mm(i,j)
            write (text,'(A)') text2(7:10)
            read (text,*) yy(j)
        else
            read(unit_cli,*,iostat=ios) dd(i,j),mm(i,j),yy(j),&
                                    tp(i,j),hm(i,j),prc(i,j),prs(i,j),rd(i,j),wd(i,j), tx(i,j),tn(i,j)
       endif
        i_exit = i
        if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
            recs(j) = i
            write (*,*) 'Year ',j, yy(j)
            realrec = j
            if (j .eq. year) ios = -10 
            exit
        endif
        if (ios .ne. 0) exit
      end do

    case (4)   ! suffix 'txt'
      if (j .eq. 1 .and. testtext.eq.0) then
          read(unit_cli,*) text
          testtext = 1
      end if    
      do i=1,366
        read(unit_cli,*,iostat=ios) dd(i,j),mm(i,j),yy(j),&
                                     tx(i,j),tp(i,j),tn(i,j),prc(i,j),hm(i,j),prs(i,j),rd(i,j),wd(i,j)
        i_exit = i
        if ((dd(i,j) .eq. 31) .and. (mm(i,j) .eq. 12)) then
            recs(j) = i
            write (*,*) 'Year ',j, yy(j)
            realrec = j
            if (j .eq. year) ios = -10 
            exit
        endif
        if (ios .ne. 0) exit
      end do
    
    case (5 )
      call read_dwd  
      exit

    end select

    IF (realrec .lt. year .and. ios .ne. 0 .and. c .eq. 'n') then
        if (dd(i_exit,j) .gt. 0)  i_exit = i_exit+1
        if (i_exit .ge. 365) i_exit = 1
        repflag = 1
        if (.not.flag_mult8910) then
            print *, ' >>>FORESEE message: Not enough data records in file'
            print *, ' IOSTAT = ', ios

            WRITE (*,'(A,I5)') ' >>>FORESEE message: Fill next values with same data from day number', i_exit
            CALL error_mess(0,'read_cli: Not enough data records in meteorology file; number of complete years: ',real(realrec))
            write(unit_err,'(A,I5)')' read_cli: Fill next values with same from first year, day: ',i_exit
            write(unit_err,'(A,I5)')' read_cli: Fill next values with same data up to years:     ',year
            exit
        endif

    end if
  if (ios .ne. 0) exit
  if (yy(j) .ge. time_b) then
      if ((j .eq. 1) .and. (yy(j) .gt. time_b)) then
        CALL error_mess(0,'read_cli: No climate data in meteorology file for year ',real(time_b))
        flag_end = 6
        return
      endif
      j = j+1
  endif
end do

end subroutine read_cli

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

subroutine climfill

integer istart

istart = i_exit 
if(istart.eq.0) istart =istart +1
do j=realrec+1,year
   print *,'Year ',j
   yy(j)=yy(j-realrec)
   recs(j)=recs(j-realrec)
   do i=istart,366
      dd(i,j) = dd(i,j-realrec)
      mm(i,j) = mm(i,j-realrec)
      tp(i,j) = tp(i,j-realrec)
      hm(i,j)  = hm(i,j-realrec)
      prc(i,j) = prc(i,j-realrec)
      prs(i,j) = prs(i,j-realrec)
      rd(i,j)  = rd(i,j-realrec)
      wd(i,j)  = wd(i,j-realrec)
      tx(i,j) = tx(i,j-realrec)
      tn(i,j) = tn(i,j-realrec)
   end do
end do

end subroutine climfill

END subroutine prepare_climate

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

SUBROUTINE store_para(hpara, simpara, parerr)

use data_simul
use data_out
use data_par
use data_species
use data_soil_cn
use data_stand
use data_tsort
implicit none

integer  inum
real  hpara, parerr
character(100):: simpara, hchar1 
integer, external :: array_num

if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' store_para'

parerr = 0.
if (trim(simpara) .eq. 'year') then
   year=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'time_b') then
    time_b=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'kpatchsize') then
    kpatchsize=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'dz') then
    dz=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'ns_pro') then
    ns_pro=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_mort') then
    flag_mort=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_reg') then
    flag_reg=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_stand') then
    flag_stand=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_sveg') then
    flag_sveg=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_mg') then
    flag_mg=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_dis') then
    flag_dis=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_light') then
    flag_light=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_folhei') then
    flag_folhei=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_volfunc') then
    flag_volfunc=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_resp') then
    flag_resp=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_limi') then
    flag_limi=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_sign') then
    flag_sign=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_decomp') then
    flag_decomp=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_wred') then
    flag_wred=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_wurz') then
    flag_wurz=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_cond') then
    flag_cond=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_int') then
    flag_int=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_eva') then
    flag_eva=hpara
   parerr = 1.
   return
endif
if ((trim(simpara) .eq. 'flag_co2') .or.(trim(simpara) .eq. 'flag_CO2')) then
    flag_co2=hpara
   parerr = 1.
   return
endif
if (adjustl(trim(simpara)) .eq. 'flag_sort') then
    flag_sort  = hpara
    parerr = 1.
    return
endif
if (adjustl(trim(simpara)) .eq. 'flag_wpm') then
    flag_wpm  = hpara
    parerr = 1.
    return
endif
if (trim(simpara) .eq. 'time_out') then
    time_out=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_dayout') then
    flag_dayout=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_cohout') then
    flag_cohout=hpara
   parerr = 1.
   return
endif
if (trim(simpara) .eq. 'flag_sum') then
    flag_sum=hpara
   parerr = 1.
   return
endif

if (trim(simpara) .eq. 'k_hum') then
    k_hum=hpara
    parerr = 1.
    return
endif

if (trim(simpara) .eq. 'k_hum_r') then
    k_hum_r=hpara
    parerr = 1.
    return
endif

if (trim(simpara) .eq. 'k_nit') then
    k_nit=hpara
    parerr = 1.
    return
endif
if (adjustl(trim(simpara)) .eq. 'alfm') then
    alfm   = hpara
    parerr = 1.
    return
endif
if (adjustl(trim(simpara)) .eq. 'gpmax') then
    gpmax  = hpara
    parerr = 1.
    return
endif
if (adjustl(trim(simpara)) .eq. 'alfm') then
    alfm  = hpara
    parerr = 1.
    return
endif

! Species parameter
hchar1 = adjustl(simpara)
inum = array_num(hchar1)
if (hchar1(1:9) .eq. 'k_opm_fol') then
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%k_opm_fol  = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:9) .eq. 'k_opm_frt') then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%k_opm_frt  = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:9) .eq. 'k_syn_fol') then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%k_syn_fol  = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:9) .eq. 'k_syn_frt') then
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%k_syn_frt  = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:3) .eq. 'psf') then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%psf  = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:7) .eq. 'Phmodel') then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%Phmodel  = hpara
        parerr = 1.
        return
	endif
endif
if ((hchar1(1:4) .eq. 'pnus') .or. (hchar1(1:4) .eq. 'Pnus')) then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%pnus = hpara
        parerr = 1.
        return
	endif
endif
if ((hchar1(1:4) .eq. 'RedN') .or. (hchar1(1:4) .eq. 'redn')) then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    svar(inum)%RedN = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:4) .eq. 'prms') then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%prms = hpara
        parerr = 1.
        return
	endif
endif
if (hchar1(1:4) .eq. 'prmr') then
    inum = array_num(hchar1)
	if (inum .gt. 0 .and. inum .le. nspecies) then
	    spar(inum)%prmr = hpara
        parerr = 1.
        return
	endif
endif

END subroutine store_para

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

integer FUNCTION array_num(string)

! reads the field numbre out of an array and hands it back as integer

implicit none

integer ipos1, ipos2, inum
character (100) string
character (10) help, hchar

    ipos1 = scan(string, '(' )
    ipos2 = scan(string, ')' )
	ipos1 = ipos1+1
	ipos2 = ipos2-1
	hchar = string(ipos1:ipos2)
    write(help,'(A3)') hchar
    read(help,*) inum
    array_num = inum

end function array_num
