!*****************************************************************!
!*                                                               *!
!*              4C (FORESEE) Simulation Model                    *!
!*                                                               *!
!*                                                               *!
!*                    Subroutines for:                           *!
!*                                                               *!
!*           preparation of statistical analysis                 *!
!*                                                               *!
!*                    Author: F. Suckow                          *!
!*                                                               *!
!*   contains:                                                   *!
!*   mess                                                        *!
!*   prep_mw                                                     *!
!*   prep_simout                                                 *!
!*   kind_pos                                                    *!
!*   store_sim_kind                                              *!
!*   prep_stat_out                                               *!
!*   read_simout                                                 *!
!*   open_sfile                                                  *!
!*                                                               *!
!*                  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 mess

use data_mess
use data_out
use data_simul

implicit none

integer i, j, k
integer :: hd = -99
real ::    hv = -9999.0
real ::    helpn, totm1, totm2, totm3   ! total match as average from several values
integer maxmess
logical ex
character(10)  :: helpsim
character(150) :: filename

allocate (app(site_nr))
if (unit_mess .lt. 0) then
  do
      inquire (File = mesfile(1), exist = ex)
      if(ex .eqv. .false.) then
        write (*, '(A)') ' >>>foresee message: File ',trim(mesfile(1)),' not exists !'
        write (*, '(A)', advance='no') ' please write full name of measurement file: ' 
        read(*,'(A)') mesfile(1)    
        cycle
      else
        exit
      endif
  enddo
endif

! error.log schreiben
write(unit_err,'(A)')
write(unit_err,'(A)')
write(unit_err,'(A)') ' * * * * *  Statistics  * * * * *'
write(unit_err,'(A)') 

fkind = 0
call prep_mw
if (tkind .eq. 1) call stat_mon
call prep_simout
if (.not. flag_mess) return
call prep_stat_out
do i = 1,site_nr
    ip = i
    app(i) = i
    nme_av   = 0.
    nmae_av  = 0.
    nrmse_av = 0.
    pme_av   = 0.
    prmse_av = 0.
    tic_av   = 0.
    meff_av  = 0.
    rsq_av   = 0.
    totm1    = 0.
    totm2    = 0.
    imk_nme  = imkind
    imk_nmae = imkind
    imk_nrmse= imkind
    imk_rsq  = imkind

    call read_simout
    call residuen(i)
    call statistik

! Mittelwert berechnen und ausdrucken   
    helpn = imkind - fkind
    nme_av   = nme_av / (imk_nme - fkind)
    nmae_av  = nmae_av / (imk_nmae - fkind)
    nrmse_av = nrmse_av / (imk_nrmse - fkind)
    pme_av   = pme_av / helpn
    prmse_av = prmse_av / helpn
    tic_av   = tic_av / helpn
    meff_av  = meff_av / helpn
    rsq_av   = rsq_av/(imk_rsq - fkind)
    
! Calculation of total match without missing values
    helpn = 2.
    totm1 = tic_av + (1.-meff_av)    
    totm2 = totm1
    totm3 = totm1   
    totm1 = totm1/helpn    
    if (rsq_av .ge. 0.) then
        helpn = helpn + 1.
        totm2 = totm2 + (1-rsq_av)  
        totm3 = totm2
        totm2 = totm2 / helpn
    endif
    if (nrmse_av .lt. -9000.) then
        helpn = helpn + 1.
        totm3 = (totm2 + nrmse_av) / helpn 
    endif
    write (unit_stat, '(I5,2X, A20,1X,A10,I8,1X,33E13.5)') ip, site_name(ip), 'average', hd, &
                        hv, hv, hv, hv, hv, hv, nme_av, hv, nmae_av, hv, hv, nrmse_av, pme_av, prmse_av, tic_av, meff_av, hv, rsq_av, &
                        hv, hv, hv, hv, hv, hv, hv, hv, hv, hv, hv, hv, totm1, totm2, totm3
    write (unit_stat,*)

