subroutine barotp(m,n)
      use mod_xc    ! HYCOM communication interface
      use mod_pipe  ! HYCOM debugging interface
c
c --- micom version 2.8
#if defined(STOKES)
      use mod_stokes  !    HYCOM Stokes Drift
#endif      
      implicit none
c
      include 'common_blocks.h'
c
      integer m,n
c
c --- ------------------------------------------------------------------------
c --- advance barotropic equations.
c ---   on entry: -n- is time t-dt, -m- is time t
c ---   on exit:                    -m- is time t, -n- is time t+dt
c ---   time level 3 is only used internally (n and m are always 1 or 2).
c
c --- LeapFrog version based on:
c ---   Y. Morel, Baraille, R., Pichon A. (2008) "Time splitting and
c ---   linear stability of the slow part of the barotropic component", 
c ---   Ocean Modeling, 23, pp 73-81.
c --- ------------------------------------------------------------------------
c
      logical    lpipe_barotp
      parameter (lpipe_barotp=.false.)
      logical    ldebug_barotp
      parameter (ldebug_barotp=.false.)
c
      real    q,pbudel,pbvdel,utndcy,vtndcy,wblpf
      real*8  sump
      integer i,j,l,lll,ml,nl,mn,lstep1,mbdy
c
      real, save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) ::
     &        pbaros,ubaros,vbaros
c
      mbdy = 6
c#if defined (BARO_SLP) 
CC Fanf: Code provided by SHOM to account for the barotropic response from the
CC SLP
C         margin = 1
C         do j = 1-margin, jj+margin
C         do l = 1, isp(j)
C         do i = max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l))
C            util4(i,j) = slp(i,j  ,l0)*w0+slp(i,j  ,l1)*w1+
C     .                   slp(i,j  ,l2)*w2+slp(i,j  ,l3)*w3
C         enddo
C         enddo
C         enddo
C         margin = 0
C         do j = 1-margin, jj+margin
C         do l = 1, isu(j)
C         do i = max(1-margin,ifu(j,l)), min(ii+margin,ilu(j,l))
C            ia = i-1
C            if (iu(ia,j).ne.1) ia = i
C            utotn(i,j) = utotn(i,j) 
C     .                 - thref*(util4(i,j)-util4(ia,j))*scuxi(i,j)   
C         enddo
C         enddo
C         enddo
C         do j = 1-margin, jj+margin
C         do l = 1, isv(j)
C         do i = max(1-margin,ifv(j,l)), min(ii+margin,ilv(j,l))
C            ja = j-1
C            if (iv(i,ja).ne.1) ja = j
C            vtotn(i,j) = vtotn(i,j) 
C     .                 - thref*(util4(i,j)-util4(i,ja))*scvyi(i,j)   
C         enddo
C         enddo
C         enddo
C#endif
c    
C$RMY c   >>>>>>>>>>>>>>>>wave induced pressure      
C$RMY #if defined(STOKES_3D)
C$RMY       if (stpres) then
C$RMY          margin = 0
C$RMY !$OMP PARALLEL DO PRIVATE(j,i)
C$RMY !$OMP&         SCHEDULE(STATIC,jblk)
C$RMY          do j = 1-margin, jj+margin
C$RMY          do i = 1-margin, ii+margin
C$RMY             if (SEA_U) then
C$RMY                utotn(i,j) = utotn(i,j)
C$RMY      .                 - (ws0*(sj(i,j,ls0)-sj(i-1,j,ls0)) +
C$RMY      .                    ws1*(sj(i,j,ls1)-sj(i-1,j,ls1)))*scuxi(i,j)
C$RMY             endif !iu
C$RMY             if (SEA_V) then
C$RMY                vtotn(i,j) = vtotn(i,j)
C$RMY      .                 - (ws0*(sj(i,j,ls0)-sj(i,j-1,ls0)) +
C$RMY      .                    ws1*(sj(i,j,ls1)-sj(i,j-1,ls1)))*scvyi(i,j)
C$RMY             endif !iv
C$RMY          enddo !i
C$RMY          enddo !j
C$RMY       endif
C$RMY #endif
C$RMY c
C$RMY c   >>>>>>>>>>>>>>>>wave induced radiation stress
C$RMY #if defined(STOKES_2D)
C$RMY       margin = 0
C$RMY !$OMP PARALLEL DO PRIVATE(j,i)
C$RMY !$OMP&         SCHEDULE(STATIC,jblk)
C$RMY       do j = 1-margin, jj+margin
C$RMY       do i = 1-margin, ii+margin
C$RMY          if (SEA_U) then
C$RMY             utotn(i,j) = utotn(i,j)
C$RMY      .              - ( ws0*d_sx(i,j,ls0) + ws1*d_sx(i,j,ls1) )
C$RMY      .               / max(onem,
C$RMY      .                 0.5*(pbot(i-1,j)*onetai(i-1,j)+pbavg(i-1,j,m) + 
C$RMY      .                      pbot(i,j)*onetai(i,j)    +pbavg(i,j,m)))
C$RMY          endif !iu
C$RMY          if (SEA_V) then
C$RMY             vtotn(i,j) = vtotn(i,j)
C$RMY      .              - ( ws0*d_sy(i,j,ls0) + ws1*d_sy(i,j,ls1) )
C$RMY      .               / max(onem,
C$RMY      .                 0.5*(pbot(i-1,j)*onetai(i,j-1)+pbavg(i,j-1,m) + 
C$RMY      .                      pbot(i,j)*onetai(i,j)    +pbavg(i,j,m)))
C$RMY          endif !iv
C$RMY       enddo !i
C$RMY       enddo !j
C$RMY #endif
c
c --- utotn,vtotn from momtum is time step t-1 to t+1 barotropic tendency      
      call xctilr(utotn(  1-nbdy,1-nbdy    ),1, 1, 6,6, halo_uv)
      call xctilr(vtotn(  1-nbdy,1-nbdy    ),1, 1, 6,6, halo_vv)
