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

      module lcdmodule

      use global, only: ncol, nrow

      implicit none

      logical, save :: lcdinit = .true.
      logical, save :: lqd
      real, save :: simcsize, xmin, ymin, xmax, ymax

      integer, dimension(:), allocatable, save :: genip
      integer, dimension(:,:), allocatable, save :: genpos

      end module lcdmodule

      module rdlcd_interface

      implicit none

      interface

      subroutine rdlcd(nlist,rlist,lstbeg,ldim,in,iout,label,ncol,nrow,n
     1lay)

      implicit none

c arguments
      integer, intent(inout) :: nlist
      integer, intent(in) :: ldim
      real, dimension(:,:), pointer,intent(inout) :: rlist
      real, dimension(:,:), pointer :: rlisttmp
      integer, intent(in) :: lstbeg
      integer, intent(in) :: in
      integer, intent(in) :: iout
      character(len=*), dimension(*), intent(in) :: label
      integer, intent(in) :: ncol
      integer, intent(in) :: nrow
      integer, intent(in) :: nlay

      end subroutine

      end interface

      end module rdlcd_interface

      subroutine rdlcd(mxlist,rlist,lstbeg,ldim,in,iout,label,ncol,nrow,
     1nlay)

c description:
c ------------------------------------------------------------------------------
c read data block of Line Column Data type.
c

c declaration section
c ------------------------------------------------------------------------------
      use imod_utl, only: imod_utl_openasc
      use gwfmetmodule, only: cdelr,cdelc
      use lcdmodule, only: genip, genpos

      implicit none

c arguments
      integer, intent(in) :: ldim
      integer, intent(inout) :: mxlist
      real, dimension(:,:), pointer,intent(inout) :: rlist
      real, dimension(:,:), pointer :: rlisttmp
      integer, intent(in) :: lstbeg
      integer, intent(in) :: in
      integer, intent(in) :: iout
      character(len=*), dimension(*), intent(in) :: label
      integer, intent(in) :: ncol
      integer, intent(in) :: nrow
      integer, intent(in) :: nlay

c local variables
      integer :: icol, irow, jcol, jrow, igen, iact, ngen, ii, kk,
     1           ip1, ip2,ic1, ic2, ir1, ir2, j, jj, nlist
      integer ::  il, is, ie, iline, jline, nline
      integer :: ilay,il1,il2,jlay
      real :: fct,HFB1EXPORT_GETDZ,c,z,c1,c2,zz
      real, dimension(:,:), allocatable :: tmp,tf,bf,res,fdz
      integer(kind=1), dimension(:,:,:), allocatable :: ipc
      character(len=1024) :: fname
      integer,dimension(:),allocatable :: lun,dun,nlun
      character(len=1024) :: str

c parameters
      character(len=24) :: aname(1)
      data aname(1) /'                     LCD'/
      integer,parameter :: ineighbours=2
      REAL :: ZF,NODATA

c program section
c ------------------------------------------------------------------------------
c init
      call initlcd()

      NODATA=HUGE(1.0)

c read number of hfb layers
      read(in,*) ngen

c allocate ipc
      allocate(tmp(ncol,nrow))
      
      IF(ALLOCATED(IPC))DEALLOCATE(IPC); ALLOCATE(IPC(NCOL,NROW,2))
      IF(ALLOCATED(TF))DEALLOCATE(TF); IF(ALLOCATED(BF))DEALLOCATE(BF)
      ALLOCATE(TF(NCOL,NROW),BF(NCOL,NROW))

      allocate(lun(nlay),dun(nlay),nlun(nlay)); lun=0; dun=0; nlun=0

      do ilay=1,nlay
       write(fname,'(a,i3.3,a)') 'hfb_l',ilay,'.gen'
       call imod_utl_openasc(lun(ilay),fname,'w')
       write(fname,'(a,i3.3,a)') 'hfb_l',ilay,'.dat'
       call imod_utl_openasc(dun(ilay),fname,'w')
       WRITE(dun(ilay),'(A)') 'no,net_resis,tot_resis,fraction'
      enddo
      
