Forked from
4C / FORESEE
182 commits behind the upstream repository.
-
Petra Lasch-Born authoredPetra Lasch-Born authored
statistik.f 15.15 KiB
!*****************************************************************!
!* *!
!* 4C (FORESEE) Simulation Model *!
!* *!
!* *!
!* Subroutines for: *!
!* *!
!* statistical analysis of model quality *!
!* *!
!* Author: F. Suckow *!
!* *!
!* contains: *!
!* residuen *!
!* statistik *!
!* mean (n, arr) *!
!* variance (n, meanv, arr) *!
!* correl (n, meanv1, arr1, meanv2, arr2) *!
!* sumsq (n, arr) *!
!* stat_mon *!
!* *!
!* 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 residuen (ip)
use data_mess
implicit none
integer i,j, ires, ip
! calculate and save residues, with date, simulation as well as measurement value
! Residuen berechnen, mit Datum, Sim.- und Messwert speichern
if (ip .eq. 1) then
allocate (val(imkind))
do i=1,imkind
ires = 0
val(i)%tkind = tkind
allocate (val(i)%day(1:anz_val))
allocate (val(i)%year(1:anz_val))
allocate (val(i)%resid(1:anz_val))
allocate (val(i)%sim(1:anz_val))
allocate (val(i)%mess(1:anz_val))
val(i)%day = -99
val(i)%year = -99
val(i)%resid = -9999.0
val(i)%mkind = sim_kind(i)
do j = 1,anz_val
if (mess2(j,i) .gt. -9000.0 .and. sim1(j,i) .gt. -9000.0) then
ires = ires + 1
val(i)%day(ires) = stz(1,j)
val(i)%year(ires) = stz(2,j)
val(i)%resid(ires)= sim1(j,i) - mess2(j,i)
val(i)%sim(ires) = sim1(j,i)
val(i)%mess(ires) = mess2(j,i)
else
endif
enddo
val(i)%imes = ires
enddo
else do i=1,imkind
ires = 0
val(i)%resid = -9999.9
val(i)%sim = -9999.9
do j = 1,anz_val
if (mess2(j,i) .gt. -9000.0 .and. sim1(j,i) .gt. -9000.0) then
ires = ires + 1
val(i)%resid(ires)= sim1(j,i) - mess2(j,i)
val(i)%sim(ires) = sim1(j,i)
else
endif
enddo
enddo
endif
END SUBROUTINE residuen
!**************************************************************
SUBROUTINE statistik
use data_mess
use data_simul
implicit none
integer imt ! aktueller Messwert-Typ
real, external :: mean, variance, correl, sumsq
integer i, n, nhelp
real help, h1, h2
real, allocatable, dimension(:):: arr, arrs, arrm, harr
real:: avs, & ! mean value simulation; Mittelwert Simulation
mins, & ! Minimum Simulation
maxs, & ! Maximum Simulation
stdevs, & ! standard deviation simulation; Standardabweichung Simulation
varis, & ! scattering of simulation; Streuung Simulation
varcos, & ! coefficient of variation for simulation; Variationskoeffizient Simulation
avm, & ! mean value measurements; Mittelwert Messwerte
minm, & ! minimum value measurements; Minimum Messwerte
maxm, & ! maximum value measurements; Maximum Messwerte
stdevm, & ! standard deviation measurements; Standardabweichung Messwerte
varim, & ! scattering of measurements; Streuung Messwerte
varcom, & ! coefficient of variation of measurements; Variationskoeffizient Messwerte
corrco, & ! coefficient of correlation; Korrelationskoeffizient
rsq, & ! coefficient of determination; Bestimmtheitsmass
avr, & ! mean error residues; Mittlerer Fehler Residuen
minr, & ! minimum residues; Minimum Residuen
maxr, & ! maximum residues; Maximum Residuen
stdevr, & ! standard deviation residues; Standardabweichung Residuen
varir, & ! scattering of residues; Streuung Residuen
varcor, & ! coefficient of variation residues; Variationskoeffizient Residuen
nme, & ! normalised mean error; Normalisierter mittlerer Fehler
mae, & ! mean absolute error of residues; Mittlerer absoluter Fehler Residuen
nmae, & ! normalised mean absolute error; Normalisierter mittlerer absoluter Fehler
sse , & ! sum of squared errors; Fehlerquadratsumme
rmse, & ! Root mean square error
nrmse, & ! Normalised root mean square error
pme, & ! mean procental error; Mittlerer prozentualer Fehler
prmse, & ! mean squared procentual error; Mittlerer quadratischer prozentualer Fehler
tic, & ! Theilsch imbalance coefficient; Theilscher Ungleichheitskoeffizient
meff ! Model efficiency (Medlyn et al. 2005)
do imt = 1, imkind
n = val(imt)%imes
if (n .gt. 0) then
allocate (arr(n))
allocate (arrs(n))
allocate (arrm(n)) allocate (harr(n))
! simulation
arrs = val(imt)%sim(1:n)
avs = mean(n, arrs)
mins = minval(arrs)
maxs = maxval(arrs)
varis = variance(n, avs, arrs)
if (varis .ge. 0.) then
stdevs = sqrt(varis)
else
stdevs = 0.
endif
if (avs .ne. 0.) then
varcos = stdevs / avs
else
varcos = -9999.0
endif
! observed
arrm = val(imt)%mess(1:n)
avm = mean(n, arrm)
minm = minval(arrm)
maxm = maxval(arrm)
varim = variance(n, avm, arrm)
if (varim .ge. 0.) then
stdevm = sqrt(varim)
else
stdevm = 0.
endif
! residuals
arr = val(imt)%resid(1:n)
avr = mean(n, arr)
minr = minval(arr)
maxr = maxval(arr)
varir = variance(n, avr, arr)
if (varir .ge. 0.) then
stdevr = sqrt(varir)
else
stdevr = 0.
endif
if (avr .ne. 0.) then
varcor = stdevr / avr
else
varcor = -9999.0
endif
corrco = correl(n, avs, arrs, avm, arrm)
if (corrco .ge. -1.) then
rsq = corrco * corrco
rsq_av = rsq_av + rsq
else
imk_rsq = imk_rsq - 1
rsq = -9999.0
endif
if (avs .ne. 0.) then
nme = (avm - avs) / avs
nme_av = nme_av + nme
else
imk_nme = imk_nme - 1
nme = -9999.0
endif
mae = mean(n, abs(arr))
sse = sumsq(n, arr)
rmse = sqrt(sse / n) if (avm .ne. 0.) then
varcom = stdevm / avm
nrmse = rmse / abs(avm)
nrmse_av = nrmse_av + nrmse
nmae = mae / abs(avm)
nmae_av = nmae_av + nmae
else
imk_nrmse = imk_nrmse - 1
imk_nmae = imk_nmae - 1
varcom = -9999.0
nrmse = -9999.0
nmae = -9999.0
endif
nhelp = n
do i = 1, n
if (arrm(i) .ne. 0.) then
harr(i) = abs(arr(i)/arrm(i))
else
nhelp = nhelp -1
harr(i) = 0
endif
enddo
pme = mean(nhelp, harr)
prmse = sumsq(nhelp, harr)
prmse = sqrt(prmse / nhelp)
tic = sse / (sumsq(n, arrs) + sumsq(n, arrm))
h1=sumsq(n, arr)
harr = arrm-avm
h2=sumsq(n, harr)
if (h2.ne.0) then
meff = 1. - (h1 / h2)
else
meff=1
end if
! Mittelwert
pme_av = pme_av + pme
prmse_av = prmse_av + prmse
tic_av = tic_av + tic
meff_av = meff_av + meff
deallocate (arr)
deallocate (arrm)
deallocate (arrs)
deallocate (harr)
write (unit_stat, '(I5,2X, A20,1X,A10,I8,1X,30E13.5)') ip, site_name(ip), val(imt)%mkind, val(imt)%imes, &
avr, minr, maxr, stdevr, varir, varcor, nme, mae, nmae, sse, rmse, nrmse, pme, prmse, tic,meff, corrco, rsq, &
avs, mins, maxs, stdevs, varis, varcos, avm, minm, maxm, stdevm, varim, varcom
endif
enddo
END SUBROUTINE statistik
!**************************************************************
REAL FUNCTION mean (n, arr)
integer n, i
real, dimension(n):: arr
real help
help = 0.
do i = 1,n
help = help + arr(i)
enddo
mean = help / n
END FUNCTION mean
!**************************************************************
REAL FUNCTION variance (n, meanv, arr)
integer n, i
real, dimension(n):: arr
real meanv, help, xx
help = 0.
if (n .gt. 1) then
do i = 1,n
xx = arr(i) - meanv
help = help + xx * xx
enddo
variance = help / (n -1)
else
variance = -9999.0
endif
END FUNCTION variance
!**************************************************************
REAL FUNCTION correl (n, meanv1, arr1, meanv2, arr2)
integer n, i
real, dimension(n):: arr1, arr2
real meanv1, meanv2, help1, help2, help3, xx1, xx2
help1 = 0.
help2 = 0.
help3 = 0.
do i = 1,n
xx1 = arr1(i) - meanv1
xx2 = arr2(i) - meanv2
help1 = help1 + xx1 * xx2
help2 = help2 + xx1 * xx1
help3 = help3 + xx2 * xx2
enddo
if ((help2 .gt. 1.E-06) .and. (help3 .gt. 1.E-06)) then
correl = help1 / sqrt(help2*help3)
else
correl = -9999.0
endif
END FUNCTION correl
!**************************************************************
REAL FUNCTION sumsq (n, arr)
integer n, i
real, dimension(n):: arr
real help
help = 0.
do i = 1,n
help = help + arr(i) * arr(i)
enddo
sumsq = help
END FUNCTION sumsq
!**************************************************************
Subroutine stat_mon
! Statistics of monthly values, derived from daily observed values
use data_mess
use data_out
use data_simul
implicit none
integer i, j, k, l
integer dd, mm, yy, doy, yanz, arranz
integer :: outunit ! output unit
character(250) text, filename
character(20) idtext, datei, vunit, obskind
character(150) htext
real, allocatable, dimension(:):: helparr ! help array with montly values of one month for all years
real, allocatable, dimension(:,:):: help_mon ! array with monthly values for each year, year
real, allocatable, dimension(:,:):: help_day ! array with mean daily values per month for each year, year
integer, allocatable, dimension(:,:):: help_num ! array with number of measurement values for each year, year
yanz = mtz(2,imess) - mtz(2,1) + 1
if (.not. allocated(unit_mon)) then
allocate(unit_mon(imkind))
allocate(unit_mon_stat(imkind))
allocate(helparr(yanz))
endif
if (.not. allocated(help_mon)) then
allocate(help_mon(12,yanz))
allocate(help_day(12,yanz))
allocate(help_num(12,yanz))
endif
do k = 1, imkind
help_mon = 0.0
help_num = 0
obskind = sim_kind(k)
filename = trim(dirout)//trim(site_name(ip))//'_'//trim(obskind)//'_mon_obs'//'.out'
unit_mon(k) = getunit()
open(unit_mon(k),file=filename,status='replace')
! Calculate mmonthly sums
do j = 1, imess
doy = mtz(1,j)
yy = mtz(2,j)
call TZINDA(dd,mm,yy,doy)
yy = mtz(2,j) - mtz(2,1) + 1
if (mess1(j,k) .Gt. -9990.) then
if (sim_kind(k) .eq. 'AET') then
if (mess1(j,k) .lt. 0.) then
mess_info = '# negative AET set to zero'
mess1(j,k) = 0. ! avoid negative AET
endif
endif
help_mon(mm,yy) = help_mon(mm,yy) + mess1(j,k)
help_num(mm,yy) = help_num(mm,yy) + 1
endif
enddo ! j
do j = 1, yanz
do i = 1,12
if (help_num(i,j) .gt. 0) then
help_day(i,j) = help_mon(i,j) / help_num(i,j)
else
help_mon(i,j) = -9999.
help_day(i,j) = -9999.
endif
enddo
enddo
! Output of monthly sums
select case (trim(obskind))
case ('AET')
vunit = 'mm' case ('GPP', 'NPP', 'TER')
vunit = 'g C/m'
case ('Snow')
return
case default
vunit = ' '
end select
write (unit_mon(k), '(A)') '# Monthly sum, daily mean of month and number of values per month of observed '//trim(obskind)
write (unit_mon(k), '(A)') '# '//trim(vunit)
write (unit_mon(k), '(A)', advance='no') '# Year'
do i = 1,12
write (unit_mon(k), '(A8,I2)', advance='no') trim(obskind)//'_',i
enddo
write (unit_mon(k), '(A)')
l = 0
do j = mtz(2,1), mtz(2,imess)
l = l + 1
write (unit_mon(k), '(A,I6,12F10.2)') 'sum ', j, (help_mon(i,l), i=1,12)
write (unit_mon(k), '(A,I6,12F10.2)') 'daily mean ', j, (help_day(i,l), i=1,12)
write (unit_mon(k), '(A,I6,12I10)') 'number ', j, (help_num(i,l), i=1,12)
enddo
! Statistics
filename = trim(dirout)//trim(site_name(ip))//'_'//trim(obskind)//'_mon_obs_stat'//'.res'
outunit = getunit()
open(outunit,file=filename,status='replace')
write (outunit, '(A)') '# Statistics over all years for each monthly sum and daily mean per month of '//trim(obskind)
write (outunit, '(A)') '# '//trim(vunit)
write (outunit, '(A, I6)') '# Simulation period (years): ', year
write (outunit, '(A)') '# site_id Month number Mean Minimum Maximum Variance Var.Coeff. Std.Dev. Skewness Excess 0.05-Quant. 0.95-Quant. Median'
write (outunit, '(A)') 'monthly sum'
do i = 1,12
arranz = 0
do j = 1,yanz
if (help_mon(i,j) .gt. -9990.) then
arranz = arranz + 1
helparr(arranz) = help_mon(i,j)
endif
enddo ! j
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! nur letzte 20 Zeichen schreiben
write (outunit, '(A20,I5,I8)', advance = 'no') idtext, i, arranz
call calc_stat(arranz, helparr, outunit)
enddo ! i
write (outunit, '(A)') ' '
write (outunit, '(A)') 'daily mean per month'
do i = 1,12
arranz = 0
do j = 1,yanz
if (help_day(i,j) .gt. -9990.) then
arranz = arranz + 1
helparr(arranz) = help_day(i,j)
endif
enddo ! j
htext = adjustr(site_name(ip))
idtext = adjustl(htext (131:150)) ! nur letzte 20 Zeichen schreiben
write (outunit, '(A20,I5,I8)', advance = 'no') idtext, i, arranz
call calc_stat(arranz, helparr, outunit)
enddo ! i
enddo ! k
continue
End Subroutine stat_mon