!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*    - WRITESIM:   Write simulation options into file           *!
!*                                                               *!
!*                  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/XXXXXXXXXXXXXXXXXXXXX     *!
!*                                                               *!
!*****************************************************************!

SUBROUTINE writesim(simfile_new)

! read simulation options from file

use data_climate
use data_mess
use data_out
use data_simul
use data_stand
use data_site
use data_tsort

implicit none

logical ex
integer i, ios, ios1, nowunit, nowunit1, k, anzclim, j
real help
character a
character (150) tspec, tname, tval, tsite, tman, ttree, tdepo, tred, tlit,  &
                    pathdir1, pathdir2,pathdir3, pathdir4, pathdir5, climszen, siteall, climall,site_name_all
character(3), dimension(100) :: clim_nam
character(150), dimension(:), allocatable:: site_name_ad
character(150), dimension(:), allocatable:: climfile_ad
character(150), dimension(:), allocatable:: manfile_ad
character(150), dimension(:), allocatable:: treefile_ad
character(150), dimension(:), allocatable:: depofile_ad
character(150), dimension(:), allocatable:: wpmfile_ad
character(10), dimension(1000)   :: climnum
character(50)  istand
character(150)  simfile_new

    nowunit = getunit()
    ios = 0
    
!!! set Filename 

    write (*, *) ' Input name of simfile'
    write (*, *) simfile_new 

    open(nowunit,file=simfile_new,iostat=ios, status='replace')

    write(nowunit,'(I6,A)',iostat=ios) flag_multi	, '  ! Run option'
    write(nowunit,'(I6,A)',iostat=ios) site_nr		, '  ! number of runs'
    write(nowunit,'(A)',iostat=ios) '! *** simulation specifications  **************************************'
    write(nowunit,'(I6,A)',iostat=ios) year			, '  ! number of simulation years'
    write(nowunit,'(I6,A)',iostat=ios) time_b		, '   ! start year for simulation'
    write(nowunit,'(F7.0,A)',iostat=ios) kpatchsize	, '  ! patch size [m²]'
    write(nowunit,'(F7.1,A)',iostat=ios) dz			, '  ! thickness of foliage layers [cm]'
    write(nowunit,'(I6,A)',iostat=ios) ns_pro		, '  ! time step photosynthesis calculations [d]'
    write(nowunit,'(A)',iostat=ios) '! *** choice of model options *****************************************'
    write(nowunit,'(I6,A)',iostat=ios) flag_mort	, '  ! mortality flag (flag_mort)'
    write(nowunit,'(I6,A)',iostat=ios) flag_reg		, '  ! regeneration flag (flag_reg)'
    write(nowunit,'(I6,A)',iostat=ios) flag_forska	, '  ! use FORSKA environmental factors and regeneration (flag_forska)'
    write(nowunit,'(I6,A)',iostat=ios) flag_stand	, '  ! initialization flag (flag_stand)'
    write(nowunit,'(I6,A)',iostat=ios) flag_sveg	, '  ! soil vegetation flag (flag_sveg)    !!! new !!!'
    write(nowunit,'(I6,A)',iostat=ios) flag_mg		, '  ! management flag (flag_mg)'
    write(nowunit,'(I6,A)',iostat=ios) flag_dis		, '  ! disturbance flag (flag_dis)'
    write(nowunit,'(I6,A)',iostat=ios) flag_light	, '  ! ligth algorithm number (flag_light)'
    write(nowunit,'(I6,A)',iostat=ios) flag_folhei	, '  ! foliage-height relationship (flag_folhei)'
    write(nowunit,'(I6,A)',iostat=ios) flag_volfunc	, '  ! volume function (flag_volfunc)'
    write(nowunit,'(I6,A)',iostat=ios) flag_resp	, '  ! respiration flag (flag_resp)'
    write(nowunit,'(I6,A)',iostat=ios) flag_limi	, '  ! limitation flag (flag_limi)'
    write(nowunit,'(I6,A)',iostat=ios) flag_decomp	, '  ! decomposition model (flag_decomp)'
    write(nowunit,'(I6,A)',iostat=ios) flag_sign	, '  ! root activity function flag (flag_sign)'
    write(nowunit,'(I6,A)',iostat=ios) flag_wred	, '  ! soil water uptake flag (flag_wred)'
    write(nowunit,'(I6,A)',iostat=ios) flag_wurz	, '  ! root distribution flag (flag_wurz)'
    write(nowunit,'(I6,A)',iostat=ios) flag_cond	, '  ! heat conductance flag (flag_cond)'
    write(nowunit,'(I6,A)',iostat=ios) flag_int		, '  ! interception flag (flag_int)'
    write(nowunit,'(I6,A)',iostat=ios) flag_eva		, '  ! evapotranspiration flag (flag_eva)'
    write(nowunit,'(I6,A)',iostat=ios) flag_co2		, '  ! CO2 flag (flag_CO2)'
    write(nowunit,'(I6,A)',iostat=ios) flag_sort	, '  ! sort flag (flag_sort)'
    write(nowunit,'(I6,A)',iostat=ios) flag_wpm		, '  ! wpm flag (flag_wpm)'
    write(nowunit,'(I6,A)',iostat=ios) flag_stat	, '  ! comparison with measurements (flag_stat)'
    write(nowunit,'(A)',iostat=ios) '! *** output specifications *******************************************'
    write(nowunit,'(I6,A)',iostat=ios) time_out
