source_code/version_2.3_windows/readsim.f

+!*****************************************************************!
+!*                                                               *!
+!*              4C (FORESEE) Simulation Model                    *!
+!*                                                               *!
+!*                    Subroutines for:                           *!
+!*    - READSIM:   Read simulation options from file             *!
+!*    - ALLOFILE:  Allocate simulation files                     *!
+!*    - READCON                                                  *!
+!*                                                               *!
+!*                  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 readsim
+
+! read simulation options from file
+
+use data_mess
+use data_out
+use data_par
+use data_simul
+use data_species
+use data_stand
+use data_site
+use data_tsort
+use data_climate
+
+
+implicit none
+
+logical ex
+integer i, ios, ios1, nowunit, nowunit1, k, anzclim, j, l, helpi, helpw, helpy, ihelp, ilen
+integer ltsunit
+character:: a, ttext
+character (150) tspec, tname, tclim, tval, tsite, tman, ttree, tdepo, tred, tlit, tsoilid, &
+                    pathdir1, pathdir2,pathdir3, pathdir4, pathdir5, pathdir6, pathdir7, &
+                    climszen, siteall, climall,site_name_all
+character(50), dimension(:), allocatable:: site_name_ad
+character(50), dimension(:), allocatable:: climfile_ad
+character(50), dimension(:), allocatable:: sitefile_ad
+character(50), dimension(:), allocatable:: manfile_ad
+character(50), dimension(:), allocatable:: treefile_ad
+character(50), dimension(:), allocatable:: wpmfile_ad
+character(50), dimension(:), allocatable:: depofile_ad
+character(50), dimension(:), allocatable:: redfile_ad
+character(50), dimension(:), allocatable:: litfile_ad
+character(150):: text 
+character(50) :: istand 
+character(10) :: helpsim, text4 
+integer :: spec1, spec2, tm
+real  :: h1, h2
+ 
+real, dimension(:), allocatable:: clim_long, clim_lat, clim_height    ! coordinates and height of climate stations
+character(10), dimension(:), allocatable:: climnum
+character(50), dimension(:), allocatable:: clim_nam
+
+    nowunit = getunit()
+    ios  = 0
+    nvar = 0
+
+    call testfile(simfile,ex)
+    if(ex .eqv. .false.) return
+    open(nowunit,file=simfile,iostat=ios,status='old',action='read')
+
+    read(nowunit,*,iostat=ios) flag_multi
+    
+    if(flag_multi .ge. 1) then      
+        read(nowunit,*,iostat=ios) site_nr
+
+        if(flag_multi .eq. 9 .or. flag_multi .eq. 10) then
+            flag_mult910 = .True.
+        else
+            flag_mult910 = .False.
+        endif
+        
+        if((flag_mult910 .or. flag_multi .eq. 8) .and. (site_nr .gt.1)) then
+            flag_mult8910 = .True.
+        else
+            flag_mult8910 = .False.
+        endif
+
+		repeat_number = site_nr
+              allocate(sitenum(site_nr))
+              allocate(clim_id(site_nr))
+              allocate(soilid(site_nr))
+              allocate(gwtable(site_nr))
+              allocate(NOdep(site_nr))
+              allocate(NHdep(site_nr))
+              clim_id = "xxx"
+              NOdep = 0.
+              NHdep = 0.
+    endif
+
+    select case (flag_multi)
+    case (1, 4)
+        flag_clim = 1
+    case (7, 8, 9, 10)
+        flag_clim = 1
+        flag_trace = .FALSE.
+    case default
+        flag_clim = 0
+    end select
+
+    read(nowunit,*,iostat=ios) ! skip comment line 'simulation specifications'
+    read(nowunit,*,iostat=ios) year
+    read(nowunit,*,iostat=ios) time_b
+    read(nowunit,*,iostat=ios) kpatchsize
+    read(nowunit,*,iostat=ios) dz
+    read(nowunit,*,iostat=ios) ns_pro
+    read(nowunit,*,iostat=ios) ! skip comment line 'choice of model options'
+    read(nowunit,*,iostat=ios) flag_mort
+    read(nowunit,*,iostat=ios) flag_reg
+    read(nowunit,*,iostat=ios) flag_lambda
+    read(nowunit,*,iostat=ios) flag_stand
+    read(nowunit,*,iostat=ios) flag_sveg
+    read(nowunit,*,iostat=ios) flag_mg
+    read(nowunit,*,iostat=ios) flag_dis
+    read(nowunit,*,iostat=ios) flag_light
+    read(nowunit,*,iostat=ios) flag_folhei
+    read(nowunit,*,iostat=ios) flag_volfunc
+    read(nowunit,*,iostat=ios) flag_resp
+    read(nowunit,*,iostat=ios) flag_limi
+    read(nowunit,*,iostat=ios) flag_decomp
+    read(nowunit,*,iostat=ios) flag_sign
+    read(nowunit,*,iostat=ios) flag_wred
+    read(nowunit,*,iostat=ios) flag_wurz
+    read(nowunit,*,iostat=ios) flag_cond
+    read(nowunit,*,iostat=ios) flag_int
+    read(nowunit,*,iostat=ios) flag_eva
+    read(nowunit,*,iostat=ios) flag_co2
+    read(nowunit,*,iostat=ios) flag_sort
+    read(nowunit,*,iostat=ios) flag_wpm
+    read(nowunit,*,iostat=ios) flag_stat
+    read(nowunit,*,iostat=ios) ! skip comment line 'output specifications'
+    read(nowunit,*,iostat=ios) time_out
+
+    ! test lamda_ts
+
+   if(flag_lambda.eq.1) then
+       allocate(lambda_ts(168,3))
+       ltsunit=getunit()
+       open (ltsunit,file='input/lambdats_oak_pine.par', IOSTAT=ios,status='old')
+       read (ltsunit,*), text, spec1, spec2
+       !write(4567,*)text, spec1,spec2
+     
+    do j=1,168
+       read(ltsunit,*) tm, h1, h2
+       
+  
+             lambda_ts(j,1)= tm
+             lambda_ts(j,2) = h1
+             lambda_ts(j,3) = h2
+!             write(4567,*) lambda_ts(j,1), lambda_ts(j,2), lambda_ts(j,3)
+                
+       end do
+       
+   end if   
+!     define name of yearly output variables
+      nyvar = 1
+      read(nowunit,*,iostat=ios) outy_file(nyvar)
+      do while (trim(outy_file(nyvar)) .ne. 'end')
+         nyvar = nyvar + 1
+         read(nowunit,*) outy_file(nyvar)
+      enddo
+
+     read(nowunit,*,iostat=ios) flag_dayout
+!     define name of daily output variables
+      ndvar = 1
+      read(nowunit,*) outd_file(ndvar)
+      do while (trim(outd_file(ndvar)) .ne. 'end')
+          ndvar = ndvar + 1
+          read(nowunit,*) outd_file(ndvar)
+      enddo
+
+    read(nowunit,*,iostat=ios) flag_cohout
+!     define name of cohort output variables
+      ncvar = 1
+      read(nowunit,*) outc_file(ncvar)
+      do while (trim(outc_file(ncvar)) .ne. 'end')
+         ncvar = ncvar + 1
+         read(nowunit,*) outc_file(ncvar)
+      enddo
+      
+    read(nowunit,*,iostat=ios) flag_sum
+    read(nowunit,*,iostat=ios) ! skip comment line 'input'
+
+    if (.not.flag_mult910) call allofile
+
+ SELECT CASE(flag_multi)
+ CASE (0,1,2,3,6)
+      jpar = 0
+      DO i=1,site_nr
+        if(i .gt. 1)then
+          read(nowunit,*,iostat=ios) ! skip comment line 'run number'
+          do
+             jpar = jpar + 1
+             read(nowunit,*) vpar(jpar), simpar(jpar)
+             if (vpar(jpar) .lt. -90.0) exit
+          enddo
+        endif
+
+        read(nowunit,'(A)',iostat=ios) specfile(i)
+        read(nowunit,'(A)') site_name(i)
+        read(nowunit,'(A)') climfile(i)
+        read(nowunit,'(A)') sitefile(i)
+        read(nowunit,'(A)') valfile(i)
+        read(nowunit,'(A)') treefile(i)
+        read(nowunit,'(A)') standid(i)
+        read(nowunit,'(A)') manfile(i)
+        read(nowunit,'(A)') depofile(i)
+        read(nowunit,'(A)') redfile(i)
+        read(nowunit,'(A)',iostat=ios) litfile(i)
+        
+      ! fill clim_id   
+        clim_id(i) = climfile(i)
+      ios1=-1
+      ! measurements
+	   if(flag_multi.ne.2) then
+        if (ios .eq. 0) read(nowunit,'(A)',iostat=ios1) text
+        if (ios1 .eq. 0) then
+            if (flag_stat .gt. 0 .and. i .eq. 1) then
+               allocate (mesfile(anz_mesf)) 
+               mesfile(1) = text
+               ttext = adjustl(text)
+               if (ttext .eq. '!' .or. ttext .eq. '*') then
+                  backspace (nowunit)
+               else
+                  if (.not.flag_mult8910) write (*, '(A, I3,A,A)')' >>>foresee message: site_nr ',i,'; filename of measurements:  ', trim(mesfile(1))
+               endif 
+            else
+               ttext = adjustl(text)
+               if (ttext .eq. '!' .or. ttext .eq. '*') backspace (nowunit)
+            endif
+        endif
+        end if
+     if (.not.flag_mult8910) print *, ' >>>foresee message: site_nr ',i,'; input of filenames completed'
+      end DO
+
+ CASE (4, 5, 8)
+      allocate(latitude(site_nr))
+      allocate(RedN_list(15, site_nr))
+      RedN_list = -99.9
+      read(nowunit,'(A)',iostat=ios) specfile(1)
+      read(nowunit,'(A)') site_name(1)
+      read(nowunit,'(A)') treefile(1)
+      read(nowunit,'(A)') manfile(1)
+      read(nowunit,'(A)') siteall       !   control xxx.con
+      read(nowunit,'(A)') climall       ! climate stations with coordination
+      read(nowunit,'(A)') pathdir1      ! path for climate scenarios 
+      read(nowunit,'(A)') pathdir2      ! path for soil file xxx.sop or name of total soil file (flag_multi=8)
+      read(nowunit,'(A)') climszen      ! labeling climate scenarios
+      if (flag_multi .eq. 8.or.flag_multi.eq.5) read(nowunit,*) text   ! BRB / BAWUE / DEU
+      if (.not.flag_mult8910) print *, ' >>>foresee message: Input of filenames completed'
+
+      site_name1   = site_name(1)
+!  define name of output variables
+      nvar = 1
+      read(nowunit,*) outvar(nvar)
+      do while (trim(outvar(nvar)) .ne. 'end')
+         nvar = nvar + 1
+         read(nowunit,*) outvar(nvar)
+      enddo
+      if (nvar .gt. 1) allocate(output_var(nvar,site_nr,0:year))
+
+          helpw = 0
+          helpi = 0
+          do i = 1, nvar-1
+             select case (trim(outvar(i)))
+
+              case ('AET_mon','AETmon','aetmon','aet_mon','cwb_mon','cwbmon','PET_mon','PETmon','petmon','pet_mon', &
+                    'GPP_mon','GPPmon','gppmon','gpp_mon','NEP_mon','NEPmon','nepmon','nep_mon','NPP_mon','NPPmon','nppmon','npp_mon', &
+                    'perc_mon','percmon','temp_mon','tempmon','prec_mon','precmon', 'resps_mon','respsmon','TER_mon','TERmon','ter_mon','termon')  
+                flag_cum = 1
+                helpi    = helpi + 1
+                output_var(i,1,0) = 1.*helpi              ! field numbre of monthly value
+
+              case ('AET_week','AETweek','aetweek','aet_week','cwb_week','cwbweek','PET_week','PETweek','petweek','pet_week', &
+                    'GPP_week','GPPweek','gppweek','gpp_week','NEP_week','NEPweek','nepweek','nep_week','NPP_week','NPPweek','nppweek','npp_week', &
+                    'perc_week','percweek','temp_week','tempweek','prec_week','precweek', 'resps_week','respsweek', 'TER_week','TERweek','ter_week','terweek')  
+                flag_cum = 1
+                helpw    = helpw + 1
+                output_var(i,1,0) = 1.*helpw              ! field numbre of weekly values
+
+             end select   ! outvar
+      
+          enddo
+          if (helpi .gt. 0) then
+             allocate(output_varm(helpi,site_nr,year,12))
+          endif
+          if (helpw .gt. 0) then
+             allocate(output_varw(helpw,site_nr,year,52))
+          endif
+
+      call errorfile(simfile, ios, nowunit)
+
+!  reading file with description of climate stations used
+      allocate(climnum(3000))
+      allocate(clim_long(3000))
+      allocate(clim_lat(3000))
+      allocate(clim_height(3000))
+      allocate(clim_nam(3000))
+
+      call testfile(climall,ex)
+      if (ex .eqv. .false.) return
+      nowunit1 = getunit()
+      ios1 = 0
+      open(nowunit1,file=climall,iostat=ios,status='old',action='read')
+      k=1
+      do
+         READ(nowunit1,'(A)',iostat=ios1) a
+         IF (a .ne. '!') exit
+
+      end do
+      backspace nowunit1
+
+      do
+          read(nowunit1,*,iostat=ios1) climnum(k), clim_long(k),clim_lat(k),    &
+                                    clim_height(k)  
+          if(ios1 .lt. 0) exit
+          k = k+1
+      end do
+      anzclim = k-1
+      ios1 = 0
+
+      call errorfile(climall, ios1, nowunit1)
+
+! reading control file with site-id, climate-id, soil-id, gwtabe-id
+      call testfile(siteall,ex)
+      if (ex .eqv. .false.) return
+      nowunit1 = getunit()
+      open(nowunit1,file=siteall,iostat=ios1,status='old',action='read')
+      do
+         READ(nowunit1,'(A)',iostat=ios1) a
+         IF (a .ne. '!') exit
+
+      end do
+      backspace nowunit1
+!      if (flag_multi .eq. 8) read(nowunit1,*) text   ! BRB / BAWUE / DEU
+ 
+      select case (trim(text))
+      case ('BRB')
+        flag_climnam = 1
+      case ('BAWUE')
+        flag_climnam = 2
+      case ('DEU')
+        flag_climnam = 3
+	  case ('REMO')
+	    flag_climnam = 4
+	  case('WETTREG')
+	    flag_climnam =5
+      end select
+
+      do i=1,site_nr
+          select case (flag_multi)
+          case (4)
+             read(nowunit1,*,iostat=ios1) sitenum(i), clim_id(i), soilid(i), gwtable(i)
+             flag_climnam = 1
+            
+             sitefile(i) =trim(pathdir2)//'wbuek'//trim(soilid(i))//'.sop'
+             valfile(i)  =trim(pathdir2)//'wbuek'//trim(soilid(i))//'.soi'
+             standid(i)   = sitenum(i)
+
+          case (5,8)
+             call readcon(i, nowunit1)  
+             soilid(i) = adjustl(soilid(i))
+             ihelp = len(trim(soilid(i)))
+             sitefile(i) = trim(pathdir2)
+               if( flag_climnam.eq.3) then
+			      climfile(i) = trim(pathdir1)//trim(clim_id(i))//trim(climszen)//'.dat'
+			   end if
+               if(flag_climnam.eq.4) then
+			      climfile(i) = trim(pathdir1)//'gp_'//trim(clim_id(i))//'_'//trim(climszen)//'.txt'
+			   end if
+
+               if(flag_climnam.eq.5) then
+			      climfile(i) = trim(pathdir1)//trim(clim_id(i))//'_'//trim(climszen)//'.dat'
+               end if
+		  end select
+
+
+         do j = 1,anzclim
+           if(clim_id(i).eq.climnum(j)) then
+             select case (flag_climnam)
+ 
+             case (1)    ! WK
+                if(flag_climtyp .eq. 5) then
+                   climfile(i) = trim(pathdir1)//trim(clim_nam(j))//trim(climszen)//'.dat'
+                else
+                   climfile(i) = trim(pathdir1)//trim(clim_nam(j))//trim(climszen)//'.cli'
+                end if
+
+             case (2)    ! Klara
+                climfile(i) = trim(pathdir1)//trim(climnum(j))//trim(climszen)//'.dat'
+             end select
+             latitude(i) = clim_lat(j)
+             exit
+           end if
+           if (j .eq. anzclim) then
+              write (unit_err,*) '***  4C-error - searching in file:', trim(climall)
+              write (unit_err,*) '                no climate station found for climate id: ', clim_id(i)
+              write (unit_err,*) 
+           endif
+         end do
+
+!          fill in sitefile
+         site_name(i) = site_name(1)
+         specfile(i)  = specfile(1)
+         treefile(i)  = treefile(1)
+         manfile(i)   = manfile(1)
+         depofile(i)  = 'dummy.dep'
+         redfile      = 'dummy.red'
+         litfile      = 'dummy.lit'
+      enddo
+
+      if ((.not.flag_mult8910) .and. (ios1 .lt. 0)) print *, 'no information for site number ', i
+      call errorfile(siteall, ios1, nowunit1)
+
+      deallocate(climnum)
+      deallocate(clim_long)
+      deallocate(clim_lat)
+      deallocate(clim_height)
+      deallocate(clim_nam)
+
+      close(nowunit1)
+
+!  variation of flag_multi= 5, especially for SILVISTRAT
+ CASE (7)
+
+      allocate(site_name_ad(site_nr))
+      allocate(climfile_ad(site_nr))
+      allocate(sitefile_ad(site_nr))
+      allocate(manfile_ad(site_nr))
+      allocate(treefile_ad(site_nr))
+      allocate(depofile_ad(site_nr))
+      allocate(redfile_ad(site_nr))
+      allocate(litfile_ad(site_nr))
+      
+      allocate(fl_co2(site_nr))
+
+      read(nowunit,'(A)',iostat=ios) specfile(1)
+      read(nowunit,'(A)') site_name_all
+      read(nowunit,'(A)') siteall
+      read(nowunit,'(A)') pathdir1    ! path climate file
+      read(nowunit,'(A)') pathdir2    ! path soil file
+      read(nowunit,'(A)') pathdir3    ! path treeini file
+      read(nowunit,'(A)') pathdir4    ! path management file
+      read(nowunit,'(A)') pathdir5    ! path deposition file
+      read(nowunit,'(A)') pathdir6    ! path RedN file
+      read(nowunit,'(A)') pathdir7    ! path litter file
+
+      call errorfile(simfile, ios, nowunit)
+
+! reading control file with site-id,name, climate scenario, soil-id, man-file, treeini-file, dep-file
+
+      call testfile(siteall,ex)
+      if (ex .eqv. .false.) return
+
+      nowunit1 = getunit()
+
+      open(nowunit1,file=siteall,iostat=ios1,status='old',action='read')
+      do
+         READ(nowunit1,'(A)',iostat=ios1) a
+         IF (a .ne. '!') exit
+
+      end do
+      backspace nowunit1
+
+      do i=1,site_nr
+         read(nowunit1,*,iostat=ios1) sitenum(i),site_name_ad(i), climfile_ad(i),sitefile_ad(i),treefile_ad(i), &
+                                      manfile_ad(i),depofile_ad(i),redfile_ad(i),litfile_ad(i), fl_co2(i) 
+         specfile(i) = specfile(1)
+         standid(i)  = sitenum(i)
+         site_name(i)= trim(site_name_all)//trim(site_name_ad(i))
+         climfile(i) = trim(pathdir1)//trim(climfile_ad(i))
+         sitefile(i) = trim(pathdir2)//trim(sitefile_ad(i))
+         treefile(i) = trim(pathdir3)//trim(treefile_ad(i))
+         manfile(i)  = trim(pathdir4)//trim(manfile_ad(i))
+         depofile(i) = trim(pathdir5)//trim(depofile_ad(i))
+         redfile(i)  = trim(pathdir6)//trim(redfile_ad(i))
+         litfile(i)  = trim(pathdir7)//trim(litfile_ad(i))
+
+      enddo
+       call errorfile(siteall, ios1, nowunit1)
+
+       deallocate(site_name_ad)
+       deallocate(climfile_ad)
+       deallocate(sitefile_ad)
+       deallocate(manfile_ad)
+       deallocate(treefile_ad)
+       deallocate(depofile_ad)
+       deallocate(redfile_ad)
+       deallocate(litfile_ad)
+	   if (allocated(wpmfile_ad)) deallocate(wpmfile_ad)
+
+      close(nowunit1)
+
+ CASE (9, 10)
+
+ ! read once only per climate station 
+    allocate(sitefile(site_nr))
+    allocate(climfile(site_nr))
+    allocate(treefile(site_nr))
+    allocate(manfile(site_nr))
+    allocate(standid(site_nr))
+    allocate(latitude(site_nr))
+    allocate(site_name(site_nr))
+    allocate(RedN_list(15, site_nr))
+    RedN_list = -99.9
+
+ ! read once only
+    allocate(specfile(1))
+    allocate(depofile(1))
+    allocate(redfile(1))
+    allocate(litfile(1))
+    allocate(valfile(1))
+
+      read(nowunit,'(A)',iostat=ios) specfile(1)
+      read(nowunit,'(A)') site_name(1)
+      read(nowunit,'(A)') treefile(1)
+      read(nowunit,'(A)') manfile(1)
+      read(nowunit,'(A)') siteall       ! control file xxx.con
+      read(nowunit,'(A)') climall       ! climate station with coordiantes
+      read(nowunit,'(A)') pathdir1      ! path of climate scenarios
+      read(nowunit,'(A)') pathdir2      ! path of soil file xxx.sop or name of total soil file (flag_multi=8)
+      read(nowunit,'(A)') climszen      ! labeling climate scenarios
+      read(nowunit,'(A)') text          ! degree of climate scenarios
+      read(nowunit,*) nrreal            ! amount of realisations/implementations
+
+      if (.not.flag_mult8910) print *, ' >>>foresee message: Input of filenames completed'
+
+      depofile(1)  = 'dummy.dep'
+      redfile(1)   = 'dummy.red'
+      litfile(1)   = 'dummy.lit'
+      site_name    = site_name(1)
+      site_name1   = site_name(1)
+
+      ilen   = len(trim(text))
+      text   = adjustl(text)
+      nrclim = 0                     
+      do while (ilen .gt. 0)
+          nrclim = nrclim + 1
+          ihelp = scan(text, ' ')
+          typeclim(nrclim) = adjustl(text(1:ihelp-1))
+          text = adjustl(text(ihelp:))
+          ilen = len(trim(text))
+      enddo
+  
+  ! processing of nrreal realisations/implementations of climate scenarios
+      site_anz = nrreal * nrclim * site_nr
+      allocate(climszenfile(site_nr, nrclim, nrreal))
+
+!  define name of output variables
+      nvar = 1
+      read(nowunit,*) outvar(nvar)
+      do while (trim(outvar(nvar)) .ne. 'end')
+         nvar = nvar + 1
+         read(nowunit,*) outvar(nvar)
+      enddo
+
+      if (nvar .gt. 1)  then
+          allocate(output_var(nvar-1,1,0:year))
+          allocate(output_unit(nvar-1))
+          allocate(climszenres(nvar-1,site_nr,nrclim,nrreal))
+          output_unit     = -99
+          output_unit_all = -99
+
+          helpy = 0
+          helpi = 0
+          helpw = 0
+          do i = 1, nvar-1
+
+             select case (trim(outvar(i)))
+
+              case ('AET_year','AETyear','aetyear','aet_year','cwb_year','cwbyear','PET_year','PETyear','petyear','pet_year', &
+                    'GPP_year','GPPyear','gppyear','gpp_year','NEP_year','NEPyear','nepyear','nep_year','NPP_year','NPPyear','nppyear','npp_year', &
+                    'perc_year','percyear','temp_year','tempyear','prec_year','precyear', 'resps_year','respsyear','TER_year','TERyear','ter_year','teryear')  
+                flag_cum = 1
+                helpy    = helpy + 1
+                output_var(i,1,0) = 1.*helpy              ! field numbre of yearly values
+
+              case ('AET_mon','AETmon','aetmon','aet_mon','cwb_mon','cwbmon','PET_mon','PETmon','petmon','pet_mon', &
+                    'GPP_mon','GPPmon','gppmon','gpp_mon','NEP_mon','NEPmon','nepmon','nep_mon','NPP_mon','NPPmon','nppmon','npp_mon', &
+                    'perc_mon','percmon','temp_mon','tempmon','prec_mon','precmon', 'resps_mon','respsmon','TER_mon','TERmon','ter_mon','termon')  
+                flag_cum = 1
+                helpi    = helpi + 1
+                output_var(i,1,0) = 1.*helpi              ! field numbre of monthly values
+
+              case ('AET_week','AETweek','aetweek','aet_week','cwb_week','cwbweek','PET_week','PETweek','petweek','pet_week', &
+                    'GPP_week','GPPweek','gppweek','gpp_week','NEP_week','NEPweek','nepweek','nep_week','NPP_week','NPPweek','nppweek','npp_week', &
+                    'perc_week','percweek','temp_week','tempweek','prec_week','precweek', 'resps_week','respsweek', 'TER_week','TERweek','ter_week','terweek')  
+                flag_cum = 1
+                helpw    = helpw + 1
+                output_var(i,1,0) = 1.*helpw              ! field numbre of weekly values
+
+             end select   ! outvar
+      
+          enddo
+          if (helpy .gt. 0) then
+             allocate(climszenyear(helpy,site_nr,nrclim,nrreal,year))
+          endif
+          if (helpi .gt. 0) then
+             allocate(climszenmon(helpi,site_nr,nrclim,nrreal,12))
+             allocate(output_varm(helpi,1,year,12))
+          endif
+          if (helpw .gt. 0) then
+             allocate(climszenweek(helpw,site_nr,nrclim,nrreal,52))
+             allocate(output_varw(helpw,1,year,52))
+          endif
+      endif   ! nvar
+
+      call errorfile(simfile, ios, nowunit)
+
+!  reading file with description of climate stations used
+      allocate(climnum(3000))
+      allocate(clim_long(3000))
+      allocate(clim_lat(3000))
+      allocate(clim_height(3000))
+      allocate(clim_nam(3000))
+
+      call testfile(climall,ex)
+      if (ex .eqv. .false.) return
+      nowunit1 = getunit()
+      ios1 = 0
+      open(nowunit1,file=climall,iostat=ios,status='old',action='read')
+      k=1
+      do
+         READ(nowunit1,'(A)',iostat=ios1) a
+         IF (a .ne. '!') exit
+
+      end do
+      backspace nowunit1
+
+      do
+          read(nowunit1,*,iostat=ios1) climnum(k), clim_long(k),clim_lat(k), clim_height(k)
+          if(ios1 .lt. 0) exit
+          k = k+1
+      end do
+      anzclim = k-1
+      ios1 = 0
+
+      call errorfile(climall, ios1, nowunit1)
+      
+! reading control file with site-id, climate-id, soil-id, gwtabe-id
+
+      call testfile(siteall,ex)
+      if (ex .eqv. .false.) return
+      nowunit1 = getunit()
+      open(nowunit1,file=siteall,iostat=ios1,status='old',action='read')
+      do
+         READ(nowunit1,'(A)',iostat=ios1) a
+         IF (a .ne. '!') exit
+      end do
+      backspace nowunit1
+
+      do i=1,site_nr
+         call readcon(i, nowunit1)  
+         
+         sitefile(i) = trim(pathdir2)
+		 if(i.gt.1) treefile(i)= treefile(1)
+		 if(i.gt.1) manfile(i) = manfile(1)
+         k = 1
+         do while (clim_id(i).ne.climnum(k)) 
+             k = k + 1
+             if (k .gt. anzclim) then
+                write (unit_err,*)
+                write (unit_err,*) ' >>>foresee message: Climate ID ', trim(clim_id(i)), ' not in file ',trim(climall)
+                write (unit_err,*) '                     Site number ',sitenum(i)
+                write (*,*)
+                write (*,*) ' >>>foresee message: Climate ID ', trim(clim_id(i)), ' not in file ',trim(climall)
+                write (*,*) '                     Site number ',sitenum(i)
+                print *,' Program will stop!'
+                flag_end = 4
+                return
+             endif
+         enddo
+      latitude(i) = clim_lat(k)
+         do l = 1, nrclim
+            do j = 1, nrreal
+                write (helpsim,'(I5)') j
+                read (helpsim,*) text4
+			    select case (flag_multi)
+			    case (9)
+			        climszenfile(i,l,j) = trim(pathdir1)//trim(typeclim(l))//'/real_'//trim(text4)//'/'//trim(clim_id(i))//trim(climszen)//'.dat'
+			    case (10)
+			        if (j .lt. 10) then
+			            text4 = '00'//text4
+			        else if (j .lt. 100)   then
+			            text4 = '0'//text4
+			        endif
+			        climszenfile(i,l,j) = trim(pathdir1)//'/q'//trim(text4)//'/'//trim(clim_id(i))//trim(climszen)//'.dat'
+                end select
+            enddo   !j
+         end do   !l
+      enddo
+
+      if ((.not.flag_mult8910) .and. (ios1 .lt. 0)) print *, 'no information for site number ', i
+      call errorfile(siteall, ios1, nowunit1)
+
+      deallocate(climnum)
+      deallocate(clim_long)
+      deallocate(clim_lat)
+      deallocate(clim_height)
+      deallocate(clim_nam)
+
+      close(nowunit1)
+
+ END SELECT
+
+    jpar = 0  ! reset jpar for restore
+
+    if(flag_multi .eq. 2)then
+              read (nowunit,*) step_sum_T,n_T_downsteps,n_T_upsteps
+              read (nowunit,*) step_fac_P,n_P_downsteps,n_P_upsteps
+
+              site_nr = (1+n_T_downsteps+n_T_upsteps) * (1+n_P_downsteps+n_P_upsteps)
+		      repeat_number = site_nr
+
+              tspec   = specfile(1)
+              tname   = site_name(1)
+              tclim   = climfile(1)
+              tsite   = sitefile(1)
+              tval    = valfile(1)
+              ttree   = treefile(1)
+              tman    = manfile(1)
+              tdepo   = depofile(1)
+              tred    = redfile(1)
+              tlit    = litfile(1)
+              istand  = standid(1)
+              tsoilid = soilid(1)
+
+              deallocate (specfile)
+              deallocate (site_name)
+              deallocate (climfile)
+              deallocate (clim_id)
+              deallocate (sitefile)
+              deallocate (valfile)
+              deallocate (treefile)
+              deallocate (manfile)
+              deallocate (depofile)
+              deallocate (redfile)
+              deallocate (litfile)
+			  deallocate (wpmfile)
+              deallocate (standid)
+              deallocate (soilid)
+              allocate (specfile(site_nr))
+              allocate (site_name(site_nr))
+              allocate (climfile(site_nr))
+              allocate (clim_id(site_nr))
+              allocate (sitefile(site_nr))
+              allocate (valfile(site_nr))
+              allocate (treefile(site_nr))
+              allocate (manfile(site_nr))
+              allocate (depofile(site_nr))
+              allocate (standid(site_nr))
+              allocate (soilid(site_nr))
+              allocate (redfile(site_nr))
+              allocate (litfile(site_nr))
+			  allocate (wpmfile(site_nr))
+
+              specfile  = tspec
+              site_name = tname
+              climfile  = tclim
+              sitefile  = tsite
+              valfile   = tval
+              treefile  = ttree
+              manfile   = tman
+              depofile  = tdepo
+              redfile   = tred
+              litfile   = tlit
+              standid   = istand
+              soilid    = tsoilid
+
+       call errorfile(simfile, ios, nowunit)
+
+    endif   ! flag_multi = 2
+close(nowunit)
+
+END subroutine readsim
+
+!**************************************************************
+
+SUBROUTINE allofile
+
+use data_simul
+
+implicit none
+
+    allocate(site_name(site_nr))
+    allocate(climfile(repeat_number))
+    allocate(sitefile(site_nr))
+    allocate(valfile(site_nr))
+    allocate(treefile(repeat_number))
+    allocate(standid(repeat_number))
+    allocate(manfile(repeat_number))
+    allocate(depofile(repeat_number))
+    allocate(redfile(repeat_number))
+    allocate(litfile(repeat_number))
+	allocate(wpmfile(repeat_number))
+    allocate(specfile(repeat_number))
+
+end subroutine allofile
+
+!**************************************************************
+
+SUBROUTINE readcon (ii, unitnum)
+
+use data_depo
+use data_out
+use data_par
+use data_simul
+use data_site
+
+implicit none
+
+integer ii, ihelp, unitnum, ios1, ilen, helpi
+character(150):: text  
+character(10):: helpsim, text4 
+
+         read(unitnum,'(A)',iostat=ios1) text     
+        ! text disassemble
+        ! sitenum
+          ilen  = len(trim(text))
+          text  = adjustl(text)
+          ihelp = verify(text, charset)      
+          text4 = adjustl(text(1:ihelp-1))
+          sitenum(ii) = text4
+          text  = adjustl(text(ihelp+1:))
+          ilen  = len(trim(text))
+          ihelp = scan(text, charset)
+          text  = text(ihelp:)
+          ihelp = verify(text, charset)
+          clim_id(ii) = adjustl(text(1:ihelp-1))
+          text  = adjustl(text(ihelp+1:))
+          ilen  = len(trim(text))
+          ihelp = scan(text, charset)
+          text  = text(ihelp:)
+          ihelp = verify(text, charset)
+          soilid(ii) = adjustl(text(1:ihelp-1))
+        ! gwtable
+          text  = adjustl(text(ihelp+1:))
+          ilen  = len(trim(text))
+          ihelp = scan(text, charset)
+          text  = text(ihelp:)
+          ihelp = verify(text, charset)
+          text4 = adjustl(text(1:ihelp-1))
+          write (helpsim,'(A)') text4
+          read (helpsim,*) gwtable(ii)
+        ! standid 
+          text  = adjustl(text(ihelp+1:))
+          ilen  = len(trim(text))
+          ihelp = scan(text, charset)
+          text  = text(ihelp:)
+          ihelp = verify(text, charset)
+          text4 = adjustl(text(1:ihelp-1))
+          standid(ii) = text4
+        ! deposition
+          text = adjustl(text(ihelp+1:))
+          ilen  = len(trim(text))
+          if (ilen .gt. 0) then
+              text  = adjustl(text)
+              ihelp = verify(text, charset)
+              text4 = adjustl(text(1:ihelp-1))
+              write (helpsim,'(A)') text4
+              read (helpsim,*) NOdep(ii)     ! hand over in readdepo as concentration
+              text  = adjustl(text(ihelp+1:))
+              ilen  = len(trim(text))
+              ihelp = scan(text, charset)
+              text  = text(ihelp:)
+              ihelp = verify(text, charset)
+              text4 = adjustl(text(1:ihelp-1))
+              write (helpsim,'(A)') text4
+              read (helpsim,*) NHdep(ii)     ! hand over in readdepo as concentration
+            ! RedN 
+              text = adjustl(text(ihelp+1:))
+              ilen = len(trim(text))
+              do while (ilen .gt. 0)
+                  ihelp = verify(text, charset)
+                  text4 = adjustl(text(1:ihelp-1))
+                  write (helpsim,'(A)') text4
+                  read (helpsim,*) helpi
+                  text = adjustl(text(ihelp+1:))
+                  ihelp = verify(text, charset)
+                  text4 = adjustl(text(1:ihelp-1))
+                  write (helpsim,'(A)') text4
+                  read (helpsim,*) RedN_list(helpi, ii)
+                  text = adjustl(text(ihelp+1:))
+                  ilen = len(trim(text))
+              enddo
+          else
+              NOdep(ii) = 0.
+              NHdep(ii) = 0.
+          endif
+
+End SUBROUTINE readcon

