      subroutine bigrid(depth)
c
c --- set loop bounds for irregular basin in c-grid configuration
c --- q,u,v,p are vorticity, u-velocity, v-velocity, and mass points, resp.
c --- 'depth' = basin depth array, zero values indicate land
c
      include 'dimension_raw.h'
c
c --- mpd = maximum number of points on a diagonal
      parameter (mpd=jdm)        !  mpd=min(idm,jdm) where allowed by compiler
      dimension depth(idm,jdm)
      dimension id(mpd,idm+jdm),jd(mpd,idm+jdm),npts(idm+jdm)
c
      do 17 i=1,ii
      do 17 j=1,jj
      ip(i,j)=0
      iq(i,j)=0
      iu(i,j)=0
 17   iv(i,j)=0
c
c --- fill single-width inlets
 16   nfill=0
      do 15 j=1,jj1
      ja=2*max(  1,j-1)-(j-1)
      jb=2*min(jj1,j+1)-(j+1)
      do 15 i=1,ii1
      ia=2*max(  1,i-1)-(i-1)
      ib=2*min(ii1,i+1)-(i+1)
      nzero=0
      if (depth(i,j).gt.0.) then
      if (depth(ia,j).le.0.) nzero=nzero+1
      if (depth(ib,j).le.0.) nzero=nzero+1
      if (depth(i,ja).le.0.) nzero=nzero+1
      if (depth(i,jb).le.0.) nzero=nzero+1
      if (nzero.ge.3) then
      write (*,'('' depth(''i4'',''i4'') has been set to zero'')') i,j
      depth(i,j)=0.
      nfill=nfill+1
      end if
      end if
 15   continue
      if (nfill.gt.0) go to 16
c
c --- mass points are defined where water depth is greater than zero
      do 2 i=1,ii1
      do 2 j=1,jj1
      if (depth(i,j).gt.0.) ip(i,j)=1
  2   continue
c
c --- u,v points are located halfway between any 2 adjoining mass points
      do 3 j=1,jj1
      do 3 i=2,ii1
      if (ip(i-1,j).gt.0.and.ip(i,j).gt.0) iu(i,j)=1
  3   continue
      do 4 j=2,jj1
      do 4 i=1,ii1
      if (ip(i,j-1).gt.0.and.ip(i,j).gt.0) iv(i,j)=1
  4   continue
c
c --- 'interior' q points require water on all 4 sides.
      do 5 j=2,jj1
      do 5 i=2,ii1
      if (min(ip(i,j),ip(i-1,j),ip(i,j-1),ip(i-1,j-1)).gt.0) iq(i,j)=1
  5   continue
c
c --- 'promontory' q points require water on 3 (or at least 2 diametrically 
c --- opposed) sides
      do 10 j=2,jj1
      do 10 i=2,ii1
ccc      if (iq(i,j).ne.0) go to 10
ccc      if (ip(i,j).gt.0.and.ip(i-1,j-1).gt.0) iq(i,j)=6
ccc      if (ip(i-1,j).gt.0.and.ip(i,j-1).gt.0) iq(i,j)=7
      if ((ip(i  ,j).gt.0.and.ip(i-1,j-1).gt.0).or.
     .    (ip(i-1,j).gt.0.and.ip(i  ,j-1).gt.0)) iq(i,j)=1
 10   continue
c
c --- determine loop bounds for vorticity points, including interior and
c --- promontory points
      call indxi(iq,ifq,ilq,isq)
      call indxj(iq,jfq,jlq,jsq)
c
c --- vorticity points on borders
c --- 'left' border (j increasing to right)
ccc      do 6 j=1,jj1
ccc      do 6 i=2,ii1
ccc      if (iq(i,j).ne.0) go to 6
ccc      if (ip(i,j).gt.0.and.ip(i-1,j).gt.0) iq(i,j)=2
ccc  6   continue
c --- 'right' border
ccc      do 7 j=2,jj
ccc      do 7 i=2,ii1
ccc      if (iq(i,j).ne.0) go to 7
ccc      if (ip(i,j-1).gt.0.and.ip(i-1,j-1).gt.0) iq(i,j)=3
ccc  7   continue
c --- 'top' border (i increasing downward)
ccc      do 8 j=2,jj1
ccc      do 8 i=1,ii1
ccc      if (iq(i,j).ne.0) go to 8
ccc      if (ip(i,j).gt.0.and.ip(i,j-1).gt.0) iq(i,j)=4
ccc  8   continue
c --- 'bottom' border
ccc      do 9 j=2,jj1
ccc      do 9 i=2,ii
ccc      if (iq(i,j).ne.0) go to 9
ccc      if (ip(i-1,j).gt.0.and.ip(i-1,j-1).gt.0) iq(i,j)=5
ccc  9   continue
c
c --- determine loop indices for mass and velocity points
      call indxi(ip,ifp,ilp,isp)
      call indxj(ip,jfp,jlp,jsp)
      call indxi(iu,ifu,ilu,isu)
      call indxj(iu,jfu,jlu,jsu)
      call indxi(iv,ifv,ilv,isv)
      call indxj(iv,jfv,jlv,jsv)
c
c --- velocity points on borders
c --- 'left' and 'right' border points for v (y pointing right)
ccc      do 13 i=1,ii1
ccc      do 13 j=1,jj1
ccc      if (ip(i,j).ne.1) go to 13
ccc      if (iv(i,j).eq.0) iv(i,j)=2
ccc      if (iv(i,j+1).eq.0) iv(i,j+1)=2
ccc 13   continue
c --- 'top' and 'bottom' border points for u (x pointing down)
ccc      do 14 i=1,ii1
ccc      do 14 j=1,jj1
ccc      if (ip(i,j).ne.1) go to 14
ccc      if (iu(i,j).eq.0) iu(i,j)=2
ccc      if (iu(i+1,j).eq.0) iu(i+1,j)=2
ccc 14   continue
      return
      end
