!!  Copyright (C) Stichting Deltares, 2005-2014.
!!
!!  This file is part of iMOD.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!  Contact: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.

subroutine pest1log()

! modules
use global, only: lipest

#ifdef IPEST
use imod_utl, only: imod_utl_printtext,imod_utl_itos
use pestvar, only: pest_iter,lgrad,llnsrch,pest_igrad,iupestout
#endif

implicit none

! locals
character(len=256) :: line

!###======================================================================

#ifndef IPEST
lipest = .false.; return
#endif

if (.not.lipest) return

#ifdef IPEST

call imod_utl_printtext('',-1,iupestout)
if(lgrad)then
   if(pest_iter.eq.0)then
      call imod_utl_printtext('Initial simulation',-1,iupestout)
   else
      line='Cycle '//TRIM(imod_utl_itos(pest_igrad))//': Derivatives by Forward Finite-Difference'
      call imod_utl_printtext(trim(line),-1,iupestout)
   ENDIF
ENDIF
if(llnsrch)then
   call imod_utl_printtext('Line Search',-1,iupestout)
endif

#endif

end subroutine pest1log

subroutine pest1alpha_grid(ptype,a,nrow,ncol,nlay,a2)

! modules
use imod_utl, only: imod_utl_printtext,imod_utl_itos,imod_utl_rtos
use global, only: lipest

#ifdef IPEST
use pestvar, only: param, pest_iter,lgrad,llnsrch,pest_igrad,iupestout
#endif

implicit none

! arguments
integer, intent(in) :: nrow, ncol, nlay
real, dimension(ncol,nrow,nlay), intent(inout) :: a
real, dimension(ncol,nrow,nlay), intent(in), optional :: a2
character(len=2), intent(in) :: ptype

! parameters
real, parameter :: tiny=1.0e-20

! locals
character(len=256) :: line
integer :: i, j, k, ils, irow, icol
real :: c, fct, ppart

!###======================================================================

#ifndef IPEST
lipest = .false.; return
#endif

if (.not.lipest) return

#ifdef IPEST

!## initialize parameters
if(pest_iter.eq.0)then

   do i=1,size(param)

      if (trim(param(i)%ptype).ne.trim(ptype)) cycle

      ils=param(i)%ils ! layer

      if(.not.associated(param(i)%x))allocate(param(i)%x(param(i)%nodes))
      select case (trim(ptype))
      case('KD','KH','KV','VA','SC','AF')
         do j=1,param(i)%nodes
            irow=param(i)%irow(j); icol=param(i)%icol(j)
            param(i)%x(j)=a(icol,irow,ils)*param(i)%f(j)
         enddo
      case('VC')   !## vertical c values
         do j=1,param(i)%nodes
            irow=param(i)%irow(j); icol=param(i)%icol(j)
            c = 1/(a(icol,irow,ils)+tiny)
            param(i)%x(j)=c*param(i)%f(j)
         enddo
      case('AA')   !## anisotropy angle
         if (.not.present(a2)) then
            call imod_utl_printtext('Cannot apply PEST scaling factor for anisotrophy angle',2)
         end if
         do j=1,param(i)%nodes
            irow=param(i)%irow(j); icol=param(i)%icol(j)
            if (a2(icol,irow,ils).lt.1.0) then
               param(i)%x(j)=a(icol,irow,ils)*param(i)%f(j)
            endif
         enddo
      case default
         call imod_utl_printtext('Something went wrong for PEST',2)
      end select
   end do
end if

do i=1,size(param)

   if (trim(param(i)%ptype).ne.trim(ptype)) cycle

   fct=param(i)%alpha(1); ils=param(i)%ils ! layer
   fct=exp(fct) ! scaling

   select case (trim(ptype))
   case('KD','KH','KV','VA','SC','AF')
      do j=1,param(i)%nodes
         irow=param(i)%irow(j); icol=param(i)%icol(j)
         ppart             =a(icol,irow,ils)-param(i)%x(j)
         a(icol,irow,ils)=ppart+param(i)%x(j)*fct
      enddo
   case('VC')   !## vertical c values
      do j=1,param(i)%nodes
         irow=param(i)%irow(j); icol=param(i)%icol(j)
         c = 1/(a(icol,irow,ils)+tiny)
         ppart = c-param(i)%x(j)
         c = ppart + param(i)%x(j)*fct
         a(icol,irow,ils)=1/(c+tiny)
      enddo
   case('AA')   !## anisotropy angle
      if (.not.present(a2)) then
         call imod_utl_printtext('Cannot apply PEST scaling factor for anisotrophy angle',2)
      end if
      do j=1,param(i)%nodes
         irow=param(i)%irow(j); icol=param(i)%icol(j)
         if (a2(icol,irow,ils).lt.1.0) then
            ppart           =a(icol,irow,ils)-param(i)%x(j)
            a(icol,irow,ils)=ppart+param(i)%x(j)*fct
         endif
      enddo
   case default
      call imod_utl_printtext('Something went wrong for PEST',2)
   end select
   line=' * '//param(i)%ptype//' adjusted ('//trim(imod_utl_itos(param(i)%nodes))//')with alpha='//trim(imod_utl_rtos(fct,'f',7))
   call imod_utl_printtext(trim(line),-1,iupestout)
end do

#endif
end subroutine

