      subroutine BathRead(d,q,h1,h2,n1,n2,m,n,grnd,indir,latimg)
c.....Reads bathymetry data in MOST or ASCII Grid format
c.....initializes variables and BC
      real*8 d(m,n),q(m,n),grnd
      real*8 h1(m),h2(n)
      real*4 val
      character*120 fname,fnameout,indir,bathdimnames(4)
      character*200 str, str2
      real*8 xllcorner, yllcorner, cellsize, nodata
      real*8 xoffsetcell, yoffsetcell,dd(m)
      integer headerlines
      integer*4  iret, ncid, latid, lonid, bathid
      integer*4 lattype, lontype, bathtype, cnt(2), strt(2)
      integer latdims(4), londims(4), bathdims(4), lonlen, latlen
      logical lon_neg
      integer*4 count
c..... switch when finding image-ordered bathymetry files
      logical latimg

      common /nc/ ns1,ne1,ns2,ne2,istep,jstep,scl_m2cm

      include 'netcdf.inc'

      read(2,'(a)') fname
      write(3,*) ' Bathymetry: ', fname(1:lnblnk(fname))
      iret = nf_open(fname,0,ncid)
c.....if not successfully opened as NetCDF, go to ASCII-format readers
      if (iret .ne. NF_NOERR) then
         goto 805
      endif
      write (3,*) 'Detected netCDF formatted bathymetry file'      
      if ( nf_inq_varid(ncid,'lon', lonid) .ne. NF_NOERR) then
        if ( nf_inq_varid(ncid,'LON', lonid) .ne. NF_NOERR) then
          call check_err(nf_inq_varid(ncid,'x1', lonid))
        endif
      endif
      if ( nf_inq_varid(ncid,'lat', latid) .ne. NF_NOERR) then
        if ( nf_inq_varid(ncid,'LAT', latid) .ne. NF_NOERR) then
          call check_err(nf_inq_varid(ncid,'x2', latid))
        endif
      endif     
      if ( nf_inq_varid(ncid,'bathy', bathid) .ne. NF_NOERR) then
        if ( nf_inq_varid(ncid,'BATHY', bathid) .ne. NF_NOERR ) then
          if ( nf_inq_varid(ncid,'depth', bathid) .ne. NF_NOERR) then
            call check_err(nf_inq_varid(ncid,'DEPTH', bathid))
          endif
        endif
      endif
      call check_err(nf_inq_vardimid(ncid, bathid, bathdims))
      call check_err(nf_inq_dimname(ncid, bathdims(1), bathdimnames(1)))
      call check_err(nf_inq_dimname(ncid, bathdims(2), bathdimnames(2)))
      call check_err(nf_inq_dimlen(ncid, lonid, lonlen))
      call check_err(nf_inq_dimlen(ncid, latid, latlen))
      n1 = lonlen
      n2 = latlen
      call bath_size_check(n1,n2,m,n)
      call check_err(nf_inq_vartype(ncid, latid, lattype))
      call check_err(nf_inq_vartype(ncid, lonid, lontype))
      call check_err(nf_inq_vartype(ncid, bathid, bathtype))
      if (lattype .eq. NF_FLOAT) then
          iret = nf_get_var_real(ncid, latid, h2)! latimg)
      else if (lattype .eq. NF_DOUBLE) then
          iret = nf_get_var_double(ncid, latid, h2)! latimg)
      else 
          goto 820
      endif

      call check_lat_image(h2,n2,n,latimg)

      if (lontype .eq. NF_FLOAT) then
          iret = nf_get_var_real(ncid, lonid, h1)
      else if (lontype .eq. NF_DOUBLE) then
          iret = nf_get_var_double(ncid, lonid, h1)
      else 
          goto 820
      endif
      
       if (bathtype .eq. NF_FLOAT) then
          ibtype=1
      else if (bathtype .eq. NF_DOUBLE) then
          ibtype=2
      else 
          goto 820
      endif
      strt(1)=1
      cnt(1)=n1
      cnt(2)=1
      do j=1,n2
        strt(2)=j
        if (ibtype .eq. 1) then
          iret = nf_get_vara_real(ncid,bathid,strt,cnt,dd)
        else
      	  iret = nf_get_vara_double(ncid,bathid,strt,cnt,dd)
        endif
        if (latimg) then
          do i=1,n1
            d(i,n2-j+1)=dd(i)
          enddo
        else
          do i=1,n1
            d(i,j)=dd(i)
          enddo
        endif
      enddo
	            
      
      iret = nf_close(ncid)
      goto 840
