subroutine bracket_z(i,j,z,keytv,levup,levdn,dzup,dzdn,zup,zdn)
c
c=======================================================================
c                                                                    ===
c  This routine determines the model levels and corresponding depths ===
c  which bracket the given depth.                                    ===
c                                                                    ===
c     On Input:                                                      ===
c                                                                    ===
c     I, J   horizontal position at which to search column. (integer)===
c     KEYTV  flag to indicate whether considering velocity or        ===
c            tracer grid.   (integer)                                ===
c     Z      depth to be bracketted (cm; real)                       ===
c                                                                    ===
c     On Output:                                                     ===
c                                                                    ===
c     DZUP   linear interpolation weight for data at LEVUP (real)    ===
c     DZDN   linear interpolation weight for data at LEVDN (real)    ===
c     LEVUP  number of the level immediately above Z (integer)       ===
c     LEVDN  number of the level immediately below Z (integer)       ===
c     ZUP    depth of level immediately above Z (cm; real)           ===
c     ZDN    depth of level immediately below Z (cm; real)           ===
c                                                                    ===
c     Calls:    DEPTHSLAB                                            ===
c                                                                    ===
c=======================================================================
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer levup,levdn,i,j,keytv
      FLOAT
     *      dzdn,dzup,z,zdn,zup
      FLOAT
     *      zslab(imt,km)
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c  Get a slab of depths at the desired location.
c-----------------------------------------------------------------------
c
      call depthslab (j,keytv,zslab)
c
c-----------------------------------------------------------------------
c  Find the levels bracketting the given depth.
c-----------------------------------------------------------------------
c
      levdn=2
      do 10 while ((z.gt.zslab(i,levdn)).and.(levdn.lt.km))
        levdn=levdn+1
  10  continue
      levup=levdn-1
c
c-----------------------------------------------------------------------
c  Get the depths at the bracketting levels.
c-----------------------------------------------------------------------
c
      zup=zslab(i,levup)
      zdn=zslab(i,levdn)
c
c-----------------------------------------------------------------------
c  Get the weights for linear interpolation.
c-----------------------------------------------------------------------
c
      if((z.ge.zslab(i,1)).and.(z.le.zslab(i,km))) then
c
c-----------------------------------------------------------------------
c  An interior value, interpolate.
c-----------------------------------------------------------------------
c
        dzup=(zdn-z)/(zdn-zup)
        dzdn=c1-dzup
      elseif(z.lt.zslab(i,1)) then
c
c-----------------------------------------------------------------------
c  Too shallow, use nearest value.
c-----------------------------------------------------------------------
c
        dzup=c1
        dzdn=c0
      else
c
c-----------------------------------------------------------------------
c  Too deep, use nearest value.
c-----------------------------------------------------------------------
c
        dzup=c0
        dzdn=c1
      endif
      return
      end