c
      if     (lpipe .and. lpipe_barotp) then
c ---   compare two model runs.
        call pipe_compare_sym2(utotn, iu,'barotp:utotn',
     &                         vtotn, iv,'barotp:vtotn')
        call pipe_compare_sym1(pvtrop,iq,'barotp:pvtrp')
      endif

    !   write(lp,*)'from initit barooo. maval ubavg=',
    ! &     maxval(ubavg(:,:))
c
c --- explicit time integration of barotropic flow (forward-backward scheme)
c --- in order to combine forward-backward scheme with leapfrog treatment of
c --- coriolis term, v-eqn must be solved before u-eqn every other time step
c
      if     (btrlfr) then
        if     (delt1.ne.baclin) then  !not on very 1st time step
C ---     start at time level t-dt and go to t+dt.
          lstep1 = lstep + lstep  !more stable, but also more expensive
!$OMP     PARALLEL DO PRIVATE(j,i)
!$OMP&             SCHEDULE(STATIC,jblk)
          do j=1,jj
            do i=1,ii
              pbavg(i,j,3) = pbavg(i,j,n)
              ubavg(i,j,3) = ubavg(i,j,n)
              vbavg(i,j,3) = vbavg(i,j,n)
c
              pbaros(i,j) = (1.-2.*wts2)*pbavg(i,j,m)+wts2*pbavg(i,j,n)
              ubaros(i,j) = (1.-2.*wts2)*ubavg(i,j,m)+wts2*ubavg(i,j,n)
              vbaros(i,j) = (1.-2.*wts2)*vbavg(i,j,m)+wts2*vbavg(i,j,n)
            enddo !i
          enddo !j
        else !1st time step
C ---     start at time level t and go to t+dt.
          lstep1 = lstep