c.....ASCII format handling; examine first line to differentiate MOST
c.....format from ESRI ASCII grid format
 805  open(unit=1,file=fname,status='old',form='formatted')
      read (1,'(a)', err=820, end=820) str
      call to_lower(str, str)
      if (str(1:5).eq.'ncols') goto 810
c.....Read MOST format
      read (str,*, err=820) n1,n2
      call bath_size_check(n1,n2,m,n)
      write (3,*) 'Reading MOST formatted bathymetry file'

      read (1,*) (h1(i), i=1,n1)
      read (1,*) (h2(i), i=1,n2)
      call check_lat_image(h2,n2,n,latimg)
      if (latimg) then
        do j=n2,1,-1
          read (1,*) (d(i,j), i=1,n1)
        end do
      else
        do j=1,n2
          read (1,*) (d(i,j), i=1,n1)
        end do
      endif
      goto 840
c.....Read ASCII Grid format
 810  read (str,*, err=820, end=820) str2,n1
      write (3,*) 'Reading ASCII Grid formatted bathymetry file'
      latimg = .true.
      read (1,*, err=820, end=820) str2,n2
      call to_lower(str2, str2)
      if (str2.ne.'nrows') then
         write(3,*) 'Error reading ASCII Grid file: expecting nrows'
         goto 820
      endif
      call bath_size_check(n1,n2,m,n)
      read (1,'(a)', err=820, end=820) str
      call to_lower(str, str)
      if (str(1:9).eq.'xllcorner') then
         xoffsetcell=0.5
      else if (str(1:9).eq.'xllcenter') then
         xoffsetcell=0
      else
         write(3,*) 'Error reading ASCII Grid file:'
         write(3,*) 'expecting xllcorner or xllcenter'
         goto 820
      endif
      read (str,*, err=820, end=820) str2,xllcorner
      if (xllcorner.gt.360.or.xllcorner.lt.-360) then
         write(3,*) 'xllcorner > 360 or xllcorner < -360'
         write(3,*) 'Grid file must be geodetic coord system (lat/lon)'
         goto 820
      endif
      read (1,'(a)', err=820, end=820) str
      call to_lower(str, str)
      if (str(1:9).eq.'yllcorner') then
         yoffsetcell=0.5
      else if (str(1:9).eq.'yllcenter') then
         yoffsetcell=0
      else
         write(3,*) 'Error reading ASCII Grid file: '
         write(3,*) 'expecting yllcorner or yllcenter'
         goto 820
      endif
      read (str,*, err=820, end=820) str2,yllcorner
      read (1,'(a)', err=820, end=820) str
      call to_lower(str, str)
      if (str(1:8).ne.'cellsize') then
         write(3,*) 'Error reading ASCII Grid file: expecting cellsize'
         goto 820
      endif
      read (str,*, err=820, end=820) str2,cellsize
      read (1,'(a)', err=820, end=820) str
      call to_lower(str, str)
c     nodata_value line is optional
      if (str(1:6).eq.'nodata') then
         headerlines=6
      else
         headerlines=5
      endif
c      write(3,*) 'ncols: ', n1
c      write(3,*) 'nrows: ', n2
c      write(3,*) 'xll: ', xllcorner
c      write(3,*) 'yll: ', yllcorner
c      write(3,*) 'cellsize: ', cellsize
c      write(3,*) 'nodata: ', nodata
      do i=1,n1
         h1(i)=xllcorner + xoffsetcell*cellsize + i*cellsize
      end do
      do i=n2,1,-1
         h2(i)=yllcorner + yoffsetcell*cellsize + i*cellsize
      end do
