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