c count number of hfb and fill
      do iact = 1, 2
         if (iact == 2) then
            do while(.true.)
               backspace(in)
               read(in,'(a)') str
               str = adjustl(str)
               call upcase(str)
               if (str(1:3).eq.'LCD') then
                  exit
               else
                  backspace(in)
               end if
            end do
            read(in,*) ngen
         end if
         ii = lstbeg-1
         !## process each genfile
         do igen = 1, ngen

            read(in,*) ilay, fct

            kk = ilay
            !## fill genpos list
            call u2drel(tmp,aname(1),nrow,ncol,kk,in,iout)

            nline = size(genip)-1

            do j = 1, nline 
               
             ipc = int(0,1)
             is = genip(j-1)+1; ie = genip(j)

             !## line not in current model dimensions
             if(ie.eq.0)cycle
             if(is.gt.size(genpos,1).or.ie.gt.size(genpos,1))cycle
             
             IF(ilay.EQ.0)THEN; TF=NODATA; BF=NODATA; ENDIF
 
             jcol = genpos(is,1); jrow = genpos(is,2);

             !## startpoint
             JLINE=GENPOS(is,3)
 
             !## get top/bottom elevations on grid
             IF(ilay.EQ.0)THEN
              !## process line
              do il = is+1, ie
               ILINE=GENPOS(il,3)
               IF(ILINE.EQ.JLINE)THEN
                IC1=GENPOS(il-1,1); IC2=GENPOS(il  ,1)
                IR1=GENPOS(il-1,2); IR2=GENPOS(il  ,2)
                IP1=GENPOS(il-1,4); IP2=GENPOS(il  ,4)
                !## fill in z-coordinates
                ZF=(genpos(IL-1,5)+genpos(IL,5))/2.0
                DO IROW=IR1,IR2; DO ICOL=IC1,IC2
                 IF(TF(ICOL,IROW).EQ.NODATA)THEN
                  TF(ICOL,IROW)=ZF/100.0
                  BF(ICOL,IROW)=ZF/100.0
                 ELSE
                  TF(ICOL,IROW)=MAX(TF(ICOL,IROW),ZF/100.0)
                  BF(ICOL,IROW)=MIN(BF(ICOL,IROW),ZF/100.0)
                 ENDIF
                ENDDO; ENDDO
               ENDIF
               IF(ILINE.NE.JLINE.OR.il.EQ.ie)JLINE=GENPOS(J,3)
              ENDDO
             ENDIF

             jcol = genpos(is,1); jrow = genpos(is,2);

             !## startpoint
             JLINE=GENPOS(is,3)

             do il = is+1, ie

              ILINE=GENPOS(il,3)
              IF(ILINE.EQ.JLINE)THEN
               IC1=GENPOS(il-1,1); IC2=GENPOS(il  ,1)
               IR1=GENPOS(il-1,2); IR2=GENPOS(il  ,2)
               IP1=GENPOS(il-1,4); IP2=GENPOS(il  ,4)
               CALL HFBGETFACES(IC1,IC2,IR1,IR2,IP1,IP2,IPC,NROW,NCOL)
              ENDIF
              
              IF(ILINE.NE.JLINE.OR.IL.EQ.IE)THEN
                      
               !## determine what layer(s)
               IF(ilay.EQ.0)THEN
                IL1=1; IL2=NLAY
               ELSE
                IL1=ilay; IL2=IL1
               ENDIF

               do irow = 1, nrow
                do icol = 1, ncol
                 !## place horizontal wall
                 if (irow.lt.nrow) then
                  if(ipc(icol,irow,2).eq.int(1,1)) then
                   
                   !## x-direction
                   DO jLAY=IL1,IL2
                    Z=1.0
                    IF(ilay.EQ.0)THEN   !ICOL,IROW,ICOL,IROW+1
                     Z=HFB1EXPORT_GETDZ(TF,BF,ICOL,IROW,ICOL,IROW+1,NODA
     1TA,jLAY,NCOL,NROW)                  
                    ENDIF
                    
                    !## skip fault on side of model or less than 0.0 fraction
                    IF(Z.LE.0.0.OR.IROW+1.GT.NROW)CYCLE
                    
                    ii = ii + 1
                    !## store into memory second cycli
                    if (iact.eq.2) then

                     rlisttmp(1,ii) = jlay
                     rlisttmp(2,ii) = irow
                     rlisttmp(3,ii) = icol
                     rlisttmp(4,ii) = irow+1
                     rlisttmp(5,ii) = icol                  
                     rlisttmp(6,ii) = FCT
                     rlisttmp(7,ii) = Z 
                    
                    endif
                   enddo
                  end if
                 end if
                 
                 !## place vertical wall
                 if (icol.lt.ncol) then
                  if(ipc(icol,irow,1).eq.int(1,1)) then
                   
                   !## y-direction
                   DO jLAY=IL1,IL2
 
                    Z=1.0
                    IF(ilay.EQ.0)THEN  !ICOL,IROW,ICOL+1,IROW
                     Z=HFB1EXPORT_GETDZ(TF,BF,ICOL,IROW,ICOL+1,IROW,NODA
     1TA,jLAY,NCOL,NROW)                  
                    ENDIF

                    !## skip fault on side of model or less than 0.0 fraction
                    IF(Z.LE.0.0.OR.ICOL+1.GT.NCOL)CYCLE

                    ii = ii + 1
                    if (iact.eq.2) then
                     rlisttmp(1,ii) = jlay
                     rlisttmp(2,ii) = irow
                     rlisttmp(3,ii) = icol
                     rlisttmp(4,ii) = irow
                     rlisttmp(5,ii) = icol+1
                     rlisttmp(6,ii) = FCT
                     rlisttmp(7,ii) = Z 
                    end if
                   enddo
                  end if
                 end if
                end do ! icol
               end do ! irow

               !## reset for the next line   
               IPC=INT(0,1); JLINE=GENPOS(is,3)            

              endif  !if(iline.eq.jline)then
             end do ! do il = is, ie
            
            end do !  iline = 1, nline ! iline

         end do ! igen

         nlist = ii-lstbeg+1
         if (iact.eq.1) then
            if (.not.associated(rlisttmp)) then
               allocate(rlisttmp(ldim,nlist))
            end if
         end if ! iact = 1
      end do ! iact