! File mit Residuen schreiben
    if (flag_stat .ge. 2) then
        write (helpsim,'(I4)') ip
        read (helpsim,*) anh
        filename = trim(dirout)//trim(site_name(ip))//'_resid'//'.res'//trim(anh)
        unit_mout = getunit() 
        open(unit_mout,file=filename,status='replace')

        write (unit_mout, '(A)') '#    Residuals  etc.'    
        write (unit_mout, '(A)') '#       Number      kind   '    
        do j = 1, imkind
           write (unit_mout, '(I14,3X,A10,26X)', advance='no') val(j)%imes, val(j)%mkind
        enddo 
        write (unit_mout, '(A)') ' '
        do j = 1, imkind
           write (unit_mout, '(A)', advance='no') '      day year     residual   simulation  measurement'
        enddo 
        write (unit_mout, '(A)') ' '
    
        maxmess = maxval(val%imes)
        do k = 1, maxmess
           do j = 1, imkind
              if (val(j)%imes .ge. k) then
                 write (unit_mout, '(4X,2I5,3E13.5)', advance='no') val(j)%day(k), val(j)%year(k), val(j)%resid(k), val(j)%sim(k), val(j)%mess(k)
              else
                 write (unit_mout, '(4X,2I5,3E13.5)', advance='no') hd, hd, hv,hv,hv
              endif
           enddo 
           write (unit_mout, '(A)') ' '
        enddo

        close(unit_mout) 
    endif
enddo

write (*,*)
write (*, '(A)') ' Statistical analysis completed'
write (*,*)

END SUBROUTINE mess

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

SUBROUTINE prep_mw

use data_mess
use data_simul

implicit none

INTERFACE
  SUBROUTINE kind_pos(pos1, pos2, ikind, imkind, vkind, text)
  ! assumed shape arrays
      integer :: ikind, imkind
      character(150) text
      character(10), dimension(ikind):: vkind  
      integer, dimension(:):: pos1, pos2   ! Position of variables in input file
  END SUBROUTINE
END INTERFACE

integer i, j, k, ios
integer id, im, iy, itz
integer idate
character(3) ttext 
character(250) text, filename

idate = 10
allocate (mtz(2,idate))
  unit_cons = getunit()
  open(unit_cons,file='con')
if (unit_mess .lt. 0) then
  filename = mesfile(1)
  unit_mess = getunit() 
  open(unit_mess,file=filename,iostat=ios,status='old',action='read')
endif

    do
      read(unit_mess,*) text
      ios = scan(text, '!')
      IF (ios .eq. 0) then
         backspace(unit_mess)
         exit
      endif
    enddo

! determin kind of measurement values; read 1. line
  read (unit_mess, '(A)') text

  ttext = adjustl(text)
  if (ttext.eq.'dat' .or. ttext.eq.'Dat' .or. ttext.eq.'DAT') then
     tkind = 1   ! day
  else
     tkind = 2   ! year  
  endif
  call store_sim_kind(imkind, sim_kind, text)

! convert measurement values to daily counter 
  select case (tkind)
  case (1)     ! daily values
      imess = 0
      do 
        read (unit_mess, '(2(I2,1X),I4)',iostat=ios) id, im, iy
        if (ios .lt. 0) exit
        call daintz(id,im,iy,itz)
        imess = imess + 1
        if (imess .gt. idate) then
           allocate (help1(2,idate))
           help1 = mtz
           deallocate (mtz)
           idate = idate + 10
           allocate (mtz(2,idate))
           do j= 1,idate - 10
              mtz(1,j) = help1(1,j)
              mtz(2,j) = help1(2,j)
           enddo
           deallocate (help1)
        endif
        mtz(1,imess) = itz
        mtz(2,imess) = iy
      enddo

    !read meassurement values
      rewind (unit_mess)
      allocate (mess1 (imess, imkind))
      mess1 = -9999.0

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

      do j = 1,imess
        read (unit_mess, *,iostat=ios) text, (mess1(j,k), k=1,imkind)
      enddo

    case (2)    ! yearly values
      imess = 0
      if(allocated(mess1)) then
         write (*,'(A)') ' Feld mess1 bereits allokiert'
         STOP
      endif
      allocate (mess1(idate, imkind))
      mess1 = -9999.0
      do 
        imess = imess + 1
        mtz(1,imess) = 0
        read (unit_mess, *,iostat=ios) mtz(2,imess), (mess1(imess,k), k=1,imkind)
        mtz(1,imess) = 0
        if (ios .lt. 0) exit
        if (imess .gt. idate-1) then
           allocate (help1(2,idate))
           allocate (help2(idate, imkind))
           help1 = mtz
           help2 = mess1
           deallocate (mtz)
           deallocate (mess1)
           idate = idate + 10
           allocate (mtz(2,idate))
           allocate (mess1(idate, imkind))
           mess1 = -9999.9
           do j= 1,idate - 10
              mtz(1,j) = 0
              mtz(2,j) = help1(2,j)
              do k=1,imkind
                 mess1(j,k) = help2(j,k)
              enddo
           enddo
           deallocate (help1)
           deallocate (help2)
        endif
      enddo
      imess = imess - 1
  end select

