c**********************************************************************
c     subroutine sbgrid
c**********************************************************************
      subroutine sbgrid
c
c     +++ PURPOSE +++
c     to calculate grid size and spacing for EROSION.
c     grid size assumes outer points are outside field boundary
c      a distance ix/2
c     to calculate number of grid points for EROSION.
c     A max 'interior' square grid of 29X29 is assigned-no barriers
c     A max 'interior' rectangular grid of 59X59 is assinged barriers
c     to assign subregion index no. to each grid point.
c
c     +++ ARGUMENT DECLARATION +++

c

c     +++ LOCAL DEFINITIONS +++
c     imax  - no. grid intervals in x-direction
c     jmax  - no. grid intervals in y-direction.
c     ix    - grid interval in x-direction (m)
c     jy    - grid interval in y-direction (m)
c     dxmin - minimum grid interval (m)
c     csr   - current subr. index at grid point i,j
c     icsr  - same as csr but not an array
c     i,j   - do loop indexes
c
c     + + + GLOBAL COMMON BLOCKS + + +
*$noreference
      include  'p1werm.inc'
      include  'm1geo.inc'
      include  'm1subr.inc'
c
c     + + + LOCAL COMMON BLOCKS + + +
      include  'erosion/m2geo.inc'
      include  'erosion/e2grid.inc'
c
*$reference
c
c     +++ LOCAL VARIABLES +++
      integer  ngdpt
      integer  icsr, i, j
      real*4   dxmin, lx, ly
c
c     +++ END SPECIFICATIONS +++
c
c     set min grid spacing
       dxmin = 7.0
c     set max no. of grid points with no barrier

       ngdpt = 30
c     barriers?
       if (nbr .gt. 0) then
c        find shortest barrier to determine dxmin
         do 5 i=1,nbr
            dxmin = min(dxmin, 5.0*amzbr(i))
    5    continue
         ngdpt = mngdpt
       endif

c     calculate max grid intervals
c       calc. lx and ly sides of field
      lx = amxsim (1,2)-amxsim(1,1)
      ly = amxsim (2,2)-amxsim(2,1)
c
c^tmp out
c      write(*,*) 'tmp out from sbgrid, line 69'
c     write (*,*) 'lx=', lx, 'ly=',ly

c     increase grid points on large field
c      if((lx .gt. 200) .or. (ly  .gt. 200)) then
c         ngdpt = mngdpt
c      endif
c
c        case where lx > ly
      if ( lx .gt. ly)then
        imax  = ifix ( lx / dxmin)
        imax = min(imax,ngdpt)
        imax = max(imax,2)
c     calculate spacing for square or with barriers a rectangular grid
        ix  = lx / (imax - 1)

         if (nbr .gt. 0) then
           jmax  = ifix (ly / dxmin)
           jmax  = min(jmax, ngdpt)
         else
           jmax = anint(ly/ix) + 1
         endif

        jmax = max(jmax,2)
        jy   = ly/(jmax - 1)

c        case where lx = ly or lx < ly
      else
        jmax  = ifix (ly / dxmin)
        jmax = min(jmax,ngdpt)
        jmax = max(jmax,2)
        jy   = ly / (jmax - 1)

        if (nbr .gt. 0) then
           imax  = ifix (lx / dxmin)
           imax  = min(imax,ngdpt)
        else
           imax = anint(lx/jy) + 1
        endif
        imax = max(imax,2)
        ix = lx/(imax-1)

      endif
c
c     determine subregion of each grid point
c     for a single subregion now
      icsr = 1
      do 20 j = 0, jmax
      do 10 i = 0, imax
c
c     for multiple subregions
c     do 5  icsr = 1, nsubr
c     for multiple subregions
c     if (i*ix .lt. amxsr(1,1,icsr) .or. i*ix .gt. amxsr(1,2,icsr))
c    & go to 10
c     if (j*jy .lt. amxsr(2,1,icsr) .or. j*jy .gt. amxsr(2,2,icsr))
c    & go to 10
c   5 continue
c
        csr(i,j) = icsr
   10 continue
   20 continue
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++