!$OMP     PARALLEL DO PRIVATE(j,i)
!$OMP&             SCHEDULE(STATIC,jblk)
          do j=1,jj
            do i=1,ii
              pbavg(i,j,n) = pbavg(i,j,m)
              ubavg(i,j,n) = ubavg(i,j,m)
              vbavg(i,j,n) = vbavg(i,j,m)
              pbavg(i,j,3) = pbavg(i,j,m)
              ubavg(i,j,3) = ubavg(i,j,m)
              vbavg(i,j,3) = vbavg(i,j,m)
            enddo !i
          enddo !j
        endif !usual:1st time step
      else
C ---   start at time level t    and go to t+dt.
        lstep1 = lstep          !original, less stable, method
!$OMP   PARALLEL DO PRIVATE(j,i)
!$OMP&           SCHEDULE(STATIC,jblk)
         do j=1,jj
           do i=1,ii
             pbavg(i,j,n) = pbavg(i,j,m)
             ubavg(i,j,n) = ubavg(i,j,m)
             vbavg(i,j,n) = vbavg(i,j,m)
           enddo !i
         enddo !j
      endif !btrlfr
c
      if     (btrlfr) then
        wblpf = 0.0   !1st minor time step, lll=1, only
      else
        wblpf = wbaro
      endif
c
      do 840 lll=1,lstep1,2
c
      call xctilr(pbavg(  1-nbdy,1-nbdy,1  ),1, 3, 6,6, halo_ps)
      call xctilr(ubavg(  1-nbdy,1-nbdy,1  ),1, 3, 6,6, halo_uv)
      call xctilr(vbavg(  1-nbdy,1-nbdy,1  ),1, 3, 6,6, halo_vv)
c
             
      if     (lpipe .and. lpipe_barotp) then
        call pipe_compare_sym1(
     &    pbavg(1-nbdy,1-nbdy,nl),ip,'barot+:pbavn')
        call pipe_compare_sym2(
     &    ubavg(1-nbdy,1-nbdy,nl),iu,'barot+:ubavn',
     &    vbavg(1-nbdy,1-nbdy,nl),iv,'barot+:vbavn')
        call pipe_compare_sym1(
     &    pbavg(1-nbdy,1-nbdy,ml),ip,'barot+:pbavm')
        call pipe_compare_sym2(
     &    ubavg(1-nbdy,1-nbdy,ml),iu,'barot+:ubavm',
     &    vbavg(1-nbdy,1-nbdy,ml),iv,'barot+:vbavm')
      endif
c
c --- odd minor time step.
c
      ml=n
      nl=3
c
c --- continuity equation
c
c --- rhs: pbavg, ubavg+, vbavg+
c --- lhs: pbavg
c
      margin = mbdy - 1
c
!$OMP PARALLEL DO PRIVATE(j,l,i,pbudel,pbvdel)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do l=1,isp(j)
          do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l))
#if defined(STOKES)
C
C   Barotropic Stokes flow included here
C
            pbudel = (ubavg(i+1,j,ml)+usdbavg(i+1,j))*
     &                    (depthu(i+1,j)*scuy(i+1,j))
     &              -(ubavg(i,  j,ml)+usdbavg(i,  j))*
     &                    (depthu(i,  j)*scuy(i,  j))
            pbvdel = (vbavg(i,j+1,ml)+vsdbavg(i,j+1))*
     &                    (depthv(i,j+1)*scvx(i,j+1))
     &              -(vbavg(i,j,  ml)+vsdbavg(i,j  ))*
     &                    (depthv(i,j  )*scvx(i,j  ))
#else
            pbudel =  ubavg(i+1,j,ml)*(depthu(i+1,j)*scuy(i+1,j))
     &               -ubavg(i  ,j,ml)*(depthu(i  ,j)*scuy(i  ,j))
            pbvdel =  vbavg(i,j+1,ml)*(depthv(i,j+1)*scvx(i,j+1))
     &               -vbavg(i,j  ,ml)*(depthv(i,j  )*scvx(i,j  ))
#endif
            pbavg(i,j,nl)=
     &        ((1.-wblpf)*pbavg(i,j,ml)+
     &             wblpf *pbavg(i,j,nl) )-
     &         (1.+wblpf)*dlt*(pbudel + pbvdel)*scp2i(i,j)
          enddo
        enddo
      enddo