END SUBROUTINE prep_mw

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

SUBROUTINE prep_simout

use data_mess
use data_out
use data_simul

implicit none

INTERFACE
  SUBROUTINE kind_pos(pos1, pos2, ikind, imkind, vkind, text)
  ! assumed shape arrays
      integer :: ikind, imkind
      character(150) text
      character(10), dimension(ikind):: vkind  
      integer, dimension(:):: pos1, pos2   ! position of variablen in input file
  END SUBROUTINE
END INTERFACE

integer i, ii, ik, j, k, year1
integer, allocatable, dimension(:):: yd, yy
character(150) :: filename

flag_mess = .FALSE.
year1 = year

! Create complete array of measurements
select case (tkind)
case (1)
    anz_val = 0
    allocate (yd(year1))
    allocate (yy(year1))
    do i=1,year1
       yy(i) = time_b + i - 1
       if (mod(yy(i),4) .eq. 0 .and. yy(i) .ne. 1900) then
          yd(i) = 366
       else
          yd(i) = 365
       endif
       anz_val = anz_val + yd(i)
    enddo

    allocate (mess2(anz_val, imkind))
    allocate (help1(2,anz_val))
    mess2 = -9999.0
    j = 1
    k = 0
       do while (mtz(2,j) .lt. time_b)
          j = j+1
       enddo
    do ii = 1, year1
       do i = 1, yd(ii)
          k = k + 1
          help1(1,k) = i
          help1(2,k) = yy(ii)
          if ((mtz(1,j) .eq. help1(1,k)) .and. (mtz(2,j) .eq. help1(2,k))) then
             do ik = 1, imkind
                mess2(k,ik) = mess1(j,ik)
                flag_mess = .TRUE.    
             enddo   ! ik
                j = j+1
          else
             do ik = 1, imkind
                mess2(k,ik) = -9999.9    
             enddo   ! ik
          endif
       enddo   ! i
    enddo   ! ii

case (2)
    allocate (yy(year1))
    anz_val = year1
    do i=1,year1
       yy(i) = time_b + i - 1
    enddo

    allocate (mess2(anz_val, imkind))
    allocate (help1(2,anz_val))
    mess2 = -9999.9
    j = 1
    do while (mtz(2,j) .lt. time_b)
      j = j+1
    enddo
    do ii = 1, year1
          help1(2,ii) = yy(ii)
          help1(1,ii) = 0
          if (mtz(2,j) .eq. help1(2,ii)) then
             do ik = 1, imkind
                mess2(ii,ik) = mess1(j,ik)  
                flag_mess = .TRUE.    
             enddo   ! ik
                j = j+1
          else
             do ik = 1, imkind
                mess2(ii,ik) = -9999.9    
             enddo   ! ik
          endif
    enddo   ! ii

end select

if (.not. flag_mess) then
    write (*,*)
    write (*, '(A)') ' Statistical analysis:'
    write (*, '(A)') ' No measurements within the simulation period'
    write (*,*)
    return
endif

! write file with complete set of meassurement values
    if (flag_stat .eq. 3) then
        filename = trim(dirout)//trim(site_name(1))//'_mess'//'.mes'
        unit_mout = getunit() 
        open(unit_mout,file=filename,status='replace')

        write (unit_mout, '(A)') '# Measurements '    
        write (unit_mout, '(A)') mess_info   
        write (unit_mout, '(A)', advance='no') '# day year'    
        do i=1,imkind
            write (unit_mout, '(A13)', advance='no') sim_kind(i)
        enddo
        write (unit_mout, '(A)') ' '
    
        do i = 1, anz_val
          write (unit_mout, '(2I5)', advance='no') help1(1,i), help1(2,i)
          do j = 1, imkind
             write (unit_mout, '(E13.5)', advance='no') mess2(i,j)
          enddo 
           write (unit_mout, '(A)') ' '
        enddo

        close(unit_mout) 
    endif

