subroutine vavgv (ipt,jpt,z,hd,u,ub,v,vb,uz,uzb,vz,vzb)
c
c=======================================================================
c                                                                    ===
c  This routine performs Vertical Averaging of horizontal Velocity.  ===
c  This particular routine averages the velocity based upon the      ===
c  fraction of a submerged sphere in a given level.                  ===
c                                                                    ===
c  ------                                                            ===
c  Input:                                                            ===
c  ------                                                            ===
c                                                                    ===
c     IPT   X-index of desired point.                    (integer)   ===
c     JPT   Y-index of desired point.                    (integer)   ===
c     Z     Desired depth of interpolation.              (real)      ===
c     HD    Diameter (height) of the drifter.            (real)      ===
c     U     Internal zonal velocities, current.          (real array)===
c     UB    Internal zonal velocities, previous.         (real array)===
c     V     Internal meridional velocities, current.     (real array)===
c     VB    Internal meridional velocities, previous.    (real array)===
c                                                                    ===
c  Common Blocks:                                                    ===
c                                                                    ===
c  /VERTICAL/                                                        ===
c                                                                    ===
c     HVZ   depth at the bottom of the velocity boxes.   (real array)===
c                                                                    ===
c  -------                                                           ===
c  Output:                                                           ===
c  -------                                                           ===
c                                                                    ===
c     UZ    Interpolated zonal velocity, current.        (real)      ===
c     UZB   Interpolated zonal velocity, previous.       (real)      ===
c     VZ    Interpolated meridional velocity, current.   (real)      ===
c     VZB   Interpolated meridional velocity, previous.  (real)      ===
c                                                                    ===
c  Calls:  BRACKET_Z                                                 ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <vertical.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer ipt,jpt,k,knd,kst,levdn,levup
      FLOAT
     *      arg1,arg2,dzd,dzu,hd,totarea,uz,uzb,vz,vzb,wt,z,zdn,zup,znd,
     *      zst
      FLOAT
     *      u(imt,km),ub(imt,km),v(imt,km),vb(imt,km)
      FLOAT
     *      wtfn
c
      wtfn(arg1) = arg1*sqrt(c1-arg1*arg1) + asin(arg1)
c
c=======================================================================
c  Begin executable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Determine averaging range.
c-----------------------------------------------------------------------
c
c  Find starting level.
c
      zst = max( (z-p5*hd) , c0 )
c
      call bracket_z (ipt,jpt,zst,vgrid,levup,levdn,dzu,dzd,zup,zdn)
c
      if (hvz(ipt,levup,jpt).gt.zst) then
         kst = levup
       else
         kst = levdn
      end if
c
c  Find ending level.
c
      znd = min( (z+p5*hd) , hvz(ipt,km,jpt) )
c
      call bracket_z (ipt,jpt,znd,vgrid,levup,levdn,dzu,dzd,zup,zdn)
c
      if (hvz(ipt,levup,jpt).lt.znd) then
         knd = levdn
       else
         knd = levup
      end if
c
c  Find total cross-sectional area of drifter within water.
c
      arg1 = (znd-z)*c2/hd
      arg2 = (zst-z)*c2/hd
      totarea = wtfn(arg1) - wtfn(arg2)
c
c-----------------------------------------------------------------------
c  Average internal mode horizontal velocities.
c-----------------------------------------------------------------------
c
      if (znd.gt.hvz(ipt,kst,jpt)) then
c
c  --- Drifter extends over several levels.
c
         arg1 = (hvz(ipt,kst,jpt)-z)*c2/hd
         arg2 = (zst-z)*c2/hd
         wt   = (wtfn(arg1)-wtfn(arg2))/totarea
         uz   = u (ipt,kst)*wt
         uzb  = ub(ipt,kst)*wt
         vz   = v (ipt,kst)*wt
         vzb  = vb(ipt,kst)*wt
c
         do 10 k = kst+1, knd-1
            arg1 = (hvz(ipt,k  ,jpt)-z)*c2/hd
            arg2 = (hvz(ipt,k-1,jpt)-z)*c2/hd
            wt   = (wtfn(arg1)-wtfn(arg2))/totarea
            uz   = uz  + u (ipt,k)*wt
            uzb  = uzb + ub(ipt,k)*wt
            vz   = vz  + v (ipt,k)*wt
            vzb  = vzb + vb(ipt,k)*wt
  10     continue
c
         arg1 = (znd-z)*c2/hd
         arg2 = (hvz(ipt,knd-1,jpt)-z)*c2/hd
         wt   = (wtfn(arg1)-wtfn(arg2))/totarea
         uz   = uz  + u (ipt,knd)*wt
         uzb  = uzb + ub(ipt,knd)*wt
         vz   = vz  + v (ipt,knd)*wt
         vzb  = vzb + vb(ipt,knd)*wt
c
       else
c
c  --- Drifter is contained within one level.
c
         uz   = u (ipt,kst)
         uzb  = ub(ipt,kst)
         vz   = v (ipt,kst)
         vzb  = vb(ipt,kst)
c
      end if
c
      return
      end