c---------------------------------------------------------
      subroutine baro_sum (npt, nz, nzi, uc, vc, ubar, vbar, dept)
c---------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension uc(npt,nz),vc(npt,nz),ubar(npt),vbar(npt),nzi(npt),dept(npt)

      do i = 1, npt
         ubar(i) = uc(i,1)
         vbar(i) = vc(i,1)
      enddo
      
      do i = 1, npt
         do k = 2, nzi(i)
            ubar(i) = ubar(i) + uc(i,k)
            vbar(i) = vbar(i) + vc(i,k)
         enddo
      enddo

      do i = 1, npt
         depi = 1./dept(i)
         ubar(i) = depi*ubar(i)
         vbar(i) = depi*vbar(i)
      enddo

      return
      end

      subroutine baro_comp (npt,nz,nzi,uforc,vforc,tfu,tfv
     *     ,corx,cory,pgfx,pgfy,cor,pgf,gradE)
c-----------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension uforc(npt),vforc(npt),tfu(npt),tfv(npt),cor(npt,2)
      dimension gradE(npt,2),pgf(npt,2)
      dimension nzi(npt),pgfx(npt,nz),pgfy(npt,nz),corx(npt,nz),cory(npt,nz)

c       the Coriolis terms are done implicitly in baro_solv
c          and the pressure gradient terms are done separately
c          so we must accumulate them in cor and gradE and subtract from them
c          from the vertically integrated forcings, uforc and vforc
      
      do i = 1, npt
         mz = nzi(i)
         do k = 1, mz
            cor(i,1) = cor(i,1) + corx(i,k)
            cor(i,2) = cor(i,2) + cory(i,k)
            pgf(i,1) = pgf(i,1) + pgfx(i,k)
            pgf(i,2) = pgf(i,2) + pgfy(i,k)
         enddo
      enddo

      do i = 1, npt
         uforc(i) = tfu(i) - cor(i,1) - pgf(i,1)
         vforc(i) = tfv(i) - cor(i,2) - pgf(i,2)
      enddo

      return
      end

c------------------------------------------------------------------
      subroutine baro_rhs (nx,ny,npt,isk,uforc,vforc,rhs,tp,
     *     gradE,gradH,binv,dt_sub,cor,pgf,dx2i,dy2i,ycos,inx,iny,lxx,lyy)
c------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)

      common/grid/nxp,nyp,nxyc,nz,nbx,nby,ncs,land,nlo,npbc
      dimension isk(npt),inx(nx),iny(ny),lxx(1),lyy(1)

      dimension uforc(npt),vforc(npt),rhs(npt),tp(npt),
     *          pgf(npt,2),gradE(npt,2),gradH(npt,2),cor(npt,2),
     *          dx2i(1), dy2i(1), ycos(1)

      dtfac = binv/dt_sub

c.................rhs = curl(forcing-gradE) + (forcing-gradE) X grad(H)/H

      call curlz (nx,ny,npt,npt,uforc,vforc,dx2i,dy2i,ycos,
     *                       isk,inx,iny,tp,nbx,nby,lxx,lyy)

      do i = 1, npt
         uforc(i) = uforc(i) - gradE(i,1)
         vforc(i) = vforc(i) - gradE(i,2)

         rhs(i) = tp(i) + uforc(i)*gradH(i,2)- vforc(i)*gradH(i,1) 

c  scale forcing (barotropic takes a dt_sub time step) 

         rhs(i) = dtfac*rhs(i)

         uforc(i) = dtfac*uforc(i)
         vforc(i) = dtfac*vforc(i)
      enddo
      
      return
      end

c---------------------------------------------------------------------------
      subroutine baro_updat1(npt,nz,nzi,h,binv,abinv,uc,vc,u,v,ubar,vbar)