!     write name of yearly output variables
    do i = 1, outy_n
       if (outy(i)%out_flag .gt. 0) write(nowunit,'(A)',iostat=ios) outy(i)%kind_name
    enddo
    write(nowunit,'(A)',iostat=ios)  'end'

    write(nowunit,'(I6,A)',iostat=ios) flag_dayout
!     write name of daily output variables
    do i = 1, outd_n
       if (outd(i)%out_flag .gt. 0) write(nowunit,'(A)',iostat=ios) outd(i)%kind_name
    enddo
    write(nowunit,'(A)',iostat=ios)  'end'

    if(flag_cohoutd .gt. 0 .or. flag_cohouty .gt. 0) then 
       flag_cohout = 1
    else
       flag_cohout = 0
    endif
    write(nowunit,'(I6,A)',iostat=ios) flag_cohout
!     define name of cohort output variables
    ncvar = ncvar + ncdvar
    do i = 1, outcy_n
       if (outcy(i)%out_flag .gt. 0) write(nowunit,'(A)',iostat=ios) outcy(i)%kind_name
    enddo
    do i = 1, outcd_n
       if (outcd(i)%out_flag .gt. 0) write(nowunit,'(A)',iostat=ios) outcd(i)%kind_name
    enddo
    write(nowunit,'(A)',iostat=ios)  'end'

    write(nowunit,'(I6,A)',iostat=ios) flag_sum

    write(nowunit,'(A)',iostat=ios) '! *** input files *****************************************************'

 SELECT CASE(flag_multi)
 CASE (0,1,2,3,6)
      jpar = 1
      DO i=1,site_nr
        if(i .gt. 1) then
        write(nowunit,'(A,I2,A)',iostat=ios) '! ******************* run ',i,' *******************************************'

          do while (vpar(jpar) .gt. -99.0)
             write(nowunit,'(F7.1, A)') vpar(jpar), '  '//simpar(jpar)
             jpar = jpar + 1
          enddo
          help = -99.0
          write(nowunit,'(F7.1, A)') help, '  end'
        endif

        write(nowunit,'(A)',iostat=ios) specfile(i)
        write(nowunit,'(A)') site_name(i)
        write(nowunit,'(A)') climfile(i)
        write(nowunit,'(A)') sitefile(i)
        write(nowunit,'(A)') valfile(i)
        write(nowunit,'(A)') treefile(i)
        write(nowunit,*) standid(i)
        write(nowunit,'(A)') manfile(i)
        write(nowunit,'(A)') depofile(i)
        write(nowunit,'(A)') redfile(i)
        write(nowunit,'(A)') litfile(i)
        if(i .eq. 1 .and. flag_stat .gt. 0) write(nowunit,'(A)') mesfile(1)

        print *, ' >>>foresee message: site_nr ',i,'; input of filenames completed'

      end DO

      if(flag_multi .ne. 2) call errorfile(simfile, ios, nowunit)

 CASE (4,5)
      write(nowunit,'(A)',iostat=ios) specfile(1)
      write(nowunit,'(A)') site_name(1)
      write(nowunit,'(A)') treefile(1)
      write(nowunit,'(A)') manfile(1)
      write(nowunit,'(A)') siteall
      write(nowunit,'(A)') climall
      write(nowunit,'(A)') pathdir1
      write(nowunit,'(A)') pathdir2
      write(nowunit,'(A)') climszen

      print *, ' >>>foresee message: Input of filenames completed'

!  define name of output variables
      nvar = 1
      write(nowunit,*) outvar(nvar)
      do while (trim(outvar(nvar)) .ne. 'end')
         nvar = nvar + 1
         write(nowunit,*) outvar(nvar)
      enddo

      if (nvar .gt. 1) allocate(output_var(nvar-1,site_nr,year))

      call errorfile(simfile, ios, nowunit)

