!
! --------------------------------------------------------------
! Module of Storage Function Model
! --------------------------------------------------------------
!
module model_sfunc
  implicit none
  public :: init_sfunc         ! initialize model
  public :: run_sfunc          ! run model
  public :: check_sfunc        ! check constraint conditions
contains

!----------------------------------------------------------------------
! Initialize storage function model
!----------------------------------------------------------------------
subroutine init_sfunc(n, xmin, xmax, xinit, xname)
  implicit none
  integer,          intent(out) :: n
  real,             intent(out) :: xinit(1:*), xmin(1:*), xmax(1:*)
  character(len=*), intent(out) :: xname(1:*)
  integer :: np
  real    :: xmi(1:5), xma(1:5), xin(1:5)
  character(len=4) :: xnm(1:5)
  data np / 5 /
  data xnm(1),xmi(1),xma(1),xin(1)  / 'k   ',  5.0,  100.0,  50.0 /
  data xnm(2),xmi(2),xma(2),xin(2)  / 'p   ',  0.1,    1.0,   0.5 /
  data xnm(3),xmi(3),xma(3),xin(3)  / 'f   ',  0.2,    1.0,   0.7 /
  data xnm(4),xmi(4),xma(4),xin(4)  / 'bf  ',  0.0,    5.0,   0.0 /
  data xnm(5),xmi(5),xma(5),xin(5)  / 's0  ',  0.0, 1000.0,   0.0 /
  n = np
  xname(1:n) = xnm(1:n)
  xmin(1:n)  = xmi(1:n)
  xmax(1:n)  = xma(1:n)
  xinit(1:n) = xin(1:n)
end subroutine init_sfunc

!----------------------------------------------------------------------
! Run storage function model
!----------------------------------------------------------------------
subroutine run_sfunc(x, nd, p, e, q)
  implicit none
  real,    intent(in)  :: x(1:8)         ! parameters
  integer, intent(in)  :: nd             ! number of data
  real,    intent(in)  :: p(1:nd)        ! rainfall (mm/day)
  real,    intent(in)  :: e(1:nd)        ! evapotranspiration (mm/day)
  real,    intent(out) :: q(1:nd)        ! simulated discharge (mm/day)
  real    :: k, pp, f, bf, s0
  real    :: dt, s, qq
  integer :: nt, i, j

  k = x(1); pp = x(2); f = x(3); bf = x(4); s0 = x(5)
  s = s0

  nt = 3
  dt = 1. / real(nt)

  do i=1, nd
    q(i) = 0.0
    do j=1, nt
       s = s + f * p(i) * dt
       qq = (s / k)**(1. / pp)
       s = s - (e(i) + qq) * dt
       s = max(s, 0.0)
       q(i) = q(i) + qq
    end do
    q(i) = q(i) + bf
  end do

end subroutine run_sfunc

!----------------------------------------------------------------------
! Check constraint conditions of storage function model
!----------------------------------------------------------------------
function check_sfunc(x) result(r)
  implicit none
  real, intent(in) :: x(1:8)
  integer          :: r              ! number of constrained condition
  r = 0
end function check_sfunc

end module model_sfunc