module vtuMod
implicit none

! parameters
integer, parameter :: mxlen = 1024

logical :: lfirst = .true.
real, dimension(:,:), allocatable :: ggxnod
integer :: ngeom

integer, parameter :: iprconn = 0
integer, parameter :: iout = 0
integer :: nodes, neqs, nja
integer, dimension(:), allocatable :: ia, ja

integer, parameter :: maxnpvd = 100
integer :: npvd = 0
integer, dimension(maxnpvd) :: pvdun

integer :: vtknlay = 0

save

end module

subroutine met1ulasavvtk(text,buff,ncol,nrow,ilay,isplit,totim)

use imod_utl, only:imod_utl_createdir
use vtuMod, only: lfirst, ggxnod, mxlen, ia, ja, ngeom, npvd, pvdun, vtknlay

implicit none

! arguments
character(len=16), intent(in)          :: text
integer, intent(in)                    :: ncol, nrow
real, dimension(ncol,nrow), intent(in) :: buff
integer, intent(in)                    :: ilay
integer, intent(in)                    :: isplit
real, intent(in)                       :: totim

! functions
character(len=1024) :: met1fnamevtk, met1fnamepvd
integer :: cfn_getlun

! locals
character(len=mxlen) :: vtufile, pvdfile, s1, s2
integer :: irow, icol, n, nrc, iu, ios, i
logical :: lex

! body
if (lfirst) then
   call vtkinit()    
   lfirst = .false.
end if
if (ilay.eq.1) then ! reset
   ggxnod(11,:) = 0.       
end if
nrc = ncol*nrow
do irow = 1, nrow
   do icol = 1, ncol
      n = icol+(irow-1)*ncol+(ilay-1)*nrc    
      ggxnod(11,n) = buff(icol,irow)
   end do         
end do
if (ilay.eq.vtknlay) then
   vtufile = met1fnamevtk(isplit,text) 
   call writevtu(text,vtufile,ngeom,ia,ja,ggxnod) 
   pvdfile = met1fnamepvd(isplit,text)
   inquire(file=pvdfile,exist=lex,iostat=ios)
   if (ios.ne.0)then
      write(*,*) 'Error, inquiring pvd file'
      call ustop(' ')
   end if   
   if(.not.lex)then
      i = index(pvdfile,char(47),.true.)
      if(i.le.0) i = index(pvdfile,char(92),.true.)
      if(i.gt.0) call imod_utl_createdir(pvdfile(:i))
      iu = cfn_getlun(10,99)
      open(file=pvdfile,unit=iu)
      npvd = npvd + 1
      pvdun(npvd) = iu
      write(iu,'(a)') '<?xml version="1.0"?>'
      write(iu,'(a)') '<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">'
      write(iu,'(a)') '<Collection>'
   else
      inquire(file=pvdfile,number=iu,iostat=ios)
      if (ios.ne.0)then
         write(*,*) 'Error, inquiring pvd file'
         call ustop(' ')
      end if   
      if (iu.le.0)then ! file exists, but is not opened
         iu = cfn_getlun(10,99)
         open(file=pvdfile,unit=iu)
         npvd = npvd + 1
         pvdun(npvd) = iu
         write(iu,'(a)') '<?xml version="1.0"?>'
         write(iu,'(a)') '<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">'
         write(iu,'(a)') '<Collection>'
      end if
   end if
   write(s1,*) totim
   write(s2,'(5a)') '<DataSet timestep="',trim(adjustl(s1)),'" part="1" group="" file="',trim(vtufile),'"/>'
   write(iu,'(a)') trim(s2)
else
   return    
end if

end subroutine

subroutine vtkinit()

! modules
use vtuMod
use global, only: nrow, ncol, delr, delc, ibound, botm, lbotm   
use gwfmetmodule, only: coord_xll, coord_yll

implicit none

! locals
integer :: il,ir, ic, nrc, m, nlay
real :: dx, dy, dz, xm, ym, zm, t, b

! body
nlay = vtknlay

allocate(ggxnod(11,ncol*nrow*nlay+1))
ggxnod = 0.
dx = delr(1); dy = delc(1) ! UNIFORM ONLY
nrc = nrow*ncol
do il = 1, nlay
   do ir = 1, nrow
      do ic = 1, ncol
         m = ic+(ir-1)*ncol+(il-1)*nrc
         xm = ic*dx-dx/2; ym = (nrow-ir)*dy+dy/2
         t = botm(ic,ir,lbotm(il)-1)
         b = botm(ic,ir,lbotm(il))
         zm = (t+b)/2
         dz = t - b                    
         ! nod dum xm ym dz dx dy dz
         ggxnod(1,m) = real(m)
         ggxnod(5,m) = xm+coord_xll
         ggxnod(6,m) = ym+coord_yll
         ggxnod(7,m) = zm
         ggxnod(8,m) = dx
         ggxnod(9,m) = dy
         ggxnod(10,m) = dz
      end do    
   end do       
end do 

! unstructured grid data
nodes = nlay*nrow*ncol
neqs = nodes
ngeom = nodes
allocate(ia(neqs+1))
ia = 0   
call fillia
allocate(ja(nja))
call fillja