c
      mn=ml
c
c --- u momentum equation, 1st
c
c --- rhs: pbavg+, vbavg+, pvtrop+
c --- lhs: ubavg
c
      margin = margin - 1
c
!$OMP PARALLEL DO PRIVATE(j,l,i,utndcy)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do l=1,isu(j)
          do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l))
            utndcy=-thref*(pbavg(i,j,nl)-pbavg(i-1,j,nl))*scuxi(i,j)+
     &       ((vbavg(i  ,j,  mn)*depthv(i  ,j)
     &        +vbavg(i  ,j+1,mn)*depthv(i  ,j+1))+
     &        (vbavg(i-1,j,  mn)*depthv(i-1,j)
     &        +vbavg(i-1,j+1,mn)*depthv(i-1,j+1)))*
     &       (0.125*(pvtrop(i,j)+pvtrop(i,j+1)))
c
            ubavg(i,j,nl)=
     &        ((1.-wblpf)*ubavg(i,j,ml)+
     &             wblpf *ubavg(i,j,nl))+
     &         (1.+wblpf)*dlt*(utndcy+utotn(i,j))
c
*           if (ldebug_barotp .and. i.eq.itest.and.j.eq.jtest) then
*             write (lp,'(i9,2i5,i3,3x,a,5f7.3)')
*    &          nstep,i+i0,j+j0,lll,
*    &          'u_old,u_new,p_grad,corio,u_star =',
*    &          ubavg(i,j,ml),ubavg(i,j,nl),
*    &           -thref*(pbavg(i,j,nl)-pbavg(i-1,j,nl))*scuxi(i,j)*dlt,
*    &          (vbavg(i  ,j,  mn)*depthv(i  ,j)
*    &          +vbavg(i  ,j+1,mn)*depthv(i  ,j+1)
*    &          +vbavg(i-1,j,  mn)*depthv(i-1,j)
*    &          +vbavg(i-1,j+1,mn)*depthv(i-1,j+1))
*    &          *(pvtrop(i,j)+pvtrop(i,j+1))
*    &          *.125 * dlt,utotn(i,j) * dlt
*           endif
          enddo
        enddo
      enddo
c
      mn = nl
c
c --- v momentum equation, 2nd
c --- rhs: pbavg+, ubavg+, pvtrop+
c --- lhs: vbavg
c
      margin = margin - 1
c
!$OMP PARALLEL DO PRIVATE(j,l,i,vtndcy)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do l=1,isv(j)
          do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l))
            vtndcy=-thref*(pbavg(i,j,nl)-pbavg(i,j-1,nl))*scvyi(i,j)-
     &       ((ubavg(i,  j  ,mn)*depthu(i,  j  )
     &        +ubavg(i+1,j  ,mn)*depthu(i+1,j  ))+
     &        (ubavg(i,  j-1,mn)*depthu(i,  j-1)
     &        +ubavg(i+1,j-1,mn)*depthu(i+1,j-1)))*
     &       (0.125*(pvtrop(i,j)+pvtrop(i+1,j)))
c
            vbavg(i,j,nl)=
     &        ((1.-wblpf)*vbavg(i,j,ml)+
     &             wblpf *vbavg(i,j,nl))+
     &         (1.+wblpf)*dlt*(vtndcy+vtotn(i,j))
c
*           if (ldebug_barotp .and. i.eq.itest.and.j.eq.jtest) then
*             write (lp,'(i9,2i5,i3,3x,a,5f7.3)')
*    &          nstep,i+i0,j+j0,lll,
*    &          'v_old,v_new,p_grad,corio,v_star =',
*    &          vbavg(i,j,ml),vbavg(i,j,nl),
*    &          -thref*(pbavg(i,j,nl)-pbavg(i,j-1,nl))*scvyi(i,j)*dlt,
*    &          -(ubavg(i,  j  ,mn)*depthu(i,j  )
*    &           +ubavg(i+1,j  ,mn)*depthu(i+1,j  )
*    &           +ubavg(i,  j-1,mn)*depthu(i,j-1)
*    &           +ubavg(i+1,j-1,mn)*depthu(i+1,j-1))
*    &          *(pvtrop(i,j)+pvtrop(i+1,j))
*    &          *.125 * dlt, vtotn(i,j) * dlt
*           endif
          enddo
        enddo
      enddo