! Read data
allocate (sim1(anz_val, imkind))
allocate (stz(2,anz_val))

END SUBROUTINE prep_simout

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

SUBROUTINE kind_pos(pos1, pos2, ikind, imkind, vkind, text)

implicit none

integer imkind, &  ! amount of read kinds of measurment value 
        ikind,  &  ! amount of allowed kinds of measurement value
        j    
character(10), dimension(ikind):: vkind  
character(150) text
integer, dimension(:):: pos1, pos2   ! position of variable in input file
  
  pos1 = 9999
  imkind = 0
  do j = 1,ikind
    pos1(j) = index (text, trim(vkind(j)))
    pos2(j) = j
    if (pos1(j) .eq. 0) then
       pos1(j) = 9999
    else
       imkind = imkind +1
    endif
  enddo  ! j
 call sort_index(ikind, pos1, pos2)

END SUBROUTINE kind_pos

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

SUBROUTINE store_sim_kind(imkind, vkind, text)

implicit none

integer imkind, &  ! amount of read kinds of measurement values
        ipos, &    ! position of space character/sign
        i, j    
character(10), dimension(30):: vkind  
character(250) text, text1, text2
character(1):: setleer = ''
character(75):: setascii
  
  setascii = ''
  do i = 48,122    
    j = i-47
    setascii(j:j) = ACHAR(i)   ! fill in with ASCII-character, no space character/signs
  enddo
  imkind = 0
  ipos = verify(adjustl(text), setascii)  ! first non-ASCII-character
  text1 = ' '
  text2 = adjustl(text)
  text1 = text2(ipos:250)         ! delete date/year
  text2 = text1
  ipos = scan(text2, setascii)   ! first ASCII-character
  text1 = text2(ipos:250)        ! delete non-ASCII-characters
  text2 = text1
  do j = 1,30
      ipos = verify(text2, setascii) ! first non_ASCII-character
      vkind(j) = text2(1:ipos-1)     ! save name of measurement value
      imkind = imkind +1
      text1 = text2(ipos:250)          ! delete saved measurment value
      text2 = text1
      ipos = scan(text2, setascii)    ! first ASCII-character
      if (ipos .eq. 0) exit
      text1 = text2(ipos:250) 
      text2 = text1
  enddo  ! j

END SUBROUTINE store_sim_kind

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

SUBROUTINE prep_stat_out

use data_mess
use data_out
use data_simul

implicit none

