module outmodule
contains
   subroutine output(q,nx,ny)
      real*8, intent(in) :: q(:)
      integer, intent(in) :: nx,ny
      integer l,j
      l = 1
      do j = 1,ny+1
         write(2,'(1p," ETA=",i5,5e26.17,/,("          ",5e26.17))') j,q(l:l+nx)
         l = l+nx+1
      enddo

   end subroutine output

end module outmodule

program xyz2grd
   use outmodule
   implicit none
   integer nx,ny,i,j,l
   character(100):: infile,grdfile,depfile
   real*8,allocatable:: x(:),y(:),z(:)
   infile = 'island.xyz'
   grdfile = 'island.grd'
   depfile = 'island.dep'

   nx = 50
   ny = 125
   allocate(x((nx+1)*(ny+1)))
   allocate(y((nx+1)*(ny+1)))
   allocate(z((nx+1)*(ny+1)))

   open(1,file=infile)

   do i=1,(nx+1)*(ny+1)
      read (1,*) x(i),y(i),z(i)
   enddo

   close(1)

   open(2,file=grdfile)
   open(3,file=depfile)

   write(2,'("*")')
   write(2,'("* Deltares, RGFGRID Version 4.20.00.34496, Jun 17 2014, 15:15:54")')
   write(2,'("* File creation date: 2015-11-19, 02:01:16")')
   write(2,'("*")')
   write(2,'("Coordinate System = Cartesian")')
   write(2,'("Missing Value     =    0.00000000000000000E+00")')
   write(2,'(2i8)') nx+1,ny+1
   write(2,'(3i2)') 0,0,0

   call output(x,nx,ny)
   call output(y,nx,ny)

   close(2)

   l=1
   do j = 1,ny+1
      write(3,'(1p,12e16.7)') z(l:l+nx),-999.0
      l = l+nx+1
   enddo

   close(3)

end program xyz2grd