c
      if     (ldebug_barotp) then
        call xcsync(flush_lp)
      endif
c
      if     (lbflag.eq.1) then
        call latbdp(nl)
      elseif (lbflag.eq.2) then
        call latbdt(nl,lll)
      elseif (lbflag.eq.3) then
        call latbdf(nl,lll)
      endif
c
      if     (lpipe .and. lpipe_barotp) then
        call pipe_compare_sym1(
     &    pbavg(1-nbdy,1-nbdy,nl),ip,'barot+:pbavn')
        call pipe_compare_sym2(
     &    ubavg(1-nbdy,1-nbdy,nl),iu,'barot+:ubavn',
     &    vbavg(1-nbdy,1-nbdy,nl),iv,'barot+:vbavn')
        call pipe_compare_sym1(
     &    pbavg(1-nbdy,1-nbdy,ml),ip,'barot+:pbavm')
        call pipe_compare_sym2(
     &    ubavg(1-nbdy,1-nbdy,ml),iu,'barot+:ubavm',
     &    vbavg(1-nbdy,1-nbdy,ml),iv,'barot+:vbavm')
      endif
c
c --- even minor time step.
c
      ml=3
      nl=n
      wblpf = wbaro  !used for all subsequent time steps: lll=2,lstep1
c
c --- continuity equation
c
c --- rhs: pbavg, ubavg+, vbavg+
c --- lhs: pbavg
c
      margin = mbdy - 1
c
!$OMP PARALLEL DO PRIVATE(j,l,i,pbudel,pbvdel)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do l=1,isp(j)
          do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l))
C=================================================================
CDAN  Barotropic Stokes Drift Velocities included here
CDAN
CDAN  August 2013
CDAN_______________________________________________________________
#if defined(STOKES)
C
C   Barotropic Stokes flow included here
C
      pbudel=(ubavg(i+1,j,ml)+usdbavg(i+1,j))*
     &                            (depthu(i+1,j)*scuy(i+1,j))
     & -(ubavg(i,j,ml)+usdbavg(i,j))*(depthu(i ,j)*scuy(i,j))
c 

c            pbvdel =  vbavg(i,j+1,ml)*(depthv(i,j+1)*scvx(i,j+1))
c     &               -vbavg(i,j  ,ml)*(depthv(i,j  )*scvx(i,j  ))     
      pbvdel=(vbavg(i,j+1,ml)+vsdbavg(i,j+1))*
     &                       (depthv(i,j+1)*scvx(i,j+1))
     & -(vbavg(i,j,ml)+vsdbavg(i,j))*(depthv(i,j)*scvx(i,j))
#else
            pbudel =  ubavg(i+1,j,ml)*(depthu(i+1,j)*scuy(i+1,j))
     &               -ubavg(i  ,j,ml)*(depthu(i  ,j)*scuy(i  ,j))
            pbvdel =  vbavg(i,j+1,ml)*(depthv(i,j+1)*scvx(i,j+1))
     &               -vbavg(i,j  ,ml)*(depthv(i,j  )*scvx(i,j  ))
  
#endif
            pbavg(i,j,nl)=
     &        ((1.-wblpf)*pbavg(i,j,ml)+
     &             wblpf *pbavg(i,j,nl) )-
     &         (1.+wblpf)*dlt*(pbudel + pbvdel)*scp2i(i,j)
          enddo
        enddo
      enddo
c
      mn=ml
c
c --- v momentum equation, 1st
c
c --- rhs: pbavg+, ubavg+, pvtrop+
c --- lhs: vbavg
c
      margin = margin - 1
