!!  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.

module gwfdgfmodule
   integer, pointer, save :: ic1
   integer, pointer, save :: ic2
   integer, pointer, save :: ir1
   integer, pointer, save :: ir2
   integer, pointer, save :: il1
   integer, pointer, save :: il2
   real,    pointer, save :: ts
   integer, pointer, save :: ndgbnd
   integer, pointer, save :: ndgcells
   real, dimension(:,:), pointer, save :: dgbnd  
   integer, dimension(:,:), pointer, save :: dgcells  
   type gwfdgftype
      integer, pointer :: ic1
      integer, pointer :: ic2
      integer, pointer :: ir1
      integer, pointer :: ir2
      integer, pointer :: il1
      integer, pointer :: il2
      real,    pointer :: ts
      integer, pointer :: ndgbnd
      integer, pointer :: ndgcells
      real, dimension(:,:), pointer :: dgbnd
      integer, dimension(:,:), pointer :: dgcells  
   end type

   type(gwfdgftype), save:: gwfdgfdat(10)
end module gwfdgfmodule

subroutine gwf2dgf1ar(in,igrid,iout)
 use gwfdgfmodule

 implicit none

! arguments
 integer, intent(in) :: in          !> unit number of input file
 integer, intent(in) :: igrid       !> grid number
 integer, intent(in) :: iout

! locals
 integer :: i

! ------------------------------------------------------------------------------
! set pointer

! package id
 write(iout,'(/,9x,a,i4)') ' DGF -- DGFLOW Exchange PACKAGE, VERSION 1'// &
                   ', 5 April 2011 INPUT READ FROM UNIT',in

! allocate scalars
 allocate(ic1,ic2,ir1,ir2,il1,il2,ts)
 read(in,*) ic1,ic2,ir1,ir2,il1,il2,ts

 call sgwf2dgf1psv(igrid)

 return
end

subroutine gwf2dgf1init(ibound,ncol,nrow,nlay,igrid)

 use gwfdgfmodule

 implicit none

! arguments 
 integer,intent(in) :: ncol,nrow,nlay,igrid
 integer, dimension(ncol,nrow,nlay), intent(inout) :: ibound
 
! locals 
 logical :: flg
 integer :: icol, irow, ilay, iact, n
 integer, dimension(:,:,:), allocatable :: iwrk
! ------------------------------------------------------------------------------

 call sgwf2dgf1pnt(igrid)
 
 allocate(ndgbnd,ndgcells)

 ! set ibound
 do ilay = il1+1,il2-1
    do irow = ir1+1,ir2-1
       do icol = ic1+1, ic2-1
          ibound(icol,irow,ilay) = 0
        end do   
    end do      
 end do
 
 allocate(iwrk(ncol,nrow,nlay))
 iwrk = 0
 do ilay = il1+1,min(il2,nlay)
    do irow = ir1,ir2
       do icol = ic1, ic2
          iwrk(icol,irow,ilay) = 1 ! label DGFLOW cells
        end do   
    end do      
 end do
 
 do iact = 1, 2
    ndgcells = 0
    do ilay = 1, nlay
       do irow = 1, nrow
          do icol = 1, ncol
             if (iwrk(icol,irow,ilay).eq.1) then
                ndgcells = ndgcells + 1
                if (iact.eq.2) then
                   dgcells(1,ndgcells) = ilay 
                   dgcells(2,ndgcells) = irow 
                   dgcells(3,ndgcells) = icol  
                   n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
                   !write(*,*) 'n dgcells=',n
                end if    
             end if             
          end do       
       end do          
    end do 
    if (iact.eq.1) then
       allocate(dgcells(3,max(ndgcells,1)))
    end if   
 end do
 
 do ilay = il1+1,il2-1
    do irow = ir1+1,ir2-1
       do icol = ic1+1, ic2-1
           iwrk(icol,irow,ilay) = 0
       end do   
    end do      
 end do

 ! count and fill
 do iact = 1, 2
    ndgbnd = 0
    do ilay = 1, nlay
       flg = .false. 
       do irow = 1, nrow
          do icol = 1, ncol
             if (iwrk(icol,irow,ilay).eq.1) then
                flg = .true. 
                ndgbnd = ndgbnd + 1      
                if (iact.eq.2) then
                   dgbnd(1,ndgbnd) = ilay 
                   dgbnd(2,ndgbnd) = irow 
                   dgbnd(3,ndgbnd) = icol 
                   dgbnd(4,ndgbnd) = 1 
                   n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
                   !write(*,*) 'n dgcells=',n
                   dgbnd(5,ndgbnd) = 0.0 ! init Q 
                end if
             end if   
          end do    
       end do
       if (flg.and.ilay.ne.il2) then
          do irow = 1, nrow
             do icol = 1, ncol
                if (iwrk(icol,irow,ilay).eq.1) then
                   flg = .true. 
                   ndgbnd = ndgbnd + 1      
                   if (iact.eq.2) then
                      dgbnd(1,ndgbnd) = ilay 
                      dgbnd(2,ndgbnd) = irow 
                      dgbnd(3,ndgbnd) = icol 
                      dgbnd(4,ndgbnd) = -1 
                      n = icol + (irow-1)*ncol + (ilay-1)*nrow*ncol
                      dgbnd(5,ndgbnd) = 0.0 ! init Q 
                   end if
                end if   
             end do    
          end do
       end if
    end do 
    if (iact.eq.1) then
       allocate(dgbnd(5,max(ndgbnd,1)))   
    end if    
 end do

 call sgwf2bas7psv(igrid) 
 call sgwf2dgf1psv(igrid)
 
 deallocate(iwrk)
 
 return
