subroutine fsigma(track)
c
c=======================================================================
c                                                                    ===
c  This routine determines the density at the initial location       ===
c  of the given Lagrangian drifter.                                  ===
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(3) -- Z coord. of current drifter position        ===
c                                                                    ===
c     On Output:                                                     ===
c                                                                    ===
c            TRACK(8) -- density at initial position                 ===
c                                                                    ===
c     Calls:  BES1D, BRACKET_Z                                       ===
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,levdn,levup
#ifdef cyclic
     *        ,ncyc
#endif
      FLOAT
     *      bes1d,dzd,dzu,x,xmax,xmix,y,ymax,ymjy,z,zdn,zup
      FLOAT
     *      track(8),sigwk1(4),sigwk2(4)
      save first,xmax,ymax
      data first/.true./
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c  On first time through, set bounds for horizontal interpolations.
c-----------------------------------------------------------------------
c
      if(first) then
        xmax=FLoaT(imtm1)
        ymax=FLoaT(jmtm1)
        first=.false.
      endif
c
c-----------------------------------------------------------------------
c  Set up parameters for interpolations.
c-----------------------------------------------------------------------
c
      x=track(1)
      y=track(2)
      z=max(track(3)*m2cm,c0)
c
      ix=int(x)
      xmix=x-FLoaT(ix)
      ibry=0
      if(x.lt.c2) then
        ibry=1
      elseif(x.gt.xmax) then
        ibry=2
      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  Interpolate density to drifter depth and zonal position.
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 bracket_z(iwk,jwk,z,tgrid,levup,levdn,dzu,dzd,zup,zdn)
          sigwk1(i)=dzu*sig_vol(iwk,jwk,levup)+
     *              dzd*sig_vol(iwk,jwk,levdn)
  10     continue
         sigwk2(j)=bes1d(xmix,sigwk1(1),sigwk1(2),sigwk1(3),sigwk1(4),
     *                   ibry)
  20  continue
c
c-----------------------------------------------------------------------
c  Interpolate density to meridional drifter position.
c-----------------------------------------------------------------------
c
      track(8)=bes1d(ymjy,sigwk2(1),sigwk2(2),sigwk2(3),sigwk2(4),jbry)
      return
      end