      subroutine bdryindx(fkmp,fkmq)
c
c=======================================================================
c                                                                    ===
c  This routine sets boundary indicators auxiliary arrays.           ===
c                                                                    ===
c  On output:                                                        ===
c                                                                    ===
c     FKMP    number of vertical levels of ocean at T points.        ===
c     FKMQ    number of vertical levels of ocean at U,V points.      ===
c                                                                    ===
c  Calls:     EXITUS                                                 ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <fullwd.h>
#include <onedim.h>
#include <fields.h>
#include <iounits.h>
#include <vertslabs.h>
#if defined coast & defined coastedge
#  include <edges.h>
#endif
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,j,l,lsegp,n
      FLOAT
     *      fkmp(imt,jmt),fkmq(imt,jmt),fkmz(imt,jmt)
#if defined coast & defined coastedge
      FLOAT
     *     fx
      integer m
      integer iadd(8),iedge,ipt,jadd(8),jedge,jpt
#endif
      equivalence (pb,fkmz)
#if defined pindex & defined coast
      integer on
#endif

#if defined coast & defined coastedge
      save iadd,jadd
      data iadd /1,-1,0,0,1,-1,1,-1/
      data jadd /0,0,1,-1,1,1,-1,-1/
#endif
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c=======================================================================
c  Set maximum level indicators for topography.
c=======================================================================
c
      do 100 j=1,jmt
      do 100 i=1,imt
        fkmp(i,j)=0
        fkmq(i,j)=0
        fkmz(i,j)=0
 100  continue
#ifdef coast
c
c  1st, set number of vertical levels for T points.  At input the
c  land/sea descriptor array LANDT has value 1 at sea points.
c
      do 110 j=1,jmt
      do 110 i=1,imt
        if(landt(i,j).eq.1) fkmp(i,j)=km
 110  continue
#else
c
c  1st, set number of vertical levels for T points.
c
      do 110 j=1,jmt
      do 110 i=1,imt
        fkmp(i,j)=km
 110  continue
#endif
c
c  2nd, compute number of vertical levels at each U,V point.
c
      do 120 j=1,jmtm1
      do 120 i=1,imtm1
        fkmq(i,j)=min(fkmp(i,j),fkmp(i+1,j),fkmp(i,j+1),fkmp(i+1,j+1))
 120  continue
      do 130 j=1,jmt
        fkmq(imt,j)=fkmq(imtm1,j)
 130  continue
      do 140 i=1,imt
        fkmq(i,jmt)=fkmq(i,jmtm1)
 140  continue
c
c  3rd, compute an array to indicate "interior" grid boxes.
c
#if defined coast & defined coastedge
      nedget=0
      do 150 j=2,jmt
      do 150 i=2,imt
        fkmz(i,j)=min(fkmq(i-1,j-1),fkmq(i,j-1),fkmq(i-1,j),fkmq(i,j))
        if(fkmp(i,j).gt.c0 .and. fkmz(i,j).eq.c0) then
           nedget=nedget+1
           ipdxt(nedget)=i+(j-1)*imt
        endif
 150  continue
      do 160 j=1,jmt
        fkmz(1,j)=fkmz(2,j)
        if(fkmp(1,j).gt.c0 .and. fkmz(1,j).eq.c0) then
           nedget=nedget+1
           ipdxt(nedget)=1+(j-1)*imt
        endif
 160  continue
      do 170 i=1,imt
        fkmz(i,1)=fkmz(i,2)
        if(fkmp(i,1).gt.c0 .and. fkmz(i,1).eq.c0) then
           nedget=nedget+1
           ipdxt(nedget)=i
        endif
 170  continue
c
      do 175 i=1,nedget
         do 171 l=1,8
            boxt(l,i)=0
 171     continue
         l=1
         call ip2ij (ipdxt(i),imt,ipt,jpt)
         do 172 m=1,8
            iedge=ipt+iadd(m)
            jedge=jpt+jadd(m)
            if ((iedge.ge.1).and.(iedge.le.imt).and.
     &                           (jedge.ge.1).and.(jedge.le.jmtm1)) then
               j=iedge+(jedge-1)*imt
               call upbox(l,j,ipdxt,nedget,fkmz,imtjmt,boxt(l,i))
            end if
 172     continue
 175  continue
