!!  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
use imod_utl, only: imod_utl_printtext,imod_utl_itos,imod_utl_dtos
use pestvar, only: pest_iter,lgrad,llnsrch,pest_igrad,iupestout,pest_ktype

implicit none

! locals
character(len=256) :: line

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

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

end subroutine pest1log

!###====================================================================
subroutine pest1_deallocate()
!###====================================================================
use pestvar, only: param, pest_single
use global, only: lipest 
implicit none
integer :: i

if(.not.lipest)return
if(pest_single.eq.0)return

write(*,*) 'release memory'
do i=1,size(param)
 if(associated(param(i)%irow))deallocate(param(i)%irow)
 if(associated(param(i)%icol))deallocate(param(i)%icol)
 if(associated(param(i)%f))   deallocate(param(i)%f)
 if(associated(param(i)%x))   deallocate(param(i)%x)
 if(associated(param(i)%xy))  deallocate(param(i)%xy)
enddo

end subroutine pest1_deallocate

!###====================================================================
subroutine pest1alpha_grid(ptype,a,nrow,ncol,nlay,iout,a2)
!###====================================================================
use imod_utl, only: imod_utl_printtext,imod_utl_itos,imod_utl_dtos,imod_utl_createdir
use gwfmetmodule, only: cdelr, cdelc, resultdir
use m_mf2005_main, only: kper
use global, only: lipest, ibound ,npertxt
use pestvar, only: param, pest_iter,lgrad,llnsrch,pest_igrad,iupestout,pest_ktype,blnkout

implicit none

! arguments
integer, intent(in) :: nrow, ncol, nlay, iout
real, dimension(ncol,nrow,nlay), intent(inout) :: a
!real(kind=8),dimension(:,:,:),allocatable :: a_dbl
real, dimension(ncol,nrow,nlay), intent(in), optional :: a2 !ncol,nrow,nlay
character(len=2), intent(in) :: ptype
real(kind=8),dimension(:,:),allocatable :: xyz,xpp
integer :: nxyz,ipp,ilay
character(len=1024) :: fname,dir
logical :: PESTALPHA_PPOINT
! parameters
real(kind=8), parameter :: tiny=1.0d-20
integer,dimension(4) :: ios
integer(kind=8) :: sdate_par,edate_par,sdate_mod,edate_mod

real,allocatable,dimension(:,:,:) :: cnt

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

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

if(.not.lipest)return

if(size(a,1).ne.ncol.or.size(a,2).ne.nrow.or.size(a,3).ne.nlay)then
 write(*,'(6i10)') size(a,1),ncol,size(a,2),nrow,size(a,3),nlay
 write(*,*) 'array a dimensions incorrect'; pause; stop
endif

call sgwf2met1pnt(1) !igrid)

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

   !## check whether too many parameters effect a location
   allocate(cnt(ncol,nrow,nlay)); cnt=0.0
   
   do i=1,size(param)

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

      ils=param(i)%ils !## layer
      
      IF(PARAM(I)%ZTYPE.EQ.0)THEN
       
       if(.not.associated(param(i)%x))allocate(param(i)%x(param(i)%nodes))
       select case (trim(ptype))
       case('KD','KH','KV','VA','SC','AF','EP','SY')
          do j=1,param(i)%nodes
             irow=param(i)%irow(j); icol=param(i)%icol(j)
!             if(irow.eq.1.and.icol.eq.499)then
!              write(*,*) param(i)%ptype,param(i)%ils; pause
!             endif
             !## only modify active/constant head nodes
             if(ibound(icol,irow,ils).ne.0)then
              param(i)%x(j)=a(icol,irow,ils)*param(i)%f(j)
              cnt(icol,irow,ils)=cnt(icol,irow,ils)+param(i)%f(j)
             endif
          enddo
       !## can become unconfined and then it is not layer 1
       case('RE','ER','ED')
          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)
             cnt(icol,irow,ils)=cnt(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)