!    construct final fault list per model layer
      allocate(fdz(ncol,nrow),res(ncol,nrow))

      do iact=1,2      
       
       !## number of adjusted hfb-elements
       jj=0
       
       !## process each layer
       DO ILAY=1,NLAY
  
        IPC=INT(0,1)
        RES=0.0
        FDZ=0.0
      
        DO ii=1,nlist
         jlay=int(rlisttmp(1,ii))
         !# not current modellayer
         if(jlay.ne.ilay)cycle
         ir1=int(rlisttmp(2,ii))
         ic1=int(rlisttmp(3,ii))
         ir2=int(rlisttmp(4,ii))
         ic2=int(rlisttmp(5,ii))
         c=      rlisttmp(6,ii)
         z=      rlisttmp(7,ii)
         IF(IC1.EQ.IC2)THEN
          IPC(IC1,IR1,2)=INT(1,1)
         ELSE
          IPC(IC1,IR1,1)=INT(1,1)
         ENDIF
         
         !## still some space left in modellayer for an additional fault
         if(fdz(ic1,ir1).lt.1.0)then
          !## available space
          zz=1.0-fdz(ic1,ir1)
          !## net available space
          zz=min(zz,z)
          !## resistance, sum conductances - ignore resistance of zero days
          IF(C.GT.0.0)RES(IC1,IR1)=RES(IC1,IR1)+(1.0/C)*zz
          !## occupation fraction
          FDZ(IC1,IR1)=FDZ(IC1,IR1)+Z
         endif
         
        ENDDO
      
        DO IROW=1,NROW; DO ICOL=1,NCOL

         !## place vertical wall (block in y-direction)
         IF(IPC(ICOL,IROW,1).EQ.INT(1,1))THEN
          IF(ICOL.LT.NCOL)THEN

           !## transform conductances to resistance
           C1=1.0/RES(ICOL,IROW)
         
           !## get total resistance related to thickness of model layer
           C2=C1*FDZ(ICOL,IROW)**4.0

           !## write fault to gen- and datfile
           call HFB1EXPORT_WRITEGEN(2,lun(ilay),dun(ilay),nlun(ilay),ico
     1l,irow,C1,C2,FDZ(ICOL,IROW))
           
           jj=jj+1
           if(iact.eq.2)then
            rlist(1,jj)=ilay
            rlist(2,jj)=irow
            rlist(3,jj)=icol
            rlist(4,jj)=irow
            rlist(5,jj)=icol+1
            rlist(6,jj)=c2
            rlist(7,jj)=0.0
           endif
           
          ENDIF 
         ENDIF

         !## place horizontal wall (block in x-direction)
         IF(IPC(ICOL,IROW,2).EQ.INT(1,1))THEN
          IF(IROW.LT.NROW)THEN

           !## transform conductances to resistance
           C1=1.0/RES(ICOL,IROW)
         
           !## get total resistance related to thickness of model layer
           C2=C1*FDZ(ICOL,IROW)**4.0

           call HFB1EXPORT_WRITEGEN(1,lun(ilay),dun(ilay),nlun(ilay),ico
     1l,irow,C1,C2,FDZ(ICOL,IROW))
        
           jj=jj+1
           if(iact.eq.2)then
            rlist(1,jj)=ilay
            rlist(2,jj)=irow
            rlist(3,jj)=icol
            rlist(4,jj)=irow+1
            rlist(5,jj)=icol
            rlist(6,jj)=c2
            rlist(7,jj)=0.0
           endif
           
          ENDIF
         ENDIF
     
        ENDDO; ENDDO
       ENDDO
      
       if(iact.eq.1)then
        mxlist=jj
        deallocate(rlist); allocate(rlist(ldim,mxlist))
       endif
       
      enddo
    
      !## close files
      do ilay=1,nlay
       if(nlun(ilay).eq.0)then
        close(lun(ilay),status='delete')
        close(dun(ilay),status='delete')
       else
        write(lun(ilay),'(a)') 'end'; close(lun(ilay)); close(dun(ilay))
       endif
      enddo
      
