!
! --------------------------------------------------------------
! Module of 2-series tank model
! --------------------------------------------------------------
!
module model_tank2
  implicit none
  public :: init_tank2         ! initialize model
  public :: run_tank2          ! run model
  public :: check_tank2        ! check constraint conditions
contains

!----------------------------------------------------------------------
! Initialize 2-seriese Tank Model
!----------------------------------------------------------------------
subroutine init_tank2(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:8), xma(1:8), xin(1:8)
  character(len=4) :: xnm(1:8)
  data np / 8 /
  data xnm(1),xmi(1),xma(1),xin(1)  / 'a11 ',  0.,    1.,   0.1 /
  data xnm(2),xmi(2),xma(2),xin(2)  / 'a12 ',  0.,    1.,  0.01 /
  data xnm(3),xmi(3),xma(3),xin(3)  / 'a2  ',  0.,    1.,  0.01 /
  data xnm(4),xmi(4),xma(4),xin(4)  / 'b   ',  0.,    1.,  0.01 /
  data xnm(5),xmi(5),xma(5),xin(5)  / 'z11 ',  0.,  400.,   50. /
  data xnm(6),xmi(6),xma(6),xin(6)  / 'z12 ',  0.,  200.,   10. /
  data xnm(7),xmi(7),xma(7),xin(7)  / 'h1  ',  0.,  500.,    0. /
  data xnm(8),xmi(8),xma(8),xin(8)  / 'h2  ',  0., 5000.,  300. /
  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_tank2

!----------------------------------------------------------------------
! Run 2-seriese Tank Model
!----------------------------------------------------------------------
subroutine run_tank2(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    :: a11, a12, a2, b, z11, z12, h1, h2
  real    :: dt, q11, q12, q2, r
  integer :: nt, i, j

  a11 = x(1); a12 = x(2); a2  = x(3); b   = x(4)
  z11 = x(5); z12 = x(6); h1  = x(7); h2  = x(8)

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

  do i=1, nd
    q(i) = 0.0
    do j=1, nt
       q11 = a11 * max(h1 - z11, 0.0)                ! surface flow
       q12 = a12 * max(h1 - z12, 0.0)                ! sub-surface flow
       r   = b * h1                                  ! infiltration
       q2  = a2 * h2                                 ! base flow
       h1 = h1 + (p(i) - e(i) - q11 - q12 - r) * dt
       h2 = h2 + (r - q2) * dt
       h2 = h2 + min(h1, 0.0)        ! evaporate from h2 if h1 is empty
       h1 = max(h1, 0.0)
       h2 = max(h2, 0.0)
       q(i) = q(i) + q11 + q12 + q2
    end do
  end do

end subroutine run_tank2

!----------------------------------------------------------------------
! Check constraint conditions of 2-seriese Tank Model
!----------------------------------------------------------------------
function check_tank2(x) result(r)
  implicit none
  real, intent(in) :: x(1:8)
  integer          :: r              ! number of constrained condition

  !====== Check constraint conditions ======
  r = 0
  if(x(1) < x(2))              r = r + 1        ! a11 < a12
  if(x(5) < x(6))              r = r + 1        ! z11 < z12
  if(x(1) + x(2) + x(4) > 1.0) r = r + 1        ! a11+a12+b > 1

end function check_tank2

end module model_tank2