subroutine badvec(Hpst,Ut,Vt,Wt)

c  Advection of biological particles, zonal (u), meridional (v),
c  and vertical (w).  Velocities used are computed in zonadv.F,
c  meradv.F, and veradv.F, and represent V*H*Area.  The mean V
c  between the 2 grid cells is used, to reduce diffusivity.
c  Uses either 1st upwind differencing method on biological 
c  particles to improve stability (UPWIND = 1) or centered
c  difference if no stability problem (UPWIND = 0).  (A stability
c  problem results if the sinking rate is high).

#define UPWIND    1

#include "gloparam.F.h"
#include "mapping.h"
#include "ostate.F.h"
#include "definebio.h"
#include "combarr.h"
#include "ogrid.F.h"
      real Wk(len)
      real rmass(len)
      real Ut(len,km),Vt(len,km),Wt(len,km)
      real Hpst(len,km)

      rmass = Area_h/3600.0   !convert to mg/m/s

#if UPWIND
      m = indext2  !index of "past"
#else
      m = indext1  !index of "present"
      m2 = indext2 !index of "past"
#endif

c  Convert all variables to mg/s
      do n = 1,ntyp
       do k = 1,km
        do i = 1,nwater
         P_tend(i,k,n) = P_tend(i,k,n)*rmass(i)*Hpst(i,k)
        enddo
       enddo
      enddo
      do n = nnut+1,nde
       do k = 1,km
        do i = 1,nwater
         ws(i,k,n) = ws(i,k,n)*rmass(i)
        enddo
       enddo
      enddo

c  Zonal advection
c      write(6,*)'Zonal'
      do i = 1,len
       Wk(i) = 0.0
      enddo
      do n = 1,ntyp
       do k = 1,km
        do i = 1,nwater
#if UPWIND
         if (Ut(i,k) .gt. 0.0)then
          Wk(i) = Ut(i,k)*P(i,k,m,n)
         else
          Wk(i) = Ut(i,k)*P(ieast(i),k,m,n)
         endif
#else
         Wk(i) = Ut(i,k)*(P(i,k,m,n)+P(ieast(i),k,m,n))*0.5
#endif
        enddo
        do i = 1,nwater
         trnd = Wk(iwest(i)) - Wk(i)
         P_tend(i,k,n) = P_tend(i,k,n) + trnd
        enddo
       enddo   !k
      enddo    !n
      do i = 1,len
       Wk(i) = 0.0
      enddo

c  Meridional advection
c      write(6,*)'Meridional'
      do i = 1,len
       Wk(i) = 0.0
      enddo
      do n = 1,ntyp
       do k = 1,km
        do i = 1,nwater
#if UPWIND
         if (Vt(i,k) .gt. 0.0)then
          Wk(i) = Vt(i,k)*P(i,k,m,n)
         else
          Wk(i) = Vt(i,k)*P(inorth(i),k,m,n)
         endif
#else
         Wk(i) = Vt(i,k)*(P(i,k,m,n)+P(inorth(i),k,m,n))*0.5
#endif
        enddo
        do i = 1,nwater
         trnd = Wk(isouth(i)) - Wk(i)
         P_tend(i,k,n) = P_tend(i,k,n) + trnd
        enddo
       enddo  !k
      enddo   !n
      do i = 1,len
       Wk(i) = 0.0
      enddo

c  Vertical advection      
c      write(6,*)'Vertical'
      do n = 1,ntyp
#if UPWIND
       do k = 1,km-1
        do i = 1,nwater
         if (Wt(i,k) .gt. 0.0) then
          trnd = P(i,k+1,m,n)*Wt(i,k)
         else
          trnd = P(i,k,m,n)*Wt(i,k)
         endif
         P_tend(i,k,n)   = P_tend(i,k,n)   + trnd
         P_tend(i,k+1,n) = P_tend(i,k+1,n) - trnd
        enddo ! i
       enddo  ! k
 
       k = km
       do i = 1,nwater
        if (Wt(i,k) .gt. 0.0) then
         trnd = Pdeep(n)*Wt(i,k)
        else
         trnd = P(i,k,m,n)*Wt(i,k)
        endif
        P_tend(i,k,n) = P_tend(i,k,n) + trnd
       enddo ! i
#else
      k = 1
      do i = 1,nwater
        if (Wt(i,k) .gt. 0.0) then
         trnd = 0.5*(P(i,k,m,n)+P(i,k+1,m,n))*Wt(i,k)
        else
         trnd = P(i,k,m2,n)*Wt(i,k)
        endif
        P_tend(i,k,n) = P_tend(i,k,n) + trnd
        P_tend(i,k+1,n) = P_tend(i,k+1,n) - trnd
      enddo ! i
      do k = 2,km-1
       do i = 1,nwater
          trnd = (P(i,k,m,n)+P(i,k+1,m,n))*0.5*Wt(i,k)
          P_tend(i,k,n) = P_tend(i,k,n) + trnd
          P_tend(i,k+1,n) = P_tend(i,k+1,n) - trnd
       enddo ! i
      enddo ! k
      k = km
      do i = 1,nwater
        if (Wt(i,k) .gt. 0.0) then
         trnd = Pdeep(n)*Wt(i,k)
         P_tend(i,k,n) = P_tend(i,k,n) + trnd
        else
         trnd = (P(i,k,m2,n)+Pdeep(n))*Wt(i,k)*0.5
         P_tend(i,k,n) = P_tend(i,k,n) + trnd
        endif
      enddo ! i
#endif
      enddo ! n


c  Sinking 
c      write(6,*)'Sinking'
      do n = nnut+1,nde
#if UPWIND
       do k = 1,km-1
        do i = 1,nwater
         trnd = P(i,k,m,n)*ws(i,k,n)
         P_tend(i,k,n)   = P_tend(i,k,n)   - trnd
         P_tend(i,k+1,n) = P_tend(i,k+1,n) + trnd
        enddo ! i
       enddo  ! k
       k = km
       do i = 1,nwater
        trnd = P(i,k,m,n)*ws(i,k,n)
        P_tend(i,k,n)   = P_tend(i,k,n)   - trnd
       enddo ! i
 
#else
       k = 1
       do i = 1,nwater
        trnd = P(i,k,m2,n)*ws(i,k,n)
        P_tend(i,k,n) = P_tend(i,k,n) - trnd
        P_tend(i,k+1,n) = P_tend(i,k+1,n) + trnd
       enddo ! i
       do k = 2,km-1
        do i = 1,nwater
         trnd = (P(i,k,m,n)+P(i,k+1,m,n))*0.5*ws(i,k,n)
         P_tend(i,k,n) = P_tend(i,k,n) - trnd
         P_tend(i,k+1,n) = P_tend(i,k+1,n) + trnd
        enddo ! i
       enddo ! k
       k = km
       do i = 1,nwater
        trnd = (P(i,k,m,n)+Pdeep(n))*0.5*ws(i,k,n)
        P_tend(i,k,n) = P_tend(i,k,n) - trnd
       enddo ! i
#endif
      enddo ! n

      return
      end