c
!$OMP PARALLEL DO PRIVATE(j,l,i,vtndcy)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do l=1,isv(j)
          do i=max(1-margin,ifv(j,l)),min(ii+margin,ilv(j,l))
            vtndcy=-thref*(pbavg(i,j,nl)-pbavg(i,j-1,nl))*scvyi(i,j)-
     &       ((ubavg(i,  j  ,mn)*depthu(i,  j  )
     &        +ubavg(i+1,j  ,mn)*depthu(i+1,j  ))+
     &        (ubavg(i,  j-1,mn)*depthu(i,  j-1)
     &        +ubavg(i+1,j-1,mn)*depthu(i+1,j-1)))*
     &       (0.125*(pvtrop(i,j)+pvtrop(i+1,j)))
c
            vbavg(i,j,nl)=
     &        ((1.-wblpf)*vbavg(i,j,ml)+
     &             wblpf *vbavg(i,j,nl))+
     &         (1.+wblpf)*dlt*(vtndcy+vtotn(i,j))
c
*           if (ldebug_barotp .and. i.eq.itest.and.j.eq.jtest) then
*             write (lp,'(i9,2i5,i3,3x,a,5f7.3)')
*    &          nstep,i+i0,j+j0,lll+1,
*    &          'v_old,v_new,p_grad,corio,v_star =',
*    &          vbavg(i,j,ml),vbavg(i,j,nl),
*    &          -thref*(pbavg(i,j,nl)-pbavg(i,j-1,nl))*scvyi(i,j)*dlt,
*    &          -(ubavg(i,  j  ,mn)*depthu(i,j  )
*    &           +ubavg(i+1,j  ,mn)*depthu(i+1,j  )
*    &           +ubavg(i,  j-1,mn)*depthu(i,j-1)
*    &           +ubavg(i+1,j-1,mn)*depthu(i+1,j-1))
*    &          *(pvtrop(i,j)+pvtrop(i+1,j))
*    &          *.125 * dlt, vtotn(i,j) * dlt
*           endif
          enddo
        enddo
      enddo
c
      mn=nl
c
c --- u momentum equation, 2nd
c
c --- rhs: pbavg+, vbavg+, pvtrop+
c --- lhs: ubavg
c
      margin = margin - 1
c
!$OMP PARALLEL DO PRIVATE(j,l,i,utndcy)
!$OMP&         SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do l=1,isu(j)
          do i=max(1-margin,ifu(j,l)),min(ii+margin,ilu(j,l))
            utndcy=-thref*(pbavg(i,j,nl)-pbavg(i-1,j,nl))*scuxi(i,j)+
     &       ((vbavg(i  ,j,  mn)*depthv(i  ,j)
     &        +vbavg(i  ,j+1,mn)*depthv(i  ,j+1))+
     &        (vbavg(i-1,j,  mn)*depthv(i-1,j)
     &        +vbavg(i-1,j+1,mn)*depthv(i-1,j+1)))*
     &       (0.125*(pvtrop(i,j)+pvtrop(i,j+1)))
c
            ubavg(i,j,nl)=
     &        ((1.-wblpf)*ubavg(i,j,ml)+
     &             wblpf *ubavg(i,j,nl))+
     &         (1.+wblpf)*dlt*(utndcy+utotn(i,j))
c
*           if (ldebug_barotp .and. i.eq.itest.and.j.eq.jtest) then
*             write (lp,'(i9,2i5,i3,3x,a,5f7.3)')
*    &          nstep,i+i0,j+j0,lll+1,
*    &          'u_old,u_new,p_grad,corio,u_star =',
*    &          ubavg(i,j,ml),ubavg(i,j,nl),
*    &          -thref*(pbavg(i,j,nl)-pbavg(i-1,j,nl))*scuxi(i,j)*dlt,
*    &          (vbavg(i  ,j,  mn)*depthv(i  ,j)
*    &          +vbavg(i  ,j+1,mn)*depthv(i  ,j+1)
*    &          +vbavg(i-1,j,  mn)*depthv(i-1,j)
*    &          +vbavg(i-1,j+1,mn)*depthv(i-1,j+1))
*    &          *(pvtrop(i,j)+pvtrop(i,j+1))
*    &          *.125 * dlt,utotn(i,j) * dlt
*           endif
          enddo
        enddo
      enddo