character(70) :: filename 
character(8) actdate
character(10) acttime

  filename = trim(site_name(1))//'_stat'//'.res'

  call date_and_time(actdate, acttime)
  unit_stat = getunit()
  open(unit_stat,file=trim(dirout)//filename,status='replace')

write (unit_stat, '(A)') '#  Comparison of simulated and observed values' 
write (unit_stat, '(10A)') '#  Date: ',actdate(7:8),'.',actdate(5:6),'.',actdate(1:4), &
                            '  Time: ',acttime(1:2),':',acttime(3:4)
write (unit_stat, 1000)
write (unit_stat, 2000)

1000 format('#                                            |--------   residuals   .......        ',       15('             '), &                                                                                                                                                                                     
            '|-----------------------------        simulation      -----------------------||-------------------------------      observed     ---------------------------|' )
2000 format( '# ipnr site_id              kind      number         mean          min          max      stand_dev     variance    var_coeff          NME          MAE         NMAE', &
            '          SSE         RMSE        NRMSE          PME        PRMSE          TIC         MEFF    cor_coeff      rsquare', &
            '         mean          min          max    stand_dev     variance    var_coeff         mean          min          max    stand_dev     variance    var_coeff   tot_match1   tot_match2   tot_match3')
                          
END SUBROUTINE prep_stat_out

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

SUBROUTINE read_simout

use data_mess
use data_out
use data_simul
use data_soil

implicit none

integer i,j, ios
character(150) :: text
character(50)  :: message
character(10)  :: helpsim
character(10)  :: styp, skind
character      :: text1
character(2)   :: text2
character(3)   :: text3
logical ex
integer        :: year1, unithelp
real, dimension(26):: help_day
real, dimension(13):: help_sum   ! size is adjusted to amount of elements in ...sum.out 
real, dimension(27):: help_veg
real, dimension(28):: help_veg_spec
real, dimension(8):: help_lit
real, dimension(33):: help_soil
real, dimension(50):: tief
real, allocatable, dimension(:) :: help_temp, help_water
real htief, hnlay

sim1        = -9999.9
unitday     = -99
unitcbal     = -99
unitlit     = -99
unittemp    = -99
unitsum     = -99
unitveg     = -99
unitveg_pi  = -99
unitveg_sp  = -99
unitveg_bi  = -99
unitsoil    = -99
unitsoilini = -99
unitwater   = -99
anz_sim     = ip

year1 = year

do i=1,imkind
    select case (sim_kind(i))
    case ('AET')
      if (tkind .eq. 1) then    ! daily values
          skind = 'day'
          styp  = 'out'
          if (unitday .lt. 0) call open_sfile (skind, styp, unitday)
          opos2(i) = 7
      else 
          skind = 'soil'
          styp  = 'out'
          if (unitsoil .lt. 0) call open_sfile (skind, styp, unitsoil)
          opos2(i) = 10
      endif

    case ('BIOM', 'STVOL')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 14

    case ('STVOL_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 14

    case ('STVOL_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 14

    case ('STVOL_bi')
      skind = 'veg_bi'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 14

    case ('DG')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 7

    case ('DG_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 7

    case ('DG_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 7

    case ('DG_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 7

    case ('DBH')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 23

    case ('DBH_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 24

    case ('DBH_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 24

    case ('DBH_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 24

    case ('Fol')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 9

    case ('Fol_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 9

    case ('Fol_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 9

    case ('Fol_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 9

    case ('GPP')
      if (tkind .eq. 1) then    ! daily values
          skind = 'sum'
          styp  = 'out'
          if (unitsum .lt. 0) call open_sfile (skind, styp, unitsum)
          opos2(i) = 11
      else 
          skind = 'c_bal'
          styp  = 'out'
          if (unitcbal .lt. 0) call open_sfile (skind, styp, unitsum)
          opos2(i) = 1
      endif

    case ('HO')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 8

    case ('HO_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 8

    case ('HO_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 8

    case ('HO_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 8

    case ('LAI')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 4

    case ('LAI_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 4

    case ('LAI_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 4

    case ('LAI_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 4

    case ('MH')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 24

    case ('MH_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 25

    case ('MH_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 25

    case ('MH_bi')
      skind = 'veg_bi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 25

    case ('NTREE')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 3

    case ('NTREE_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 3

    case ('NTREE_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 3

    case ('NTREE_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 3

    case ('NEE')
      skind = 'sum'
      styp  = 'out'
      if (unitsum .lt. 0) call open_sfile (skind, styp, unitsum)
      opos2(i) = 6

    case ('NEP')
      skind = 'c_bal'
      styp  = 'out'
      if (unitcbal .lt. 0) call open_sfile (skind, styp, unitcbal)
      opos2(i) = 3

    case ('Litter')
      skind = 'litter'
      styp  = 'out'
      if (unitlit .lt. 0) call open_sfile (skind, styp, unitlit)
      opos2(i) = 1

    case ('prec_stand')
      skind = 'soil'
      styp  = 'out'
      if (unitsoil .lt. 0) call open_sfile (skind, styp, unitsoil)
      opos2(i) = 2

    case ('prec_st_d')
      skind = 'day'
      styp  = 'out'
      if (unitday .lt. 0) call open_sfile (skind, styp, unitday)
      opos2(i) = 4

    case ('s_resp')
      skind = 'day'
      styp  = 'out'
      if (unitday .lt. 0) call open_sfile (skind, styp, unitday)
      opos2(i) = 12

    case ('Snow')
      skind = 'day'
      styp  = 'out'
      if (unitday .lt. 0) call open_sfile (skind, styp, unitday)
      opos2(i) = 5

    case ('STBIOM')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 10

    case ('STBIOM_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 10

    case ('STBIOM_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 10

    case ('STBIOM_bi')
      skind = 'veg_bi'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 10

    case ('Stem_inc')
      skind = 'veg'
      styp  = 'out'
      if (unitveg .lt. 0) call open_sfile (skind, styp, unitveg)
      opos2(i) = 13

    case ('Stem_inc_pi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_pi .lt. 0) call open_sfile (skind, styp, unitveg_pi)
      opos2(i) = 13

    case ('Stem_inc_sp')
      skind = 'veg_sp'
      styp  = 'out'
      if (unitveg_sp .lt. 0) call open_sfile (skind, styp, unitveg_sp)
      opos2(i) = 13

    case ('Stem_inc_bi')
      skind = 'veg_pi'
      styp  = 'out'
      if (unitveg_bi .lt. 0) call open_sfile (skind, styp, unitveg_bi)
      opos2(i) = 13

    case ('TER')
      if (tkind .eq. 1) then    ! daily values
          skind = 'sum'
          styp  = 'out'
          if (unitsum .lt. 0) call open_sfile (skind, styp, unitsum)
          opos2(i) = 12
      else 
          skind = 'c_bal'
          styp  = 'out'
          if (unitcbal .lt. 0) call open_sfile (skind, styp, unitsum)
          opos2(i) = 6
      endif

    case ('transtree')
      skind = 'day'
      styp  = 'out'
      if (unitday .lt. 0) call open_sfile (skind, styp, unitday)
      opos2(i) = 9

    case ('WC_002')
      skind = 'watvol'
      styp  = 'out'
      if (unitwater .lt. 0) call open_sfile (skind, styp, unitwater)
      opos2(i) = 1

    case ('TS_002')
      skind = 'temp'
      styp  = 'out'
      if (unittemp .lt. 0) call open_sfile (skind, styp, unittemp)
      opos2(i) = 2
    
    case default

      text2 = sim_kind(i) (1:2)     
      if ((text2 .eq. 'TS') .or. (text2 .eq. 'WC')) then       
          skind = 'soil'
          styp  = 'ini'
          if (unitsoilini .lt. 0) then
            call open_sfile (skind, styp, unitsoilini)
            read (unitsoilini, *) text
            read (unitsoilini, *) text
            do j=1, 50
                read (unitsoilini, *,iostat=ios) hnlay, tief(j)
                if (hnlay .eq. 0) then 
                    exit
                else
                    nlay = hnlay
                endif
                if (ios .ne. 0) exit
            enddo
          endif

          select case  (text2)
          case ('TS')
              skind = 'temp'
              styp  = 'out'
              if (unittemp .lt. 0) call open_sfile (skind, styp, unittemp)
      
              text3 = sim_kind(i) (4:6)
              write (helpsim, *) text3
              read (helpsim,*) htief
        !      htief = 5.
              do j=2,nlay
                  if ((tief(j)-tief(1)) .ge. htief) then 
                      opos2(i) = j+1
                      exit
                  endif
              enddo
              if (opos2(i) .le.0) then
                 message = "no simulation values of "//text2//" for depth "
                  opos2(i) = nlay
                  write(unit_err,'(A)',advance='no') trim(message)
                  write(unit_err,'(F5.0,A)') htief, " cm"
              else
                  message = "simulation values of "//text2//" for depth "
                  write(unit_err,'(A)',advance='no') trim(message)
                  write(unit_err,'(F5.0,A)') htief, " cm"
                  message = "                        selected layer: "
                  write(unit_err,'(A)',advance='no') trim(message)
                  write(unit_err,'(I3)') j
              endif

          case ('WC')
              skind = 'watvol'
              styp  = 'out'
              if (unitwater .lt. 0) call open_sfile (skind, styp, unitwater)
      
              text3 = sim_kind(i) (4:6)
              write (helpsim, *) text3
              read (helpsim,*) htief
              do j=2,nlay
                  if ((tief(j)-tief(1)) .ge. htief) then 
                      opos2(i) = j
                      exit
                  endif
              enddo
              if (opos2(i) .le.0) then
                 message = "no simulation values of "//text2//" for depth "
                  opos2(i) = nlay
                  write(unit_err,'(A)',advance='no') trim(message)
                  write(unit_err,'(F5.0,A)') htief, " cm"
              else
                  message = "simulation values of "//text2//" for depth "
                  write(unit_err,'(A)',advance='no') trim(message)
                  write(unit_err,'(F5.0,A)') htief, " cm"
                  message = "                        selected layer: "
                  write(unit_err,'(A)',advance='no') trim(message)
                  write(unit_err,'(I3)') j
              endif

          end select  ! text2
      else
          fkind = fkind + 1
          write (unit_err, *) 
          write (unit_err, '(A)') 'Statistics - Undefined kind of measurement  '//sim_kind(i)
      endif

    end select
enddo   ! i - imkind

! read in results file

! read day-file 
  if (unitday .ge. 0) then
    do
      read(unitday,*) text
      IF (adjustl(text) .ne. '#') then
         backspace(unitday)
         exit
      endif
    enddo

    do j = 1,anz_val
       read (unitday, *) stz(1,j), stz(2,j), help_day
       do i=1,imkind
         select case (sim_kind(i))
         case ('AET','Snow','prec_st_d','s_resp','transtree')
            sim1(j,i) = help_day(opos2(i))   
         end select
       enddo
    enddo
  endif   ! unitday

! read temp-file 
  if (unittemp .ge. 0) then
    do
      read(unittemp,*) text
      IF (adjustl(text) .ne. '#') then
         backspace(unittemp)
         exit
      endif
    enddo
    allocate (help_temp(nlay))

    do j = 1,anz_val
       read (unittemp, *) stz(1,j), stz(2,j), help_temp
       do i=1,imkind
         if (opos2(i) .gt. 0) then
             select case (sim_kind(i) (1:2))
             case ('TS')
                sim1(j,i) = help_temp(opos2(i))    
             end select
         endif
       enddo
    enddo
    deallocate (help_temp)
  endif   ! unittemp

! read water-file 
  if (unitwater .ge. 0) then
    do
      read(unitwater,*) text
      IF (adjustl(text) .ne. '#') then
         backspace(unitwater)
         exit
      endif
    enddo
    allocate (help_water(nlay))

    do j = 1,anz_val
       read (unitwater, *) stz(1,j), stz(2,j), help_water
       do i=1,imkind
         if (opos2(i) .gt. 0) then
         select case (sim_kind(i) (1:2))
             case ('WC')
                sim1(j,i) = help_water(opos2(i))    
             end select
         endif
       enddo
    enddo
    deallocate (help_water)
  endif   ! unitwater

! read sum-file 
  if (unitsum .ge. 0) then
    do
      read(unitsum,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         backspace(unitsum)
         exit
      endif
    enddo

    do j = 1,anz_val
       read (unitsum, *) stz(1,j), stz(2,j), help_sum
       do i=1,imkind
         select case (sim_kind(i))
         case ('NEE','GPP','TER')
            sim1(j,i) = help_sum(opos2(i))    
         end select
       enddo
    enddo
  endif   ! unitsum

! read c_bal-file 
  if (unitcbal .ge. 0) then
    do
      read(unitcbal,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit      ! 1. line for standard values is skiped
      endif
    enddo
    do j = 1,year1
       read (unitcbal, *) stz(2,j), help_veg
       do i=1,imkind
         select case (sim_kind(i))
         case ('NEP','GPP','TER')
            sim1(j,i) = help_veg(opos2(i))   

         end select
       enddo
    enddo
  endif   ! unitcbal

! read litter-file 
  if (unitlit .ge. 0) then
    do
      read(unitlit,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit
      endif
    enddo

    do j = 1,year1
       read (unitlit, *) stz(2,j), help_lit
       do i=1,imkind
         select case (sim_kind(i))

         case ('Litter')
            sim1(j,i) = help_lit(opos2(i))    

         end select
       enddo
    enddo
  endif   ! unitlit

! read soil-file 
  if (unitsoil .ge. 0) then
    do
      read(unitsoil,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit      ! 1. line of standard values is skiped
      endif
    enddo
    do j = 1,year1
       read (unitsoil, *) stz(2,j), help_soil
       do i=1,imkind
         select case (sim_kind(i))
         case ('prec_stand')
            sim1(j,i) = help_soil(opos2(i)) - help_soil(opos2(i)+1)    

         case ('AET')
            sim1(j,i) = help_soil(opos2(i))    
         end select
       enddo
    enddo
  endif   ! unitsoil

! read veg-file 
  if (unitveg .ge. 0) then
    do
      read(unitveg,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit
      endif
    enddo
    do j = 1,year1
       read (unitveg, *) stz(2,j), help_veg
       do i=1,imkind
         select case (sim_kind(i))
         case ('STBIOM')
            sim1(j,i) = (help_veg(opos2(i)) + help_veg(opos2(i)+2))    

         case ('BIOM','DG','DBH','Fol','LAI','NTREE','Stem_inc')
            sim1(j,i) = help_veg(opos2(i))    

         case ('HO','MH')
            sim1(j,i) = help_veg(opos2(i)) / 100.    

         end select
       enddo
    enddo
  endif   ! unitveg

! read veg_pi-file 
  if (unitveg_pi .ge. 0) then
    do
      read(unitveg_pi,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit
      endif
    enddo
    do j = 1,year1
       read (unitveg_pi, *) stz(2,j), help_veg_spec
       do i=1,imkind
         select case (sim_kind(i))
         case ('STBIOM_pi')
            sim1(j,i) = (help_veg_spec(opos2(i)) + help_veg_spec(opos2(i)+2))    

         case ('BIOM_pi','DG_pi','DBH_pi','Fol_pi','LAI_pi','NTREE_pi','Stem_inc_pi')
            sim1(j,i) = help_veg_spec(opos2(i))    

         case ('HO_pi','MH_pi')
            sim1(j,i) = help_veg_spec(opos2(i)) / 100.    

         end select
       enddo
    enddo
  endif   ! unitveg_pi

! read veg_sp-file 
  if (unitveg_sp .ge. 0) then
    do
      read(unitveg_sp,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit
      endif
    enddo
    do j = 1,year1
       read (unitveg_sp, *) stz(2,j), help_veg_spec
       do i=1,imkind
         select case (sim_kind(i))
         case ('STBIOM_sp')
            sim1(j,i) = (help_veg_spec(opos2(i)) + help_veg_spec(opos2(i)+2))    

         case ('BIOM_sp','DG_sp','DBH_sp','Fol_sp','LAI_sp','NTREE_sp','Stem_inc_sp')
            sim1(j,i) = help_veg_spec(opos2(i))    

         case ('HO_sp','MH_sp')
            sim1(j,i) = help_veg_spec(opos2(i)) / 100.    

         end select
       enddo
    enddo
  endif   ! unitveg_sp

! read veg_bi-file 
  if (unitveg_bi .ge. 0) then
    do
      read(unitveg_bi,*) text
      text1 = adjustl(text)
      IF (text1 .ne. '#') then
         exit
      endif
    enddo
    do j = 1,year1
       read (unitveg_bi, *) stz(2,j), help_veg_spec
       do i=1,imkind
         select case (sim_kind(i))
         case ('STBIOM_bi')
            sim1(j,i) = (help_veg_spec(opos2(i)) + help_veg_spec(opos2(i)+2))    

         case ('BIOM_bi','DG_bi','DBH_bi','Fol_bi','LAI_bi','NTREE_bi','Stem_inc_bi')
            sim1(j,i) = help_veg_spec(opos2(i))    

         case ('HO_bi','MH_bi')
            sim1(j,i) = help_veg_spec(opos2(i)) / 100.    

         end select
       enddo
    enddo
  endif   ! unitveg_bi
        
END SUBROUTINE read_simout


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

SUBROUTINE open_sfile (okind, otyp, unitnr)

use data_mess
use data_out
use data_simul

implicit none

integer unitnr
character(150) :: simsumfile  ! simulation output sum-file
character(150) :: simoutfile  ! simulation output file
character(10)  :: helpsim
character(10)  :: otyp, okind
logical ex

      WRITE(helpsim,'(I2)') app(ip)
      read(helpsim,*) anh
      simoutfile = trim(dirout)//trim(site_name(ip))//'_'//trim(okind)//'.'//trim(otyp)//trim(anh)
      inquire (File = simoutfile, exist = ex)
      if(ex .eqv. .false.) then
        write (*, '(A)') ' >>>foresee message: no such file  ', adjustl(simoutfile)    
        return
      else
         write (*, '(A)')  ' >>>foresee message: Filetest - file exists ',trim(simoutfile)
      endif
      unitnr   = getunit()
      open(unitnr,file=simoutfile,status='old')

END SUBROUTINE open_sfile