!             if(irow.eq.1.and.icol.eq.499)then
!              write(*,*) i,param(i)%ptype,param(i)%ils,param(i)%izone; pause
!             endif
             if(ibound(icol,irow,ils).ne.0)then
              c = 1.0/(a(icol,irow,ils)+tiny)
              param(i)%x(j)=c*param(i)%f(j)
              cnt(icol,irow,ils)=cnt(icol,irow,ils)+param(i)%f(j)
             endif
          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(ibound(icol,irow,ils).ne.0)then
              if (a2(icol,irow,ils).lt.1.0) then
                param(i)%x(j)=a(icol,irow,ils)*param(i)%f(j)
                cnt(icol,irow,ils)=cnt(icol,irow,ils)+param(i)%f(j)
              endif
             endif
          enddo
       case default
          call imod_utl_printtext('Something went wrong for PEST, missing parameter '//trim(ptype),2)
       end select
      !## pilot points, save initial values for f -- hoeft toch niet???
      ELSEIF(PARAM(I)%ZTYPE.EQ.1)THEN

      end if
   end do
    
   !## check whether too much for a single location
   i=0
   do ils=1,nlay
    do irow=1,nrow
     do icol=1,ncol
      !## accept small round-off errors
      if(cnt(icol,irow,ils).gt.1.01)then
       if(i.eq.0)THEN
        write(*,'(/a/)') 'Multiple definitions (fraction > 1.01) on these locations for '//trim(ptype)
        write(*,'(5a10)') 'i','icol','irow','ilay','count'
       endif
       i=i+1; if(i.le.10)write(*,'(4i10,f15.7)') i,icol,irow,ils,cnt(icol,irow,ils)
      endif
     enddo
    enddo
   enddo
   if(i.gt.10)write(*,'(A,I10)') '11 --> ',i
   if(i.gt.0)then
    write(*,'(/A/)') 'Check your zones as similar zone numbers appear on identical locations'
    pause; stop 99
   endif
   deallocate(cnt)
end if

DIR=trim(resultdir)//'\pest\parameters_cycle'//trim(imod_utl_itos(pest_iter))
CALL IMOD_UTL_CREATEDIR(DIR)

do i=1,size(param)

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

   !## skip pilot-points
   IF(PARAM(I)%ZTYPE.EQ.1)CYCLE

   !## check date if needed
   if(param(i)%sdate.ne.''.and.param(i)%edate.ne.'')then
    if(trim(param(i)%sdate).eq.'STEADY-STATE'.and.trim(param(i)%edate).eq.'STEADY-STATE')then
     if(npertxt(kper,1).ne.'STEADY-STATE')cycle
    else
     read(npertxt(kper,2),*,iostat=ios(1)) sdate_mod
     read(npertxt(kper,3),*,iostat=ios(2)) edate_mod
     read(param(i)%sdate ,*,iostat=ios(3)) sdate_par
     read(param(i)%edate ,*,iostat=ios(4)) edate_par
     !## ready to evaluate
     if(sum(ios).eq.0)then
      !## make sure all are in same dimensions
      if(sdate_mod.lt.99999999)sdate_mod=sdate_mod*1e6
      if(edate_mod.lt.99999999)edate_mod=edate_mod*1e6
      if(sdate_par.lt.99999999)sdate_par=sdate_par*1e6
      if(edate_par.lt.99999999)edate_par=edate_par*1e6
      if(edate_par.le.sdate_mod)cycle
      if(sdate_par.gt.edate_mod)cycle
     else
      cycle
     endif
    endif
   endif

   fct=param(i)%alpha(1); ils=param(i)%ils !## layer
   IF(PARAM(I)%ILOG.EQ.1)FCT=EXP(FCT)
   IF(PARAM(I)%ILOG.EQ.2)FCT=10.0**FCT

   select case (trim(ptype))
   case('KD','KH','KV','VA','SC','AF','EP','SY')
      do j=1,param(i)%nodes
         irow=param(i)%irow(j); icol=param(i)%icol(j)
         if(ibound(icol,irow,ils).ne.0)then
           ppart           =a(icol,irow,ils)-param(i)%x(j)
           a(icol,irow,ils)=ppart+param(i)%x(j)*fct
         endif
      enddo
   !## can become unconfined and then it is not layer 1
   case('RE','ER','ED')
      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)
         if(ibound(icol,irow,ils).ne.0)then
          c = 1.0/(a(icol,irow,ils)+tiny)
          ppart = c-param(i)%x(j)
          c = ppart + param(i)%x(j)*fct
          a(icol,irow,ils)=1.0/(c+tiny)
         endif
      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(ibound(icol,irow,ils).ne.0)then
          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
         endif
      enddo

   case default
      call imod_utl_printtext('Something went wrong for PEST, missing parameter '//trim(ptype),2)

   end select
   
   if(param(i)%acronym.eq.'')then
    line='Param. '//trim(imod_utl_itos(i))//'['//param(i)%ptype//',n='//trim(imod_utl_itos(param(i)%nodes))// &
       ',f='//trim(imod_utl_dtos(fct,'f',7))//',ils='//trim(imod_utl_itos(ils))//']'
   else
    line='Param. '//trim(imod_utl_itos(i))//'['//trim(param(i)%acronym)//','//param(i)%ptype//',n='//trim(imod_utl_itos(param(i)%nodes))// &
       ',f='//trim(imod_utl_dtos(fct,'f',7))//',ils='//trim(imod_utl_itos(ils))//']'
   endif
   write(*,'(a)') trim(line)
   if(iupestout.ne.0)write(iupestout,'(a)') trim(line)
   
 end do

 !## process any pilotpoints per modellayer/ adjustable parameter
 ALLOCATE(XPP(NCOL,NROW))
 !## make sure all are positive again
 do i=1,size(param); param(i)%ztype=abs(param(i)%ztype); enddo
 DO ILS=1,NLAY

  IF(.NOT.PESTALPHA_PPOINT(ILS,PTYPE,XPP,NCOL,NROW,DIR,IOUT,1))CYCLE
  
  DO IROW=1,NROW; DO ICOL=1,NCOL
   IF(IBOUND(ICOL,IROW,ILS).EQ.0)CYCLE
   A(ICOL,IROW,ILS)=A(ICOL,IROW,ILS)*XPP(ICOL,IROW)
  ENDDO; ENDDO
 
 ENDDO
 !## make sure all are positive again
 do i=1,size(param); param(i)%ztype=abs(param(i)%ztype); enddo

 DEALLOCATE(XPP) 

 !## export only initially
 if(PEST_IGRAD.eq.0)then
  allocate(xpp(ncol,nrow))
  do i=1,size(param)
   if (trim(param(i)%ptype).ne.trim(ptype)) cycle

   do ils=1,nlay
    fname=trim(dir)//'\'//trim(ptype)//'_l'//trim(imod_utl_itos(ils))//'.idf'
    n=len(fname); call pks7mpifname(fname,n)
!    if(ils.eq.nlay)then
!     if(trim(ptype).eq.'VC')cycle
!     if(trim(ptype).eq.'KV')cycle
!    endif
    xpp=a(:,:,ils)
    call met1wrtidf(fname,xpp,ncol,nrow,-999.0d0,iout)
   enddo
   exit
  enddo
  deallocate(xpp)
 endif


 end subroutine

 !###====================================================================
 LOGICAL FUNCTION PESTALPHA_PPOINT(ILS,PTYPE,XPP,NCOL,NROW,DIR,IOUT,IGLTYPE)
 !###====================================================================
 use imod_utl, only: imod_utl_printtext,imod_utl_itos,imod_utl_dtos,imod_utl_createdir, &
   utl_kriging_range,utl_kriging_main
 use gwfmetmodule, only: cdelr, cdelc, resultdir   
 use global, only: lipest, ibound 
 use pestvar, only: param, pest_iter,lgrad,llnsrch,pest_igrad,iupestout,pest_ktype,pest_krange,blnkout,pest_maxpnt,pest_ksettings
 IMPLICIT NONE
 INTEGER,INTENT(IN) :: NCOL,NROW,ILS,IOUT,IGLTYPE
 CHARACTER(LEN=2),INTENT(IN) :: PTYPE
 CHARACTER(LEN=*),INTENT(IN) :: DIR
 INTEGER :: I,J,K,NXYZ,IROW,ICOL,ILOG,N,ICYCLE,KTYPE
 REAL(KIND=8) :: FCT,NODATA,RANGE,SILL,NUGGET
 CHARACTER(LEN=256) :: FNAME,LINE
 REAL(KIND=8),DIMENSION(NCOL,NROW),INTENT(OUT) :: XPP
 REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: XYZ,XST
 
 PESTALPHA_PPOINT=.FALSE.
 
 SILL=-100.0D0
 
 DO K=1,2
  NXYZ=0
  DO I=1,SIZE(PARAM)
   
   IF(TRIM(PARAM(I)%PTYPE).NE.TRIM(PTYPE)) CYCLE
   !## skip NONE pilot-points
   IF(PARAM(I)%ZTYPE.NE.1)CYCLE
   !## not correct modellayer
   IF(PARAM(I)%ILS.NE.ILS)CYCLE

   !## always krige in log-space
   FCT=PARAM(I)%ALPHA(1)
   !## set log-type (be aware not possible to use mix of them)
   ILOG=PARAM(I)%ILOG
   
   !## convert standard deviation to variance
   SILL=PARAM(I)%STDEV**2.0D0
   
   !## make it negative
   IF(K.EQ.2)PARAM(I)%ZTYPE=-1*ABS(PARAM(I)%ZTYPE)
 
   IF(K.EQ.2)THEN
    IF(PARAM(I)%ILOG.EQ.1)THEN
     LINE=' * Module '//PARAM(I)%PTYPE//' adjusted ('//TRIM(IMOD_UTL_ITOS(SIZE(PARAM(I)%XY,1)))// &
      ') location(s) as PILOTPOINT with (LN) alpha='//TRIM(IMOD_UTL_dTOS(EXP(FCT),'F',7))
    ELSEIF(PARAM(I)%ILOG.EQ.2)THEN
     LINE=' * Module '//PARAM(I)%PTYPE//' adjusted ('//TRIM(IMOD_UTL_ITOS(SIZE(PARAM(I)%XY,1)))// &
        ') location(s) as PILOTPOINT with (LOG10) alpha='//TRIM(IMOD_UTL_dTOS(10.0**FCT,'F',7))
    ENDIF
    CALL imod_utl_printtext(TRIM(LINE),1) 
   ENDIF
    
   DO J=1,PARAM(I)%NODES !SIZE(PARAM(I)%XY,1)
    NXYZ=NXYZ+1
    IF(K.EQ.2)THEN
     !## kriging on log-values
     XYZ(NXYZ,1)=PARAM(I)%XY(J,1); XYZ(NXYZ,2)=PARAM(I)%XY(J,2); XYZ(NXYZ,3)=FCT 
    ENDIF
   ENDDO
  ENDDO
  IF(NXYZ.EQ.0)EXIT; IF(K.EQ.1)ALLOCATE(XYZ(NXYZ,3))
 ENDDO
 
 !# next parameter, this one is not used
 IF(NXYZ.EQ.0)RETURN
 
 NODATA=-999.99D0; XPP=NODATA
 RANGE =PEST_KRANGE
 KTYPE =PEST_KTYPE
 NUGGET=0.0D0

 !## set kriging setting per parameter - if available
 IF(PEST_KRANGE.LE.0.0)THEN
  !## set range to be automatic by default
  RANGE=0.0D0
  DO I=1,SIZE(PEST_KSETTINGS)
   !## settings for current type
   IF(TRIM(PTYPE).NE.PEST_KSETTINGS(I)%PTYPE)CYCLE
   IF(ILS.NE.PEST_KSETTINGS(I)%ILS)CYCLE
   !## overwrite current settings
   RANGE =PEST_KSETTINGS(I)%RANGE   
   KTYPE =PEST_KSETTINGS(I)%KTYPE
   IF(PEST_KTYPE.LT.0)KTYPE=-1*KTYPE
   SILL  =PEST_KSETTINGS(I)%SILL
   NUGGET=PEST_KSETTINGS(I)%NUGGET
   EXIT
  ENDDO
 ENDIF
 
 IF(SILL.LE.0.0D0)STOP 'SILL NOT DEFINED'

 iF(ABS(PEST_KTYPE).EQ.6)THEN
  CALL imod_utl_printtext('Nearest Neighbour is applied',1)
 else
  IF(RANGE.GT.0.0)THEN
   CALL imod_utl_printtext('Ordinary Kriging applied Range:'//TRIM(IMOD_UTL_DTOS(RANGE,'F',2))//' meter',1)
  ELSE
   CALL imod_utl_printtext('Ordinary Kriging applied Range Automatically',1)
  ENDIF
 ENDIF
  
 WRITE(*,'(A3,2A5,3A15  ,A10)') 'PRM','ILS','TYPE','NUGGET','SILL','RANGE','POINTS'
 WRITE(*,'(A3,2I5,3F15.3,I10)') PTYPE,ILS,KTYPE,NUGGET,SILL,RANGE,NXYZ

 ALLOCATE(XST(NCOL,NROW)); XST=0.0D0
 
 !## apply kriging interpolation
 CALL UTL_KRIGING_MAIN(NXYZ,XYZ(:,1),XYZ(:,2),XYZ(:,3),CDELR,CDELC,NROW,NCOL,XPP,NODATA,RANGE,KTYPE, &
    BLNKOUT,PEST_MAXPNT,XST,SILL,NUGGET)

 DO IROW=1,NROW; DO ICOL=1,NCOL
  IF(IGLTYPE.EQ.1.AND.IBOUND(ICOL,IROW,ILS).EQ.0)THEN
   XPP(ICOL,IROW)=NODATA
  ELSE
   IF(BLNKOUT%X(ICOL,IROW).LE.0.0)XPP(ICOL,IROW)=0.0
   SELECT CASE (ILOG)
    CASE (1)
     XPP(ICOL,IROW)=EXP(XPP(ICOL,IROW))
!     XST(ICOL,IROW)=EXP(XST(ICOL,IROW))
    CASE (2)
     XPP(ICOL,IROW)=10.0D0**XPP(ICOL,IROW)
!     XST(ICOL,IROW)=10.0D0**XST(ICOL,IROW)
   END SELECT
   !## areas outside range do get a value of 0.0, convert to 1.0
   IF(XPP(ICOL,IROW).EQ.0.0)THEN
    XPP(ICOL,IROW)=1.0D0
    XST(ICOL,IROW)=0.0D0
   ENDIF
  ENDIF
 ENDDO; ENDDO
 
 READ(DIR(INDEX(DIR,'cycle',.TRUE.)+5:),*) ICYCLE
 FNAME=DIR(:INDEX(DIR,'\',.TRUE.)-1)//CHAR(92)//'FACTORS'//TRIM(IMOD_UTL_ITOS(ICYCLE))//'\KRIGING_'//TRIM(PTYPE)//'_L'//TRIM(IMOD_UTL_ITOS(ILS))//'.IDF'
 N=LEN(FNAME); CALL PKS7MPIFNAME(FNAME,N); CALL MET1WRTIDF(FNAME,XPP,NCOL,NROW,NODATA,IOUT)
 FNAME=DIR(:INDEX(DIR,'\',.TRUE.)-1)//CHAR(92)//'FACTORS'//TRIM(IMOD_UTL_ITOS(ICYCLE))//'\KRIGING_VARIANCE_'//TRIM(PTYPE)//'_L'//TRIM(IMOD_UTL_ITOS(ILS))//'.IDF'
 N=LEN(FNAME); CALL PKS7MPIFNAME(FNAME,N); CALL MET1WRTIDF(FNAME,XST,NCOL,NROW,NODATA,IOUT)
  
 DEALLOCATE(XYZ,XST)
 
 PESTALPHA_PPOINT=.TRUE.
  
 END FUNCTION PESTALPHA_PPOINT
  
!###====================================================================
subroutine pest1alpha_list(ptype,nlist,rlist,ldim,mxlist,iopt1,iopt2,iout)
!###====================================================================

! modules
use global, only: lipest,buff,ncol,nrow,npertxt
use imod_utl, only: imod_utl_printtext, imod_utl_itos, imod_utl_dtos, imod_utl_createdir
use m_mf2005_main, only: kper
use gwfmetmodule, only: resultdir
use pestvar, only: param, pest_iter, iupestout

implicit none

! arguments
character(len=2), intent(in) :: ptype
integer, intent(in) :: nlist, ldim, mxlist,iout
real, dimension(ldim,mxlist), intent(inout) :: rlist
integer, intent(in) :: iopt1, iopt2
logical :: pestalpha_ppoint
! locals
character(len=256) :: line,dir
character(len=1024) :: errmsg
integer :: i, j, k, ils, irow, icol, idat, irivsubsys, irivrfact,ilay, &
           idrnsubsys, iwelsubsys, ihfbfact, ighbsubsys, nadj,inode
integer,dimension(4) :: ios
integer(kind=8) :: sdate_par,edate_par,sdate_mod,edate_mod
real(kind=8) :: ppart, fct
real(kind=8),allocatable,dimension(:,:) :: xpp

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

if(.not.lipest)return

!## make sure all are positive again
do i=1,size(param)
 param(i)%ztype=abs(param(i)%ztype)
enddo

!## first, mark the cells
do i=1,size(param)

  !## not this package - skip it
  if (trim(param(i)%ptype).ne.trim(ptype)) cycle
  !## allready processed in combination with pilotpoints
  if(param(i)%ztype.eq.-1)cycle
  
  !## check date if needed
  if(param(i)%sdate.ne.''.and.param(i)%edate.ne.'')then
   if(trim(param(i)%sdate).eq.'STEADY-STATE'.and.trim(param(i)%edate).eq.'STEADY-STATE')then
    if(npertxt(kper,1).ne.'STEADY-STATE')cycle
   else
    read(npertxt(kper,2),*,iostat=ios(1)) sdate_mod
    read(npertxt(kper,3),*,iostat=ios(2)) edate_mod
    read(param(i)%sdate ,*,iostat=ios(3)) sdate_par
    read(param(i)%edate ,*,iostat=ios(4)) edate_par
    !## ready to evaluate
    if(sum(ios).eq.0)then
     !## make sure all are in same dimensions
     if(sdate_mod.lt.99999999)sdate_mod=sdate_mod*1e6
     if(edate_mod.lt.99999999)edate_mod=edate_mod*1e6
     if(sdate_par.lt.99999999)sdate_par=sdate_par*1e6
     if(edate_par.lt.99999999)edate_par=edate_par*1e6
     if(edate_par.le.sdate_mod)cycle
     if(sdate_par.gt.edate_mod)cycle
    else
     cycle
    endif
   endif
  endif
  
  !## system/layer
  ils=param(i)%ils 

  !## fill buff with location of zone
  buff=0.0

  !## get 2d-array of pilot points for similar parameter (if needed)
  if(param(i)%ztype.eq.1)then
   allocate(xpp(ncol,nrow))
   DIR=trim(resultdir)//'\pest\parameters_cycle'//trim(imod_utl_itos(pest_iter))
   CALL IMOD_UTL_CREATEDIR(DIR)
   if(.not.pestalpha_ppoint(ils,ptype,xpp,ncol,nrow,dir,iout,2))then
    if(allocated(xpp))deallocate(xpp); cycle
   endif
   buff=1.0
  else
   do j=1,param(i)%nodes
    irow=param(i)%irow(j); icol=param(i)%icol(j)
    !## save zone location number
    buff(icol,irow,1)=1.0 !real(j,8)
   enddo
   !## factor
   fct=param(i)%alpha(1)
   if(param(i)%ilog.eq.1)fct=exp(fct)
   if(param(i)%ilog.eq.2)fct=10.0**fct
  endif
  
  !## adjust
  nadj = 0

  !## modify package with location in zone
  select case (trim(ptype))
   case('RC','IC') ! river/isg conductances
     errmsg = 'Cannot apply PEST scaling factor for river/isg conductance'
     irivsubsys = iopt1
     if (irivsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)

     idat = 5
     if (trim(ptype).eq.'IC') ils=-ils
     do j = 1, nlist
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(irivsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('RL','IL') ! river/isg river levels
     errmsg = 'Cannot apply PEST scaling factor for river/isg levels'
     irivsubsys = iopt1
     if (irivsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)

     idat = 4
     if (trim(ptype).eq.'IL') ils=-ils
     do j = 1, nlist
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(irivsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('RB','IB') ! river/isg bottom levels
     errmsg = 'Cannot apply PEST scaling factor for river/isg bottom levels'
     irivsubsys = iopt1
     if (irivsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)

     idat = 6
     if (trim(ptype).eq.'IB') ils=-ils
     do j = 1, nlist
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(irivsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('RI','II') ! river/isg infiltration factors
     errmsg = 'Cannot apply PEST scaling factor for river/isg infiltration factor'
     irivsubsys = iopt1; irivrfact = iopt2
     if (irivrfact.eq.0.or.irivsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)

     idat = irivrfact
     if (trim(ptype).eq.'II') ils=-ils
     do j = 1, nlist
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(irivsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
         !## infiltration factor cannot become > 1.0
         rlist(idat,j)=min(1.0,rlist(idat,j))
        endif
     end do

   case('DC') ! drain conductances
     errmsg = 'Cannot apply PEST scaling factor for drain conductance'
     idrnsubsys = iopt1
     if (idrnsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)
     idat = 5
     do j = 1, nlist ! match sybsystem number
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(idrnsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('DL') ! drain level
     errmsg = 'Cannot apply PEST scaling factor for drain level'
!     if (.not.present(iopt1)) call imod_utl_printtext(trim(errmsg),2)
     idrnsubsys = iopt1
     if (idrnsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)
     idat = 4
     do j = 1, nlist ! match sybsystem number
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(idrnsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('GC') ! general head conductances
     errmsg = 'Cannot apply PEST scaling factor for general conductance'
     ighbsubsys = iopt1
     if (ighbsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)
     idat = 5
     do j = 1, nlist ! match sybsystem number
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(ighbsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('SF') ! sfr conductance
     errmsg = 'Cannot apply PEST scaling factor for sfr conductance'
     idat = 4
     do j = 1, nlist ! match sybsystem number
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        nadj=nadj+1
        if(abs(param(i)%ztype).eq.1)then
         rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
        else
         rlist(idat,j)=rlist(idat,j)*fct
        endif
     end do

   case('QR','FQ') ! well rates WEL and FHB package
     errmsg = 'Cannot apply PEST scaling factor for well rates'
     iwelsubsys = iopt1
     if (iwelsubsys.eq.0) call imod_utl_printtext(trim(errmsg),2)
     idat = 4
     do j = 1, nlist ! match sybsystem number
        irow=rlist(2,j); icol=rlist(3,j)
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        !## adjust for selected system
        if (ils.eq.int(rlist(iwelsubsys,j)))then
         nadj=nadj+1
         if(abs(param(i)%ztype).eq.1)then
          rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
         else
          rlist(idat,j)=rlist(idat,j)*fct
         endif
        endif
     end do

   case('MQ') ! multinode well rates
     errmsg = 'Cannot apply PEST scaling factor for multinode-well rates'
     idat = 3
     do j = 1, nlist ! match subsystem number
        irow=rlist(1,j); icol=rlist(2,j)
        !## skip as list can be longer than neccessary
        if(irow.le.0.or.icol.le.0)cycle
        !## not in current zone
        if(buff(icol,irow,1).eq.0.0d0)cycle
        nadj=nadj+1
        if(abs(param(i)%ztype).eq.1)then
         rlist(idat,j)=rlist(idat,j)*xpp(icol,irow)
        else
         rlist(idat,j)=rlist(idat,j)*fct
        endif
     end do

    case('HF')
      errmsg = 'Cannot apply PEST scaling factor for horizontal flow barrier'
      ihfbfact = iopt1
      if (ihfbfact.ne.1) call imod_utl_printtext(trim(errmsg),2)
      do j = 1, nlist
         !# from node in zone
         irow=rlist(2,j); icol=rlist(3,j); if(buff(icol,irow,1).eq.0.0d0)cycle
         !## to node  in zone
         irow=rlist(4,j); icol=rlist(5,j); if(buff(icol,irow,1).eq.0.0d0)cycle
         !## adjust for selected system
         if (int(rlist(7,j)).eq.ils) then  
            nadj = nadj + 1
            if(abs(param(i)%ztype).eq.1)then
             rlist(6,j)=rlist(6,j)*xpp(icol,irow)
            else
             rlist(6,j) = rlist(6,j)*fct
            endif     
         end if
      end do
   end select

   if(trim(param(i)%acronym).eq.'')then
     line='Param. '//trim(imod_utl_itos(i))//'['//param(i)%ptype//',n_adjusted='//trim(imod_utl_itos(nadj))// &
       ',f='//trim(imod_utl_dtos(fct,'f',7))//',ils='//trim(imod_utl_itos(ils))//']'
   else
      line='param. '//trim(imod_utl_itos(i))//'['//trim(param(i)%acronym)//','//param(i)%ptype//',n_adjusted='//trim(imod_utl_itos(nadj))// &
       ',f='//trim(imod_utl_dtos(fct,'f',7))//',ils='//trim(imod_utl_itos(ils))//']'
   endif
   write(*,'(a)') trim(line)
   if(iupestout.ne.0)write(iupestout,'(a)') trim(line)
   
   if(abs(param(i)%ztype).eq.1)then
    deallocate(xpp)
   endif

end do

!## make sure all are positive again
do i=1,size(param)
 param(i)%ztype=abs(param(i)%ztype)
enddo

end subroutine