Skip to content
Snippets Groups Projects
Select Git revision
  • 82320a3f4e348de452568c45ac02ed754c25d278
  • master default protected
  • dev-mahnken
  • cherry-pick-c40afffb
  • cherry-pick-f1d9f851
5 results

readsim.f

Blame
  • Forked from 4C / FORESEE
    206 commits behind the upstream repository.
    Petra Lasch-Born's avatar
    Petra Lasch-Born authored
    70ab13ef
    History
    readsim.f 30.66 KiB
    !*****************************************************************!
    !*                                                               *!
    !*              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
    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 
     
    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_forska
        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
    
    
    !     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