************************************************************************
*      program gridmod.f
*
*      purpose:  Modifies sducti.x.bin(grid used in NPARC example case)
*                so that block boundaries do not overlap.
*               (Gets rid of k = 1 surface in block 1.)
*
*      input files:
*
*        fort.8     sduct.i.x.bin (PLOT3D/mgrid/unformatted  binary)
*
*
*       output files:
*
*        fort.9     sduct.new.x.bin (PLOT3D/mgrid/unformatted binary)
*
*       input variables:
*
*
************************************************************************
*
      implicit none
      integer i, igrid, j, k, k1, ngrid
      integer idim(2), jdim(2), kdim(2)
      real x(123,31,22), y(123,31,22), z(123,31,22)
*
*-----Read original grid sducti.x.bin
*
      read(8) ngrid
      write(9) ngrid
      read(8) (idim(igrid), jdim(igrid), kdim(igrid), igrid = 1, ngrid)
      kdim(1) = kdim(1) - 1
      write(6,*) "kdim1 = ", kdim(1)
      write(9) (idim(igrid), jdim(igrid), kdim(igrid), igrid = 1, ngrid)
      kdim(1) = kdim(1) + 1
      write(6,*) "kdim1 = ", kdim(1)
      do 10 igrid = 1, ngrid
        read(8) 
     c  (((x(i,j,k),i=1,idim(igrid)),j=1,jdim(igrid)),k=1,kdim(igrid)),
     c  (((y(i,j,k),i=1,idim(igrid)),j=1,jdim(igrid)),k=1,kdim(igrid)),
     c  (((z(i,j,k),i=1,idim(igrid)),j=1,jdim(igrid)),k=1,kdim(igrid))
        if (igrid .eq. 1) then
           k1 = 2
        else
           k1 = 1
        endif
        write(9) 
     c  (((x(i,j,k),i=1,idim(igrid)),j=1,jdim(igrid)),k=k1,kdim(igrid)),
     c  (((y(i,j,k),i=1,idim(igrid)),j=1,jdim(igrid)),k=k1,kdim(igrid)),
     c  (((z(i,j,k),i=1,idim(igrid)),j=1,jdim(igrid)),k=k1,kdim(igrid))
10    continue
      stop
      end