c
c v-points
c
      nedgev=0
      do 176 j=1,jmtm1
      do 176 i=1,imtm1
        fx=min(fkmz(i,j),fkmz(i,j+1),fkmz(i+1,j),fkmz(i+1,j+1))
        if(fkmq(i,j).gt.c0 .and. fx.eq.c0) then
           nedgev=nedgev+1
           ipdxv(nedgev)=i+(j-1)*imt
        endif
 176  continue
      do 179 i=1,nedgev
         do 177 l=1,8
            boxv(l,i)=0
 177     continue
         l=1
         call ip2ij (ipdxv(i),imt,ipt,jpt)
         do 178 m=1,8
            iedge=ipt+iadd(m)
            jedge=jpt+jadd(m)
            if ((iedge.ge.1).and.(iedge.le.imt).and.
     &                             (jedge.ge.1).and.(jedge.le.jmt)) then
               j=iedge+(jedge-1)*imt
               call upbox(l,j,ipdxv,nedgev,fkmq,imtjmt,boxv(l,i))
            end if
 178     continue
 179  continue
#else
      do 150 j=2,jmt
      do 150 i=2,imt
        fkmz(i,j)=min(fkmq(i-1,j-1),fkmq(i,j-1),fkmq(i-1,j),fkmq(i,j))
 150  continue
      do 160 j=1,jmt
        fkmz(1,j)=fkmz(2,j)
 160  continue
      do 170 i=1,imt
        fkmz(i,1)=fkmz(i,2)
 170  continue
#endif
c
c-----------------------------------------------------------------------
c  Compute start and end indices at tracer points for stream function
c  calculations and filters.
c-----------------------------------------------------------------------
c
#if defined pindex & defined coast
c
      do 1171 j=1,jmt
        do 1171 l=1,lseg
          isp(j,l)=0
          iep(j,l)=0
1171  continue
c
      do 1172 j=2,jmt-1
        l=0
        on=0
        do 1173 i=2,imt-1
c
          if(on.eq.0.and.landt(i,j).gt.0) then
            on=1
            l=l+1
            if(l.gt.lseg) then
              write(stdout,900) 'LSEG = ',lseg
              call exitus('BDRYINDX')
            endif
            isp(j,l)=i
          endif
c
          if(on.eq.1.and.landt(i,j).eq.0) then
            on=0
            iep(j,l)=i-1
          endif
c
          if(on.eq.1.and.i.eq.imt-1) then
            iep(j,l)=i
          endif
c
1173    continue
1172  continue
c
c
      do 1176 i=1,imt
        do 1176 l=1,lseg
          jsp(i,l)=0
          jep(i,l)=0
1176   continue
c
      do 1174 i=2,imt-1
        l=0
        on=0
        do 1175 j=2,jmt-1
c
          if(on.eq.0.and.landt(i,j).gt.0) then
            on=1
            l=l+1
            if(l.gt.lseg) then
              write(stdout,900) 'LSEG = ',lseg
              call exitus('BDRYINDX')
            endif
            jsp(i,l)=j
          endif
c
          if(on.eq.1.and.landt(i,j).eq.0) then
            on=0
            jep(i,l)=j-1
          endif
c
          if(on.eq.1.and.j.eq.jmt-1) then
            jep(i,l)=j
          endif
c
1175   continue
1174  continue
c
#endif

      do 180 n=1,lseg
      do 180 j=1,jmt
        isz(j,n)=0
        iez(j,n)=0
 180  continue
      lsegp=lseg+1
      do 190 j=2,jmtm1
        l=1
      do 190 i=2,imtm1
        if((fkmz(i-1,j).eq.c0).and.(fkmz(i,j).ne.c0)) isz(j,l)=i
#ifndef coast
        if(i.eq.2) isz(j,l)=2
#else
        if((i.eq.2).and.(fkmz(i-1,j).ne.c0).and.(fkmz(i,j).ne.c0))
     *    isz(j,l)=i
#endif
        if((fkmz(i,j).ne.c0).and.(fkmz(i+1,j).eq.c0)) then
          iez(j,l)=i
          l=l+1
        endif
#ifndef coast
        if((i.eq.imtm1).and.(fkmz(i+1,j).ne.c0)) iez(j,l)=i
#else
        if((i.eq.imtm1).and.(fkmz(i+1,j).ne.c0).and.(fkmz(i,j).ne.c0))
     *    iez(j,l)=i
#endif
        if(l.ge.lsegp) then
          write(stdout,900) 'JSEG = ',lsegp
          call exitus('BDRYINDX')
        endif
 190  continue
c
#ifdef coast
      do 200 n=1,lseg
      do 200 i=1,imt
        jsz(i,n)=0
        jez(i,n)=0
 200  continue
      lsegp=lseg+1
      do 210 i=2,imtm1
        l=1
      do 210 j=2,jmtm1
        if((fkmz(i,j-1).eq.c0).and.(fkmz(i,j).ne.c0)) jsz(i,l)=j
        if((j.eq.2).and.(fkmz(i,j-1).ne.c0).and.(fkmz(i,j).ne.c0))
     *    jsz(i,l)=2
        if((fkmz(i,j).ne.c0).and.(fkmz(i,j+1).eq.c0)) then
          jez(i,l)=j
          l=l+1
        endif
        if((j.eq.jmtm1).and.(fkmz(i,j+1).ne.c0).and.(fkmz(i,j).ne.c0))
     *    jez(i,l)=j
        if(l.ge.lsegp) then
          write(stdout,900) 'ISEG = ',lsegp
          call exitus('BDRYINDX')
        endif
 210  continue