!  writeing file with desription of climate stations used
      nowunit1 = getunit()
      ios1 = 0
      open(nowunit1,file=climall,iostat=ios,status='old',action='write')
      k=1
      do
         write(nowunit1,'(A)',iostat=ios1) a
         IF (a .ne. '!') exit

      end do
      backspace nowunit1

      do

          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

      nowunit1 = getunit()

      open(nowunit1,file=siteall,iostat=ios1,status='old',action='read')
      do
         write(nowunit1,'(A)',iostat=ios1) a
         IF (a .ne. '!') exit

      end do
      backspace nowunit1

      do i=1,site_nr
         write(nowunit1,*,iostat=ios1) sitenum(i), clim_id(i), soilid(i), gwtable(i)
!          Fuellen der sitefile
         standid(i) = sitenum(i)
         site_name(i) = site_name(1)
         specfile(i) = specfile(1)
         treefile(i) = treefile(1)
         manfile(i)  = manfile(1)
         do j = 1,anzclim
           if(clim_id(i).eq.climnum(j)) then
                if(flag_climtyp.ne.0) 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
                exit
           end if
         end do
         sitefile(i) =trim(pathdir2)//'wbuek'//trim(soilid(i))//'.sop'
         valfile(i)  =trim(pathdir2)//'wbuek'//trim(soilid(i))//'.soi'
         depofile(i) ='dummy.dep'
         redfile = 'dummy.red'
         litfile = 'dummy.lit'
      enddo

       call errorfile(siteall, ios1, nowunit1)

!  variation of flag_multi= 5, especially for SILVISTRAT

 CASE (7)

      allocate(site_name_ad(site_nr))
      allocate(climfile_ad(site_nr))
      allocate(manfile_ad(site_nr))
      allocate(treefile_ad(site_nr))
	  allocate(wpmfile_ad(site_nr))
      allocate(depofile_ad(site_nr))

      allocate(fl_co2(site_nr))

      write(nowunit,'(A)',iostat=ios) specfile(1)
      write(nowunit,'(A)') site_name_all
      write(nowunit,'(A)') sitefile(1)
      write(nowunit,'(A)') valfile(1)
      write(nowunit,'(A)') siteall
      write(nowunit,'(A)') pathdir1
      write(nowunit,'(A)') pathdir2
      write(nowunit,'(A)') pathdir3
      write(nowunit,'(A)') depofile(1)
      write(nowunit,'(A)') redfile(1)
      write(nowunit,'(A)') litfile(1)

      call errorfile(simfile, ios, nowunit)

! reading control file with site-id,name, climate scenario, man-file, treeini-file, dep-file

      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),manfile_ad(i), treefile_ad(i),depofile_ad(i),fl_co2(i)
!          Fuellen der sitefile
         standid(i) = sitenum(i)
         climfile(i)= trim(pathdir1)//climfile_ad(i)
         site_name(i) = trim(site_name_all)//trim(site_name_ad(i))
         specfile(i) = specfile(1)
         sitefile(i) = sitefile(1)
         valfile(i) = valfile(1)
         treefile(i) =trim(pathdir2)//trim(treefile_ad(i))
         manfile(i)  =trim(pathdir3)//trim(manfile_ad(i))
         depofile(i)  =depofile(1)
         redfile(i)  = redfile(1)
         litfile(i)  = litfile(1)

      enddo
       flag_co2=fl_co2(1)
       call errorfile(siteall, ios1, nowunit1)

       deallocate(site_name_ad)
       deallocate(climfile_ad)
       deallocate(manfile_ad)
       deallocate(treefile_ad)
       deallocate(depofile_ad)

 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)
              tsite = sitefile(1)
              tval  = valfile(1)
              ttree = treefile(1)
              tman  = manfile(1)
              tdepo = depofile(1)
              tred  = redfile(1)
              tlit  = litfile(1)
              istand = standid(1)

              deallocate (specfile)
              deallocate (site_name)
              deallocate (sitefile)
              deallocate (valfile)
              deallocate (treefile)
              deallocate (manfile)
              deallocate (depofile)
              deallocate (redfile)
              deallocate (litfile)
              deallocate (standid)
              allocate (specfile(site_nr))
              allocate (site_name(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 (redfile(site_nr))
              allocate (litfile(site_nr))

              specfile  = tspec
              site_name = tname
              sitefile  = tsite
              valfile   = tval
              treefile  = ttree
              manfile   = tman
              depofile  = tdepo
              redfile   = tred
              litfile   = tlit
              standid   = istand

       call errorfile(simfile, ios, nowunit)

    endif   ! flag_multi = 2

END subroutine writesim

