output_nashsutcliffe_efficiency Subroutine

  subroutine output_nashsutcliffe_efficiency(qo, qs, icd, istyr, iy)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES CRITERIA OF FIT:
    !              difference in calulated water balance, relative difference**2,
    !              Nash&Sutcliffe Efficiency & LOG-Nash&Sutcliffe Efficiency
    !**** CALLED IN: MAIN
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> PARAMETERS & VARIABLES in TITLE
    !      qo(1:inn) = observed water discharge, m3/sec.
    !      qs(1:inn) = simulated water discharge, m3/sec.
    !      inn = number of days
    !      icd = code: 1 - one year, 2 - whole period

    !      >>>>> COMMON PARAMETERS
    !      istyr = starting year
    !      iy = current year
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      akk(20) = accumulated values, internal
    !      crdif = difference in calulated water balance
    !      crdifp = difference in calulated water balance, %
    !      dif = difference between qs and qo
    !      eff = Nash & Sutcliffe Efficiency for variables qo and qs
    !      efflog = Nash & Sutcliffe Efficiency for log(qo) and log(qs)
    !      f00 = local par
    !      f11 = local par
    !      ik = local par
    !      im = local par
    !      qmid = local par
    !      qmid2 = local par
    !      reldif2 = relative difference **2
    !      xdif = xqo - xqs
    !      xqo = log(qo)
    !      xqs = log(qs)
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !**** Include common parameters & descriptions

    integer, intent(in) :: istyr, icd, iy
    real(dp), dimension(:), intent(in) :: qo, qs

    integer ik, im
    real(dp) efflog, f00, f11, qmid, qmid2, reldif2, xdif, xqo, xqs, akk(20), crdif, crdifp, dif, eff

    eff = 0.
    efflog = 0.

    do im = 1, 20
      akk(im) = 0.
    end do

    do ik = 1, size(qo)
      dif = qs(ik) - qo (ik)
      qmid = 0.5 * (abs(qo(ik)) + qs(ik))
      qmid2 = qmid * qmid
      akk(1) = akk(1) + 1.
      akk(2) = akk(2) + qo(ik)
      akk(3) = akk(3) + qo(ik) * qo(ik)
      akk(4) = akk(4) + dif
      akk(5) = akk(5) + dif * dif
      if (qmid2 .gt. 0.001) then
        akk(6) = akk(6) + dif * dif / qmid2
      else
      endif

      if (qo(ik) .gt. 0. .and. qs(ik) .gt. 0.) then
        xqo = log(qo(ik))
        xqs = log(qs(ik))
        xdif = xqo - xqs
        akk(7) = akk(7) + 1.
        akk(8) = akk(8) + xqo
        akk(9) = akk(9) + xqo * xqo
        akk(10) = akk(10) + xdif * xdif
      endif
    end do

    crdif = akk(4)

    if (akk(2) .ne. 0) then
      crdifp = akk(4) * 100 / akk(2)
    else
      crdifp = 0.
    endif

    reldif2 = akk(6)

    if (akk(1) .gt. 0.) then
      f00 = akk(3) - akk(2) * akk(2) / akk(1)
      if (f00 .gt. 0.0001) eff = (f00 - akk(5)) / f00
    endif

    if (akk(7) .gt. 0.) then
      f11 = akk(9) - akk(8) * akk(8) / akk(7)
      if (f11 .gt. 0.0001) efflog = (f11 - akk(10)) / f11
    endif

    if (icd .eq. 1) then
      call log_debug("output_nashsutcliffe_efficiency", &
        'CRITERIA of fit in (year): Abs. difference, Difference %, Rel.dif**2, L-Efficiency, Efficiency:', &
        i1=istyr+iy-1, reals=(/crdif, crdifp, reldif2, efflog, eff/))
    else
      call log_debug("output_nashsutcliffe_efficiency", &
        'CRITERIA of fit TOTAL: Abs. difference, Difference %, Rel.dif**2, L-Efficiency, Efficiency:', &
        reals=(/crdif, crdifp, reldif2, efflog, eff/))
    endif

    return
  end subroutine output_nashsutcliffe_efficiency