end subroutine

subroutine writevtu(text,fname,ngeom,ia,ja,ggxnod)

! modules
USE IR_Precision
USE Lib_VTK_IO
use vtuMod, only: mxlen

implicit none

! arguments
character(len=16), intent(in)         :: text
character(len=*), intent(in)          :: fname
integer, intent(in)                   :: ngeom
integer, dimension(*), intent(in)     :: ia, ja
real, dimension(11,ngeom), intent(in) :: ggxnod

! functions
integer :: getlun

! locals
logical :: lmatch
integer, parameter :: i0=1, i1=2, i2=3, i3=4, i4=5, i5=6, i6=7, i7=8 
!character(len=mxlen) :: fname
integer :: lun, n, i, j, k, ii, jj, kk, kkk, ic0, ic1, mxnbr, np
double precision, dimension(:,:,:), allocatable :: xyzl
double precision :: xm, ym, zm, dx, dy, dz
character(len=100), dimension(10) :: sarr

! variables for lib_vtk_io
integer(I4P) :: Nn ! number of points
integer(I4P) :: Ne ! number of cells
integer(I4P), dimension(:,:), allocatable :: connect
integer(I4P), dimension(:), allocatable :: offset
integer(I1P), dimension(:), allocatable :: cell_type
integer(I4P):: E_IO
real(R4P), dimension(:), allocatable :: x, y, z

! debug
character(len=mxlen) :: f
character(len=256) :: s
double precision :: xmin, xmax, ymin, ymax, zmin, zmax
integer :: m, nlast
logical :: lskip
double precision, dimension(:,:), allocatable :: dwrk 

! body

! write the
if (.false.) then
   allocate(dwrk(6,ngeom))
   do n = 1, ngeom
      xm  = dble(ggxnod(5,n))
      ym  = dble(ggxnod(6,n))
      zm  = dble(ggxnod(7,n))
      dx  = dble(ggxnod(8,n))
      dy  = dble(ggxnod(9,n))      
      dz  = dble(ggxnod(10,n))         
      dwrk(1,n) = xm-dx/2.
      dwrk(2,n) = xm+dx/2.
      dwrk(3,n) = ym-dy/2.
      dwrk(4,n) = ym+dy/2.
      dwrk(5,n) = zm-dz/2.
      dwrk(6,n) = zm+dz/2.
   end do 
    
   ! determine bounding box
   xmin = dble(minval(dwrk(1,:))) 
   xmax = dble(maxval(dwrk(2,:))) 
   ymin = dble(minval(dwrk(3,:))) 
   ymax = dble(maxval(dwrk(4,:))) 
   zmin = dble(minval(dwrk(5,:))) 
   zmax = dble(maxval(dwrk(6,:))) 
   write(*,*) xmin,xmax,ymin,ymax,zmin,zmax
   write(f,'(2a)') trim(fname), '.txt'    
   open(unit=99,file=f)
   write(99,'(a)') '['
   nlast = 0
   do n = 1, ngeom       
      lskip = .true.
      if (dwrk(1,n).eq.xmin .or. dwrk(2,n).eq.xmax) lskip = .false.
      if (dwrk(3,n).eq.ymin .or. dwrk(4,n).eq.ymax) lskip = .false.
      if (dwrk(5,n).le.zmin+2.0 .or. dwrk(6,n).ge.zmax-5.0) lskip = .false.
      if (lskip) cycle
      nlast = n
   end do 
   m = 0
   do n = 1, ngeom
      lskip = .true.
      if (dwrk(1,n).eq.xmin .or. dwrk(2,n).eq.xmax) lskip = .false.
      if (dwrk(3,n).eq.ymin .or. dwrk(4,n).eq.ymax) lskip = .false.
      if (dwrk(5,n).le.zmin+2.0 .or. dwrk(6,n).ge.zmax-5.0) lskip = .false.
      if (lskip) cycle     
      m = m + 1
      xm  = dble(ggxnod(5,n))
      ym  = dble(ggxnod(6,n))
      zm  = dble(ggxnod(7,n))
      dx  = dble(ggxnod(8,n))
      dy  = dble(ggxnod(9,n))      
      dz  = dble(ggxnod(10,n))        
      write(sarr(1),*) (xm-xmin)/(xmax-xmin)
      write(sarr(3),*) (ym-ymin)/(ymax-ymin)
      write(sarr(2),*) (zm-zmin)/(zmax-zmin)
      write(sarr(4),*) dx/(xmax-xmin)
      write(sarr(6),*) dy/(ymax-ymin)     
      write(sarr(5),*) dz/(zmax-zmin)
      write(sarr(7),*) m
      write(sarr(8),*) ggxnod(11,n)
      do i = 1, 8
         sarr(i) = adjustl(sarr(i))
      end do   
      write(s,'(17a)') '{"xc": [',trim(sarr(1)),',',trim(sarr(2)),',',trim(sarr(3)),'], "dx": [',trim(sarr(4)),',',trim(sarr(5)),',',trim(sarr(6)),'], "nod":',trim(sarr(7)),', "val":',trim(sarr(8)),'}'
      if(n.ne.nlast) s = trim(s)//','
      write(99,'(a)') trim(s)
   end do
   write(99,'(a)') ']'
   close(99)
   deallocate(dwrk)