subroutine pest1alpha_list(ptype,nlist,rlist,ldim,mxlist,iopt1,iopt2)

! modules
use global, only: lipest

#ifdef IPEST
use global, only: buff
use imod_utl, only: imod_utl_printtext, imod_utl_itos, imod_utl_rtos
use m_mf2005_main, only: kper
use pestvar, only: param, pest_iter, iupestout
#endif

implicit none

! arguments
character(len=2), intent(in) :: ptype
integer, intent(in) :: nlist, ldim, mxlist
real, dimension(ldim,mxlist), intent(inout) :: rlist
integer, intent(in), optional :: iopt1, iopt2

! locals
character(len=256) :: line
character(len=1024) :: errmsg
integer :: i, j, k, ils, irow, icol, idat, irivsubsys, irivrfact,&
           idrnsubsys, ihfbfact, nadj
real :: ppart, fct

!###======================================================================

#ifndef IPEST
lipest = .false.; return
#endif

if (.not.lipest) return

#ifdef IPEST

! first, mark the cells
buff = 0.
do i=1,size(param)

   if (trim(param(i)%ptype).ne.trim(ptype)) cycle
   ils=param(i)%ils ! system/layer

   select case (trim(ptype))
   case('RC','IC') ! river/isg conductances
      errmsg = 'Cannot apply PEST scaling factor for river/isg conductance'
      if (.not.present(iopt1)) call imod_utl_printtext(trim(errmsg),2)
      irivsubsys = iopt1
      if (irivsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)

      if (trim(ptype).eq.'IC') ils=-ils
      do j = 1, nlist ! match sybsystem number
         irow=rlist(1,j); icol=rlist(2,j)
         if (int(ils).eq.int(rlist(irivsubsys,j))) buff(icol,irow,1) = real(j)
      end do
      idat = 4
   case('RI','II') ! river/isg infiltration factors
      errmsg = 'Cannot apply PEST scaling factor for river/isg infiltration factor'
      if (.not.present(iopt1)) call imod_utl_printtext(trim(errmsg),2)
      if (.not.present(iopt2)) call imod_utl_printtext(trim(errmsg),2)
      irivsubsys = iopt1; irivrfact = iopt2
      if (irivrfact.eq.0.or.irivsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)

      if (trim(ptype).eq.'II') ils=-ils
      do j = 1, nlist ! match sybsystem number
         irow=rlist(1,j); icol=rlist(2,j)
         if (int(ils).eq.int(rlist(irivsubsys,j))) buff(icol,irow,1) = real(j)
      end do
      idat = irivrfact
   case('DC') ! drain conductances
      errmsg = 'Cannot apply PEST scaling factor for drain conductance'
      if (.not.present(iopt1)) call imod_utl_printtext(trim(errmsg),2)
      idrnsubsys = iopt1
      if (idrnsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)
      do j = 1, nlist ! match sybsystem number
         irow=rlist(1,j); icol=rlist(2,j)
         if (int(ils).eq.int(rlist(idrnsubsys,j))) buff(icol,irow,1) = real(j)
      end do
      idat = 4
   end select
end do

! allocate
if (pest_iter.eq.0) then
   do i=1,size(param)
      if (trim(param(i)%ptype).ne.trim(ptype)) cycle
      ils=param(i)%ils ! system

      if(.not.associated(param(i)%x))allocate(param(i)%x(param(i)%nodes))
      select case (trim(ptype))
      case('RC','RI','IC','II','DC')
         do j=1,param(i)%nodes
            irow=param(i)%irow(j); icol=param(i)%icol(j)
            k = int(buff(icol,irow,1))
            if (k.gt.0) then
               param(i)%x(j)=rlist(idat,k)*param(i)%f(j)
            end if
         end do
      end select
   end do
end if

! adjust
do i=1,size(param)
   nadj = 0
   if (trim(param(i)%ptype).ne.trim(ptype)) cycle
   fct=param(i)%alpha(1); ils=param(i)%ils ! system

   select case (trim(ptype))
   case('RC','RI','IC','II','DC')
      do j=1,param(i)%nodes
         irow=param(i)%irow(j); icol=param(i)%icol(j)
         k = int(buff(icol,irow,1))
         if (k.gt.0) then
            nadj = nadj + 1
            ppart        =rlist(idat,k)-param(i)%x(j)
            rlist(idat,k)=ppart+param(i)%x(j)*fct
         end if
      end do
   case('HF')
      errmsg = 'Cannot apply PEST scaling factor for horizontal flow barrier'
      if (.not.present(iopt1)) call imod_utl_printtext(trim(errmsg),2)
      ihfbfact = iopt1
      if (idrnsubsys.ne.1) call imod_utl_printtext(trim(errmsg),2)
      fct = param(i)%alpha(1); fct = exp(fct)
      do j = 1, nlist
         if (int(rlist(1,j)).eq.ils) then
            nadj = nadj + 1
            rlist(6,j) = rlist(6,j)*fct
         end if
      end do
   end select

   line=' * '//param(i)%ptype//' adjusted ('//trim(imod_utl_itos(nadj))//') with alpha='//trim(imod_utl_rtos(fct,'f',7))
   if(kper.eq.1)then
      call imod_utl_printtext(trim(line),-1,iupestout)
   else
      call imod_utl_printtext(trim(line),1)
   endif

end do

#endif

end subroutine