!*****************************************************************!
!*                                                               *!
!*              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
