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