end if

! count maximum number of neigbors
mxnbr = 0
do n = 1, ngeom
   ic0 = ia(n); ic1 = ia(n+1)-1   
   mxnbr = max(mxnbr,ic1-ic0)
end do

allocate(xyzl(3,8,mxnbr+1))
allocate(cell_type(ngeom),offset(ngeom),connect(8,ngeom)); connect = 0
allocate(x(8*ngeom),y(8*ngeom),z(8*ngeom))
np = 0
do n = 1, ngeom
   ic0 = ia(n); ic1 = ia(n+1)-1   
   jj = 0
   do ii=ic0,ic1 ! loop over parent and neighbors
      jj = jj + 1 
      j = ja(ii)
      ! determine corner nodes
      xm  = dble(ggxnod(5,j))
      ym  = dble(ggxnod(6,j))
      zm  = dble(ggxnod(7,j))
      dx  = dble(ggxnod(8,j))
      dy  = dble(ggxnod(9,j))      
      dz  = dble(ggxnod(10,j)) 
      ! determine corner points
    ! x                        y                        z 
      xyzl(1,i0,jj) = xm - dx/2; xyzl(2,i0,jj) = ym - dy/2; xyzl(3,i0,jj) = zm - dz/2  
      xyzl(1,i1,jj) = xm + dx/2; xyzl(2,i1,jj) = ym - dy/2; xyzl(3,i1,jj) = zm - dz/2  
      xyzl(1,i2,jj) = xm - dx/2; xyzl(2,i2,jj) = ym + dy/2; xyzl(3,i2,jj) = zm - dz/2  
      xyzl(1,i3,jj) = xm + dx/2; xyzl(2,i3,jj) = ym + dy/2; xyzl(3,i3,jj) = zm - dz/2  
      
      xyzl(1,i4,jj) = xm - dx/2; xyzl(2,i4,jj) = ym - dy/2; xyzl(3,i4,jj) = zm + dz/2  
      xyzl(1,i5,jj) = xm + dx/2; xyzl(2,i5,jj) = ym - dy/2; xyzl(3,i5,jj) = zm + dz/2  
      xyzl(1,i6,jj) = xm - dx/2; xyzl(2,i6,jj) = ym + dy/2; xyzl(3,i6,jj) = zm + dz/2  
      xyzl(1,i7,jj) = xm + dx/2; xyzl(2,i7,jj) = ym + dy/2; xyzl(3,i7,jj) = zm + dz/2  
   end do
   
   ! add node numbers
   j = ja(ic0) ! parent node
   do kk=i0,i7 ! coordinates of parent node
      if (connect(kk,j) == 0) then
         np = np + 1 ! number of points
         connect(kk,j) = np
         x(np) = xyzl(1,kk,1)
         y(np) = xyzl(2,kk,1)
         z(np) = xyzl(3,kk,1)
      end if   
   end do    
   
   ! check for double nodes
   do kk = i0,i7 ! loop over parent nodes
      jj = 1
      do ii = ic0+1,ic1 ! loop over neighbor nodes
         jj = jj + 1 
         j = ja(ii)
         do kkk = i0,i7
            if (xyzl(1,kkk,jj)==xyzl(1,kk,1).and.&
                xyzl(2,kkk,jj)==xyzl(2,kk,1).and.&
                xyzl(3,kkk,jj)==xyzl(3,kk,1)) then
               connect(kkk,j) = connect(kk,ja(ic0))    
            end if   
         end do
      end do 
   end do    
end do

! write the vtu-file
cell_type = 11; Nn = np; Ne = ngeom
offset(1) = 8
do i = 2, Ne
   offset(i) = offset(i-1)+8
end do   
do i = 1, Ne
   do kk = i0,i7
      connect(kk,i) = connect(kk,i)-1
   end do
end do   
!write(fname,'(2a)') trim(vtkPref), '.vtu'
!E_IO = VTK_INI_XML(output_format = 'ascii', filename = fname, mesh_topology = 'UnstructuredGrid')
E_IO = VTK_INI_XML(output_format = 'binary', filename = fname, mesh_topology = 'UnstructuredGrid')
E_IO = VTK_GEO_XML(NN = Nn, NC = Ne, X = x, Y = y, Z = z)
E_IO = VTK_CON_XML(NC = Ne, connect = reshape(connect,(/8*Ne/)), offset = offset, cell_type = cell_type )
E_IO = VTK_DAT_XML(var_location = 'cell', var_block_action = 'opeN')
E_IO = VTK_VAR_XML(NC_NN = Ne, varname = trim(adjustl(text)), var = ggxnod(11,:))
E_IO = VTK_DAT_XML(var_location = 'cell', var_block_action = 'CLOSE')
E_IO = VTK_GEO_XML()
E_IO = VTK_END_XML()

! clean up
deallocate(xyzl,cell_type,offset,connect,x,y,z)

end subroutine writevtu