!
! ------------------------------------------
! Sample program of parameter calibration
!    (multi model, multi objective function)
! ------------------------------------------
!
!----------------------------------------------------------------------
! Module for Common Data
!----------------------------------------------------------------------
module common_data
  implicit none
  integer, parameter, public :: ndmax = 100000

  !====== Hydrological data ======
  integer, public :: nd             ! number of data
  real,    public :: p(ndmax)       ! precipitation (mm/day)
  real,    public :: e(ndmax)       ! evapotranspiration (mm/day)
  real,    public :: q(ndmax)       ! discharge (mm/day)

  !====== Control parameters ======
  integer, public :: id_model       ! model ID
  integer, public :: id_algo        ! algorithm ID
  integer, public :: id_objf        ! objective function ID

  !====== Control parameter file ======
  character(len=80), public :: pfname = 'param.txt'

end module common_data

!----------------------------------------------------------------------
! Main Program
!----------------------------------------------------------------------
program test5
  use common_data
  use util_model
  use util_eval
  implicit none
  integer, parameter :: nmax = 100
  integer          :: n, i
  real             :: xinit(1:nmax), xmin(1:nmax), xmax(1:nmax)
  character(len=4) :: xname(1:nmax)
  real             :: x(1:nmax), r
  real             :: qc(1:ndmax)

  !====== Initialize control parameters ======
  call init_system

  !====== Read input date ======
  open(1, file='indata.txt', status='old')
  do i=1, ndmax
    read(1, *, end=99) p(i), e(i), q(i)
  end do
  99 continue
  close(1)
  nd = i - 1

  !====== Initialize model ======
  call init_model(id_model, n, xmin, xmax, xinit, xname) ! util_model

  !====== Optimize ======
  call optimize(n, xmin, xmax, xinit, x, r)

  !====== Calculate discharge with optimum parameters ======
  call run_model(id_model, x, nd, p, e, qc)              ! util_model

  !====== Print result ======
  do i=1, n
    print *, xname(i), ' =', x(i)
  end do
  print *
  print *, 'RMSE =', rmse(nd, q, qc)                     ! util_eval
  print *, 'NSE  =', nse(nd, q, qc)                      ! util_eval
  print *

  !====== Write output date ======
  open(2, file='outdata.txt')
  do i=1, nd
    write(2, *) p(i), e(i), q(i), qc(i)
  end do
  close(2)

end program test5

!----------------------------------------------------------------------
! Evaluation Function
!----------------------------------------------------------------------
function evfunc(x, ev) result(r)
  use common_data
  use util_model
  use util_eval
  implicit none
  real, intent(in)  :: x(:)      ! parameter list
  real, intent(out) :: ev        ! evaluation value
  integer           :: r         ! number of constrained condition
  real :: qc(1:nd)

  !====== Check constraint conditions ======
  r = check_model(id_model, x)                        ! util_model

  !====== Run model and evaluate ======
  if(r == 0) then
    call run_model(id_model, x, nd, p, e, qc)         ! util_model
    ev = eval_objfunc(id_objf, nd, q, qc)             ! util_eval
  else
    ev = 0.0
  end if

end function evfunc

!----------------------------------------------------------------------
! Optimization
!----------------------------------------------------------------------
subroutine optimize(n, xmin, xmax, xinit, x, r)
  use common_data
  use opti_sceua
  use opti_pso
  use util_eval
  implicit none
  integer, intent(in)  :: n
  real,    intent(in)  :: xmin(1:n), xmax(1:n), xinit(1:n)
  real,    intent(out) :: x(1:n), r
  integer  :: minmax                    ! minimize-maximize flag
  external :: evfunc

  minmax = eval_minmax(id_objf)         ! util_eval
  select case (id_algo)
    case(1)
      call do_sceuaf(n, xmin, xmax, xinit, minmax, evfunc, pfname, x, r)
    case(2)
      call do_psof(n, xmin, xmax, xinit, minmax, evfunc, pfname, x, r)
    case default
      print *, 'invalid algorithm'
      stop
  end select

end subroutine optimize

!----------------------------------------------------------------------
! Initialize control parametes
!----------------------------------------------------------------------
subroutine init_system
  use common_data
  implicit none
  integer  :: model, algorithm, objfunc
  namelist /ctrl_opti/ model, algorithm, objfunc

  !====== Read control parameter file ======
  open(3, file=pfname, status='old')
  read(3, ctrl_opti)
  close(3)

  !====== Setup control parameters ======
  id_model = model
  id_algo  = algorithm
  id_objf  = objfunc

end subroutine init_system