subroutine sigma2z(track)
c
c=======================================================================
c                                                                    ===
c  This routine determines the correct depth for a Lagrangian        ===
c  drifter which follows isopycnal surfaces. If no such level        ===
c  is found, the depth remains  unchanged.   If more than one        ===
c  such level is found, the nearest is used.                         ===
c                                                                    ===
c     On Input:                                                      ===
c                                                                    ===
c     TRACK  descriptor of Lagrangian drifter (real array; input)    ===
c            TRACK(1) -- X coord. of current drifter position        ===
c            TRACK(2) -- Y coord. of current drifter position        ===
c            TRACK(8) -- density anomaly at initial position         ===
c                                                                    ===
c     On Output:                                                     ===
c                                                                    ===
c            TRACK(3) -- Z coord. of current drifter position        ===
c                                                                    ===
c     Calls:  BES1D, SIG2Z                                           ===
c                                                                    ===
c     Common blocks: (only relevant variables documented)            ===
c                                                                    ===
c     /TRKFLD/                                                       ===
c                                                                    ===
c     SIG_VOL   density at every grid point (real array; input)      ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <trkfld.h>
c
c-----------------------------------------------------------------------
c  Define local data
c-----------------------------------------------------------------------
c
      logical first
      integer i,ibry,iwk,ix,j,jbry,jwk,jy
#ifdef cyclic
     *        ,ncyc
#endif
      FLOAT
     *      bes1d,sig,x,xmix,y,ymax,ymjy,z
#ifndef cyclic
     *      ,xmax
#endif
      FLOAT
     *      track(8),zwk1(4),zwk2(4)
#ifndef cyclic
      save first,xmax,ymax
# else
      save first,ymax
#endif
      data first/.true./
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c  On first pass, determine bounds for horizontal interpolations.
c-----------------------------------------------------------------------
c
      if(first) then
#ifndef cyclic
        xmax=FLoaT(imtm1)
#endif
        ymax=FLoaT(jmtm1)
        first=.false.
      endif
c
c-----------------------------------------------------------------------
c  Determine position and desired density.  Compute parameters for
c  horizontal interpolations.
c-----------------------------------------------------------------------
c
      x=track(1)
      y=track(2)
      z=max(track(3)*m2cm,c0)
      sig=track(8)
c
      ix=int(x)
      xmix=x-FLoaT(ix)
      ibry=0
#ifndef cyclic
      if(x.lt.c2) then
        ibry=1
      elseif(x.gt.xmax) then
        ibry=2
      endif
#endif
c
      jy=int(y)
      ymjy=y-FLoaT(jy)
      jbry=0
      if(y.lt.c2) then
        jbry=1
      elseif(y.gt.ymax) then
        jbry=2
      endif
c
c-----------------------------------------------------------------------
c  Determine depths over needed "footprint" and zonally interpolate.
c-----------------------------------------------------------------------
c
      do 20 j=1,4
        jwk=max(min(jy+j-2,jmt),1)
        do 10 i=1,4
#ifndef cyclic
          iwk=max(min(ix+i-2,imt),1)
# else
          iwk=ix+i-2
          ncyc=(iwk-1)/imtm2
          if(iwk.lt.1) ncyc=ncyc-1
          iwk=iwk-ncyc*imtm2
#endif
          call sig2z(iwk,jwk,z,sig,sig_vol,zwk1(i))
 10     continue
        zwk2(j)=bes1d(xmix,zwk1(1),zwk1(2),zwk1(3),zwk1(4),ibry)
 20   continue
c
c-----------------------------------------------------------------------
c  Meridionally interpolate depths.
c-----------------------------------------------------------------------
c
      z=bes1d(ymjy,zwk2(1),zwk2(2),zwk2(3),zwk2(4),jbry)
      track(3)=max(cm2m*z,c0)
      return
      end