c
      if     (ldebug_barotp) then
        call xcsync(flush_lp)
      endif
c
      if     (lbflag.eq.1) then
        call latbdp(nl)
      elseif (lbflag.eq.2) then
        call latbdt(nl,lll+1)
      elseif (lbflag.eq.3) then
        call latbdf(nl,lll+1)
      endif
c
 840  continue  ! lll=1,lstep1,2
c
      if     (btrlfr .and. delt1.ne.baclin) then  !not on very 1st time step
c ---   Asselin filtering
!$OMP   PARALLEL DO PRIVATE(j,i)
!$OMP&           SCHEDULE(STATIC,jblk)
        do j=1,jj
          do i=1,ii
            ubavg(i,j,m) = ubaros(i,j)+wts2*ubavg(i,j,n)
            vbavg(i,j,m) = vbaros(i,j)+wts2*vbavg(i,j,n)
            pbavg(i,j,m) = pbaros(i,j)+wts2*pbavg(i,j,n)
          enddo !i
        enddo !j
      endif !btrlfr & not on very 1st time step
c
      if     (lbflag.eq.1) then
c
c ---   correct mean height.
c ---   this should not be required - so there may be a bug in the bc.
c
!$OMP   PARALLEL DO PRIVATE(j,l,i)
!$OMP&         SCHEDULE(STATIC,jblk)
        do j=1,jj
          do l=1,isp(j)
            do i=max(1,ifp(j,l)),min(ii,ilp(j,l))
              util1(i,j) = pbavg(i,j,nl)*scp2(i,j)
            enddo
          enddo
        enddo
        call xcsum(sump, util1,ip)
        q = sump/area
c
c ---   rhs: pbavg
c ---   lhs: pbavg
c
        margin = 0
c
!$OMP   PARALLEL DO PRIVATE(j,l,i)
!$OMP&         SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do l=1,isp(j)
            do i=max(1-margin,ifp(j,l)),min(ii+margin,ilp(j,l))
              pbavg(i,j,1) = pbavg(i,j,1) - q
              pbavg(i,j,2) = pbavg(i,j,2) - q
              pbavg(i,j,3) = pbavg(i,j,3) - q
            enddo
          enddo
        enddo
      endif
      if     (lpipe .and. lpipe_barotp) then
        call pipe_compare(pbavg(1-nbdy,1-nbdy,1), ip,'barotp:pbav1')
        call pipe_compare(pbavg(1-nbdy,1-nbdy,2), ip,'barotp:pbav2')
        call pipe_compare(pbavg(1-nbdy,1-nbdy,3), ip,'barotp:pbav3')
      endif
c
      return
      end subroutine barotp
c
c
c> Revision history:
c>
c> Mar. 1995 - changed vertical velocity averaging interval from 10 cm to 1 m
c>             (loops 33,35)
c> Mar. 1995 - changed order of loop nesting in loop 842
c> July 1997 - eliminated 3-D arrays -uold,vold- (used in time smoothing)
c> Aug. 1997 - transferred loops preceding loop 840 to momeq2.f
c> Jan. 2000 - added latbdp for lateral boundary ports
c> Aug. 2001 - two barotropic time steps per loop, for halo efficiency
c> Nov. 2006 - added lbflag==3 (latbdf) and thref_bt (mod_tides)
c> Nov. 2006 - removed thref_bt (and mod_tides)
c> Apr. 2007 - added btrlfr: leapfrog time step; see also momtum
c> Apr. 2010 - bugfixes for 1st time step and 1st miner time step
c> Apr. 2010 - added Asselin filtering to btrlfr