c deallocate arrays
      deallocate(rlisttmp)
      if (allocated(ipc)) deallocate(ipc)
      if (allocated(tmp)) deallocate(tmp)
      IF(ALLOCATED(TF))DEALLOCATE(TF); IF(ALLOCATED(BF))DEALLOCATE(BF)
      deallocate(lun,dun,nlun)
      
c end of program
      return
      end

      !###====================================================================
      subroutine HFB1EXPORT_WRITEGEN(it,iu,ju,n,icol,irow,C,RES,FDZ)
      !###====================================================================
      use gwfmetmodule, only: cdelr,cdelc
      implicit none
      real,intent(in) :: C,RES,FDZ
      integer,intent(IN) :: IT,iu,ju,icol,irow
      integer,intent(INOUT) :: n
      
      if(it.eq.1)THEN
       n=n+1
       write(ju,'(i10,3(1x,e15.7))') n,c,res,fdz
       write(iu,'(i10,1x,e15.7)') n
       write(iu,'(2(f10.2,a1))') cdelr(icol-1),',',cdelc(irow)
       write(iu,'(2(f10.2,a1))') cdelr(icol),',',cdelc(irow)
       write(iu,'(a)') 'end'
      endif

      if(it.eq.2)then
       n=n+1
       write(ju,'(i10,3(1x,e15.7))') n,c,res,fdz
       write(iu,'(i10,1x,e15.7)') n
       write(iu,'(2(f10.2,a1))') cdelr(icol),',',cdelc(irow-1)
       write(iu,'(2(f10.2,a1))') cdelr(icol),',',cdelc(irow)
       write(iu,'(a)') 'end'
      end if

      end subroutine HFB1EXPORT_WRITEGEN

      !###====================================================================
      REAL FUNCTION HFB1EXPORT_GETDZ(TF,BF,IC1,IR1,IC2
     1,IR2,NODATA,ILAY,NCOL,NROW)
      !###====================================================================
      USE GLOBAL,ONLY : BOTM
      IMPLICIT NONE     
      INTEGER,INTENT(IN) :: IC1,IR1,IC2,IR2,ILAY,NCOL,NROW
      REAL,INTENT(IN) :: NODATA
      REAL,INTENT(IN),DIMENSION(NCOL,NROW) :: TF,BF
      REAL :: DZ,TFV,BFV,TPV,BTV,C1,C2,CT,FFCT
      INTEGER :: IL1,IL2
      
      HFB1EXPORT_GETDZ=0.0

      !## determine values
      IF(TF(IC1,IR1).NE.NODATA.AND.TF(IC2,IR2).NE.NODATA)THEN
       TFV=(TF(IC1,IR1)+TF(IC2,IR2))/2.0
      ELSEIF(TF(IC1,IR1).NE.NODATA)THEN
       TFV=TF(IC1,IR1)
      ELSE
       TFV=TF(IC2,IR2)
      ENDIF
      IF(BF(IC1,IR1).NE.NODATA.AND.BF(IC2,IR2).NE.NODATA)THEN
       BFV=(BF(IC1,IR1)+BF(IC2,IR2))/2.0
      ELSEIF(BF(IC1,IR1).NE.NODATA)THEN
       BFV=BF(IC1,IR1)
      ELSE
       BFV=BF(IC2,IR2)
      ENDIF

      !## get internal layer number of vector of top/bottom information
      IL1=(ILAY*2)-1
      IL2=(ILAY*2)
      !## array starts at 0
      IL1=IL1-1
      IL2=IL2-1
      
      TPV=(BOTM(IC1,IR1,IL1)+BOTM(IC2,IR2,IL1))/2.0
      BTV=(BOTM(IC1,IR1,IL2)+BOTM(IC2,IR2,IL2))/2.0

      !## nett appearance of fault in modellayer
      DZ=MIN(TFV,TPV)-MAX(BFV,BTV)

      !## not in current modellayer
      IF(DZ.LT.0.0)RETURN

      IF(TPV-BTV.GT.0.0)THEN
       !## fraction of fault in modellayer
       DZ=DZ/(TPV-BTV)
      ENDIF
          
      HFB1EXPORT_GETDZ=DZ

      END FUNCTION HFB1EXPORT_GETDZ

      !###====================================================================
      REAL FUNCTION HFB1EXPORT_GETFACTOR(FCT,TF,BF,IC1,IR1,IC2,IR2,NODAT
     1A,ILAY,NCOL,NROW)
      !###====================================================================
      USE GLOBAL,ONLY : BOTM
      IMPLICIT NONE     
      INTEGER,INTENT(IN) :: IC1,IR1,IC2,IR2,ILAY,NCOL,NROW
      REAL,INTENT(IN) :: FCT,NODATA
      REAL,INTENT(IN),DIMENSION(NCOL,NROW) :: TF,BF
      REAL :: DZ,TFV,BFV,TPV,BTV,C1,C2,CT,FFCT
      INTEGER :: IL1,IL2
      
      HFB1EXPORT_GETFACTOR=0.0

      !## determine values
      IF(TF(IC1,IR1).NE.NODATA.AND.TF(IC2,IR2).NE.NODATA)THEN
       TFV=(TF(IC1,IR1)+TF(IC2,IR2))/2.0
      ELSEIF(TF(IC1,IR1).NE.NODATA)THEN
       TFV=TF(IC1,IR1)
      ELSE
       TFV=TF(IC2,IR2)
      ENDIF
      IF(BF(IC1,IR1).NE.NODATA.AND.BF(IC2,IR2).NE.NODATA)THEN
       BFV=(BF(IC1,IR1)+BF(IC2,IR2))/2.0
      ELSEIF(BF(IC1,IR1).NE.NODATA)THEN
       BFV=BF(IC1,IR1)
      ELSE
       BFV=BF(IC2,IR2)
      ENDIF

      !## get internal layer number of vector of top/bottom information
      IL1=(ILAY*2)-1
      IL2=(ILAY*2)
      !## array starts at 0
      IL1=IL1-1
      IL2=IL2-1
      
      TPV=(BOTM(IC1,IR1,IL1)+BOTM(IC2,IR2,IL1))/2.0
      BTV=(BOTM(IC1,IR1,IL2)+BOTM(IC2,IR2,IL2))/2.0

      !## nett appearance of fault in modellayer
      DZ=MIN(TFV,TPV)-MAX(BFV,BTV)

      !## not in current modellayer
      IF(DZ.LT.0.0)RETURN

      IF(TPV-BTV.GT.0.0)THEN
       !## fraction of fault in modellayer
       DZ=DZ/(TPV-BTV)
      ENDIF
      
      !## if dz.eq.0, modellayer has thickness of zero, but fault to be retained
      IF(DZ.EQ.0.0)DZ=1.0
      
      !## resistance of fault
      FFCT=FCT; IF(FCT.EQ.0.0)FFCT=10.0E10

      !## factor declines quadratically with layer occupation
      HFB1EXPORT_GETFACTOR=FFCT*DZ**4.0

      END FUNCTION HFB1EXPORT_GETFACTOR

      !###====================================================================
      SUBROUTINE HFBGETFACES(IC1,IC2,IR1,IR2,IP1,IP2,IPC,NROW,NCOL)
      !###====================================================================
      IMPLICIT NONE
      INTEGER,INTENT(IN) :: IC1,IC2,IR1,IR2,IP1,IP2
      INTEGER,INTENT(IN) :: NROW,NCOL
      INTEGER(KIND=1),INTENT(INOUT),DIMENSION(NCOL,NROW,2) :: IPC
      INTEGER,DIMENSION(2) :: JPC,JPR,JC,JR,JP
      INTEGER :: I,IC,IR 

      JC(1)=IC1; JC(2)=IC2
      JR(1)=IR1; JR(2)=IR2
      JP(1)=IP1; JP(2)=IP2

      DO I=1,2
       IF(JP(I).EQ.2.OR.JP(I).EQ.3)JPC(I)=JC(I)
       IF(JP(I).EQ.1.OR.JP(I).EQ.4)JPC(I)=JC(I)-1
       IF(JP(I).EQ.1.OR.JP(I).EQ.2)JPR(I)=JR(I)-1
       IF(JP(I).EQ.3.OR.JP(I).EQ.4)JPR(I)=JR(I)
      ENDDO

      !## do nothing, is similar point
      IF(JPR(1).EQ.JPR(2).AND.JPC(1).EQ.JPC(2))RETURN 

      !## do nothing whenever jpc.eq.0 or jpr.eq.0
      IF(JPC(1).EQ.0.OR.JPC(2).EQ.0)RETURN
      IF(JPR(1).EQ.0.OR.JPR(2).EQ.0)RETURN

      !## horizontal fault ipc(,,1)=1
      IF(JPR(1).EQ.JPR(2).AND.JPC(1).NE.JPC(2))THEN
       IC=MAX(JPC(1),JPC(2)); IR=JPR(1); IPC(IC,IR,2)=INT(1,1)
      ENDIF
      !## vertical fault ipc(,,2)=1
      IF(JPC(1).EQ.JPC(2).AND.JPR(1).NE.JPR(2))THEN
       IC=JPC(1); IR=MAX(JPR(1),JPR(2)); IPC(IC,IR,1)=INT(1,1)
      ENDIF
      !## diagonal, add two faults
      IF(JPR(1).NE.JPR(2).AND.JPC(1).NE.JPC(2))THEN
       !## goto to the west
       IF(JPC(1).GT.JPC(2))THEN
        !## goto to the north-west
        IF(JPR(1).GT.JPR(2))THEN
         IC=MIN(JPC(1),JPC(2)); IR=MAX(JPR(1),JPR(2))
         IPC(IC,IR,1)=INT(1,1) !## vertical
         IC=MAX(JPC(1),JPC(2)); IR=MAX(JPR(1),JPR(2))
         IPC(IC,IR,2)=INT(1,1) !## horizontal
        !## goto to the south-west
        ELSE
         IC=MIN(JPC(1),JPC(2)); IR=MAX(JPR(1),JPR(2))
         IPC(IC,IR,1)=INT(1,1) !## vertical
         IC=MAX(JPC(1),JPC(2)); IR=MIN(JPR(1),JPR(2))
         IPC(IC,IR,2)=INT(1,1) !## horizontal
        ENDIF
       !## goto to the east
       ELSE
        !## goto to the north-east
        IF(JPR(1).GT.JPR(2))THEN
         IC=MIN(JPC(1),JPC(2)); IR=MAX(JPR(1),JPR(2))
         IPC(IC,IR,1)=INT(1,1) !## vertical
         IC=MAX(JPC(1),JPC(2)); IR=MIN(JPR(1),JPR(2))
         IPC(IC,IR,2)=INT(1,1) !## horizontal   
        !## goto to the south-east
        ELSE
         IC=MIN(JPC(1),JPC(2)); IR=MAX(JPR(1),JPR(2))
         IPC(IC,IR,1)=INT(1,1) !## vertical
         IC=MAX(JPC(1),JPC(2)); IR=MAX(JPR(1),JPR(2))
         IPC(IC,IR,2)=INT(1,1) !## horizontal  
        ENDIF
       ENDIF
      ENDIF

      END SUBROUTINE HFBGETFACES

      subroutine initlcd()
