subroutine sig2z(i,j,zguess,sig,sig_vol,zfinal)
c
c=======================================================================
c                                                                    ===
c  This routine determines at what depth the target density occurs   ===
c  at the given horizontal position.                                 ===
c                                                                    ===
c     On Input:                                                      ===
c                                                                    ===
c     I, J      horizontal grid location of column search (integer)  ===
c     SIG       target density (real)                                ===
c     SIG_VOL   density at every grid point (real array)             ===
c     ZGUESS    initial guess for depth (real)                       ===
c                                                                    ===
c     On Output:                                                     ===
c                                                                    ===
c     ZFINAL    final result for depth (real)                        ===
c                                                                    ===
c     Calls:  BRACKET_Z, SIG_SEARCH                                  ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data
c-----------------------------------------------------------------------
c
      logical down,found_dn,found_up,up
      integer i,j,levddn,levdn,levdup,levudn,levup,levuup
      FLOAT
     *      dsig,dzdn,dzup,sig,sigmax,sigmin,small,zddn,zdn,zdup,zfinal,
     *      zguess,zudn,zup,zuup
      FLOAT
     *      sig_vol(imt,jmt,km)
      parameter (down=.false.,up=.true.,small=c1em30)
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c  Find model vertical levels which bracket current tracer vertical
c  position.
c-----------------------------------------------------------------------
c
      zfinal=zguess
      call bracket_z(i,j,zfinal,tgrid,levup,levdn,dzup,dzdn,zup,zdn)
c
c-----------------------------------------------------------------------
c  Determine minimum & maximum of density anomalies at bracketting
c  levels.
c-----------------------------------------------------------------------
c
      sigmin=sig_vol(i,j,levup)
      sigmax=sigmin
      sigmin=min(sigmin,sig_vol(i,j,levdn))
      sigmax=max(sigmax,sig_vol(i,j,levdn))
c
c-----------------------------------------------------------------------
c  Find depth corresponding to initial density anomaly.
c-----------------------------------------------------------------------
c
      if((sigmin.le.sig).and.(sig.le.sigmax)) then
c
c-----------------------------------------------------------------------
c  If desired anomaly within range described by anomalies at
c  bracketting levels, linearly interpolate depth.
c-----------------------------------------------------------------------
c
        dsig=sig_vol(i,j,levdn)-sig_vol(i,j,levup)
        if(abs(dsig).ge.small) zfinal=(sig_vol(i,j,levdn)-sig)/dsig*zup+
     *                                (sig-sig_vol(i,j,levup))/dsig*zdn
      else
c
c-----------------------------------------------------------------------
c  Search above and below current level for desired sigma.
c-----------------------------------------------------------------------
c
        call sig_search(up,i,j,sig,sig_vol,levup,levuup,levudn,zuup,
     *                  zudn,found_up)
        call sig_search(down,i,j,sig,sig_vol,levdn,levdup,levddn,zdup,
     *                  zddn,found_dn)
        if(found_up) then
c
c-----------------------------------------------------------------------
c  Compute upper candidate, if found.
c-----------------------------------------------------------------------
c
          dsig=sig_vol(i,j,levudn)-sig_vol(i,j,levuup)
          if(abs(dsig).ge.small) then
            zup=(sig_vol(i,j,levudn)-sig)/dsig*zuup+
     *          (sig-sig_vol(i,j,levuup))/dsig*zudn
          else
            zup=zudn
          endif
        endif
        if(found_dn) then
c
c-----------------------------------------------------------------------
c  Compute lower candidate, if found.
c-----------------------------------------------------------------------
c
          dsig=sig_vol(i,j,levddn)-sig_vol(i,j,levdup)
          if(abs(dsig).ge.small) then
            zdn=(sig_vol(i,j,levddn)-sig)/dsig*zdup+
     *          (sig-sig_vol(i,j,levdup))/dsig*zddn
          else
            zdn=zdup
          endif
        endif
        if(found_up.and.found_dn) then
c
c-----------------------------------------------------------------------
c  Possible Z-values both above & below, pick nearest.
c-----------------------------------------------------------------------
c
          if(abs(zguess-zup).le.abs(zdn-zguess)) then
            zfinal=zup
          else
            zfinal=zdn
          endif
        elseif(found_up) then
c
c-----------------------------------------------------------------------
c  Only possible Z value is above current level.
c-----------------------------------------------------------------------
c
          zfinal=zup
        elseif(found_dn) then
c
c-----------------------------------------------------------------------
c Only possible z-value is below current level.
c-----------------------------------------------------------------------
c
          zfinal=zdn
        endif
      endif
      return
      end