c     Close and reopen file, read all header lines now that we know how many
      close(unit=1,status='keep')
      open(unit=1,file=fname,status='old',form='formatted')
      do i=1,headerlines
         read (1,'(a)', err=820, end=820) str
      end do
c     Read the gridded data
      do j=n2,1,-1
         read (1,*) (d(i,j), i=1,n1)
      end do
c     ASCII grids usually have bathymetry negative - invert
      do i=1,n1
         do j=1,n2
            d(i,j) = -d(i,j)
         end do
      end do
      goto 840
 820  write (3,'(a,$)') 'ERROR: reading bathymetry file: '
      write (3,*) fname
      write (3,*) 'Halting execution.'
      call flush(3)
      close(unit=1,status='keep')
      stop
 840  continue
      call flush(3)
      close(unit=1,status='keep')
      lon_neg = .false.
      do i=1,n1
      	if (h1(i).LT.0.) lon_neg = .true.
      end do
      if (lon_neg) then
        do i=1,n1
          h1(i)=360.+h1(i)
        end do
      end if
c     Save original d in q to see how many points are affected
      do i=1,n1
         do j=1,n2
            q(i,j)=d(i,j)
         end do
      end do
      write(3,*) 'Eliminating one-node islands...'
      call RemoveIslands(d,q,n1,n2,m,n,grnd)
      count=0
      do i=1,n1
         do j=1,n2
           if (q(i,j).ne.d(i,j)) then
             count = count + 1
             q(i,j)=d(i,j)
           end if
         end do
      end do
      write(3,845) count
 845  format('  Number of points affected: ',i12)
      if (count.GT.0) then
        fnameout = fname(1:lnblnk(fname))//".1nod"
        write(3,*) '  Writing 1-node eliminated bathy file: ',fnameout
        call bath_write(fnameout,n1,n2,m,n,h1,h2,d)
      end if
      return
      end
      
      subroutine to_lower(instr,outstr)
        character *(*) instr, outstr
        integer i,ic,nlen
        nlen = len(instr)
        outstr = instr(1:nlen)
        do i=1,nlen
          ic = ichar(instr(i:i))
          if (ic .ge. 65 .and. ic .le. 90) outstr(i:i) = char(ic+32)
        end do
      end
      
      subroutine bath_size_check(n1,n2,m,n)
       implicit none
       integer n1,n2,m,n
       if (n1.GT.m) then
        write(3,*) 'ERROR: bathymetry grid x-dimension too large:',n1
        write(3,*) 'x-dimension limit:',m
        stop
       end if
       if (n2.GT.n) then
        write(3,*) 'ERROR: bathymetry grid y-dimension too large:',n2
        write(3,*) 'x-dimension limit:',n
        stop
       end if
      end

      subroutine check_lat_image(h2,n2,n,latimg)
        implicit none
        integer i, n2, n
        real*8 h2(n),lattmp(n)
        logical latimg

        if (h2(2).LT.h2(1)) then
          write(3,*) '...image-format bathy file (latitude decreasing)'
          latimg = .true.
          do i=1, n2
            lattmp(i) = h2(n2-i+1)
          enddo
          do i=1, n2
            h2(i) = lattmp(i)
          enddo
        else
          write(3,*) '...data-format bathy file (latitude increasing)'
          latimg = .false.
        endif
      end

      subroutine bath_write(fname,n1,n2,m,n,ax1,ax2,bath)
      implicit none
      character*(*) fname
      integer i,j,n1,n2,m,n
      real*8 ax1(m),ax2(n)
      real*8 bath(m,n)
      open (unit=12,file=fname,status='unknown',form='formatted')
      write(12,850) n1,n2
      do i=1,n1
        write(12,860) ax1(i)
      enddo
      do j=n2,1,-1
        write(12,860) ax2(j)
      enddo
      do j=n2,1,-1
        write(12,870) (bath(i,j),i=1,n1)
      enddo
 850  format(1X,I5,1X,I5)
 860  format(F13.8)  ! horizontal accuracy: 1.1cm (at equator)
 870  format(6000F12.5) ! vertical accuracy: 0.1cm

      end
