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     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  icsr, i, j
      real*4   dxmin
c
c     +++ END SPECIFICATIONS +++
c
c     set min grid spacing
       dxmin = 7.0
c     calculate max grid intervals
c        case where x > y
      if ((amxsim (1,2)-amxsim(1,1)).gt.(amxsim (2,2)-amxsim(2,1)))then
        imax  = ifix ((amxsim (1,2)-amxsim(1,1)) / dxmin)
        imax = min(imax,mngdpt)
        imax = max(imax,2)
c     calculate spacing for uniform grid.
        ix  = (amxsim (1,2)-amxsim(1,1)) / (imax - 1)
        jmax = anint((amxsim(2,2)-amxsim(2,1))/ix) + 1
        jmax = max(jmax,2)
        jy   = (amxsim(2,2)-amxsim(2,1))/(jmax - 1)
c        case where x = y or x < y
      else
        jmax  = ifix ((amxsim (2,2)-amxsim(2,1)) / dxmin)
        jmax = min(jmax,mngdpt)
        jmax = max(jmax,2)
        jy   = (amxsim (2,2)-amxsim(2,1)) / (jmax - 1)
        imax = anint((amxsim(1,2)-amxsim(1,1))/jy) + 1
        imax = max(imax,2)
        ix = (amxsim(1,2)-amxsim(1,1))/(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++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++