c description:
c ------------------------------------------------------------------------------
c read file of LCD type
c

c declaration section
c ------------------------------------------------------------------------------
      use lcdmodule
      use global, only: iunit, delc, delr
      use gwfmetmodule
      use m_mf2005_iu, only: iumet

c arguments

c local variables
      integer :: icol, irow

c parameters

c program section
c ------------------------------------------------------------------------------
      if (lcdinit) then
         lcdinit = .false.
      else
         return
      end if

c check if metadata package is activated
      if (IUNIT(IUMET).le.0) then
         write(*,*) 'Error: lcd file, please activate met package.'
         call ustop(' ')
      end if
      if (.not.associated(coord_xll) .or.
     &    .not.associated(coord_yll) ) then
         write(*,*) 'Error: lcd file can only be used with xll and yll'
         call ustop(' ')
      end if

c check if grid is uniform
      lqd = .true.
      if ((maxval(delr).ne.minval(delr)).or.
     &    (maxval(delc).ne.minval(delc))) lqd = .false.
      if (lqd) then
         simcsize = delr(1)
      else
         write(*,*) 'Error: non-uniform grids not yet supported.'
         call ustop(' ')
      end if

      if (lqd) then
         xmin=cdelr(0)
         ymax=cdelc(0)
         xmax=xmin+(simcsize*ncol)
         ymin=ymax-(simcsize*nrow)
      endif

c end of program
      return
      end