end  

subroutine gwf2dgf1fm(igrid)

 use global,       only: ibound, rhs
 use gwfdgfmodule, only: ndgbnd, dgbnd

 implicit none
 
! arguments
 integer, intent(in) :: igrid

! locals 
 integer :: i, ir, ic, il, flg
 real :: q
! ------------------------------------------------------------------------------

 call sgwf2dgf1pnt(igrid)
 
 if (ndgbnd.le.0) return
 
 do i = 1, ndgbnd
    il  = dgbnd(1,i)
    ir  = dgbnd(2,i)
    ic  = dgbnd(3,i)
    flg = dgbnd(3,i) 
    q   = dgbnd(4,i)  
    if (ibound(ic,ir,il).gt.0.and.flg.eq.1) then
       rhs(ic,ir,il) = rhs(ic,ir,il) - q    
    end if
 end do
 
 return 
end

subroutine sgwf2dgf1psv(igrid)
 use gwfdgfmodule
 implicit none
 integer, intent(in) :: igrid

 gwfdgfdat(igrid)%ic1      => ic1
 gwfdgfdat(igrid)%ic2      => ic2
 gwfdgfdat(igrid)%ir1      => ir1
 gwfdgfdat(igrid)%ir2      => ir2
 gwfdgfdat(igrid)%il1      => il1
 gwfdgfdat(igrid)%il2      => il2
 gwfdgfdat(igrid)%ts       => ts
 gwfdgfdat(igrid)%ndgbnd   => ndgbnd
 gwfdgfdat(igrid)%ndgcells => ndgcells
 gwfdgfdat(igrid)%dgbnd    => dgbnd
 gwfdgfdat(igrid)%dgcells  => dgcells
 return
end

subroutine sgwf2dgf1pnt(igrid)
 use gwfdgfmodule
 implicit none
 integer, intent(in) :: igrid

 ic1      => gwfdgfdat(igrid)%ic1
 ic2      => gwfdgfdat(igrid)%ic2
 ir1      => gwfdgfdat(igrid)%ir1
 ir2      => gwfdgfdat(igrid)%ir2
 il1      => gwfdgfdat(igrid)%il1
 il2      => gwfdgfdat(igrid)%il2
 ts       => gwfdgfdat(igrid)%ts
 ndgbnd   => gwfdgfdat(igrid)%ndgbnd
 ndgcells => gwfdgfdat(igrid)%ndgcells
 dgbnd    => gwfdgfdat(igrid)%dgbnd
 dgcells  => gwfdgfdat(igrid)%dgcells
 return
end

subroutine gwf2dgf1da(igrid)
 use gwfdgfmodule
 implicit none
 integer, intent(in) :: igrid

 call sgwf2dgf1pnt(igrid)
 deallocate(ic1)
 deallocate(ic2)
 deallocate(ir1)
 deallocate(ir2)
 deallocate(il1)
 deallocate(il2)
 deallocate(ts)
 deallocate(ndgbnd)
 deallocate(ndgcells)
 deallocate(dgbnd)
 deallocate(dgcells)
 return
end