c----------------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension h(npt,nz),uc(npt,nz),vc(npt,nz),
     *     u(npt),v(npt),ubar(npt), vbar(npt), nzi(npt)

      do i = 1, npt
         do k = 1, nzi(i)
            hik = h(i,k)
            uc(i,k) = uc(i,k) + hik*(u(i) - ubar(i))
            vc(i,k) = vc(i,k) + hik*(v(i) - vbar(i))
         enddo
      enddo

      do i = 1, npt
         ubar(i) = u(i)
         vbar(i) = v(i)
      enddo

      return
      end


c---------------------------------------------------------------------------
      subroutine baro_updat2(npt,nz,binv,abinv,cor,pgf,gradE,tfu,tfv)
c----------------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      dimension cor(npt,2),gradE(npt,2), tfu(npt), tfv(npt),pgf(npt,2)

      do i = 1, npt
         cor(i,1)   = abinv * cor(i,1)
         cor(i,2)   = abinv * cor(i,2)
         pgf(i,1)   = abinv * pgf(i,1)
         pgf(i,2)   = abinv * pgf(i,2)
         gradE(i,1) = abinv * gradE(i,1)
         gradE(i,2) = abinv * gradE(i,2)
         tfu(i) = abinv * tfu(i)
         tfv(i) = abinv * tfv(i)
      enddo

      return
      end


c---------------------------------------------------------------------------
      subroutine baro_reset (nstep,npt,nz,nzi,nzibx,nziby,dept,h,ubar,vbar,u,v,
     *                      lxxk,lyyk,uc,vc)
c-------------------------------------------------------------------------


c    This subroutine is responsible for preserving the rigid lid 
c    assumption. It must be called immediately before ddiv to
c    work properly
      dimension h(npt,nz),uc(npt,nz),vc(npt,nz),ubar(npt),vbar(npt)
      dimension u(npt,nz),v(npt,nz), dept(npt),nzi(npt),nzibx(npt),nziby(npt)
      dimension lxxk(1), lyyk(1)
      common /param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20),
     *               itherm,mlc,limp

      call bcset (mbc,lxxk,lyyk,npt,u,v)

      do i = 1, npt
         do k = 1, nzi(i)
            hi = h(i,k)
            uc(i,k) = u(i,k)*hi
            vc(i,k) = v(i,k)*hi
         enddo
      enddo

      call baro_sum (npt, nz, nzi, uc, vc, u, v, dept)
      
      do i = 1, npt
         dep = 0.
         mz = nzibx(i)
         do k = 1, mz
            dep = dep + h(i,k)
         enddo
         fdep = dept(i)/dep
         do k = 1, mz
            uc(i,k) = uc(i,k) + h(i,k)*fdep*(ubar(i) - u(i,1))
         enddo

         dep = 0.
         mz = nziby(i)
         do k = 1, mz
            dep = dep + h(i,k)
         enddo
         fdep = dept(i)/dep
         do k = 1, mz
            vc(i,k) = vc(i,k) + h(i,k)*fdep*(vbar(i) - v(i,1))
         enddo
      enddo

c      do i = 1, npt
c         mz = nzibx(i)
c         uc(i,mz) = uc(i,mz) + dept(i)*(ubar(i) - u(i,1))
c         mz = nziby(i)
c         vc(i,mz) = vc(i,mz) + dept(i)*(vbar(i) - v(i,1))
c      enddo

      call decap (npt, nz, nzi, u,v,uc,vc,h)

      return
      end

      subroutine baro_split (npt,nz,nzi,h,dept,fu,fv,u,v,tfu,tfv)
c-----------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension h(npt,1),fu(npt,1),fv(npt,1),tfu(1),tfv(1),
     *          nzi(npt),dept(npt),u(1),v(1)


      call baro_sum (npt,nz,nzi,fu,fv,u,v,dept)

      do i = 1, npt
         dep = dept(i)
         tfu(i)  = tfu(i) + dep*u(i)
         tfv(i)  = tfv(i) + dep*v(i)
         do k = 1, nzi(i)
            hk = h(i,k)
            fu(i,k) = fu(i,k) - hk* u(i)
            fv(i,k) = fv(i,k) - hk* v(i)
         enddo
      enddo

      return
      end

