-
Petra Lasch-Born authoredPetra Lasch-Born authored
out_var_stat.f 17.59 KiB
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* output of variables with statistics for climate scenarios *!
!* *!
!* contains *!
!* OUT_VAR_STAT compressing of output variables *!
!* CALC_STAT calculation of statistics *!
!* 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 out_var_stat(kind, act_real)
! compressing of output variables with statistics (multi run 9, 10)
use data_out
use data_par
use data_simul
use data_site
IMPLICIT NONE
integer kind ! 1 - aggregation per realisation (average)
! 2 - aggregation per climate scenario over all realisations with statistics
! 3 - statistics per month over all years
integer act_real ! number of actual realisation
integer i, j, k, unit_nr, ii
real varerr, help
character(50) :: filename ! complete name of output file
real, dimension(nrreal) :: helparr
real, dimension(year):: helpmon
character(30) :: helpvar
character(20) idtext, datei
character(150) htext
! mit Numerical Recipies
REAL:: adev,ave,var, &
curt=-99. , &
sdev=-99. , &
skew=0.
! Statistische Masszahlen fuer Klimaszen.-Realisierungen
real:: avcl, & ! Mittelwert
mincl, & ! Minimum
maxcl, & ! Maximum
median, & ! Median
stdevcl=-99. , & ! Standardabweichung
varicl, & ! Streuung
varcocl ! Variationskoeffizient
real quant05, quant95 ! 0.05 and 0.95 quantile
real, external :: mean, variance
if (flag_trace) write (unit_trace, '(I4,I10,A,2I5)') iday, time_cur, ' out_var_stat ',kind,act_real
select case (kind)
case (1,2)
if (output_unit_all .le.0) then
filename = trim(site_name1)//'_var_all.out'
output_unit_all = getunit()
open(output_unit_all,file=trim(dirout)//filename,status='replace') write (output_unit_all, '(A)') '# Output of mean annual values for each site and each realization of climate scenarios'
write (output_unit_all, '(A, I6)') '# Simulation period (years): ', year
write (output_unit_all, '(A, I6)') '# Number of climate scenarios: ', nrclim
write (output_unit_all, '(A, I6)') '# Number of realizations: ', nrreal
write (output_unit_all, *)
write (output_unit_all, '(A)', advance='no') '# Type_clim.scen. Site_ip Real.'
do i = 1, nvar-1
select case (trim(outvar(i)))
case ('AET_year','cwb_year','GPP_year','NEP_year','NPP_year','perc_year','PET_year','temp_year','TER_year','prec_year','resps_year')
continue
case ('AET_mon','cwb_mon','GPP_mon','NEP_mon','NPP_mon','perc_mon','PET_mon','temp_mon','TER_mon','prec_mon','resps_mon')
continue
case ('AET_week','cwb_week','GPP_week','NEP_week','NPP_week','perc_week','PET_week','temp_week','TER_week','prec_week','resps_week')
continue
case default
write (output_unit_all, '(A12)', advance='no') trim(outvar(i))
end select
enddo
write (output_unit_all, '(A)') ''
endif
case (3)
do i = 1, nvar-1
if (output_unit_mon(i) .le.0) then ! for monthly values
filename = trim(site_name1)//'_'//trim(outvar(i))//'_stat.res'
output_unit_mon(i) = getunit()
open(output_unit_mon(i),file=trim(dirout)//filename,status='replace')
write (output_unit_mon(i), '(A)') '# Output of mean monthly values for '//trim(outvar(i))
write (output_unit_mon(i), '(A, I6)') '# Simulation period (years): ', year
varerr = 0
endif
enddo
end select
select case (kind)
case (1) ! after each realisation
write (output_unit_all, '(2X, A15, 1X, A10, I5,2X)', advance = 'no') trim(typeclim(iclim)), sitenum(ip), act_real
do i = 1, nvar-1
select case (trim(outvar(i)))
case ('AET_year','cwb_year','GPP_year','NEP_year','NPP_year','perc_year','PET_year','temp_year','TER_year','prec_year','resps_year')
ii = output_var(i,1,0)
do j = 1, year
climszenyear(ii,ip,iclim,act_real,j) = output_var(i,1,j)
enddo
case ('AET_mon','cwb_mon','GPP_mon','NEP_mon','NPP_mon','perc_mon','PET_mon','temp_mon','TER_mon','prec_mon','resps_mon')
ii = output_var(i,1,0)
do k = 1,12
help = 0.
do j = 1, year
help = help + output_varm(ii,1,j,k)
enddo
help = help / year
climszenmon(ii,ip,iclim,act_real,k) = help
enddo
case ('AET_week','cwb_week','GPP_week','NEP_week','NPP_week','perc_week','PET_week','temp_week','TER_week','prec_week','resps_week')
ii = output_var(i,1,0)
do k = 1,52
help = 0.
do j = 1, year
help = help + output_varw(ii,1,j,k) enddo
help = help / year
climszenweek(ii,ip,iclim,act_real,k) = help
enddo
case default
help = 0.
do j = 1, year
help = help + output_var(i,1,j)
enddo ! j
help = help / year
climszenres(i,ip,iclim,act_real) = help
write (output_unit_all, '(E12.4)', advance = 'no') help
end select ! outvar
end do ! i
write (output_unit_all, '(A)') ''
case (2) ! am Ende der Simulation
do i = 1, nvar-1
if (output_unit(i) .lt. 0) then
helpvar = outvar(i)
call out_var_select(helpvar, varerr, unit_nr)
if (varerr .ne. 0.) then
output_unit(i) = unit_nr
write (unit_nr, '(A, I6)') '# Simulation period (years): ', year
write (unit_nr, '(A, I6)') '# Number of climate scenarios: ', nrclim
write (unit_nr, '(A, I6)') '# Number of realizations: ', nrreal
select case (trim(outvar(i)))
case ('AET_year','cwb_year','GPP_year','NEP_year','NPP_year','perc_year','PET_year','temp_year','TER_year','prec_year','resps_year')
write (unit_nr, '(A)') '# Statistics over all realizations for each year '
write (unit_nr, '(A)') '# Type_clim.scen. Site_ip Year Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
case ('AET_mon','cwb_mon','GPP_mon','NEP_mon','NPP_mon','perc_mon','PET_mon','temp_mon','TER_mon','prec_mon','resps_mon')
write (unit_nr, '(A)') '# Statistics over all realizations and all years for each month '
write (unit_nr, '(A)') '# Type_clim.scen. Site_ip Month Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
case ('AET_week','cwb_week','GPP_week','NEP_week','NPP_week','perc_week','PET_week','temp_week','TER_week','prec_week','resps_week')
write (unit_nr, '(A)') '# Statistics over all realizations and all years for each week '
write (unit_nr, '(A)') '# Type_clim.scen. Site_ip Week Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
case default
write (unit_nr, '(A)') '# Statistics over all realizations (mean of all years) '
write (unit_nr, '(A)') '# Type_clim.scen. Site_ip Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
end select
else
write (*,*)
write (*,*) '*** 4C-error - output of variables (out_var_file): ', trim(outvar(i)), ' not found'
write (*,*)
write (unit_err,*)
write (unit_err,*) '*** 4C-error - no such output variable (out_var_file): ', trim(outvar(i))
endif
endif
if (output_unit(i) .ge. 0) then
select case (trim(outvar(i)))
case ('AET_year','cwb_year','GPP_year','NEP_year','NPP_year','perc_year','PET_year','temp_year','TER_year','prec_year','resps_year')
ii = output_var(i,1,0)
do k = 1, year
write (output_unit(i), '(2X, A15, 1X, A10, I7)', advance = 'no') trim(typeclim(iclim)), sitenum(ip), k
do j = 1, nrreal
helparr(j) = climszenyear(ii,ip,iclim,j,k)
enddo
call calc_stat(nrreal, helparr, output_unit(i))
enddo
case ('AET_mon','cwb_mon','GPP_mon','NEP_mon','NPP_mon','perc_mon','PET_mon','temp_mon','TER_mon','prec_mon','resps_mon')
ii = output_var(i,1,0)
do k = 1, 12 write (output_unit(i), '(2X, A15, 1X, A10, I7)', advance = 'no') trim(typeclim(iclim)), sitenum(ip), k
do j = 1, nrreal
helparr(j) = climszenmon(ii,ip,iclim,j,k)
enddo
call calc_stat(nrreal, helparr, output_unit(i))
enddo
case ('AET_week','cwb_week','GPP_week','NEP_week','NPP_week','perc_week','PET_week','temp_week','TER_week','prec_week','resps_week')
ii = output_var(i,1,0)
do k = 1, 52
write (output_unit(i), '(2X, A15, 1X, A10, I7)', advance = 'no') trim(typeclim(iclim)), sitenum(ip), k
do j = 1, nrreal
helparr(j) = climszenweek(ii,ip,iclim,j,k)
enddo
call calc_stat(nrreal, helparr, output_unit(i))
enddo
case default
write (output_unit(i), '(2X, A15, 1X, A10)', advance = 'no') trim(typeclim(iclim)), sitenum(ip)
do j = 1, nrreal
helparr(j) = climszenres(i,ip,iclim,j)
enddo
call calc_stat(nrreal, helparr, output_unit(i))
end select
endif
enddo
case (3) ! Monthly values
do i = 1, nvar-1
helpvar = outvar(i)
select case (trim(outvar(i)))
case ('AET_year','cwb_year','GPP_year','NEP_year','NPP_year','perc_year','PET_year','temp_year','TER_year','prec_year','resps_year')
ii = output_var(i,1,0)
do j = 1, year
climszenyear(ii,ip,iclim,act_real,j) = output_var(i,1,j)
enddo
case ('GPP_mon','NPP_mon','TER_mon')
ii = output_var(i,1,0)
if (ip .eq.1) then
write (output_unit_mon(i), '(A)') '# Statistics over all years for each month '
write (output_unit_mon(i), '(A)') '# g C/m '
write (output_unit_mon(i), '(A)') '# ipnr site_id Month Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
endif
do k = 1,12
help = 0.
do j = 1, year
helpmon(j) = output_varm(ii,1,j,k) * 100. ! tC/ha --> gC/m
enddo
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! only write last 20 signs
write (output_unit_mon(i), '(I5,2X, A20,I5)', advance = 'no') ip, idtext, k
call calc_stat(year, helpmon, output_unit_mon(i))
enddo
case ('NEE_mon')
ii = output_var(i,1,0)
if (ip .eq.1) then
write (output_unit_mon(i), '(A)') '# Statistics over all years for each month '
write (output_unit_mon(i), '(A)') '# g C/m '
write (output_unit_mon(i), '(A)') '# ipnr site_id Month Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
endif
do k = 1,12
help = 0.
do j = 1, year
helpmon(j) = output_varm(ii,1,j,k) ! gC/m
enddo
htext = adjustr(site_name(ip)) idtext = adjustl(htext (131:150)) ! only write last 20 signs
write (output_unit_mon(i), '(I5,2X, A20,I5)', advance = 'no') ip, idtext, k
call calc_stat(year, helpmon, output_unit_mon(i))
enddo
case ('resps_mon')
ii = output_var(i,1,0)
if (ip .eq.1) then
write (output_unit_mon(i), '(A)') '# Statistics over all years for each month '
write (output_unit_mon(i), '(A)') '# g C/m '
write (output_unit_mon(i), '(A)') '# ipnr site_id Month Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
endif
do k = 1,12
help = 0.
do j = 1, year
helpmon(j) = output_varm(ii,1,j,k) * kgha_in_gm2 ! kgC/ha --> gC/m
enddo
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! only write last 20 signs
write (output_unit_mon(i), '(I5,2X, A20,I5)', advance = 'no') ip, idtext, k
call calc_stat(year, helpmon, output_unit_mon(i))
enddo
case ('AET_mon','cwb_mon','perc_mon','PET_mon','temp_mon','prec_mon')
ii = output_var(i,1,0)
if (ip .eq.1) then
write (output_unit_mon(i), '(A)') '# Statistics over all years for each month '
write (output_unit_mon(i), '(A)') '# '
write (output_unit_mon(i), '(A)') '# ipnr site_id Month Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
endif
do k = 1,12
help = 0.
do j = 1, year
helpmon(j) = output_varm(ii,1,j,k)
enddo
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! only write last 20 signs
write (output_unit_mon(i), '(I5,2X, A20,I5)', advance = 'no') ip, idtext, k
call calc_stat(year, helpmon, output_unit_mon(i))
enddo
case ('AET_week','cwb_week','GPP_week','NEP_week','NPP_week','perc_week','PET_week','temp_week','TER_week','prec_week','resps_week')
ii = output_var(i,1,0)
do k = 1,52
help = 0.
do j = 1, year
help = help + output_varw(ii,1,j,k)
enddo
help = help / year
climszenweek(ii,ip,iclim,act_real,k) = help
enddo
case default
help = 0.
do j = 1, year
help = help + output_var(i,1,j)
enddo ! j
help = help / year
climszenres(i,ip,iclim,act_real) = help
write (output_unit_all, '(E12.4)', advance = 'no') help
end select ! outvar
end do ! i
write (output_unit_all, '(A)') ''
end select
END SUBROUTINE out_var_stat
!**************************************************************
SUBROUTINE calc_stat(nreal, helparr, outunit)
! calculate statistics
use data_out
use data_simul
IMPLICIT NONE
integer :: outunit ! output unit
integer :: nreal ! number of elements
real, dimension(nreal) :: helparr ! input-array with dimension nreal
! with numerical recipies
REAL:: adev,ave,var, &
curt=-99. , &
sdev=-99. , &
skew=0.
! statistical measurment figures for climate scenario realisation
real:: avcl, & ! mean
mincl, & ! minimum
maxcl, & ! maximum
median, & ! median
stdevcl=-99. , & ! standard deviation
varicl, & ! dispersion
varcocl ! coefficient of variance
real quant05, quant95 ! 0.05 and 0.95 quantile
real, external :: mean, variance
avcl = mean(nreal, helparr)
mincl = minval(helparr)
maxcl = maxval(helparr)
varicl = variance(nreal, avcl, helparr)
if (varicl .ge. 0.) stdevcl = sqrt(varicl)
if (avcl .ne. 0.) then
varcocl = stdevcl / avcl
else
varcocl = -9999.0
endif
call quantile(nreal, helparr, quant05, quant95, median)
! with numerical recipies
if (nreal .gt. 1) call moment(helparr, nreal, ave,adev,sdev,var,skew,curt)
write (outunit, '(11E12.4)') avcl, mincl, maxcl, varicl, varcocl, sdev, skew, curt, quant05, quant95, median
END SUBROUTINE calc_stat