#endif
c
c-----------------------------------------------------------------------
c  Compute start and end indices at velocity points for filters.
c-----------------------------------------------------------------------
c
      do 280 n=1,lseg
      do 280 j=1,jmt
        isq(j,n)=0
        ieq(j,n)=0
 280  continue
      lsegp=lseg+1
      do 290 j=2,jmtm2
        l=1
      do 290 i=2,imtm2
        if((fkmq(i-1,j).eq.c0).and.(fkmq(i,j).ne.c0)) isq(j,l)=i
#ifndef coast
        if(i.eq.2) isq(j,l)=2
#else
        if((i.eq.2).and.(fkmq(i-1,j).ne.c0).and.(fkmq(i,j).ne.c0))
     *    isq(j,l)=i
#endif
        if((fkmq(i,j).ne.c0).and.(fkmq(i+1,j).eq.c0)) then
          ieq(j,l)=i
          l=l+1
        endif
#ifndef coast
        if((i.eq.imtm2).and.(fkmq(i+1,j).ne.c0)) ieq(j,l)=i
#else
        if((i.eq.imtm2).and.(fkmq(i+1,j).ne.c0).and.(fkmq(i,j).ne.c0))
     *    ieq(j,l)=i
#endif
        if(l.ge.lsegp) then
          write(stdout,900) 'JSEG = ',lsegp
          call exitus('BDRYINDX')
        endif
 290  continue
c
#ifdef coast
      do 300 n=1,lseg
      do 300 i=1,imt
        jsq(i,n)=0
        jeq(i,n)=0
 300  continue
      lsegp=lseg+1
      do 310 i=2,imtm2
        l=1
      do 310 j=2,jmtm2
        if((fkmq(i,j-1).eq.c0).and.(fkmq(i,j).ne.c0)) jsq(i,l)=j
        if((j.eq.2).and.(fkmq(i,j-1).ne.c0).and.(fkmq(i,j).ne.c0))
     *    jsq(i,l)=2
        if((fkmq(i,j).ne.c0).and.(fkmq(i,j+1).eq.c0)) then
          jeq(i,l)=j
          l=l+1
        endif
        if((j.eq.jmtm2).and.(fkmq(i,j+1).ne.c0).and.(fkmq(i,j).ne.c0))
     *    jeq(i,l)=j
        if(l.ge.lsegp) then
          write(stdout,900) 'ISEG = ',lsegp
          call exitus('BDRYINDX')
        endif
 310  continue
#endif
c
c-----------------------------------------------------------------------
c  Compute the surface area and volume of the ocean.
c-----------------------------------------------------------------------
c
      area=c0
      volume=c0
      do 320 j=2,jmtm1
      do 320 i=2,imtm1
        if(fkmp(i,j).gt.c0) then
          area=area+cst(j)*dxt(i)*dyt(j)
          volume=volume+cst(j)*dxt(i)*dyt(j)*hd(i,j)
        endif
 320  continue
c
c-----------------------------------------------------------------------
c  Report out area and volume of the ocean.
c-----------------------------------------------------------------------
c
      write(stdout,910) area,volume
c
#ifdef coast
c-----------------------------------------------------------------------
c  Set land-sea mask at velocity points: LANDV=KM sea point,
c                                        LANDV=0 land point.
c-----------------------------------------------------------------------
c
      do 370 j=1,jmt
      do 370 i=1,imt
        landv(i,j)=fkmq(i,j)
 370  continue
c
#endif
#if defined bndy_rlx | defined cstfrc | defined botfrc
c-----------------------------------------------------------------------
c  Set boundary relaxation coefficients.
c-----------------------------------------------------------------------
c
      call set_brlx
c
#endif
#if (defined nest2smaller | defined fltpbd) & defined coast
c-----------------------------------------------------------------------
c  Set streamfunction mask.
c-----------------------------------------------------------------------
c
      call set_pmask
c
#endif
#ifdef fltpbd
c-----------------------------------------------------------------------
c  Set boundary filtering indices.
c-----------------------------------------------------------------------
c
      call getpbn
c
#endif
#ifdef sunflush
c  Flush output buffers.
c
      call flush(stdout)
#endif
c
 900  format(/' BDRYINDX - number of segments exceeded, ',a,i4)
 910  format(/' Surface Area = ',1pe13.6,5x,'  Volume = ',1pe13.6)
c
      return
      end