source_code/version_2.3_windows/root.f

+!*****************************************************************!
+!*                                                               *!
+!*              4C (FORESEE) Simulation Model                    *!
+!*                                                               *!
+!*                                                               *!
+!*                    Subroutines for:                           *!
+!*                   Root distribution                           *!
+!*                                                               *!
+!* - ROOT_DISTR                                                  *!
+!* - ROOT_EFF                                                    *!
+!* - ROOT_DEPTH                                                  *!
+!* - ROOT_INI                                                    *!
+!* - DEALLOC_ROOT                                                *!
+!* - ROOTC_NEW (nicht benutzt wegen Problemen bei Verkettung)    *!
+!* - CR_DEPTH                                                    *!
+!*                                                               *!
+!*                  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 root_distr
+
+! Calculation of root distribution for each cohorte
+
+use data_simul
+use data_soil
+use data_stand
+use data_par
+use data_species
+
+implicit none
+
+integer specn     ! species type (number)
+integer i, j, nj, k, jlay
+integer nr        ! aux. var. for nroot (rooting depth)
+integer rkind     ! kind of calculation of root depth
+real    frtrel_1, frtrel_j  ! rel fine root fraction of previous layer
+real    frtrel_s  ! Sum of fine root fractions
+real    radius    ! radius of cylyndric space created by roots of the root length density 
+
+real beta         ! base of power
+real help
+real alpha, b	! Parameters for Arora function
+real troot2				! theoretical root biomass of population (coarse and fine roots) only for Arora funktion spereated according to cohorts [kg/m]
+real :: part_coef=0.0	! Verteilungskoeffizient um Verhltnis zwischen fr_loss und redis zu bestimmen
+real, dimension (1:nlay) :: fr_loss1, valspace, frtrelcoh	    !auxiliary vectors
+
+rkind = rdepth_kind
+
+if ((anz_tree + anz_sveg) .eq. 0) return
+select case (flag_wurz)
+
+case (0)
+   root_fr = 0.
+   zeig => pt%first
+   do while (associated(zeig))
+
+      call root_depth (rkind, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
+      zeig%coh%nroot = nr
+      zeig%coh%frtrel = thick/depth(nr)
+      specn = zeig%coh%species
+
+	  do j = 1, nr
+         root_fr(j) = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+      enddo
+      do j = nr+1, nlay
+         zeig%coh%frtrel(j) = 0.
+      enddo
+      zeig%coh%rooteff = 0.   ! zero after use
+      zeig => zeig%next
+   enddo
+
+case (1)    ! Funktion
+   root_fr = 0.
+   zeig => pt%first
+   do while (associated(zeig))
+      call root_depth (rkind, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot) !nderung MG: bergabe von Grob und Feinwurzelmasse an root_depth
+      zeig%coh%nroot = nr
+      specn = zeig%coh%species
+      if (specn .eq. 2 .or. specn .eq. 3) then
+         beta = 0.976
+      else
+         beta = 0.966
+      endif
+	  frtrel_1        = 1.
+      zeig%coh%frtrel = 0.
+      do j=1,nr
+         frtrel_j           = beta ** depth(j)
+         zeig%coh%frtrel(j) = frtrel_1 - frtrel_j
+         frtrel_1           = frtrel_j
+      enddo
+      frtrel_s = SUM(zeig%coh%frtrel)
+      frtrel_s = 1./frtrel_s
+      do j=1,nr
+ !     scaling of root distribution
+         zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+         root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+      enddo
+
+      zeig%coh%rooteff = 0.   ! zero after use
+      zeig => zeig%next
+   enddo
+
+case (2)    ! read/use default distribution; not changed 
+
+   root_fr = 0.
+   zeig => pt%first
+   do while (associated(zeig))
+      if (zeig%coh%frtrel(1) .gt. 0.) then
+          do j = 1,nroot_max
+            root_fr(j) = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+          enddo
+      else
+           root_fr = 0.
+              call root_depth (1, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot) 
+              zeig%coh%nroot = nr
+              specn = zeig%coh%species
+              if (specn .eq. 2 .or. specn .eq. 3) then
+                 beta = 0.98
+              else
+                 beta = 0.967
+              endif
+
+	          frtrel_1        = 1.
+              zeig%coh%frtrel = 0.
+              do j=1,nr
+                 frtrel_j           = beta ** depth(j)
+                 zeig%coh%frtrel(j) = frtrel_1 - frtrel_j
+                 frtrel_1           = frtrel_j
+              enddo
+              frtrel_s = SUM(zeig%coh%frtrel)
+              frtrel_s = 1./frtrel_s
+              do j=1,nr
+                 zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+                 root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+              enddo
+
+      endif
+
+      zeig%coh%rooteff = 0.   ! zero after use
+      zeig => zeig%next
+   enddo
+
+case (3)
+root_fr = 0.
+rkind=5
+
+zeig => pt%first
+   do while (associated(zeig))
+	  call root_depth (rkind, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot) !nderung MG: bergabe von Grob und Feinwurzelmasse an root_depth
+      zeig%coh%nroot = nr
+      specn = zeig%coh%species      
+	  alpha=0.7
+	  if (specn .eq. 2 .or. specn .eq. 3 .or. specn .eq. 6 .or. specn .eq. 7) then
+         b = 7.95
+      else
+         b = 10.91
+      endif
+ 
+	  frtrel_1        = 1.
+      zeig%coh%frtrel = 0.
+      do j=1,nr
+!        root distribution (Arora et al., 2003)
+         frtrel_j           = exp((-b/troot2**alpha)*(depth(j)/100))
+         zeig%coh%frtrel(j) = frtrel_1 - frtrel_j
+         frtrel_1           = frtrel_j		 
+	  enddo
+      frtrel_s = SUM(zeig%coh%frtrel)
+      frtrel_s = 1./frtrel_s
+      do j=1,nr
+ !     scaling of root distribution
+         zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+         root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA		 
+	  enddo
+
+      zeig%coh%rooteff = 0.   ! zero after use
+      zeig => zeig%next
+   enddo
+
+case(4)		! TRAP-model Rasse et al. (2001)
+root_fr = 0.
+rkind   = 6
+fr_loss1= 0
+k       = 0
+
+zeig => pt%first
+   do while (associated(zeig))
+	k=k+1
+	zeig%coh%x_rdpt=gr_depth(k)  
+    specn = zeig%coh%species
+  if (specn .eq. 12) then
+  continue
+  endif
+	call root_depth (rkind, specn, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
+	zeig%coh%nroot = nr
+	
+	frtrel_1        = 1.
+    zeig%coh%frtrel = 0.
+      do j=1,nr
+         if (j .eq. 1) then
+		  zeig%coh%frtrel(j) = (zeig%coh%x_rdpt**3-(zeig%coh%x_rdpt-depth(j))**3)/zeig%coh%x_rdpt**3
+		   elseif (j .eq. nr) then
+		    zeig%coh%frtrel(j)= frtrel_1 
+		     else
+			  zeig%coh%frtrel(j) = ((zeig%coh%x_rdpt-depth(j-1))**3-((zeig%coh%x_rdpt-depth(j))**3))/zeig%coh%x_rdpt**3
+         endif
+		 frtrel_1 = frtrel_1-zeig%coh%frtrel(j)		 
+	  enddo
+	frtrel_s = SUM(zeig%coh%frtrel)
+    frtrel_s = 1./frtrel_s
+	zeig%coh%frtrel = zeig%coh%frtrel * frtrel_s
+	
+	fr_loss1 = zeig%coh%frtrel
+	fr_loss  = zeig%coh%frtrel*svar(specn)%Smean(1:nlay)
+	fr_loss  = part_coef*(fr_loss1-fr_loss)
+	redis    = zeig%coh%frtrel*svar(specn)%Smean(1:nlay)
+	redis    = part_coef*(fr_loss1-redis)
+
+	do j=1,nr
+ !     scaling of root distribution
+       if (sum(svar(specn)%Smean(1:nr)) .lt. 0.0001) then
+       zeig%coh%frtrel(j) = 0.
+       else
+       zeig%coh%frtrel(j) = zeig%coh%frtrel(j)*svar(specn)%Smean(j)+(sum(redis)*svar(specn)%Smean(j)/sum(svar(specn)%Smean(1:nr)))
+       endif
+	enddo
+
+	  frtrel_s = SUM(zeig%coh%frtrel)
+      if (frtrel_s .lt. 1.E-6) then
+          do j=1,nr
+             zeig%coh%frtrel(j) = 0
+	      enddo
+      else
+          frtrel_s = 1./frtrel_s
+          do j=1,nr
+   !     scaling of root distribution
+             zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+             root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA		 
+	      enddo
+      endif
+			
+	zeig%coh%rooteff = 0. 
+	zeig => zeig%next
+   enddo
+
+case(5)
+root_fr = 0.
+rkind=5
+   zeig => pt%first
+   do while (associated(zeig))
+
+
+      call root_depth (rkind, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot) !nderung MG: bergabe von Grob und Feinwurzelmasse an root_depth
+      zeig%coh%nroot = nr
+      specn = zeig%coh%species
+      if (specn .eq. 2 .or. specn .eq. 3) then
+         beta = 0.98
+      else
+         beta = 0.967
+      endif
+
+ 
+	  frtrel_1        = 1.
+      zeig%coh%frtrel = 0.
+      do j=1,nr
+!       root distribution (Jackson et al., 1996): beta ** depth
+         frtrel_j           = beta ** depth(j)
+         zeig%coh%frtrel(j) = frtrel_1 - frtrel_j
+         frtrel_1           = frtrel_j
+      enddo
+      frtrel_s = SUM(zeig%coh%frtrel)
+      frtrel_s = 1./frtrel_s
+      do j=1,nr
+ !     scaling of root distribution
+         zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+         root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+      enddo
+
+      zeig%coh%rooteff = 0.   ! zero after use
+      zeig => zeig%next
+   enddo
+
+case(6)
+
+root_fr = 0.
+rkind=7
+   zeig => pt%first
+   k=1
+   do while (associated(zeig))
+	  
+	  zeig%coh%x_rdpt=gr_depth(k)
+      call root_depth (rkind, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt,zeig%coh%nroot) !nderung MG: bergabe von Grob und Feinwurzelmasse an root_depth
+      
+	  if (time .le. 1) then
+	   root_lay(k)=nr
+	  else
+	   root_lay(k)=root_lay(k)+nr
+	  endif
+	  
+	  if (root_lay(k) .gt. nroot_max) root_lay(k) = nroot_max
+	  
+	  zeig%coh%nroot=root_lay(k)
+	  nr=root_lay(k)
+	  
+      specn = zeig%coh%species
+      if (specn .eq. 2 .or. specn .eq. 3) then
+         beta = 0.98
+      else
+         beta = 0.967
+      endif
+
+ 
+	  frtrel_1        = 1.
+      zeig%coh%frtrel = 0.
+      do j=1,nr
+!       root distribution (Jackson et al., 1996): beta ** depth
+         frtrel_j           = beta ** depth(j)
+         zeig%coh%frtrel(j) = frtrel_1 - frtrel_j
+         frtrel_1           = frtrel_j
+      enddo
+      frtrel_s = SUM(zeig%coh%frtrel)
+      frtrel_s = 1./frtrel_s
+      do j=1,nr
+ !     scaling of root distribution
+         zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+         root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+      enddo
+
+      zeig%coh%rooteff = 0.   ! zero after use
+	  k=k+1
+      zeig => zeig%next
+   enddo
+
+case (7)    ! Funktion nach Jackson (1996) mit fester Tiefe
+   root_fr = 0.
+   nr = nroot_max
+   zeig => pt%first
+   do while (associated(zeig))
+
+      zeig%coh%nroot = nroot_max
+      specn = zeig%coh%species
+      if (specn .eq. 2 .or. specn .eq. 3) then
+         beta = 0.98
+      else
+         beta = 0.967
+      endif
+
+	  frtrel_1        = 1.
+      zeig%coh%frtrel = 0.
+      do j=1,nr
+!       root distribution (Jackson et al., 1996): beta ** depth
+         frtrel_j           = beta ** depth(j)
+         zeig%coh%frtrel(j) = frtrel_1 - frtrel_j
+         frtrel_1           = frtrel_j
+      enddo
+      frtrel_s = SUM(zeig%coh%frtrel)
+      frtrel_s = 1./frtrel_s
+      do j=1,nr
+ !     scaling of root distribution
+         zeig%coh%frtrel(j) = zeig%coh%frtrel(j) * frtrel_s
+         root_fr(j)         = root_fr(j) + zeig%coh%frtrel(j) * zeig%coh%ntreeA
+      enddo
+
+      zeig%coh%rooteff = 0.   ! zero after use
+
+      zeig => zeig%next
+   enddo
+
+end select
+
+root_fr = root_fr / (anz_tree + anz_sveg)   ! normieren
+zeig => pt%first
+do while (associated(zeig))
+    help = zeig%coh%x_frt * zeig%coh%ntreea
+    do jlay = 1, nroot_max
+        if (root_fr(jlay) .gt. zero) then 
+            zeig%coh%frtrelc(jlay) = zeig%coh%frtrel(jlay) * help / (root_fr(jlay) * totfrt_p)  ! mass of root part of total cohort in a layer
+        else
+            zeig%coh%frtrelc(jlay) = 0.
+        endif
+    enddo
+    zeig => zeig%next
+enddo
+
+if (flag_wred .eq. 9) then
+
+   !Calculation of root length density
+   zeig => pt%first
+   do while (associated(zeig))
+        if (specn .le. nspec_tree) then
+            radius = (zeig%coh%diam/6.)*100.     ! formula bhd [cm]/6 yield radius in [m] so *100 (aus Wagner 2005)
+            valspace = pi * radius**2 * thick
+        else
+            valspace = kpatchsize * 100*100 * thick
+        endif                                                  !circular cylinder 
+    
+        frtrelcoh = zeig%coh%frtrel * zeig%coh%x_frt * zeig%coh%ntreea          
+    
+        if (zeig%coh%ntreea .gt. 0 .AND. minval(valspace(1:nr)) .gt. 0.) then
+            zeig%coh%rld = (frtrelcoh*1000*spar(specn)%spec_rl*100)/(valspace* zeig%coh%ntreea)   !in cm root length /cm3 volume
+	    else
+            zeig%coh%rld = -99
+        endif
+       
+      zeig => zeig%next
+   enddo
+endif
+
+if (allocated(wat_root)) wat_root=0.
+
+END	subroutine root_distr
+
+!**************************************************************
+
+SUBROUTINE root_eff
+
+! Calculation of root efficiency in dependence of water and N uptake
+use data_soil
+use data_soil_cn
+use data_stand
+
+implicit none
+
+integer i,j
+integer nr    ! layer number of root depth
+real hroot    ! root depth
+real fdc      ! discounting function describing transport resistance
+real gw, gN   ! accounting functions of water resp. N uptake
+real glimit   ! limitation constant for use of rooting layer
+
+glimit = 0.   ! min. assumption
+
+   i = 1
+   zeig => pt%first
+   do while (associated(zeig))
+      nr = zeig%coh%nroot
+      do j = 1,nr
+         fdc   = 50./depth(j)
+         if (zeig%coh%supply .gt. 1e-06) then
+            gw = xwatupt(i,j)/zeig%coh%supply
+            gw = gw / thick(j)
+         else
+		    gw = 0.
+		 endif
+
+         gw = xwatupt(i,j)     
+         zeig%coh%rooteff(j) = zeig%coh%rooteff(j) + gw
+      enddo
+      zeig%coh%watuptc = zeig%coh%watuptc + zeig%coh%supply
+      i = i + 1
+      zeig => zeig%next
+   enddo
+
+END	subroutine root_eff
+
+!**************************************************************
+
+SUBROUTINE root_depth(rkind, specn, agec, heightc, froot, croot, nr, troot2, crdepth, nrooth)
+
+use data_simul
+use data_soil
+use data_soil_cn
+use data_stand
+
+implicit none
+
+! input:
+integer rkind ! kind of calculation of root depth
+integer specn ! species number
+integer agec  ! tree age
+integer nrooth ! for case(7)
+
+real heightc, froot, croot					! tree height of cohort, fine and coarse root mass[kg]/ tree 
+real troot, troot1,troot2, troot_stand		! total root mass 1./tree 2./ha according to  cohorts 3. /m according to cohorts Kohorten 4./ha of 4C
+real :: wat_demand							! query whether one cohort was unable to cover water demand with the from root penetrated soil layer 
+real rootingdepth, crdepth				    ! rooting depth nach Arora function in [m]
+real alpha, b								! parameter for Arorafunction	
+! output:
+integer nr                                  ! last root layer
+
+integer i,j
+real hc, wtiefe
+real, dimension(4,3):: rdepth  ! effective rooting depth depending on tree age and soil texture
+! data from Raissi et al. (2001)
+data rdepth /85, 130, 175, 95, 140, 185, 135, 180, 225, 90, 110, 135/
+
+select case (rkind)
+
+case (1)
+ ! nroot depending on tree height and soil profile depth
+   nr = 1
+   do j=1,nlay
+      if (heightc .ge. depth(j)) nr = j
+   enddo
+   if (nr .gt. nroot_max) nr = nroot_max
+   crdepth = depth(nr)
+
+case (2)
+ ! fixed nroot for all adult cohorts
+   if (agec .lt. 10) then
+     nr = 1
+     wtiefe=depth(nroot_max)/(1+exp(1.5-0.55*real(agec)))		! logicla function to determin root depth [cm] until age 10	
+		do j=1,nlay
+			if (wtiefe .ge. depth(j)) nr = j
+		enddo
+     if (nr .gt. nroot_max) nr = nroot_max
+   else
+     nr = nroot_max
+   endif
+    crdepth = depth(nr)
+
+case (3)
+ ! nroot depending on root efficiency
+   nr = nlay  
+   crdepth = depth(nr)
+
+case (4)
+ ! nroot depending on soil texture and age
+   if (agec .lt. 15) then
+      i = 1
+   else if (agec .gt. 45) then
+      i = 3
+   else
+      i = 2
+   endif
+
+   nr = 1
+   if (heightc .gt. rdepth(s_typen,i)) then
+      hc = rdepth(s_typen,i)
+   else
+      hc = heightc
+   endif
+   do j=1,nlay
+      if (hc .ge. depth(j)) nr = j
+   enddo
+   if (nr .gt. nroot_max) nr = nroot_max
+
+case (5)
+	alpha=0.7
+	if (specn .eq. 2 .or. specn .eq. 3 .or. specn .eq. 6 .or. specn .eq. 7) then
+         b = 7.95
+      else
+         b = 10.91
+      endif
+	
+	troot=froot+croot
+	troot1=troot*anz_tree_ha        ! total root biomass per ha if population of a cohort is soley comprised of trees
+	troot_stand=totfrt+totcrt		! total root biomass per ha calculated by 4C
+	troot2=troot1/10000				! conversion to m
+	rootingdepth=(3*troot2**alpha)/b	!Arora function
+		nr = 1
+		do j=1,nlay
+			if (rootingdepth*100 .ge. depth(j)) nr = j
+		enddo
+   if (nr .gt. nroot_max) nr = nroot_max
+    crdepth = depth(nr)
+
+case (6)	!Calculation in soil.f in cr_depth
+
+     if (crdepth .eq.0) then
+	    ! nroot depending on soil texture and age
+       if (agec .lt. 15) then
+          i = 1
+       else if (agec .gt. 45) then
+          i = 3
+       else
+          i = 2
+       endif
+
+       nr = 1
+       if (heightc .gt. rdepth(s_typen,i)) then
+          crdepth = rdepth(s_typen,i)
+       else
+          crdepth = heightc
+       endif
+   
+     endif
+
+     do j=1,nlay
+	    if (depth(j) .le. crdepth) nr=j
+     enddo
+     if (nr .gt. nroot_max) nr = nroot_max
+
+case (7)	!further growth only if next layer bears water
+	wat_demand=maxval(wat_root)
+	if (time .le. 1) then
+	 crdepth=30.0
+	 do j=1,nlay
+	  if (depth(j) .le. 30.) nr=j
+	 enddo
+	else
+      if (wat_demand .gt. 0) then
+	   nr=1
+	  else
+	   nr=0
+	  endif
+	endif
+
+   if (nr .gt. nroot_max) nr = nroot_max
+    crdepth = depth(nr)
+end select
+
+if (crdepth < 0.) then
+continue
+endif
+
+END	subroutine root_depth
+
+!**************************************************************
+
+SUBROUTINE root_ini
+
+! Allocation and initialisation of root distribution
+
+use data_simul
+use data_soil
+use data_species
+use data_stand
+
+implicit none
+
+integer i, j, nj, rkind, hspec, ios
+integer unit_root
+integer nr        ! aux. var. for nroot (rooting depth)
+real    frtrel_j, frtrel_1  ! rel fine root fraction of previous layer
+real    frtrel_s  ! Sum of fine root fractions
+real hfrt, help, troot2
+real, allocatable, dimension(:,:):: hd,hr
+integer, allocatable, dimension(:):: nlspec
+character text
+character (150) file_root
+
+logical :: pruefer=.false.
+
+root_fr = 0.
+  if (wlam(3) .gt. 0.4) then
+     s_typen = 1     ! sand
+  else if (wlam(3) .le. 0.15) then
+     s_typen = 4     ! clay
+  else if (wlam(3) .gt. 0.25) then
+     s_typen = 3     ! silt
+  else
+     s_typen = 2     ! loam
+  endif
+
+  if (nroot_max .lt. 0) then
+     nroot_max = 1
+     rkind = 4
+  else
+     rkind = 2
+  endif
+  rdepth_kind = rkind
+
+select case (flag_wurz)
+
+case (0,1,5)
+   if (anz_tree .gt. 0 .or. (anz_tree.eq.0 .and. flag_sveg .eq.1)) call root_distr
+
+case (3,4,6)
+    !intercept the case that the ground vegetatuin is already initialised but no trees have been initialised so cohorts are not finalised
+    if (anz_tree.eq.0 .and. flag_sveg .eq.1) then
+	    if (.not. allocated(wat_root)) then
+	     allocate(wat_root(anz_coh))
+	     wat_root=0.
+	     allocate(root_lay(anz_coh))
+         root_lay=0
+         allocate(gr_depth(anz_coh))
+         gr_depth=0.
+	     Pruefer=.true.
+	    endif
+    else
+        if (Pruefer .OR. (.not. allocated(wat_root))) then
+	     if (Pruefer) deallocate(wat_root)
+	     allocate(wat_root(anz_coh))
+	     wat_root=0.
+	     if (Pruefer) deallocate(root_lay)
+	     allocate(root_lay(anz_coh))
+         root_lay=0
+         if (Pruefer) deallocate(gr_depth)
+	     allocate(gr_depth(anz_coh))
+         gr_depth=0.
+         Pruefer=.false.
+	    endif
+    endif
+    if (anz_tree .gt. 0 .or. (anz_tree.eq.0 .and. flag_sveg .eq.1)) call root_distr
+
+case (2)
+! read root distribution once in the beginning alone
+   write (*,*)
+   write (*,'(A)', advance='no') 'Define root distribution, name of input file: '
+   read (*,'(A)') file_root
+   unit_root = getunit()
+   open (unit_root, file=trim(file_root), status='unknown')
+   allocate (hd(0:40, 1:nspecies))
+   allocate (hr(0:40, 1:nspecies))
+   allocate (nlspec(nspecies))
+
+   do
+       read (unit_root,'(A)') text
+       if (text .ne. '!') then
+            backspace(unit_root);exit
+       endif
+   enddo
+
+   ios    = 0
+   hd     = 0.
+   hr     = 0.
+   nlspec = 0
+   do while (ios .ge. 0)
+       j = 1
+       read (unit_root, *, iostat=ios) hspec
+       if (ios .lt. 0) exit
+       read (unit_root, *, iostat=ios) hd(1,hspec), hr(1,hspec)
+       do while (hd(j,hspec) .ge. 0.)
+          nlspec(hspec) = j
+          j = j+1
+          read (unit_root, *, iostat=ios) hd(j,hspec), hr(j,hspec)
+       enddo
+       if (hd(j,hspec) .lt. depth(nlay)) hd(j,hspec) = depth(nlay)
+   enddo
+   close (unit_root)
+
+   zeig => pt%first
+   do while (associated(zeig))
+
+      ns = zeig%coh%species
+      zeig%coh%frtrel = 0.
+
+  ! rel. root distribution of cohorts to species allocated
+     if (nlspec(ns) .gt. 0) then
+
+         frtrel_j = 0.
+         hfrt = 0.
+         j= 1
+         do while (hd(j,ns) .lt. depth(1))
+            hfrt = hfrt + hr(j,ns)
+            j = j+1
+         enddo
+            frtrel_j           = hr(j,ns) / (hd(j,ns)-hd(j-1,ns))
+            hfrt               = hfrt + frtrel_j * (depth(1)-hd(j-1,ns))
+            zeig%coh%frtrel(1) = hfrt
+            nj = j
+
+         do i=2,nlay
+             hfrt = 0.
+             do j = nj,nlspec(ns)+1
+                if (hd(j,ns) .lt. depth(i)) then
+                    frtrel_j = hr(j,ns) / (hd(j,ns)-hd(j-1,ns))
+                    hfrt     = hfrt + frtrel_j * (hd(j,ns)-depth(i-1))
+                else
+                    if (depth(i-1) .gt. hd(j-1,ns)) then
+                        help = depth(i)-depth(i-1)
+                    else
+                        help = depth(i)-hd(j-1,ns)
+                    endif
+                    frtrel_j = hr(j,ns) / (hd(j,ns)-hd(j-1,ns))
+                    hfrt     = hfrt + frtrel_j * help
+                    nj = j
+                    exit
+                endif
+             enddo
+
+             zeig%coh%frtrel(i) = hfrt
+          enddo
+     else
+
+            continue
+     endif
+
+      frtrel_s  = SUM(zeig%coh%frtrel)
+      zeig%coh%rooteff   = 0.
+      zeig => zeig%next
+   enddo
+   rdepth_kind = 2
+end select
+
+END	subroutine root_ini
+
+!**************************************************************
+
+SUBROUTINE dealloc_root
+
+use data_simul
+use data_stand
+
+if (flag_wurz .eq. 1) then
+   zeig => pt%first
+   do while (associated(zeig))
+
+      deallocate (zeig%coh%frtrel)
+      deallocate (zeig%coh%rooteff)
+
+      zeig => zeig%next
+   enddo
+endif
+
+END	subroutine dealloc_root
+
+!**************************************************************
+
+SUBROUTINE rootc_new (zeig1)
+
+! root initialisation of a new cohort
+
+use data_stand
+use data_soil
+
+implicit none
+
+type(coh_obj), pointer   :: zeig1     ! pointer variable for cohorts
+real troot2
+integer j, nr
+
+        allocate (zeig1%coh%frtrel(nlay))
+        allocate (zeig1%coh%rooteff(nlay))
+        zeig1%coh%frtrel  = 0.      ! initialisation
+        call root_depth (1, zeig1%coh%species, zeig1%coh%x_age, zeig1%coh%height, zeig1%coh%x_frt, zeig1%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
+        zeig1%coh%nroot = nr
+        do j=1,nr
+           zeig1%coh%rooteff = 1.   ! assumption for the first use
+        enddo
+        do j=nr+1, nlay
+           zeig1%coh%rooteff = 0.   ! layers with no roots
+        enddo
+
+END	subroutine rootc_new
+
+!**************************************************************
+
+SUBROUTINE cr_depth 
+
+! Calculation of the rooting depth after Rasse et al. 2001 
+
+use data_soil
+use data_stand
+use data_simul
+use data_climate
+use data_species
+
+implicit none
+real :: vcr		! growth rate rootdepth [cm]
+integer :: j,k
+
+vcr=0.
+
+select case (flag_wurz)
+
+case(4,6)
+zeig => pt%first
+ k=1
+ do while (associated(zeig))
+  	do j=1,nlay
+	 if (zeig%coh%x_rdpt .lt. depth(j)) then
+	  if (zeig%coh%x_age .le. 100) then
+		if (j .eq. 1) then
+        vcr=spar(zeig%coh%species)%v_growth*((100-real(zeig%coh%x_age))/100)*svar(zeig%coh%species)%Rstress(j)
+	    zeig%coh%x_rdpt=zeig%coh%x_rdpt+(vcr/recs(time))
+		gr_depth(k)=zeig%coh%x_rdpt
+		exit
+	    else
+	    vcr=spar(zeig%coh%species)%v_growth*((100-real(zeig%coh%x_age))/100)*svar(zeig%coh%species)%Rstress(j)
+	    zeig%coh%x_rdpt=zeig%coh%x_rdpt+(vcr/recs(time))
+	    gr_depth(k)=zeig%coh%x_rdpt
+		exit
+	    endif
+	   endif
+	  endif
+	enddo
+  
+  if (zeig%coh%x_rdpt .gt. depth(nroot_max)) zeig%coh%x_rdpt = depth(nroot_max)
+  k=k+1	  
+  zeig => zeig%next
+ enddo
+end select
+
+END subroutine cr_depth
+
+!*******************************************************************************
+

source_code/version_2.3_windows/script1.rc

+// Generated by ResEdit 1.5.9
+// Copyright (C) 2006-2011
+// http://www.resedit.net
+
+#include <windows.h>
+#include <commctrl.h>
+#include <richedit.h>
+#include "4C_dialogs.h"
+
+
+
+
+//
+// Bitmap resources
+//
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDB_BITMAP1        BITMAP         "4c_logo_klein.bmp"
+
+
+
+//
+// Dialog resources
+//
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_coh_daily DIALOG 100, 0, 540, 213
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_DAILYCOH_BUTTON_OK, 380, 190, 50, 14
+    PUSHBUTTON      "Cancel", ID_CANCEL_OUTF, 439, 190, 50, 14
+    GROUPBOX        "", IDC_STATIC_year, 10, 10, 510, 170
+    AUTOCHECKBOX    "", IDC_CHECK_y1, 26, 30, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y2, 26, 45, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y3, 26, 60, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y4, 26, 75, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y5, 25, 90, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y6, 26, 105, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y7, 26, 120, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y8, 26, 135, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y9, 26, 150, 271, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y10, 26, 165, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y11, 295, 30, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y12, 295, 45, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y13, 295, 60, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y14, 295, 75, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y15, 295, 90, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y16, 295, 105, 217, 8
+    LTEXT           "  Choose daily cohort output files  ", IDC_STATIC_yfile, 25, 10, 107, 8, SS_LEFT
+    AUTOCHECKBOX    "", IDC_CHECK_y17, 295, 120, 217, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y18, 295, 135, 217, 8
+    PUSHBUTTON      "Select all", ID_DAILYCOH_BUTTON_SELECT, 321, 190, 50, 14
+    PUSHBUTTON      "Deselect all", ID_DAILYCOH_BUTTON_DESELECT, 262, 190, 50, 14
+    AUTOCHECKBOX    "", IDC_CHECK_y19, 295, 150, 217, 8
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_coh_yearly DIALOG 0, 0, 718, 318
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_YEARLYCOH_BUTTON_OK, 559, 290, 50, 14
+    PUSHBUTTON      "Cancel", ID_CANCEL_OUTF, 625, 290, 50, 14
+    GROUPBOX        "", IDC_STATIC_year, 15, 10, 683, 275
+    AUTOCHECKBOX    "", IDC_CHECK_y1, 21, 30, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y2, 20, 45, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y3, 21, 60, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y4, 21, 75, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y5, 21, 90, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y6, 21, 105, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y7, 21, 120, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y8, 21, 135, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y9, 21, 150, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y10, 21, 165, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y11, 21, 180, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y12, 21, 195, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y13, 21, 210, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y14, 21, 225, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y15, 21, 240, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y16, 21, 255, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y17, 21, 270, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y18, 244, 30, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y19, 244, 45, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y20, 244, 60, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y21, 244, 75, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y22, 244, 90, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y23, 244, 105, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y24, 244, 120, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y25, 244, 135, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y26, 244, 150, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y27, 244, 165, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y28, 244, 180, 222, 8
+    LTEXT           "  Choose yearly cohort output files  ", IDC_STATIC_yfile, 30, 10, 111, 8, SS_LEFT
+    AUTOCHECKBOX    "", IDC_CHECK_y29, 244, 195, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y30, 244, 210, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y31, 244, 225, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y32, 244, 240, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y33, 244, 255, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y34, 244, 270, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y35, 465, 30, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y36, 465, 45, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y37, 465, 60, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y38, 465, 75, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y39, 465, 90, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y40, 465, 105, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y41, 465, 120, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y42, 465, 135, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y43, 465, 150, 222, 8
+    PUSHBUTTON      "Select all", ID_YEARLYCOH_BUTTON_SELECT, 493, 290, 50, 14
+    PUSHBUTTON      "Deselect all", ID_YEARLYCOH_BUTTON_DESELECT, 427, 290, 50, 14
+    AUTOCHECKBOX    "", IDC_CHECK_y44, 465, 165, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y45, 465, 178, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y46, 465, 193, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y47, 465, 210, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y48, 465, 225, 222, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y49, 465, 240, 222, 13
+    AUTOCHECKBOX    "", IDC_CHECK_y50, 465, 255, 222, 13
+    AUTOCHECKBOX    "", IDC_CHECK_y51, 465, 270, 222, 13
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_ctr DIALOG 100, 0, 280, 349
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_CTR_BUTTON_OK, 147, 319, 50, 14
+    PUSHBUTTON      "Cancel", IDCANCEL, 207, 319, 50, 14
+    LTEXT           "Edit simulation control file", IDC_STATIC_simul, 15, 14, 91, 11, SS_LEFT
+    GROUPBOX        "Run option control", IDC_STATIC_runo, 15, 30, 241, 173
+    LTEXT           "Run option", IDC_STATIC_runv, 23, 47, 36, 8, SS_LEFT
+    LTEXT           "Number of runs", IDC_STATIC_runnr, 23, 68, 49, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_runv, 69, 47, 175, 98, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    EDITTEXT        IDC_EDIT_runnr, 162, 67, 40, 12, ES_AUTOHSCROLL
+    LTEXT           "Number of simulation years", IDC_STATIC_yearn, 23, 89, 85, 8, SS_LEFT
+    LTEXT           "Start year", IDC_STATIC_start, 23, 111, 31, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_yearn, 162, 89, 40, 14, ES_AUTOHSCROLL
+    EDITTEXT        IDC_EDIT_start, 162, 111, 40, 14, ES_AUTOHSCROLL
+    LTEXT           "Patch size [m]", IDC_STATIC_patch, 23, 134, 47, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_patch, 162, 134, 40, 14, ES_AUTOHSCROLL
+    LTEXT           "Thickness of foliage layers [cm]", IDC_STATIC_thickf, 23, 158, 100, 8, SS_LEFT
+    LTEXT           "Time step photosynthesis calculations [d]", IDC_STATIC_timeph, 23, 184, 130, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_thickf, 162, 157, 40, 14, ES_AUTOHSCROLL
+    EDITTEXT        IDC_EDIT_timeph, 162, 180, 40, 14, ES_AUTOHSCROLL
+    GROUPBOX        "Run flags control", IDC_STATIC_runo3, 15, 214, 240, 93
+    COMBOBOX        IDC_COMBO_runv3, 65, 235, 101, 98, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "Run number", IDC_STATIC_runv3, 23, 238, 40, 8, SS_LEFT
+    DEFPUSHBUTTON   "Set Flags", ID_CTR_BUTTON_FLAGS, 177, 281, 65, 14
+    DEFPUSHBUTTON   "Set Files", ID_CTR_BUTTON_FILES, 177, 258, 65, 14
+    DEFPUSHBUTTON   "Apply", ID_CTR_BUTTON_RUNNR, 213, 66, 30, 14
+    DEFPUSHBUTTON   "Change Output Id", ID_CTR_BUTTON_IDS, 178, 235, 64, 14
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_daily DIALOG 0, 0, 568, 221
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_DAILY_BUTTON_OK, 408, 195, 50, 14
+    PUSHBUTTON      "Cancel", ID_CANCEL_OUTF, 471, 195, 50, 14
+    GROUPBOX        "", IDC_STATIC_day, 11, 14, 534, 171
+    AUTOCHECKBOX    "", IDC_CHECK_y1, 37, 30, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y2, 37, 45, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y3, 37, 60, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y4, 37, 75, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y5, 37, 90, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y6, 37, 105, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y7, 37, 120, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y8, 37, 135, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y9, 37, 150, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y10, 37, 165, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y11, 284, 30, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y12, 284, 45, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y13, 284, 60, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y14, 284, 75, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y15, 284, 90, 240, 8
+    LTEXT           "  Choose daily output files  ", IDC_STATIC_yfile, 22, 15, 85, 8, SS_LEFT
+    AUTOCHECKBOX    "", IDC_CHECK_y16, 284, 105, 240, 8
+    PUSHBUTTON      "Select all", ID_DAILY_BUTTON_SELECT, 345, 195, 50, 14
+    PUSHBUTTON      "Deselect all", ID_DAILY_BUTTON_DESELECT, 282, 195, 50, 14
+    AUTOCHECKBOX    "", IDC_CHECK_y17, 284, 120, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y18, 284, 135, 240, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y19, 284, 150, 240, 8
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_default_dir DIALOG 10, 100, 320, 99
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_DEFAULT_DIR_BUTTON_OK, 138, 79, 50, 14
+    LTEXT           "Default input and output directories", IDC_STATIC_4C, 15, 4, 216, 9, SS_LEFT
+    LTEXT           "Input directory", IDC_STATIC_dirin, 16, 31, 51, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_DIR_IN, 81, 27, 210, 12, ES_AUTOHSCROLL | ES_MULTILINE
+    LTEXT           "Output directory", IDC_STATIC_dirout, 16, 56, 51, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_DIR_OUT, 81, 52, 210, 12, ES_AUTOHSCROLL | ES_MULTILINE
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_files DIALOG 100, 0, 310, 330
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_FILES_BUTTON_OK, 189, 307, 50, 14
+    PUSHBUTTON      "Cancel", ID_CANCEL_FILES, 251, 307, 50, 14
+    GROUPBOX        "Model input files", IDC_STATIC_files, 7, 10, 295, 209
+    LTEXT           "Species parameter file", IDC_STATIC_specpar, 16, 51, 70, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_specpar, 109, 53, 146, 12, ES_AUTOHSCROLL | ES_MULTILINE
+    PUSHBUTTON      "Browse", IDC_BUTTON_specpar, 259, 54, 35, 12, WS_GROUP
+    LTEXT           "Soil parameter file", IDC_STATIC_sop, 16, 69, 56, 8, SS_LEFT
+    LTEXT           "Soil initialisation file", IDC_STATIC_soi, 16, 87, 60, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_sop, 109, 70, 146, 12, ES_AUTOHSCROLL
+    EDITTEXT        IDC_EDIT_soi, 109, 87, 146, 12, ES_AUTOHSCROLL
+    LTEXT           "Stand initialisation file", IDC_STATIC_ini, 13, 238, 68, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_ini, 107, 235, 145, 14, ES_AUTOHSCROLL
+    GROUPBOX        "", IDC_STATIC_ini1, 7, 222, 295, 73
+    LTEXT           "Stand identifier", IDC_STATIC_standid, 15, 256, 48, 8, SS_LEFT
+    PUSHBUTTON      "Browse", IDC_BUTTON_ini, 259, 235, 35, 12, WS_GROUP
+    LTEXT           "Management file", IDC_STATIC_man, 16, 107, 53, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_man, 109, 104, 146, 12, ES_AUTOHSCROLL
+    LTEXT           "Deposition data file", IDC_STATIC_dep, 16, 123, 61, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_dep, 109, 121, 146, 12, ES_AUTOHSCROLL
+    LTEXT           " N-reduction ( RedN) file", IDC_STATIC_red, 16, 141, 77, 8, SS_LEFT
+    LTEXT           "Litter intilisation file", IDC_STATIC_lit, 16, 158, 59, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_red, 109, 139, 146, 12, ES_AUTOHSCROLL
+    EDITTEXT        IDC_EDIT_lit, 109, 155, 146, 12, ES_AUTOHSCROLL
+    PUSHBUTTON      "Browse", IDC_BUTTON_sop, 259, 71, 35, 12, WS_GROUP
+    PUSHBUTTON      "Browse", IDC_BUTTON_soi, 259, 88, 35, 12, WS_GROUP
+    PUSHBUTTON      "Browse", IDC_BUTTON_man, 259, 105, 35, 12, WS_GROUP
+    PUSHBUTTON      "Browse", IDC_BUTTON_dep, 259, 123, 35, 12, WS_GROUP
+    PUSHBUTTON      "Browse", IDC_BUTTON_red, 259, 139, 35, 12, WS_GROUP
+    PUSHBUTTON      "Browse", IDC_BUTTON_lit, 259, 156, 35, 12, WS_GROUP
+    LTEXT           "Climate data file", IDC_STATIC_cli, 17, 32, 70, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_cli, 110, 34, 146, 12, ES_AUTOHSCROLL | ES_MULTILINE
+    PUSHBUTTON      "Browse", IDC_BUTTON_cli, 260, 35, 35, 12, WS_GROUP
+    COMBOBOX        IDC_COMBO_standid, 106, 254, 85, 77, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "WPM spinup intilisation file", IDC_STATIC_spinup, 15, 177, 84, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_wpm, 109, 174, 146, 12, ES_AUTOHSCROLL
+    PUSHBUTTON      "Browse", IDC_BUTTON_wpm, 259, 175, 35, 12, WS_GROUP
+    EDITTEXT        IDC_EDIT_standid, 106, 273, 76, 14, ES_AUTOHSCROLL
+    LTEXT           "Measurement file", IDC_STATIC_mes, 16, 194, 54, 8, SS_LEFT
+    EDITTEXT        IDC_EDIT_mes, 110, 191, 146, 12, ES_AUTOHSCROLL
+    PUSHBUTTON      "Browse", IDC_BUTTON_mes, 260, 191, 35, 12, WS_GROUP
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_flags DIALOGEX 100, 0, 580, 297
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+EXSTYLE WS_EX_TOOLWINDOW
+CAPTION "4C"
+FONT 8, "MS Sans Serif", 0, 0, 1
+{
+    DEFPUSHBUTTON   "OK", ID_FLAGS_BUTTON_OK, 443, 263, 50, 14
+    PUSHBUTTON      "Cancel", ID_CANCEL_FLAGS, 509, 262, 50, 14
+    LTEXT           "Mortality flag (flag_mort) ", IDC_STATIC_mort, 18, 33, 77, 8, SS_LEFT
+    LTEXT           "Regeneration flag (flag_reg) ", IDC_STATIC_reg, 18, 50, 90, 8, SS_LEFT
+    LTEXT           "Use FORSKA factors (flag_forska)", IDC_STATIC_forska, 18, 67, 113, 13, WS_TABSTOP | NOT WS_GROUP | SS_LEFT
+    LTEXT           "Stand initialization flag (flag_stand)", IDC_STATIC_stand, 18, 84, 109, 8, WS_TABSTOP | NOT WS_GROUP | SS_LEFT
+    LTEXT           "Soil vegetation flag (flag_sveg)", IDC_STATIC_sveg, 18, 101, 98, 8, SS_LEFT
+    LTEXT           "Management flag (flag_mg)", IDC_STATIC_mg, 18, 118, 86, 8, SS_LEFT
+    LTEXT           "Disturbance flag (flag_dis)", IDC_STATIC_dis, 18, 135, 82, 8, SS_LEFT
+    LTEXT           "Ligth algorithm number (flag_light)", IDC_STATIC_light, 18, 152, 106, 8, SS_LEFT
+    LTEXT           "Foliage-height relationship (flag_folhei)", IDC_STATIC_folhei, 18, 169, 120, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_mort, 143, 33, 133, 42, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_reg, 143, 50, 133, 116, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL
+    COMBOBOX        IDC_COMBO_forska, 143, 67, 133, 25, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "Volume function (flag_volfunc)", IDC_STATIC_volfunc, 18, 186, 96, 8, SS_LEFT
+    LTEXT           "Respiration flag (flag_resp)", IDC_STATIC_resp, 18, 203, 84, 8, SS_LEFT
+    LTEXT           "Limitation flag (flag_limi)", IDC_STATIC_limi, 18, 220, 74, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_stand, 143, 84, 133, 54, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_sveg, 143, 101, 133, 48, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_mg, 143, 118, 133, 78, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_dis, 143, 135, 133, 32, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_light, 143, 152, 133, 55, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_folhei, 143, 169, 133, 53, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_volfunc, 143, 186, 133, 63, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "Decomposition model (flag_decomp) ", IDC_STATIC_decomp, 308, 33, 117, 8, SS_LEFT
+    LTEXT           "Root activity function flag (flag_sign)", IDC_STATIC_sign, 308, 50, 115, 8, SS_LEFT
+    LTEXT           "Soil water uptake flag (flag_wred)", IDC_STATIC_wred, 308, 67, 106, 8, SS_LEFT
+    LTEXT           "Root distribution flag (flag_wurz)", IDC_STATIC_wurz, 308, 84, 101, 8, SS_LEFT
+    LTEXT           "Heat conductance flag (flag_cond)", IDC_STATIC_cond, 308, 101, 111, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_resp, 143, 203, 133, 64, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_limi, 143, 220, 133, 65, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_decomp, 431, 33, 133, 54, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_sign, 431, 50, 133, 67, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_wred, 431, 67, 133, 61, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_wurz, 431, 84, 133, 64, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_cond, 431, 101, 133, 78, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_int, 431, 118, 133, 75, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_eva, 431, 135, 133, 78, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "Interception flag (flag_int)", IDC_STATIC_int, 308, 118, 80, 8, SS_LEFT
+    LTEXT           "Evapotranspiration flag (flag_eva)", IDC_STATIC_eva, 308, 135, 106, 8, SS_LEFT
+    LTEXT           "Assortment flag (flag_sort)", IDC_STATIC_sort, 308, 169, 82, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_CO2, 431, 152, 133, 51, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_sort, 431, 169, 133, 48, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    GROUPBOX        "Model specification", IDC_STATIC_model, 0, 15, 564, 239
+    LTEXT           "wpm flag (flag_wpm)", IDC_STATIC_wpm, 307, 186, 65, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_wpm, 431, 186, 133, 64, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "Statistical analysis flag (flag_stat)", IDC_STATIC_stat, 307, 203, 104, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_stat, 431, 203, 133, 48, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    LTEXT           "CO2 flag (flag_CO2)", IDC_STATIC_CO2, 308, 152, 64, 8, SS_LEFT
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_ids DIALOG 100, 0, 130, 78
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_IDS_BUTTON_OK, 10, 55, 50, 14
+    GROUPBOX        "Output File Identifiers", IDC_STATIC_model, 6, 14, 111, 35
+    EDITTEXT        IDC_EDIT_ID, 14, 29, 95, 14, ES_AUTOHSCROLL
+    PUSHBUTTON      "Cancel", ID_CANCEL_IDS, 67, 55, 50, 14
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_main DIALOGEX 300, 50, 245, 268
+STYLE DS_3DLOOK | DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_VISIBLE | WS_POPUP | WS_SYSMENU
+EXSTYLE WS_EX_CLIENTEDGE | WS_EX_STATICEDGE
+CAPTION "4C"
+FONT 8, "MS Sans Serif", 0, 0, 1
+{
+    GROUPBOX        "", IDC_STATIC_year, 5, 6, 233, 258, 0, WS_EX_STATICEDGE
+    CTEXT           "Forest Growth Model 4C   PIK   2016", IDC_STATIC_4C, 25, 53, 82, 23, SS_CENTER
+    DEFPUSHBUTTON   "Start 4C", ID_START_4C, 103, 233, 56, 14
+    PUSHBUTTON      "Exit 4C", IDSTOP, 174, 233, 50, 14
+    GROUPBOX        "Simulation Control", IDC_STATIC_Control, 20, 126, 207, 68
+    AUTORADIOBUTTON "Start simulation", IDC_RADIO_start, 35, 148, 121, 10
+    AUTORADIOBUTTON "Edit control file", IDC_RADIO_edit, 35, 166, 62, 10
+    CONTROL         "", IDC_REBAR1, REBARCLASSNAME, WS_TABSTOP | 0x00000401, 4294967294, 4294967295, 249, 4, WS_EX_DLGMODALFRAME | WS_EX_ACCEPTFILES | WS_EX_STATICEDGE
+    CONTROL         IDB_BITMAP1, IDC_STATIC_pic, WC_STATIC, SS_BITMAP, 123, 14, 100, 104
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_out DIALOG 100, 0, 293, 212
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "Start simulation", ID_START, 89, 182, 56, 14
+    PUSHBUTTON      "Exit 4C", IDCANCEL, 226, 182, 50, 14
+    LTEXT           "Choice of output files", IDC_STATIC, 15, 18, 67, 8, SS_LEFT
+    PUSHBUTTON      "Yearly output", IDC_BUTTON_yearly, 207, 37, 50, 14
+    PUSHBUTTON      "Daily output", IDC_BUTTON_daily, 207, 60, 50, 14
+    PUSHBUTTON      "Cohorts yearly", IDC_BUTTON_coh_yearly, 207, 84, 50, 14
+    PUSHBUTTON      "Cohorts daily", IDC_BUTTON_coh_daily, 207, 108, 50, 14
+    GROUPBOX        "", IDC_STATIC_choice_out, 8, 7, 268, 164
+    LTEXT           "Summation output", IDC_STATIC_SUM, 19, 145, 58, 8, SS_LEFT
+    COMBOBOX        IDC_COMBO_sum, 88, 142, 165, 74, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_SORT
+    DEFPUSHBUTTON   "Back", ID_OUT_BUTTON_BACK, 161, 182, 50, 14
+    DEFPUSHBUTTON   "Save changes", ID_SAVE, 17, 182, 56, 14
+    COMBOBOX        IDC_COMBO_yearly, 14, 37, 181, 77, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_daily, 13, 60, 181, 77, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_coh_yearly, 13, 85, 181, 77, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+    COMBOBOX        IDC_COMBO_coh_daily, 13, 109, 181, 77, WS_TABSTOP | WS_VSCROLL | CBS_DROPDOWN | CBS_AUTOHSCROLL | CBS_SORT
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_4C_yearly DIALOG 100, 0, 606, 352
+STYLE DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_POPUP | WS_SYSMENU
+CAPTION "4C"
+FONT 8, "MS Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", ID_YEARLY_BUTTON_OK, 455, 327, 50, 14
+    PUSHBUTTON      "Cancel", ID_CANCEL_OUTF, 519, 327, 50, 14
+    GROUPBOX        "", IDC_STATIC_year, 7, 15, 573, 302
+    AUTOCHECKBOX    "", IDC_CHECK_y1, 21, 30, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y2, 21, 45, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y3, 21, 60, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y4, 21, 75, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y5, 21, 90, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y6, 21, 105, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y7, 21, 120, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y8, 21, 135, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y9, 21, 150, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y10, 21, 165, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y11, 21, 180, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y12, 21, 195, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y13, 21, 210, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y14, 21, 225, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y15, 21, 240, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y16, 21, 255, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y17, 21, 270, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y18, 21, 285, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y19, 21, 300, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y20, 294, 30, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y21, 294, 45, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y22, 294, 60, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y23, 294, 75, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y24, 294, 90, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y25, 294, 105, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y26, 294, 120, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y27, 294, 135, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y28, 294, 150, 270, 8
+    LTEXT           "  Choose yearly output files  ", IDC_STATIC_yfile, 21, 14, 89, 8, SS_LEFT
+    AUTOCHECKBOX    "", IDC_CHECK_y29, 294, 165, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y30, 294, 180, 270, 8
+    PUSHBUTTON      "Select all", ID_YEARLY_BUTTON_SELECT, 391, 327, 50, 14
+    PUSHBUTTON      "Deselect all", ID_YEARLY_BUTTON_DESELECT, 327, 327, 50, 14
+    AUTOCHECKBOX    "", IDC_CHECK_y31, 294, 195, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y32, 294, 210, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y33, 294, 225, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y34, 294, 240, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y35, 294, 255, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y36, 294, 270, 270, 8
+    AUTOCHECKBOX    "", IDC_CHECK_y37, 294, 285, 270, 8
+}
+
+
+
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+IDD_DIALOG1 DIALOG 0, 0, 186, 95
+STYLE DS_3DLOOK | DS_CENTER | DS_MODALFRAME | DS_SETFONT | WS_CAPTION | WS_VISIBLE | WS_POPUP | WS_SYSMENU
+CAPTION "Dialog"
+FONT 8, "Microsoft Sans Serif"
+{
+    DEFPUSHBUTTON   "OK", IDOK, 129, 7, 50, 14
+    PUSHBUTTON      "Cancel", IDCANCEL, 129, 24, 50, 14
+}
+
+
+
+//
+// Version Information resources
+//
+LANGUAGE LANG_GERMAN, SUBLANG_GERMAN
+VS_VERSION_INFO VERSIONINFO
+    FILEVERSION     1,0,0,0
+    PRODUCTVERSION  0,0,0,0
+    FILEOS          VOS_NT_WINDOWS32
+    FILETYPE        VFT_APP
+    FILESUBTYPE     VFT2_UNKNOWN
+    FILEFLAGSMASK   0x0000003F
+    FILEFLAGS       0x00000000
+{
+    BLOCK "StringFileInfo"
+    {
+        BLOCK "080904B0"
+        {
+            VALUE "Comments", "\0"
+            VALUE "CompanyName", "Potsdam-Institut fr Klimafolgenforschung\0"
+            VALUE "FileDescription", "FORESEE - Forest Ecosystems in a changing Environment\0"
+            VALUE "FileVersion", "0.99e\0"
+            VALUE "InternalName", "4C\0"
+            VALUE "LegalCopyright", "Copyright  2004\0"
+            VALUE "LegalTrademarks", "\0"
+            VALUE "OriginalFilename", "4C.exe\0"
+            VALUE "PrivateBuild", "\0"
+            VALUE "ProductName", "4C - PIK\0"
+            VALUE "ProductVersion", "0.99e\0"
+            VALUE "SpecialBuild", "\0"
+        }
+    }
+    BLOCK "VarFileInfo"
+    {
+        VALUE "Translation", 0x0809, 0x04B0
+    }
+}

source_code/version_2.3_windows/seed_multi.f

+!*****************************************************************!
+!*                                                               *!
+!*              4C (FORESEE) Simulation Model                    *!
+!*                                                               *!
+!*                                                               *!
+!*         SR   SEED_multi                                       *!
+!*                                                               *!
+!*        including SR/Function                                  *!
+!*         function rtflsp    (regula falsi solving equation)    *!
+!*         function weight                                       *!
+!*         function weight1                                      *!
+!*                                                               *!
+!*   generation of a variety of seedling cohorts for             *!
+!*   one seed number according to seedmass distribution          *!
+!*   (for given mean value and standard deviation)               *!
+!*                                                               *!
+!*                  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 seed_multi(nseed,nsp)
+
+USE data_species
+use data_stand
+use data_help
+use data_par
+use data_soil
+use data_simul
+
+IMPLICIT NONE
+integer    :: nseed, nseedha, nsclass , k, j, nr
+integer,dimension(:),allocatable  :: nsc
+
+real, dimension(:), allocatable   :: msc,      &
+                                     shooth,   &
+                                     nschelp
+integer    :: nsp
+REAL       :: shoot
+REAL       :: ms, msclass, x1,x2,xacc,shelp, nshelp,ntot,help
+REAL	   :: troot2
+
+real       :: standdev
+real       :: rtflsp, weight
+
+TYPE(cohort)    ::tree_ini
+
+external weight
+external rtflsp
+
+if(nseed.eq.0) return
+ standdev = spar(nsp)%seedsd*1000.
+ hnspec = nsp
+ ms = spar(nsp)%seedmass *1000.        ! g ---> mg
+ nseedha = nseed
+ nshelp = nseedha/10000.
+
+! calculation of seed class number
+ if(flag_reg.eq.3) then
+     nsclass = int(100.*nshelp**0.6)
+ else if(flag_reg.eq.30) then
+  nsclass = int(10.*nshelp**0.6)+1
+ end if
+allocate(nsc(nsclass))
+allocate(nschelp(nsclass))
+allocate(msc(nsclass))
+allocate(shooth(nsclass))
+
+! seed weight and number of seeds per class
+  msclass = 6.*standdev/nsclass
+  ntot = 0
+  help =  (1/(sqrt(2*pi)*standdev))
+  do k=1, nsclass
+
+        msc(k) = (ms - 3.*standdev) + msclass*(k-1)
+        nschelp(k) = help*exp(-((msc(k)-ms)**2)/(2*(standdev)**2))
+        ntot = ntot + nschelp(k)
+
+  end do
+
+  do k= 1,nsclass
+
+       nsc(k) = nint((nschelp(k)*nseedha/ntot) + 0.5)
+
+  end do
+! calculation of shoot weight per seed class and initilization
+
+  do k = 1,nsclass
+
+    mschelp = msc(k)/1000000.      ! mg ---> kg
+    x1 = 0.
+    x2 = 0.1
+    xacc=(1.0e-10)*(x1+x2)/2
+    
+! solve mass equation; determine root
+    shelp=rtflsp(weight,x1,x2,xacc)
+    shooth(k)= shelp
+    max_coh = max_coh + 1
+
+    call coh_initial (tree_ini)
+
+    tree_ini%ident =  max_coh
+    tree_ini%species = nsp
+    tree_ini%ntreea = nsc(k)
+    tree_ini%nta = nsc(k)
+    shoot = shooth(k)
+    tree_ini%x_sap = shoot                                ! [kg]
+    shoot = shoot * 1000.                                  ! [g]
+    tree_ini%x_fol= (spar(nsp)%seeda*(tree_ini%x_sap** spar(nsp)%seedb))   ![kg]
+    tree_ini%x_frt = tree_ini%x_fol                                          ! [kg]
+! Leder
+    tree_ini%x_hrt = 0.
+    tree_ini%med_sla = spar(nsp)%psla_min + spar(nsp)%psla_a*0.5
+    tree_ini%t_leaf = tree_ini%med_sla* tree_ini%x_fol                              ! [m-2]
+    tree_ini%ca_ini = tree_ini%t_leaf
+    tree_ini%crown_area = tree_ini%ca_ini
+    tree_ini%underst = 1
+
+! tranformation of shoot biomass kg --> mg
+
+    if(nsp.ne.2)tree_ini%height = spar(nsp)%pheight1*(shoot*1000.)**spar(nsp)%pheight2      ! [cm] berechnet aus shoot biomass (mg)
+! Leder
+
+    if(nsp.eq.2) tree_ini%height = 10**(spar(nsp)%pheight1+ spar(nsp)%pheight2*LOG10(shoot*1000.)+ &
+                              spar(nsp)%pheight3*(LOG10(shoot*1000.))**2)
+       IF(nsc(k).ne.0.) then
+       IF (.not. associated(pt%first)) THEN
+          ALLOCATE (pt%first)
+          pt%first%coh = tree_ini
+          NULLIFY(pt%first%next)
+!     root distribution
+        call root_depth (1, pt%first%coh%species, pt%first%coh%x_age, pt%first%coh%height, pt%first%coh%x_frt, pt%first%coh%x_crt, nr, troot2, pt%first%coh%x_rdpt, pt%first%coh%nroot)
+        pt%first%coh%nroot = nr
+        do j=1,nr
+           pt%first%coh%rooteff = 1.   ! assumption for the first use
+        enddo
+        do j=nr+1, nlay
+           pt%first%coh%rooteff = 0.   ! layers with no roots
+        enddo
+       ELSE
+           
+          ALLOCATE(zeig)
+          zeig%coh = tree_ini
+          zeig%next => pt%first
+          pt%first => zeig
+
+          call root_depth (1, zeig%coh%species, zeig%coh%x_age, zeig%coh%height, zeig%coh%x_frt, zeig%coh%x_crt, nr, troot2, zeig%coh%x_rdpt, zeig%coh%nroot)
+        zeig%coh%nroot = nr
+        do j=1,nr
+           zeig%coh%rooteff = 1.   ! assumption for the first use
+        enddo
+        do j=nr+1, nlay
+           zeig%coh%rooteff = 0.   ! layers with no roots
+        enddo
+
+      END IF
+      anz_coh=anz_coh+1
+    END IF
+  end do
+
+deallocate(nsc)
+deallocate(nschelp)
+deallocate(msc)
+deallocate(shooth)
+
+END SUBROUTINE  seed_multi
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!    weight: seed mass function
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function weight (x)
+
+use data_help
+use data_species
+
+implicit none
+
+real     :: x
+real     :: p1,p2,  weight
+
+p1 = spar(hnspec)%seeda
+p2 = spar(hnspec)%seedb
+
+weight = p1*2*(x**p2) + x - 0.7*mschelp
+
+end function weight
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!    weight1: coarse root mass function
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+function weight1 (x)
+
+use data_help
+use data_species
+
+real     :: x
+real     :: p1,p2
+
+p1 = spar(hnspec)%seeda
+p2 = spar(hnspec)%seedb
+
+weight1 = p1*(x**p2) + x - mschelp
+
+end function weight1
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!      rtflsp: regula falsi solving euation
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+FUNCTION rtflsp(func,x1,x2,xacc)
+  INTEGER MAXIT
+  REAL rtflsp,x1,x2,xacc,func
+  EXTERNAL func
+  PARAMETER (MAXIT=30)
+  INTEGER j
+  REAL del,dx,f,fh,fl,swap,xh,xl
+  fl=func(x1)
+  fh=func(x2)
+  if(fl.lt.0.)then
+        xl=x1
+        xh=x2
+  else
+        xl=x2
+        xh=x1
+        swap=fl
+        fl=fh
+        fh=swap
+   endif
+   dx=xh-xl
+   do j=1,MAXIT
+        rtflsp=xl+dx*fl/(fl-fh)
+        f=func(rtflsp)
+
+        if(f.lt.0.) then
+          del=xl-rtflsp
+          xl=rtflsp
+          fl=f
+        else
+          del=xh-rtflsp
+          xh=rtflsp
+          fh=f
+        endif
+        dx=xh-xl
+        if(abs(del).lt.xacc.or.f.eq.0.)return
+   end do
+END function rtflsp
\ No newline at end of file

source_code/version_2.3_windows/sim_ini.f

+!*****************************************************************!
+!*                                                               *!
+!*              4C (FORESEE) Simulation Model                    *!
+!*                                                               *!
+!*                                                               *!
+!*                    Subroutines for:                           *!
+!*    - Simulation initialisation (SIM_INI)                      *!
+!*                                                               *!
+!*                  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 sim_ini
+
+use data_biodiv
+use data_climate
+use data_depo
+use data_evapo
+use data_inter
+use data_manag
+use data_simul
+use data_site
+use data_stand
+use data_soil
+use data_soil_cn
+use data_species
+use data_par
+use data_frost
+
+implicit none
+
+type(Coh_Obj), pointer :: p  ! pointer to cohort list
+
+anz_sim  = anz_sim + 1
+time     = 0
+time_cur = time_b - 1    ! before Sim.-Start in year_ini time_cur=time_cur+1
+iday     = 0
+
+act_thin_year = 1
+flag_cum  = 0
+flag_lit  = 0
+flag_sens = 0
+flag_redn = .FALSE.
+
+lai=0.
+gp_can    = 0.
+sumbio    = 0.
+totfrt    = 0.
+sumNPP    = 0.
+nppsum    = 0.
+gppsum    = 0.
+cum_sumNPP= 0.
+NEE_mon   = 0.
+NPP_mon   = 0.
+autresp   = 0.
+autresp_m = 0.
+anrspec   = 0
+anz_coh   = 0
+anz_spec  = 0
+anz_tree  = 0
+med_diam  = 0.
+hdom      = 0.
+mean_drIndAl = 0.
+
+med_air    = 0.
+med_rad    = 0.
+med_air_cm = 0.
+med_air_wm = 0.
+med_air_ms = 0.
+med_air_mj = 0.
+med_wind   = 0.
+temp_mon   = 0.
+prec_mon   = 0.
+sum_prec   = 0.
+sum_prec_ms= 0.;
+sum_prec_mj= 0.
+gdday = 0.
+days_summer = 0
+days_hot    = 0
+days_ice    = 0
+days_dry    = 0
+days_hrain  = 0
+days_rain   = 0
+days_rain_mj= 0
+days_snow   = 0
+days_wof = 0
+gdday_all = 0.
+med_air_all  = 0.
+sum_prec_all = 0.
+med_rad_all = 0.
+int_cum_can  = 0.
+int_cum_sveg = 0.
+interc_m_can = 0.
+interc_m_sveg= 0.
+perc_cum    = 0.
+perc_mon    = 0.
+wupt_cum    = 0.
+wupt_r_c    = 0.
+wupt_e_c    = 0.
+tra_tr_cum  = 0.
+tra_sv_cum  = 0.
+dew_m       = 0.
+aet_cum     = 0.
+pet_cum     = 0.
+pet_m		= 0.
+AET_m       = 0.
+wupt_r_m    = 0.
+perc_m      = 0.
+wat_tot     = 0.
+gp_can_mean = 0.
+gp_can_max  = 0.
+snow        = 0.
+snow_day    = 0
+totFPARcan  = 0.
+Rnet_cum    = 0.
+
+! fire index
+fire(1)%mean_m = 0
+fire(2)%mean_m = 0
+fire(3)%mean_m = 0
+fire_indb_m    = 0
+
+ind_arid_an = 0.
+ind_lang_an = 0.
+ind_cout_an = 0.
+ind_wiss_an = 0.
+ind_mart_an = 0.
+ind_mart_vp = 0.
+ind_emb		= 0.
+ind_weck	= 0.
+ind_reich	= 0.
+con_gor		= 0.
+con_cur		= 0.
+con_con		= 0.
+cwb_an		= 0.
+cwb_an_m    = 0.
+ind_bud     = 0.
+ind_shc     = 0.
+
+ind_arid_an_m = 0.
+ind_lang_an_m = 0.
+ind_cout_an_m = 0.
+ind_wiss_an_m = 0.
+ind_mart_an_m = 0.
+ind_mart_vp_m = 0.
+ind_emb_m     = 0.
+ind_weck_m    = 0.
+ind_reich_m   = 0.
+con_gor_m     = 0.
+con_cur_m     = 0.
+con_con_m     = 0.
+ind_bud_m     = 0.
+ind_shc_m     = 0.
+ntindex       = 0.
+
+tempmean_mo = 0
+
+aet_dec  = 0.
+temp_dec = 0.
+prec_dec = 0.
+rad_dec  = 0.
+hum_dec  = 0.
+
+! frost index
+if(flag_climtyp .ge. 3) then
+ ! calculation for airtemp_min > -90. 
+    tminmay=0
+    lfind=0
+    dlfabs=0.
+    tminmay_sp=0
+    dlfabs_sp=0.
+    flag_tveg=0
+else
+    tminmay=-99
+    lfind=-99
+    dlfabs=-99.
+    tminmay_sp=-99
+    dlfabs_sp=-99.
+    flag_tveg=-99
+endif
+
+
+!! initialisation of root distribution
+RedN_mean = 0.
+anz_RedN  = 0
+N_min     = 0.
+N_min_m   = 0.
+resps_c   = 0.
+resps_c_m = 0.
+resps_mon = 0.
+N_tot     = 0.
+N_an_tot  = 0.
+N_hum_tot = 0.
+C_tot     = 0.
+C_hum_tot = 0.
+N_lit     = 0.
+C_lit     = 0.
+Nupt_c    = 0.
+Nupt_m    = 0.
+C_accu      = 0.
+Nleach_c    = 0.
+Nleach_m    = 0.
+N_lit_m     = 0.
+C_lit_m     = 0.
+totfol_lit  = 0.
+totfol_lit_tree = 0.
+totfrt_lit  = 0.
+totfrt_lit_tree = 0.
+tottb_lit   = 0.
+totcrt_lit  = 0.
+totstem_lit = 0.
+C_opm_fol   = 0.
+C_opm_frt   = 0.
+C_opm_crt   = 0.
+C_opm_tb    = 0.
+C_opm_stem  = 0.
+N_opm_stem  = 0.
+N_opm_fol   = 0.
+N_opm_frt   = 0.
+N_opm_crt   = 0.
+N_opm_tb    = 0.
+Ndep_cum    = 0.
+Ndep_cum_all= 0.
+if (flag_multi .ne. 8) then
+    if ((flag_multi .ne. 2) .or. (ip .le. 1)) then
+        NOdep(ip) = 0.
+        NHdep(ip) = 0.
+    endif
+endif
+
+flag_bc = 0
+totsteminc   = 0.
+cumsteminc   = 0.
+cumsumvsdead = 0.
+cumsumvsab   = 0.
+sumvsdead    = 0.
+sumvsab      = 0.
+
+p => pt%first
+do while (associated(p))
+  p%coh%N_pool  = 0.
+
+  p => p%next
+enddo  ! p (cohorts)
+
+allocate(dayfract(ns_pro))
+
+!  fields for frost index
+allocate(dnlf(year))
+allocate(tminmay_ann(year))
+allocate(date_lf(year))
+allocate(date_lftot(year))
+allocate(dnlf_sp(year))
+allocate (anzdlf(year))
+allocate (sumtlf(year))
+
+dnlf_sp=0
+dnlf = 0
+tminmay_ann = 0.
+date_lf = 0
+date_lftot = 0
+anzdlf = 0.
+sumtlf = 0.
+
+end subroutine sim_ini
\ No newline at end of file

source_code/version_2.3_windows/simul.f

+!*****************************************************************!
+!*                                                               *!
+!*              4C (FORESEE) Simulation Model                    *!
+!*                                                               *!
+!*                                                               *!
+!*                    Subroutines for:                           *!
+!*    - Simulation control: SIM_CONTROL                          *!
+!*                          SIMULATION_4C                        *!
+!*                                                               *!
+!*                  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 sim_control
+
+use data_climate
+use data_simul
+use data_site
+use data_out
+
+implicit none
+
+integer run_nr, ipp, irl, icl, i
+character a
+character(8) actdate
+character(10) acttime, helpsim, text1, text2
+real time1, time2, time3
+
+unit_err=getunit()
+if(flag_multi.eq.5) dirout = './'
+open(unit_err,file=trim(dirout)//trim(site_name(1))//'_error.log',status='replace', position='append')
+unit_trace=getunit()
+open(unit_trace,file=trim(dirout)//trim(site_name(1))//'_trace.log',status='replace', position='append')
+write (unit_trace, *) ' Trace of calls  -  subroutines of 4C '
+write (unit_trace, *) 
+write (unit_trace, *) 'iday  time_cur    subroutine '
+write (unit_trace, '(I4,I10,A)') iday, time_cur, ' sim_control'
+
+! check daily output
+if (year > 5 .and. flag_dayout .ge. 1) then
+       write(*,*) '  Warning: Your choice of daily output is ON with a simulation time of'
+       write(*,'(I6,A,I8,A)') year,'  years. This option will create ',365*year,' data records per file!'
+       write(*,'(A)',advance='no') '  Do you really want do use daily output (y/n)? '
+       read *,a
+       IF (a .eq. 'n' .or. a .eq. 'N') then
+          flag_dayout = 0
+       ENDIF
+ENDIF
+
+! open file ycomp (yearly compressed output (multi run))
+IF (time_out .ne. -2) call prep_out_comp
+
+!call epsilon
+	IF(flag_multi.eq.1) THEN
+	   run_nr = site_nr
+	ELSE IF (flag_multi.eq.5) THEN
+	   run_nr = 1
+	ELSE
+	   run_nr = repeat_number
+	ENDIF
+
+call date_and_time(actdate, acttime)
+write (unit_err, *)
+time3 = 0.
+
+  if (.not.flag_mult910) then
+    nrreal = 1
+    nrclim = 1
+  endif
+   
+  do icl = 1, nrclim   ! climate scenarios
+    iclim = icl
+	DO ipp = 1, run_nr ! sites
+     ip = ipp
+     if (flag_trace) write (unit_trace, '(I4,I10,A,I3)') iday, time_cur, ' sim_control ip=',ip  
+      do irl = 1, nrreal  ! realization f climate scenarios
+           if (flag_mult910) then
+               climfile(ip) = climszenfile(ip, icl, irl)
+               site_name (ip) = trim(site_name1)//'_'//trim(sitenum(ip))
+               write (helpsim,'(I10)') icl
+               read (helpsim,*) text1
+               write (helpsim,'(I10)') irl
+               read (helpsim,*) text2
+               site_name (ip) = trim(site_name (ip))//'_'//trim(text1)//'_'//trim(text2)
+               write (unit_err, *)
+               write (unit_err, '(A,3I5,2X, A50)')'*  ip, cli-scenario, real., site:  ', ip, icl, irl, site_name(ip)
+               write (unit_err, '(A,A)') 'Climate file:  ', trim(climfile(ip))
+  	       else
+               write (unit_err, *)
+               write (unit_err, '(A10,I5,2X, A50)') ' ip/site ', ip, trim(site_name(ip))
+               site_name1 = trim(site_name(ip))
+           endif
+           call CPU_time (time1)
+           if(ip.ne.0) then
+	           CALL sim_ini
+	           CALL prepare_site
+	           if (flag_multi.eq.5) then
+!                  call m4c_simenv_in
+                  unit_comp2 = 6    ! standard output 
+               end if
+
+	           if(flag_end.gt.0) then
+                  select case (flag_end)
+                  case (1)
+     		        print*,ip, '  stop in prepare_stand (see error.log)'
+                  case (2)
+     		        print*, ip, 'stop in prepare_stand, stand ',   &
+     		        'identificator not found in prepare_stand'
+                  case (3)
+                    print*,ip, 'stop in canopy'
+                  case (4)
+                    print*,ip, 'stop in readsim, climate ID not found'
+                  case (5)
+                    print*,ip, ' stop in readsoil, soil ID not found ', adjustl(soilid(ip))
+                  case (6)
+                    write(*,'(A,I5)') '  >>>foresee message: stop in read_cli, no climate data for year ',time_b
+                    call finish_simul
+                    stop
+                  case (7)
+                    print*,ip, ' stop in readsoil, error during reading soil data ', adjustl(soilid(ip))
+                    call finish_simul
+                    stop
+                 case default
+                    print*,ip, 'flag_end = ', flag_end
+                  end select
+
+	              call finish_simul
+	              flag_end = 0
+	           else
+	              IF(flag_multi==2) CALL fixclimscen
+                    if (.not.flag_mult910) then
+                        write (*,*)
+                        write (*,*)  '>>> Start FORESEE-Simulation site ', ipp
+                        write (*,*)
+                    endif
+
+	              CALL simulation_4c
+	              CALL finish_simul
+	           endif
+               if (flag_mult910)  then
+                    call out_var_stat(1, irl)
+               else
+                    if ((flag_multi .ne. 8) .and. (nvar .gt. 1)) call out_var_stat(3, 1)
+               endif 
+           endif  ! flag_end
+ 	       call CPU_time (time2)
+           if (.not.flag_mult910)  then
+               print *, ' run time for simulation ',ip, time2-time1, ' sec'
+           endif
+           write (unit_err, *) ' run time for simulation ',ip, time2-time1, ' sec'
+           time3 = time3 + (time2-time1)
+	  enddo  ! irl 
+      if (flag_mult910)  call out_var_stat(2, -99) 
+      write (unit_err, *)
+      write (unit_err, *)
+      write (unit_err, *) '* * * * *  New ip/site  * * * * *'
+    ENDDO   ! ip until site_nr (page number)
+    write (unit_err, *)
+    write (unit_err, *) '************  New climate scenario  **********'
+  enddo   ! icl
+
+	if (nvar .gt. 1) then 
+        select case (flag_multi)
+        case (5, 9, 10)
+            continue
+        case (1)
+            continue
+        case default    
+            call out_var_file
+        end select
+    endif
+
+!  comparison with measurements
+    if (flag_stat .gt. 0) CALL mess	
+    
+ 	call CPU_time (time1)
+    time3 = time3 + (time1-time2)
+    write (unit_err, *)
+    write (unit_err, *) ' total run time ', time3, ' sec'
+    CALL finish_all
+	PRINT *,'All simulations finished!'
+
+END SUBROUTINE sim_control
+
+!**************************************************************
+
+SUBROUTINE simulation_4C
+
+  !*** Declaration part ***!
+
+  USE data_simul
+  USE data_species    ! species specific parameters
+  USE data_site       ! site specific data
+  USE data_climate    ! climate data
+  USE data_soil
+  USE data_soil_cn
+  USE data_stand           ! state variables of stand, cohort and cohort element
+  USE data_out
+  USE data_manag
+  USE data_plant
+  USE data_par
+  IMPLICIT NONE
+
+if (flag_trace) write (unit_trace, '(I4,I10,A)') iday, time_cur, ' simulation_4C'
+
+! allocation of environmental variable fields
+ if(flag_wpm.ne.4 .and. flag_wpm.ne.5.and.flag_wpm.ne.6) then
+  ! time loop
+  DO time = 1, year
+     iday = 1
+  ! Update population variable for new year if population is changed through interventions 
+     if (flag_standup .gt. 0 .or. flag_dis==1 .or. flag_dis==1) then
+        call stand_balance
+        call standup
+        flag_standup = 0
+     endif
+
+      CALL year_ini
+      ! Calculate RedN from soil C/N
+      ! read or create Redn for areal application 
+
+     IF (time.EQ.1 .and. flag_redn) CALL RedN_ini
+     
+     IF (flag_dis .eq. 1 .or. flag_dis .eq. 2) then
+         CALL dis_manag
+      endif
+  
+      ! simulation of processes with subannual resolution (fluxes and soil)
+      CALL stand_daily
+      if(flag_end.ge.1) exit ! exit do loop time
+
+      ! compressed output of start values
+      IF (lcomp1) THEN                           
+          CALL out_comp(unit_comp1)
+          lcomp1 = .FALSE.
+      ENDIF
+
+      ! cohort litter production
+      CALL senescence
+
+      ! calculation of stand variables over all patches
+       CALL stand_balance
+
+      ! calculation of soil variables for yearly output
+      CALL s_year
+
+      ! calculation of fire variables for yearly output
+      CALL fire_year
+
+      ! calculation of indices for yearly output
+      CALL t_indices(temp_mon)
+
+      ! summation output
+      IF(flag_sum.eq.4) THEN
+         write(unit_sum,'(I5,9F11.3)') time_cur,photsum,npppotsum,nppsum,resosum,lightsum,nee,abslightsum,precsum, tempmean
+         photsum=0.;npppotsum=0.;nppsum=0.;resosum=0.;lightsum=0.;nee=0.;abslightsum=0.;precsum=0.
+      ENDIF
+
+       totsteminc = 0.
+       totsteminc_m3 = 0.
+      ! cohort loop for change in crown dimensions, allocation and tree dimension calculations
+      zeig=>pt%first
+      DO
+          IF (.not.ASSOCIATED(zeig)) exit
+ 
+          IF (zeig%coh%height.ge.thr_height .and. zeig%coh%species.le.nspec_tree) then
+
+              ! determine crown movement
+              CALL CROWN( zeig )
+
+              ! allocate NPP to the various tree compartments
+              CALL PARTITION( zeig )
+             if(flag_end.ge.1) exit  ! do loop
+          ENDIF
+          IF (zeig%coh%species.EQ.nspec_tree+1) then            ! Ground vegetation
+              ! allocate NPP to the various ground vegetation compartments
+              CALL PARTITION_SV( zeig )
+          ENDIF
+          IF (zeig%coh%species.eq.nspec_tree+2)  then             ! Mistletoe
+              CALL PARTITION_MI( zeig )
+              if(flag_end.ge.1) exit  ! do loop
+          ENDIF
+        zeig=>zeig%next
+      END DO
+      if(flag_end.ge.1) exit ! exit do loop time
+
+ ! calculation of annual mortality
+      IF(flag_mort.ge.1) CALL stand_mort
+
+! annual growth of trees below thr_height, which are initialized by planting (not seeded!)
+! at the beginning of the simulation or during management (shelter-wood)
+        if(flag_reg.ne.2.and.flag_sprout.eq.0) CALL growth_seed
+         CALL mort_seed
+      if(flag_sprout.eq.1) flag_sprout=0
+      IF(flag_mg==1) then
+	     if(thin_year(act_thin_year)==time_cur) then
+            CALL management
+             act_thin_year = act_thin_year+1
+		 end if
+       ELSE IF((flag_mg.ge.2 .or. flag_mg.eq.3 .or. flag_mg.eq.33.or. flag_mg.eq.9 .or. flag_mg.eq.10 .or. flag_mg.eq. 333).and.anz_spec.ne.0) THEN
+         CALL management
+         if(flag_wpm.ne.0) CALL timsort
+       ENDIF
+
+! no assortment if wpm is not called
+      if(flag_mg.eq.0.and.anz_spec.ne.0) then
+         if(flag_wpm.ne.0) call timsort
+      end if
+      CALL litter
+
+! input of dead biomass into soil compartments
+      CALL cn_inp
+
+ !     if(flag_multi.eq.5) call m4c_simenv_out
+! annual establishment for all species
+      IF (flag_reg.eq.1.or.flag_reg.eq.2.or.flag_reg.eq.3.or.flag_reg.eq.30)  CALL stand_regen
+
+! cumsteminc = cumsteminc + totsteminc   
+! planting of seedlings/saplings at the beginning of simulation
+      if(flag_reg.ge.9..and. flag_reg.lt.100. .and. time.eq.1) call planting
+      if(flag_reg.ge.9..and. flag_reg.lt.100. .and. flag_mg .eq.44) call planting
+
+! Update stand variables if stand changed
+     if (flag_standup.gt.0 .or. anz_spec.eq.0) then
+        call stand_balance
+!        if (flag_standup .gt. 1) call root_distr   ! wird generell in year_ini berechnet
+     endif
+    cumsteminc = cumsteminc + totsteminc   
+
+	
+
+      ! yearly output
+      IF (time_out .gt. 0) THEN
+         IF (mod(time,time_out) .eq. 0) then
+             CALL outyear (1)
+             CALL outyear (2)
+         endif
+     ENDIF
+
+      ! store of output variables (multi run 4, 8, 9)
+      IF (nvar .gt. 1) CALL outstore
+
+    ! RedN calculation
+     if ((flag_limi .eq. 10) .or. (flag_limi .eq. 15)) call RedN_calc
+
+!     CALL list_cohort
+     CALL del_cohort
+
+     if (.not.flag_mult910) PRINT *, ' * Year ', time, time_cur,' finished... '
+
+   END DO ! time
+    ! calculation of stand variables over all patches at the end!
+      CALL stand_balance
+
+         !*****    wpm    ******
+	! check if management
+	if(flag_mg == 0) then
+		flag_wpm = 0
+	endif	
+
+	if (flag_wpm == 1 .or. flag_wpm == 21 .or. flag_wpm == 11) call wpm
+	if (flag_wpm == 2) call sea
+	if (flag_wpm == 3) then
+		call wpm
+		call sea
+	end if
+         !*** * * * * * * * ****
+else
+    call wpm
+end if
+   if (flag_wpm .gt. 0) call out_wpm(1)   
+
+   CALL out_comp(unit_comp2)
+
+   if(flag_end.eq.1) print*,ip, 'stop in partitio'
+   if(flag_end.eq.3) print*,ip, 'stop in calc_la in canopy: toplayer = 125 m'
+   flag_end = 0
+   if (.not.flag_mult910) PRINT *, ' * Simulation ',ip,' finished.'
+
+END SUBROUTINE simulation_4C
+
+
+!**************************************************************

source_code/version_2.3_windows/soil.f

+!*****************************************************************!
+!*                                                               *!
+!*              4C (FORESEE) Simulation Model                    *!
+!*                                                               *!
+!*                                                               *!
+!*                    Subroutines for:                           *!
+!*                Soil and Water - Programs                      *!
+!*                                                               *!
+!*   contains:                                                   *!
+!*   SOIL                                                        *!
+!*   SOIL_INI                                                    *!    
+!*   SOIL_WAT                                                    *!
+!*   UPT_WAT                                                     *!
+!*   FRED1 - ...11                                                *!
+!*   TAKE_WAT                                                    *!
+!*   BUCKET                                                      *!
+!*   SNOWPACK                                                    *!
+!*   HUM_ADD                                                     *!
+!*   BC_APPL: application of biochar                             *!
+!*                                                               *!
+!*                  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 soil
+
+!   Soil processes (frame)
+
+use data_climate
+use data_depo
+use data_out 
+use data_simul
+use data_soil
+use data_soil_cn
+
+implicit none
+
+call evapo
+call intercep
+call soil_wat
+call soil_temp
+
+if (flag_wurz .eq. 4 .or. flag_wurz .eq. 6) then
+    call soil_stress	!calculate ground stress factors
+    call cr_depth		!define root depth 
+endif
+call soil_cn
+call root_eff
+
+END subroutine soil
+
+!**************************************************************
+
+SUBROUTINE soil_ini
+
+! Initialisation of soil data and parameters
+
+use data_inter
+use data_evapo
+use data_out
+use data_par
+use data_simul
+use data_soil
+use data_soil_cn
+use data_species
+use data_stand
+
+implicit none
+
+integer i,j,k
+real d_0, t_0
+! Table of quarz and clay content (mass%) versus wlam
+real, dimension(17) :: xwlam  = (/1.5, 1.15, 0.9, 0.67, 0.6,  0.5, 0.38, 0.37, 0.3, 0.29, 0.27, 0.26, 0.25, 0.24, 0.23, 0.22, 0.15/),  &
+                       yquarz = (/93.0,85.0,80.0, 82.0,76.0, 64.0, 65.0, 51.0,60.0, 30.0, 14.0, 10.0, 12.0, 20.0, 30.0, 43.0, 23.0/),  &
+                       yclay  = (/3.0, 3.0,  3.0, 12.0, 6.0,  6.0, 10.0,  4.0,21.0, 12.0, 10.0, 37.0, 15.0, 40.0, 30.0, 35.0, 55.0/) 
+real value
+real, dimension(nlay):: humush(nlay)
+real, dimension(nlay):: xfcap, xwiltp, xpv   ! output of addition for water capacities
+
+! estimation of soil layer values
+d_0 = 0.
+do j = 1, nlay
+   t_0      = thick(j)
+   mid(j)   = d_0 + 0.5*t_0
+   d_0      = d_0 + t_0
+   depth(j) = d_0
+enddo
+
+perc	= 0.
+wupt_r	= 0.
+wupt_ev	= 0.
+thick_1 = thick(1)
+
+select case (flag_soilin)
+case (0,2)
+    do i=1,nlay
+        if (i .gt. 1) then
+            call tab_int(xwlam, yquarz,  17, wlam(i), value)
+            sandv(i)  = value / 100.   ! Mass% of mineral fraction
+            call tab_int(xwlam, yclay,  17, wlam(i), value)
+            clayv(i)  = value / 100.
+            siltv(i)  = 1. - clayv(i) - sandv(i)
+        else
+            sandv(1)  = 0.0
+            clayv(1)  = 0.0
+            siltv(1)  = 0.0
+        endif
+    enddo
+
+case (1,3,4)
+    clayv  = clayv / 100.
+    sandv  = sandv / 100.
+    siltv  = 1. - clayv - sandv
+    if ((sandv(1) .le. zero) .and. (clayv(1) .le. zero)) siltv(1) = 0.
+    skelv  = skelv / 100.
+    humusv = humusv / 100.
+end select
+
+! Settings for subroutine take_wat
+skelfact  = 1.
+pv        = skelfact * pv_v * thick * 0.1	       ! mm
+wilt_p    = skelfact * wilt_p_v * thick * 0.1	   ! mm
+field_cap = skelfact * f_cap_v * thick * 0.1	   ! mm
+wats      = field_cap   ! mm
+watvol    = f_cap_v                    ! vol%
+
+n_ev_d = nlay
+nlgrw  = nlay+1
+do i=1,nlay
+	if (w_ev_d .gt. depth(i)) n_ev_d = i
+	if (grwlev .gt. depth(i)) nlgrw = i+1
+    vol(i)    = thick(i) * 10000.  
+enddo
+
+! dry mass of first layer
+dmass = vol * dens
+rmass1 = dmass(1) - (C_hum(1) + C_opm(1)) / cpart    ! corection term of first layer
+
+humush = humusv
+if (2.*C_hum(1) .lt. humusv(1)*dmass(1)) then
+    humusv(1) = C_hum(1) / (dmass(1) * cpart)
+endif
+do i=2, nlay
+    humusv(i) = C_hum(i) / (dmass(i) * cpart)
+enddo
+
+if (flag_bc .gt. 0) y_bc_n = 1      ! actual number of biochar application
+
+! calculation of additions for water capacities 
+call hum_add(xfcap, xwiltp, xpv)
+
+fcaph  = f_cap_v - xfcap
+wiltph = wilt_p_v - xwiltp
+pvh    = pv_v - xpv
+
+! ground water
+do i = nlgrw, nlay
+    wats(i) = pv(i)
+enddo
+
+interc_can  = 0.
+int_st_can  = 0.
+interc_sveg = 0.
+int_st_sveg = 0.
+aev_i       = 0.
+
+
+
+
+wat_tot = SUM(wats)
+END subroutine soil_ini
+
+!**************************************************************
+
+SUBROUTINE soil_wat
+use data_out
+! soil water balance
+
+use data_climate
+use data_evapo
+use data_inter
+use data_out
+use data_par 
+use data_simul
+use data_soil
+use data_species
+use data_stand
+
+implicit none
+ 
+real :: eva_dem		! evaporation demand of soil
+real :: p_inf = 0.  ! infiltrated water
+real :: pev		    ! local: soil evaporation
+real :: watot, wetot    ! total water content at start and end
+real :: wutot           ! total water uptake from the soil
+real :: wutot_ev        ! total water uptake by soil evaporation
+real :: wutot_r         ! total water uptake by roots
+real, allocatable, dimension(:) :: upt   ! local array for uptake
+real, external :: wat_new
+
+real enr, wa, we, percolnl, snow_sm, buckdepth
+integer j
+
+allocate (upt(nlay))
+wupt_ev = 0.
+aev_s = 0.
+prec_stand = MAX(0., prec - interc_can - interc_sveg)    ! stand precipitation	
+
+if (flag_int .gt. 1000) then
+    prec_stand = prec_stand * prec_stand_red / 100.
+endif
+
+call snowpack(snow_sm, p_inf, pev)
+
+if (anz_coh .le. 0) pev = pet
+eva_dem = MAX(0., p_inf) - pev		! evaporation demand of soil
+
+!     if all stand precipitation is evaporated and there is still a demand
+!     there is an uptake from soil layers (up to an certain depth)
+if (eva_dem .lt. 0.) then
+   if (snow .le. 0) call take_wat(eva_dem, cover) 
+   aev_s = aev_s + p_inf + SUM(wupt_ev)
+else
+   aev_s = aev_s + pev
+endif	! eva_dem
+
+aet   = aev_s + aev_i
+p_inf = MAX(eva_dem, 0.)
+upt   = wupt_ev
+watot = SUM(wats)			! total initial water content
+
+do j = 1, nlgrw-1
+
+  enr   = p_inf - upt(j)
+  wa    = wats(j) - field_cap(j)
+  we    = wat_new(wa, enr, j)
+  p_inf = enr + wa - we
+
+  perc(j) = MAX(p_inf, 0.)
+  wats(j) = MAX(we+field_cap(j), wilt_p(j))
+
+enddo
+
+do j = nlgrw, nlay
+
+  enr   = p_inf - upt(j)
+  wa    = wats(j) - field_cap(j)
+  we    = pv(j) - field_cap(j)   ! ground water level is constant!
+  p_inf = enr + wa - we
+  perc(j) = MAX(p_inf, 0.)
+  wats(j) = MAX(we+field_cap(j), wilt_p(j))
+
+enddo
+
+if (flag_wred .le. 10) then
+   call upt_wat
+else
+   call upt_wat1
+endif
+
+! root uptake balanced imediate after calculation
+upt  = upt + wupt_r
+wats = wats - wupt_r
+
+watvol = 10.*wats/(thick * skelfact)   ! estimation for complete layer in Vol% without skeleton (only soil substrate)
+
+! total water quantities
+wetot    = SUM(wats)		! total final water content
+wutot_ev = SUM(wupt_ev)		! total water uptake by soil evaporation
+wutot_r  = SUM(wupt_r)		! total water uptake by roots
+wutot    = wutot_ev + wutot_r	! total water uptake
+aet      = aet + wutot_r        ! daily total aet
+
+percolnl = perc(nlay)
+
+trans_tree = 0.
+trans_sveg = 0.
+   zeig => pt%first
+   do while (associated(zeig))
+     if (zeig%coh%species .le. nspec_tree) then
+        trans_tree = trans_tree + zeig%coh%supply 
+     else
+        trans_sveg = trans_sveg + zeig%coh%supply 
+     endif
+     zeig => zeig%next
+   enddo  ! zeig (cohorts) 
+
+! cumulative water quantities
+perc_cum   = perc_cum + perc(nlay)
+wupt_r_c   = wupt_r_c + wutot_r
+wupt_e_c   = wupt_e_c + wutot_ev
+wupt_cum   = wupt_cum + wutot
+aet_cum    = aet_cum + aet
+dew_cum    = dew_cum + dew_rime
+tra_tr_cum = tra_tr_cum + trans_tree 
+tra_sv_cum = tra_sv_cum + trans_sveg
+
+call bucket(bucks_100, bucks_root, buckdepth)
+
+! number of drought days per layer
+where ((wats-0.2) .le. wilt_p) s_drought = s_drought+1
+
+ if (flag_dayout .ge. 2) then
+   write (unit_wat, '(2I5, 7F7.2, 24F8.2)') time_cur, iday, airtemp, prec, interc_can, int_st_can, &
+                    interc_sveg, int_st_sveg, snow, snow_sm,  pet, trans_dem, &
+                    pev, aev_s, aev_i, perc(nlay), watot, wetot, wutot, wutot_ev, wutot_r,&
+                    trans_tree,trans_sveg, eva_dem, gp_can, aet, ceppot_can, ceppot_sveg
+endif
+deallocate (upt)
+		    
+END	subroutine soil_wat
+
+!**************************************************************
+SUBROUTINE upt_wat
+
+! Water uptake by roots
+
+use data_simul
+use data_evapo
+use data_soil
+use data_stand 
+use data_par
+use data_species
+use data_climate
+
+implicit none
+
+real, dimension(1:anz_coh) :: tr_dem   ! auxiliary arrays for cohorts
+real wat_ava, hdem, frtrel, frtrel_1, hupt, hupt_c, totfrt_2, hv, demand_mistletoe_canopy
+real wat_at    ! total available water per layer
+real wat_ar    ! total available water per layer with uptake resistance
+real hred      ! resistance coefficient
+real, external :: fred1, fred2, fred3, fred4, fred5, fred6, fred7, fred11
+integer i, ianz, j, nroot3
+
+
+! Calculation of Water Demand
+ianz = anz_coh
+
+tr_dem = 0
+hdem   = 0
+hv     = pet-aev_i-pev_s
+if (hv .lt. 0.) hv = 0.
+
+select case (flag_eva)
+case (0,1,3)
+ if((pet .gt. 0.) .and. (hv .gt. 0.)) then
+    trans_dem = hv * alfm * (1. - exp(-gp_tot/gpmax))  ! pet (potential evapotranspiration) is reduced by intereption evaporation and potential ground evaporation
+ else
+ trans_dem = 0.0 
+ endif 
+   if (gp_tot .gt. zero) then
+      hdem   = trans_dem / gp_tot
+   else
+      hdem= 0.
+   endif
+case (8)  
+ ! potential transpiration demand of each cohort
+   if ((gp_tot .gt. zero) .and. (hv .gt. 0.)) then
+      hdem   = (pet-aev_i-aev_s) / gp_tot
+   else
+      hdem= 0.
+   endif
+
+case (2,4)
+   trans_dem = 0.
+
+case (6,16,36)
+ ! Eucalyptus
+ hv = pet
+
+ if((pet .gt. 0.) .and. (hv .gt. 0.)) then
+    trans_dem = hv * alfm * (1. - exp(-gp_tot/gpmax)) 
+ else
+    trans_dem = 0.
+ endif 
+
+ ! preparation: potential transpiration demand of each cohort
+   if (gp_tot .gt. zero) then
+      hdem   = trans_dem / gp_tot
+   else
+      hdem= 0.
+   endif
+
+case (7,17,37)
+ trans_dem = hv
+
+ ! potential transpiration demand of each cohort
+   if (gp_tot .gt. zero) then
+      hdem   = trans_dem / gp_tot
+   else
+      hdem= 0.
+   endif
+
+end select
+
+hdem = max(0., hdem)
+
+! Distribution of total Demand into Demands of Cohorts
+
+!extraction of demand of mistletoe cohort (for case flag eva = 1,3,6,7...)
+zeig => pt%first
+ do while (associated(zeig))
+  if (zeig%coh%species.eq.nspec_tree+2) then
+      demand_mistletoe_canopy=zeig%coh%gp * zeig%coh%ntreea * hdem
+  end if
+  zeig => zeig%next
+ enddo
+zeig => pt%first
+i = 1
+do while (associated(zeig))
+
+  select case (flag_eva)
+  case (0, 1, 3, 6, 7, 16, 17, 36, 37)
+
+            !uppermost tree cohort (with flag mistletoe) gets additinal demand of mistletoe
+            if (zeig%coh%mistletoe.eq.1) then
+                zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem + demand_mistletoe_canopy
+            elseif (zeig%coh%species.eq.nspec_tree+2) then                ! set demand of mistletoe to zero as it will be fullfilled by the tree
+               zeig%coh%demand=0.                                         ! set to zero because demand has been added to the infested tree cohort
+            else
+                zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem    ! all other cohorts get their demand
+            end if
+
+  case (2,4)
+ !uppermost tree cohort (with flag mistletoe) gets additinal demand, that of mistletoe
+            if (zeig%coh%mistletoe.eq.1) then
+                zeig%coh%demand = (max(0., zeig%coh%demand - zeig%coh%aev_i)  + demand_mistletoe_cohort)
+            endif
+            if (zeig%coh%species.eq.nspec_tree+2) then                ! set demand of mistletoe to zero as it will be fullfilled by the tree
+                zeig%coh%demand=0.
+            endif
+            if (zeig%coh%mistletoe.ne.1 .AND. zeig%coh%species.ne.nspec_tree+2) then
+               zeig%coh%demand = max(0., zeig%coh%demand - zeig%coh%aev_i)   
+            end if
+            trans_dem = trans_dem + zeig%coh%demand
+    end select
+
+  tr_dem(i) = zeig%coh%demand   ! demand of transpiration per cohort 
+   i = i + 1
+  zeig => zeig%next
+enddo  ! zeig (cohorts)  
+
+! Calculation of Water Supply
+frtrel_1 = 1.
+
+select case (flag_wurz) 
+case (0)
+  if (nroot_max .gt. 5) then 
+     nroot3 = 5
+  else
+     nroot3 = nroot_max
+  endif
+case default
+  nroot3=nroot_max
+end select  
+
+
+
+
+! layers with seedlings   
+do j = 1,nroot3
+   ! determination of resisctance coefficient
+   
+   select case (flag_wred)
+   case(1)
+      hred = fred1(j)
+   case(2)
+      hred = fred2(j)
+   case(3)
+      hred = fred3(j)
+   case(4)
+      hred = fred4(j)
+   case(5)
+      hred = 1.
+   case(6)
+      hred = 0.5
+   case(7)
+      hred = 0.25
+   case(8)  
+      hred = fred6(j)
+   case(10)
+      hred = fred7(j)
+   case(11)
+      hred = fred11(j)
+   end select
+   
+   wat_res(j) = hred
+
+   if (temps(j) .gt. -0.3)  then
+      wat_at = max(wats(j) - wilt_p(j), 0.)     ! total available water per layer
+      wat_ar = hred * wat_at                    ! total available water per layer with uptake resistance
+	  hupt = 0.
+   else
+      wat_ar = 0.                              ! frost
+      wat_at = 0.
+      hupt   = 0.
+   endif
+
+
+
+! Distribution of Water Supply into the Cohorts
+! Distribution of Fine Roots
+  
+   zeig => pt%first
+   i = 1   ! cohort index
+   do while (associated(zeig))
+     if (zeig%coh%species .ne. nspec_tree+2) then     ! not for mistletoe
+      frtrel = zeig%coh%frtrelc(j)
+     	
+     wat_ava = frtrel * wat_ar        ! available water per tree cohort and layer 
+
+    	
+	 if (wat_ava .ge. tr_dem(i)) then
+         hupt_c    = tr_dem(i)
+         tr_dem(i) = 0.
+     else
+         hupt_c    = wat_ava
+         tr_dem(i) = tr_dem(i) - wat_ava
+     endif
+	     
+     select case (flag_dis)    !xylem clogger disturbance 
+     case (1,2)
+     hupt_c = hupt_c * xylem_dis
+     end select
+	 
+     xwatupt(i,j) = hupt_c   ! water uptake per cohorte and layer
+     zeig%coh%supply = zeig%coh%supply + hupt_c
+     if (zeig%coh%supply .lt.0.) then
+     continue
+     endif
+     hupt = hupt + hupt_c
+  
+
+	 i = i + 1
+    end if              ! exclusion of mistletoe
+     zeig => zeig%next
+   enddo  ! zeig (cohorts) 
+
+   wupt_r(j) = hupt 
+
+enddo    ! j
+
+! layers without seedlings   
+
+if (totfrt_p.gt.(seedlfrt+zero)) then
+ totfrt_2 = 1./(totfrt_p-seedlfrt)
+
+  do j = nroot3+1, nroot_max
+   ! determination of resisctance coefficient
+   select case (flag_wred)
+   case(1)
+      hred = fred1(j)
+   case(2)
+      hred = fred2(j)
+   case(3)
+      hred = fred3(j)
+   case(4)
+      hred = fred4(j)
+   case(5)
+      hred = 1.
+   case(6)
+      hred = 0.5
+   case(7)
+      hred = 0.25
+   end select
+   
+   wat_res(j) = hred
+   
+   if (temps(j) .gt. -0.3)  then
+      wat_at = max(wats(j) - wilt_p(j), 0.)     ! total available water per layer
+      wat_ar = hred * wat_at                    ! total available water per layer with uptake resistance
+      hupt   = 0.
+   else
+      wat_ar = 0.
+   endif
+   
+   zeig => pt%first
+   i = 1   ! cohort index
+   do while (associated(zeig))
+     frtrel = zeig%coh%frtrelc(j)		
+     wat_ava = frtrel * wat_ar        ! available water per tree cohort and layer 
+     if (wat_ava .ge. tr_dem(i)) then
+         hupt_c    = tr_dem(i)
+         tr_dem(i) = 0.
+     else
+         hupt_c    = wat_ava
+         tr_dem(i) = tr_dem(i) - wat_ava
+     endif
+       
+     select case (flag_dis)    !xylem clogger disturbance 
+     case (1,2)
+     hupt_c = hupt_c * xylem_dis
+     end select
+	 
+	 xwatupt(i,j) = hupt_c
+     zeig%coh%supply = zeig%coh%supply + hupt_c
+     hupt = hupt + hupt_c
+     i = i + 1
+     zeig => zeig%next
+   enddo  ! zeig (cohorts) 
+
+   wupt_r(j) = hupt 
+enddo     ! j
+endif
+END	subroutine upt_wat
+
+!**************************************************************
+
+SUBROUTINE upt_wat1
+
+! Water uptake by roots
+! 2. Version
+
+use data_simul
+use data_evapo
+use data_soil
+use data_stand 
+use data_par
+
+implicit none
+
+real, dimension(1:anz_coh) :: tr_dem,frt_rel   ! help arrays for cohorts
+real wat_ava, hdem, frtrel, frtrel_1, hupt, hupt_c, totfrt_2
+real wat_at    ! total available water per layer
+real wat_ar    ! total available water per layer with uptake resistance
+real hred      ! resistance coefficient
+real, external :: fred1, fred2, fred3, fred4, fred5, fred6, fred7, fred11
+
+integer i, ianz, j, nroot3
+
+ianz = anz_coh
+tr_dem=0
+trans_dem = (pet-aev_i) * alfm * (1. - exp(-gp_can/gpmax))  ! pet NOT reduced by ground evaporation
+if (trans_dem .lt. 0.) trans_dem = 0.
+
+! potential transpiration demand of each cohort
+if (gp_can .gt. zero) then
+   hdem   = trans_dem / gp_can
+else
+   hdem= 0.
+endif
+
+! Estimation of transpiration demand of tree cohorts and total fine root mass
+! in layers with and without seedlings
+zeig => pt%first
+i = 1
+do while (associated(zeig))
+  select case (flag_eva)
+  case (0, 1, 3)
+    zeig%coh%demand = zeig%coh%gp * zeig%coh%ntreea * hdem  
+  case (2)
+    zeig%coh%demand = zeig%coh%demand - zeig%coh%aev_i
+  end select
+  tr_dem(i) = zeig%coh%demand
+  i = i + 1
+  zeig => zeig%next
+enddo  ! zeig (cohorts)  
+
+! uptake controlled by share of roots 
+frtrel_1 = 1.
+
+! layers with seedlings   
+do j = 1,nroot_max
+   ! determination of resisctance coefficient
+   select case (flag_wred)
+   case(1)
+      hred = fred1(j)
+   case(2)
+      hred = fred2(j)
+   case(3)
+      hred = fred3(j)
+   case(4)
+      hred = fred4(j)
+   case(5)
+      hred = 1.
+   case(6)
+      hred = 0.5
+   case(7)
+      hred = 0.25
+   case(8)  ! BKL, ArcEGMO
+      hred = fred6(j)
+   case(10)
+      hred = fred7(j)
+   end select
+
+   wat_at = max(wats(j) - wilt_p(j), 0.)     ! total available water per layer
+   wat_ar = hred * wat_at                    ! total available water per layer with uptake resistance
+   hupt   = 0.
+   
+   zeig => pt%first
+   i = 1   ! cohort index
+   do while (associated(zeig))
+     
+     frtrel = zeig%coh%frtrel(j) * zeig%coh%x_frt * zeig%coh%ntreea * totfrt_1    
+     wat_ava = frtrel * wat_ar        ! available water per tree cohort and layer 
+     if (wat_ava .ge. tr_dem(i)) then
+         hupt_c    = tr_dem(i)
+         tr_dem(i) = 0.
+     else
+         hupt_c    = wat_ava
+         tr_dem(i) = tr_dem(i) - wat_ava
+     endif
+	     
+     select case (flag_dis)    !xylem clogger disturbance 
+     case (1,2)
+     hupt_c = hupt_c * xylem_dis
+     end select
+	 
+     xwatupt(i,j) = hupt_c
+     zeig%coh%supply = zeig%coh%supply + hupt_c
+     hupt = hupt + hupt_c
+     i = i + 1
+     zeig => zeig%next
+   enddo  ! zeig (cohorts) 
+
+   wupt_r(j) = hupt 
+enddo    ! j
+END	subroutine upt_wat1
+
+!**************************************************************
+
+real FUNCTION fred1(j)
+
+! Function for calculating uptake resistance
+! from CHEN (1993)
+! empirical relation between soil water content and resistance
+! fred1=1 if (field_cap - 10%*field_cap) <= wats <= (field_cap + 10%*field_cap)
+
+use data_par
+use data_soil
+
+implicit none
+real hf, f09, f11, wc, diff
+integer j
+
+f09 = 0.9 * field_cap(j)
+f11 = 1.1 * field_cap(j)
+wc  = wats(j)
+
+if (wc .lt. wilt_p(j)) then
+   hf = 0.
+else if (wc .lt. f09) then
+   diff = f09-wilt_p(j)
+   if (diff .lt. zero) diff = 0.001
+   hf = 1. - (f09-wc) / diff
+else if (wc .gt. f11) then
+   diff = pv(j)-f11
+   if (diff .lt. zero) diff = 0.001
+   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
+   if (hf .lt. zero) hf = 0.001
+else
+   hf = 1.
+endif
+fred1 = hf
+END	function fred1
+
+!**************************************************************
+ 
+real FUNCTION fred2(j)
+
+! Function for calculating uptake resistance
+! from Aber and Federer (f=0.04 fuer Wasser in cm)
+! only 40% of total available water are plant available per day
+
+implicit none
+
+integer j
+
+fred2 = 0.05 
+
+END	function fred2
+
+!**************************************************************
+ 
+real FUNCTION fred3(j)
+
+! Function for calculating uptake resistance
+! from CHEN (1993)
+! modified to a profile defined in fred 
+! fred3 may be described by a function (old version): 
+!     fred3 = 0.0004*j*j - 0.0107*j + 0.0735
+! this case: set from a root profile, defined by input of root_fr
+
+use data_par
+use data_soil
+
+implicit none
+real hf, f09, f11, wc, diff
+! hf is a reduction factor in dependence on water content
+real fred(15)
+
+integer j
+
+! uptake reduction depending on water content
+f09 = 0.9 * field_cap(j)
+f11 = 1.1 * field_cap(j)
+wc  = wats(j)
+
+if (wc .lt. wilt_p(j)) then
+   hf = 0.
+else if (wc .lt. f09) then
+   diff = f09-wilt_p(j)
+   if (diff .lt. zero) diff = 0.001
+   hf = 1. - (f09-wc) / diff
+else if (wc .gt. f11) then
+   diff = pv(j)-f11
+   if (diff .lt. zero) diff = 0.001
+   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
+   if (hf .lt. zero) hf = 0.001
+else
+   hf = 1.
+endif
+
+fred3 = root_fr(j) * hf 
+
+END	function fred3
+
+!**************************************************************
+ 
+real FUNCTION fred4(j)
+
+! Function for calculating uptake resistance
+! modified to a profile defined in fred 
+! profile at Beerenbusch
+
+use data_soil
+
+implicit none
+real fred(15)
+
+integer j
+
+fred = (/ 0.0, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01,    &
+              0.01, 0.01, 0.01, 0.01, 0.01 /)  ! fred fuer Beerenbusch
+ 
+fred4 = fred(j) 
+
+END	function fred4
+
+!**************************************************************
+
+real FUNCTION fred6(j)
+
+! Function for calculating uptake resistance
+! from Kloecking (2006) simular to fred1
+! empirical relation between soil water content and resistance
+! fred6=1 if field_cap  <= wats <= (field_cap + 10%*field_cap)
+
+use data_soil
+
+implicit none
+real hf, f09, f11, wc
+integer j
+
+f09 = field_cap(j)
+f11 = 1.1 * field_cap(j)
+wc  = wats(j)
+
+if (wc .le. wilt_p(j)) then
+   hf = 0.
+else if (wc .lt. f09) then
+   hf = 0.1 + (0.9 *(wc-wilt_p(j)) / (f09-wilt_p(j)))
+else if (wc .gt. f11) then
+   hf = 0.3 + 0.7 * (pv(j)-wc) / (pv(j)-f11)
+   if (hf .lt. 0.) hf = 0.001
+else
+   hf = 1.
+endif
+
+fred6 = hf
+
+END	function fred6
+
+!**************************************************************
+
+real FUNCTION fred7(j)
+
+! Function for calculating uptake resistance
+! from CHEN (1993)
+! empirical relation between soil water content and resistance
+! fred1=1 if (field_cap - 10%*field_cap) <= wats <= (field_cap + 10%*field_cap)
+
+use data_par
+use data_soil
+
+implicit none
+real hf, f09, f11, wc, diff
+integer j
+
+f09 = 0.9 * field_cap(j)
+f11 = 1.1 * field_cap(j)
+wc  = wats(j)
+
+if (wc .lt. wilt_p(j)) then
+   hf = 0.
+else if (wc .lt. f09) then
+   diff = f09-wilt_p(j)
+   if (diff .lt. zero) diff = 0.001
+   hf = exp(-5.*(f09-wc) / diff)
+else if (wc .gt. f11) then
+   diff = pv(j)-f11
+   if (diff .lt. zero) diff = 0.001
+   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
+   if (hf .lt. zero) hf = 0.001
+else
+   hf = 1.
+endif
+
+fred7 = hf
+
+END	function fred7
+
+!**************************************************************
+
+real FUNCTION fred11(j)
+
+! Function for calculating uptake resistance, especially adapted for Mistletoe disturbance
+! function after van Wijk, 2000
+
+use data_par
+use data_soil
+implicit none
+real hf, S, f11, wc, diff
+integer j
+f11 = 1.1 * field_cap(j)
+wc  = wats(j)
+if (wc .lt. wilt_p(j)) then
+   hf = 0.
+else if (wc .lt. field_cap(j)) then
+   S=(field_cap(j)-wc)/(field_cap(j)-wilt_p(j))
+   hf = exp(-30*S)                                !30 = strong reduction in water avail.
+else if (wc .gt. f11) then
+   diff = pv(j)-f11
+   if (diff .lt. zero) diff = 0.001
+   hf = 0.3 + 0.7 * (pv(j)-wc) / diff
+   if (hf .lt. zero) hf = 0.001
+else
+   hf = 1.
+endif
+fred11 = hf
+END function fred11
+
+!**************************************************************
+
+ 
+
+
+
+
+
+ 
+ SUBROUTINE take_wat(eva_dem, psi)   
+      
+! Estimation of water taking out for uncovered soil 
+use data_soil
+use data_simul
+
+implicit none
+
+!input:
+real :: eva_dem	   ! evaporation demand
+real :: psi        ! covering
+
+integer i, ii, j, ntag  ! max. layer of taking out
+real, allocatable, dimension(:) :: gj
+real, external :: b_r, funcov
+real diff, gj_j, depth_j, depth_n, rij, rmax, rr, rs, sr
+
+allocate (gj(nlay))
+
+do i=1,nlay
+   wupt_ev(i)=0.0
+   gj(i)=0.0 
+enddo  
+      
+ntag    = 0
+rmax    = 0.0
+depth_n = depth(n_ev_d)
+	
+do i=1,n_ev_d
+  rij = 0.0
+  rr  = depth_n/depth(i)
+  sr  = 0.0
+  rs  = 0.0
+	    
+  do j=1,i
+!         depth for uncovered take out
+     depth_j = depth(j)
+     gj(j)   = FUNCOV(w_ev_d, rs*rr, rr*depth_j)
+     rs      = depth_j
+     sr      = sr + gj(j)
+  enddo		! i
+	  
+  if (sr.gt.1.E-7) then
+     sr = 1.0/sr
+	    
+     do j=1,i
+!   	   water take out
+!         (psi = 1.-psi) no soil evaporation in case of total covering 
+!          and maximal evaporation for uncovered soil
+        gj_j = -B_R(wats(j), field_cap(j), wilt_p(j))    &
+               * eva_dem * (1.-psi) * gj(j) * sr
+        gj_j = max(gj_j,0.0)
+        gj(j)= gj_j
+        rij  = rij + gj_j
+     enddo		! i
+
+     if (rij .gt. rmax) then
+         rmax = rij
+         ntag = i
+
+         do ii=1,ntag
+            wupt_ev(ii) = gj(ii)
+         enddo	
+
+     endif	! rij 
+   endif	! sr	  
+enddo	! n_ev_d
+
+! balance
+do i=1,nlay
+   diff = wats(i) - wilt_p(i)
+   if (wupt_ev(i) .gt. diff) then
+       wupt_ev(i) = diff
+   endif
+enddo ! nlay
+
+deallocate (gj)
+
+END  subroutine  take_wat
+
+!*******************************************************************************
+   
+real FUNCTION B_R(water, f_cap, wilting)
+
+! Reduction function for water taking out (uncovered soil)
+
+implicit none
+
+!input:
+real :: water	    ! water storage
+real :: f_cap	    ! field capacity
+real :: wilting	    ! wilting point
+
+b_r  = 1.0
+
+if (water .lt. f_cap) B_R = max((water-wilting)/(f_cap-wilting), 0.0)
+
+END  function B_R
+
+!******************************************************************************
+
+real FUNCTION funcov(wt_d, a, bb)
+
+! take out density function for uncovered soil
+
+implicit none
+
+!input:
+real :: wt_d    ! depth of water taking out by evaporation (cm)
+real :: a, bb   ! relative upper and lower depth of actual layer
+real fk, wt_5, b
+
+      fk   = .455218234
+      wt_5   = 0.05 * wt_d
+      b      = min(bb,wt_d)
+      funcov = (- b + a + 1.05*wt_d*log((b+wt_5)/(a+wt_5)))*fk/wt_d
+      
+END function funcov
+
+!******************************************************************************
+
+real FUNCTION wat_new(wat_us, wat_in, ilayer)
+!      FUNCTION WIEN(WIA,NIST,ALAM,DTI,TT,DICK)
+
+! Estimation of additional water after infiltration and percolation
+
+use data_par
+use data_soil
+
+implicit none
+
+! input:
+real :: wat_us	   ! water content in relation to field capacity
+real :: wat_in	   ! water infiltration into actual layer
+integer :: ilayer  ! number of actual layer
+real dti	 !time step
+real awi, b1, b2, la, hsqr, exphelp
+
+dti = 1.
+fakt = 0.4
+
+if (fakt .ge. 0.0) then				! percolation?
+   la = 100.0 * fakt * dti * wlam(ilayer)/thick(ilayer)**2
+   if (wat_us .le. zero) then		! water near zero?
+      if (wat_in .le. zero) then	! infiltrated water near zero?
+         wat_new = wat_us + wat_in
+      else
+         if (wat_us+wat_in .gt. zero) then
+            exphelp = sqrt(la*wat_in) * (1 + wat_us/wat_in)*1
+            if (exphelp .le.10.) then     ! avoid underflow   
+               b1 = -exp(-2. * exphelp)
+            else
+               b1 = 0.
+            endif
+            wat_new = sqrt(wat_in/la) * (1+b1)/(1-b1)
+		 else
+		    wat_new = wat_us + wat_in
+		 endif
+	  endif	! wat_in
+
+   else 
+        if (wat_in .lt. 0.) then
+            awi = abs(wat_in)
+            b1 = atan(wat_us/sqrt(awi/la)) / sqrt(la * awi)
+			if (b1 .gt. 1) then
+			   b2 = sqrt (awi * la)
+			   b1 = sin(b2) / cos(b2)
+			   b2 = sqrt(awi / la)
+			   wat_new = b2 * (wat_us - b2*b1) / (b2 + wat_us*b1)
+			else
+			   wat_new = wat_in * (1-b1)
+			endif	! b1
+        else 
+             if (wat_in .gt. 0.) then
+                b1 = sqrt(wat_in / la)
+                hsqr = sqrt(la*wat_in)
+                if (hsqr .lt. 10.) then
+                   b2 = (wat_us - b1) * exp(-2.* hsqr) / (wat_us + b1)
+                   if (b2 .ge. 1.0) then
+                       b2 = 0.99999
+                   endif
+                else
+                   b2 = 0.
+                endif
+                wat_new = b1 * (1.+b2) / (1.-b2)
+             else
+                wat_new = wat_us / (1. + la*wat_us)
+             endif	
+       endif	! wat_in		   
+   endif	! wat_us
+else
+    wat_new = wat_us
+endif		! fakt
+
+END  function wat_new
+ 
+!******************************************************************************
+
+SUBROUTINE bucket(bucksize1, bucksize2, buckdepth)
+
+! calculation of bucket size (1m; without humus layer)
+
+use data_soil
+
+implicit none
+
+real bucksize1,   &    ! bucket size of 1 m depth  (nFK)
+     bucksize2,   &    ! bucket size of rooting zone
+     buckdepth, diff
+integer j
+
+bucksize1 = 0.
+bucksize2 = 0.
+buckdepth = 0.
+do j=2,nlay
+   if ((depth(j)-depth(1)) .lt. 100.) then
+      bucksize1  = bucksize1 + wats(j) - wilt_p(j)
+	  buckdepth = depth(j) - depth(1)
+   else
+      diff      = 100. - buckdepth
+	  bucksize1  = bucksize1 + (wats(j) - wilt_p(j))*diff/thick(j)
+      buckdepth = 100.
+	  exit
+   endif
+enddo
+
+do j=2,nroot_max
+      bucksize2  = bucksize2 + wats(j) - wilt_p(j)
+enddo
+
+END  subroutine bucket
+
+!******************************************************************************
+
+SUBROUTINE snowpack(snow_sm, p_inf, pev)
+
+! properties of snow 
+! calculation of soil surface temperature under snow pack 
+
+use data_climate
+use data_evapo
+use data_inter
+use data_par
+use data_simul
+use data_soil
+use data_soil_t
+
+implicit none
+
+real p_inf		   ! infiltrated water
+real snow_sm
+real pev
+real airtemp_sm    ! melting temperature 
+real snow_old      ! old snow pack
+real tc_snow       ! thermal conductivity of snow  J/cm/s/K
+real thick_snow    ! thickness of snow
+real dens_snow     ! density of snow
+real:: dens_sn_new  = 0.1 ! density of fresh snow
+real fakta
+
+snow_old = snow
+
+!substract evaporation of snowcover from snow in both cases
+if (airtemp .lt. temp_snow) then        ! frost conditions
+    snow     = snow + prec_stand	    ! precipitation as snow
+    snow_sm  = 0.0          	        ! no snow melting
+    p_inf    = 0.0                      ! no infiltrated precipitation
+    pev	     = max((pev_s - aev_i), 0.) ! interc. evapor. reduces soil evapor.
+
+else
+
+    airtemp_sm = max(airtemp, 0.)
+    snow_sm    = airtemp*(0.45+0.2*airtemp)  ! snow melting
+    snow_sm    = MIN(snow_sm, snow)
+    snow       = snow - snow_sm
+    p_inf      = prec_stand + snow_sm	  ! infiltrated precipitation	
+    pev	       = max((pev_s - aev_i), 0.) ! interc. evapor. reduces soil evapor.
+
+end if		! airtemp
+
+if (snow .ge. zero) then
+   snow_day = snow_day + 1
+   days_snow = days_snow + 1
+   if (pev .le. zero) then
+      pev = 0.	
+   else
+    ! snow sublimation 
+      aev_s = max(min(snow, pev), 0.)   
+      snow  = snow - aev_s
+      pev   = pev - aev_s
+   endif
+
+   ! soil surface temperature under snow pack
+   ! snow hight = 0.2598 * water equivalent + 8.6851; adjustment from measurement values (see Bodentemperatur.xls)
+
+   if (snow .ge. 0.05) then
+      
+      dens_snow = dens_sn_new + snow_day*0.025
+      dens_snow = MIN(dens_snow, 1.)
+      dens_snow = 0.5*(dens_sn_new*prec_stand + dens_snow*snow_old)/snow
+      dens_snow = MIN(dens_snow, 1.)
+      tc_snow   = 0.7938*EXP(3.808*dens_snow)*0.001  ! thermal conductivity of snow  J/cm/s/K
+      thick_snow = snow / dens_snow 
+      fakta      = tc_snow * 86400. * (thick(1)/2.) / (t_cond(1) * thick_snow)   ! s --> day
+      temps_surf  = (0.5*temps(1) + fakta*airtemp) / (1. + fakta)   ! CoupModel (Jansson, 2001)
+   endif   					
+else
+   snow_day = 0
+endif
+END  subroutine snowpack
+
+!******************************************************************************
+
+SUBROUTINE soil_stress 
+
+! Calculation of the stress factors
+
+use data_soil
+use data_species
+use data_stand
+use data_par
+
+implicit none
+integer :: i, k
+
+real :: m_1, m_2, n_1, n_2
+real :: wratio, wafpo
+real, dimension (1:4) ::  allstress, xvar, yvar
+
+!temperature stress
+do i=1,nlay
+ do k=1,nspecies
+    if (temps(i) .ge. spar(k)%tbase) then
+        svar(k)%tstress(i) = sin((pi/2)*(temps(i)-spar(k)%tbase)/(spar(k)%topt-spar(k)%tbase))
+   	else
+	    svar(k)%tstress(i) = 0.
+    endif
+
+ !soil strength
+    wratio=0.    
+    if (dens(i) .le. BDopt(i)) then
+		svar(k)%BDstr(i) = 1
+		svar(k)%BDstr(i) = 1
+	elseif (dens(i) .ge. svar(k)%BDmax(i)) then
+		svar(k)%BDstr(i) = 0
+	else
+		svar(k)%BDstr(i) = (svar(k)%BDmax(i)-dens(i))/(svar(k)%BDmax(i)-BDopt(i))
+	endif
+
+	if (watvol(i) .lt. wilt_p_v(i)) then
+	wratio = 0.
+	elseif (watvol(i) .gt. f_cap_v(i)) then
+	wratio = 1.
+	else
+	wratio = (watvol(i)-wilt_p_v(i))/(f_cap_v(i)-wilt_p_v(i))
+	endif
+	
+	svar(k)%sstr(i)=svar(k)%BDstr(i)*sin(1.57*wratio)
+ 
+ !aeration
+    wafpo=watvol(i)/pv_v(i)
+	if (wafpo .ge. svar(k)%porcrit(i)) then
+		svar(k)%airstr(i) = (1.-wafpo)/(1.-svar(k)%porcrit(i))
+	else
+		svar(k)%airstr(i) = 1.
+	endif
+	
+    if (svar(k)%airstr(i) .lt. 0.) svar(k)%airstr(i) = 0.
+	
+ !soil acidity
+    xvar=(/spar(k)%ph_min, spar(k)%ph_opt_min, spar(k)%ph_opt_max, spar(k)%ph_max/)
+    yvar=(/0,1,1,0/)
+    m_1=(yvar(1)-yvar(2))/(xvar(1)-xvar(2))
+    n_1=yvar(2)-m_1*xvar(2)
+    m_2=(yvar(3)-yvar(4))/(xvar(3)-xvar(4))
+    n_2=yvar(4)-m_2*xvar(4)
+
+	if (phv(i) .gt.  spar(k)%ph_opt_max .and. phv(i) .le. spar(k)%ph_max ) then
+		svar(k)%phstr(i)=m_2*phv(i)+n_2
+		elseif (phv(i) .lt. spar(k)%ph_opt_min .and. phv(i) .ge. spar(k)%ph_min ) then
+			svar(k)%phstr(i)=m_1*phv(i)+n_1
+			elseif (phv(i) .gt. spar(k)%ph_max .or. phv(i) .lt. spar(k)%ph_min) then
+			svar(k)%phstr(i)=0.
+				else
+				svar(k)%phstr(i)=1.
+	endif
+
+ ! total stress (Rstress) is taken as the largest of the four
+    
+	allstress(1)=svar(k)%tstress(i)
+	allstress(2)=svar(k)%sstr(i)
+	allstress(3)=svar(k)%airstr(i)
+	allstress(4)=svar(k)%phstr(i)
+	
+	svar(k)%Rstress(i)= minval(allstress)
+	svar(k)%Smean(i)=svar(k)%Rstress(i)+svar(k)%Smean(i)
+	
+ enddo
+enddo
+
+END subroutine soil_stress
+
+!*******************************************************************************
+
+SUBROUTINE hum_add(xfcap, xwiltp, xpv) 
+! Soil parameter according to [Kuntze et al., Bodenkunde, 1994], S. 172 
+
+use data_simul
+use data_soil
+use data_soil_cn
+
+implicit none
+integer :: i, k
+real    :: fcapi, clayvi, siltvi, humvi, humvi2, wiltpi, pvi, nfki, hcbc 
+real, dimension(nlay):: xfcap, xwiltp, xpv   ! output of addition mm/dm
+
+    xfcap(1)  = 0.0
+    xwiltp(1) = 0.0
+    xpv(1)    = 0.0
+
+do i = 1, nlay
+    fcapi  = 0.
+    wiltpi = 0.
+    pvi    = 0.
+    clayvi = clayv(i)
+    humvi  = humusv(i)*100.
+    humvi2 = humusv(i)*humusv(i)
+    if (humvi .lt. 15.) then
+        if (clayvi .le. 0.05) then
+            wiltpi = 0.0609 * humvi2 + 0.33 * humvi 
+            pvi    = 0.0436 * humvi2 + 0.631 * humvi
+            nfki   = -0.0009 * humvi2 + 1.171 * humvi
+            fcapi  = nfki + wiltpi
+        
+        else if (clayvi .le. 0.12) then
+            wiltpi = 0.0357 * humvi2 + 0.0762 * humvi 
+            pvi    = 0.0441 * humvi2 + 0.5455 * humvi
+            nfki   = 0.0252 * humvi2 + 0.7462 * humvi
+            fcapi  = nfki + wiltpi
+            
+        else if (clayvi .le. 0.17) then
+            wiltpi = 0.0374 * humvi2 - 0.1777 * humvi 
+            pvi    = 0.0552 * humvi2 + 0.2936 * humvi
+            nfki   = 0.0324 * humvi2 + 0.6243 * humvi
+            fcapi  = nfki + wiltpi
+        
+        else if (clayvi .le. 0.35) then
+            wiltpi = 0.0179 * humvi2 - 0.0385 * humvi 
+            pvi    = 0.0681 * humvi2 + 0.0768 * humvi
+            nfki   = 0.0373 * humvi2 + 0.3617 * humvi
+            fcapi  = nfki + wiltpi
+       
+        else if (clayvi .le. 0.65) then
+            wiltpi = 0.0039 * humvi2 + 0.0254 * humvi 
+            pvi    = 0.0613 * humvi2 + 0.0947 * humvi
+            nfki   = 0.0338 * humvi2 + 0.0904 * humvi
+            fcapi  = nfki + wiltpi
+        
+        else 
+            wiltpi = 0.0 
+            pvi    = 0.0613 * humvi2 + 0.0947 * humvi
+            nfki   = 0.0104 * humvi2 + 0.2853 * humvi
+            fcapi  = nfki + wiltpi
+        
+        endif     
+    else   ! humvi > 15
+       ! organic soils
+        continue
+    endif   ! humvi
+
+    xfcap(i)  = fcapi
+    xwiltp(i) = wiltpi
+    xpv(i)    = pvi
+enddo
+    
+if (flag_bc .gt. 0) then
+    do i = 1, nlay
+        if (C_bc(i) .gt. 0.) then
+            fcapi  = f_cap_v(i)
+            clayvi = clayv(i)
+            siltvi = siltv(i)
+            humvi  = humusv(i)*100.
+            hcbc   = C_bc(i)*100.*100. / (cpart_bc(y_bc_n) * dmass(i))
+            if ((clayvi .le. 0.17) .and. (siltvi .le. 0.5)) then     ! sand
+                fcapi  = 0.0619 * hcbc
+                wiltpi = 0.0375 * hcbc
+                nfki = 7.0   
+            elseif ((clayvi .le. 0.45) .and. (siltvi .gt. 0.17)) then   ! loam
+                fcapi  = 0.015 * hcbc
+                wiltpi = 0.0157 * hcbc
+                nfki = 10.
+            else      ! clay
+                fcapi  = -0.0109 * hcbc
+                wiltpi = -0.0318 * hcbc
+                nfki = 16.
+            endif
+            xfcap(i)  = xfcap(i) + fcapi
+            xwiltp(i) = xwiltp(i) + wiltpi
+       endif       
+    
+    enddo
+endif
+
+END subroutine hum_add
+
+!*******************************************************************************
+
+SUBROUTINE bc_appl 
+
+! application of biochar
+
+use data_out
+use data_simul
+use data_soil
+use data_soil_cn
+
+implicit none
+
+character :: text
+integer :: ios, inunit, j
+logical :: ex
+real    :: hcbc 
+
+    call testfile(valfile(ip),ex)
+    IF (ex .eqv. .true.) then
+      inunit = getunit()
+      ios=0
+      open(inunit,file=valfile(ip),iostat=ios,status='old',action='read')
+      if (.not.flag_mult8910) then
+          print *,'***** Reading application values of biochar from file ',valfile(ip),'...'
+          write (unit_err, *) 'Application values of biochar from file   ',trim(valfile(ip))
+      endif
+
+      do
+         read(inunit,*) text
+         IF(text .ne. '!')then
+            backspace(inunit)
+            exit
+         endif
+      enddo
+      
+      read (inunit,*,iostat=ios) n_appl_bc
+      allocate (C_bc_appl(n_appl_bc))
+      allocate (N_bc_appl(n_appl_bc))
+      allocate (bc_appl_lay(n_appl_bc))
+      allocate (cnv_bc(n_appl_bc))
+      allocate (dens_bc(n_appl_bc))
+      allocate (cpart_bc(n_appl_bc))
+      allocate (y_bc(0 : n_appl_bc + 1))
+      y_bc = 0
+      C_bc_appl = 0.
+      N_bc_appl = 0.
+      do j = 1, n_appl_bc
+          read (inunit,*,iostat=ios) y_bc(j), cpart_bc(j), cnv_bc(j), dens_bc(j)
+          read (inunit,*,iostat=ios) bc_appl_lay(j), C_bc_appl(j)      
+      enddo
+    endif   ! ex
+
+END subroutine bc_appl
+
+!*******************************************************************************