#include "options.inc"

c----------------------------------------------------------------
c
c     BBL module for spflame
c 
c     this module implements the bottom boundary layer 
c     parameterization by Beckmann and Doescher, JPO, 97
c
c     it is linked in the main code by calls of 
c     bbl_init in setup.F, adv_vel_bbl at the
c     end of adv_vel.F and of add_bbl at the end 
c     of the tracer loop in tracer.F
c     (see also advective operators in tracer.F)
c
c     Mar 2001, C.Eden
c
c     written for MOM2.1 by:   Ralf Doescher
c     date:                    Dec, 1997
c
c     based on:
c
c !---------------------------------------------------------------
c |              Beckmann and Doescher, 1997:                    |
c |    A method for improved representation of dense water       |
c | spreading over topography in geopotential-coordinate models. |
c |      J. Phys. Oceanogr., 27, 4, 1997, 581 - 591.             |
c !---------------------------------------------------------------
c
c     FOR NON-COMMERCIAL USE ONLY
c
c----------------------------------------------------------------


      module bbl_module
      implicit none

      real, allocatable    :: h_sigma(:,:)
      real, allocatable    :: rhodiffxp(:,:),rhodiffyp(:,:)
      real, allocatable    :: ahc_north_sigma(:),ahc_south_sigma(:)

      real, allocatable    :: sigma_weight(:,:,:)
      integer, allocatable :: kt_upper_fac(:,:,:)

      real, allocatable    :: adv_vnt_s(:,:)
      real, allocatable    :: adv_vet_s(:,:)
      real, allocatable    :: adv_vbt_su(:,:,:)

      integer, allocatable :: ipdown(:,:),jpdown(:,:)
      integer, allocatable :: ifac(:,:,:),jfac(:,:,:)
      integer, allocatable :: im1fac(:,:,:),jm1fac(:,:,:)

      end module bbl_module


      subroutine bbl_init
c----------------------------------------------------------------
c     do some initializations for BBL
c----------------------------------------------------------------
      use spflame_module
      use bbl_module
      implicit none
      integer is,ie,js,je
      integer i,j,k,jm1,jp1
      integer kk,kkm1,kanf,ifacc,im1facc,jfacc,jm1facc,ibottom
      integer ijall
      integer merkupper(is_pe:ie_pe)

      if (my_pe == 0) then
       print*,''
       print*,' Initialization of BBL parameterization'
       print*,''
       print*,' ah_sigma               = ',ah_sigma
       print*,' ah_sigma_min           = ',ah_sigma_min
       print*,' enable adv.  transport = ',enable_bbl_advection
       print*,' enable diff. transport = ',enable_bbl_diffusion
       print*,''
       if (enable_bbl_only_north) then
        print*,' Note:'
        print*,' BBL is working only north of 45N and south of 45S'
       endif
      endif

#ifdef partial_cell
      if (my_pe == 0) then
       print*,''
       print*,'ERROR: option partial cell is incompatible with BBL'
       print*,''
      endif
      call halt_stop(' in bbl_init')
#endif
c
c     allocate momory for the working arrays
c
      is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

      allocate( h_sigma(is-1:ie+1,js-1:je+1) )
      h_sigma=0.
      allocate( rhodiffxp(is-1:ie,js-1:je))
      rhodiffxp=0.
      allocate( rhodiffyp(is-1:ie,js-1:je))
      rhodiffyp=0.
      allocate( ahc_north_sigma(js-1:je+1))
      ahc_north_sigma=0.
      allocate( ahc_south_sigma(js-1:je+1))
      ahc_south_sigma=0.

      if (enable_bbl_advection) then

       allocate( sigma_weight(is:ie,km,js:je) )
       sigma_weight=0.
       allocate( kt_upper_fac(is:ie,km,js:je) )
       kt_upper_fac=0
       allocate(adv_vnt_s(is-1:ie+1,      js-1:je+1) )
       adv_vnt_s=0.
       allocate(adv_vet_s(is-1:ie+1,      js-1:je+1) )
       adv_vet_s=0.
       allocate(adv_vbt_su(is-1:ie+1,0:km,js-1:je+1) )
       adv_vbt_su=0.

       allocate(ipdown(is-1:ie+1,js-1:je+1) )
       ipdown=1
       allocate(jpdown(is-1:ie+1,js-1:je+1) )
       jpdown=1

       allocate(ifac(is-1:ie+1,km,js-1:je+1) )
       ifac=1
       allocate(jfac(is-1:ie+1,km,js-1:je+1) )
       jfac=1
       allocate(im1fac(is-1:ie+1,km,js-1:je+1) )
       im1fac=1
       allocate(jm1fac(is-1:ie+1,km,js-1:je+1) )
       jm1fac=1

      endif
c
c     set up things which are done only once in the run
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      h_sigma=0.
      do j=js-1,je+1
       do i=is-1,ie+1
        if(kmt_big(i,j).ne.0) h_sigma(i,j) = dzt(kmt_big(i,j))
       enddo
      enddo

      do j=js-1,je+1
        jm1 = max(1,j-1)
        jp1 = min(jmt,j+1)
        ahc_north_sigma(j) = ah_sigma*csu(j)*dyur(j)
     &                           *cstr(j)*dytr(j)
        ahc_south_sigma(j) = ah_sigma*csu(jm1)*dyur(jm1)
     &                           *cstr(j)*dytr(j)
      enddo

      if (enable_bbl_advection) then
c
c      define weight factor for sigma advection
c
       sigma_weight(:,:,:) = 0.0
 
       do j=js,je
         do i=is,ie
           if(kmt_big(i,j).ne.0) then
            kk   = kmt_big(i,j)
            kkm1 = kmt_big(i,j)-1
            sigma_weight(i,kk  ,j)   = 1.0    ! big Ts
            sigma_weight(i,kkm1,j)   = 1.0    ! first small ts
            kanf=max(1,min(kmt_big(i,j),kmt_big(i,j+1),kmt_big(i,j-1),
     &                                 kmt_big(i+1,j),kmt_big(i-1,j)))
!           second small t:
            sigma_weight(i,kanf:kmt_big(i,j),j) =   1.0
           endif
         enddo
       enddo

       do j=js,je
         do i=is,ie
          merkupper(i)=0
         enddo
         do k=1,km
          do i=is,ie

           ifacc   = 1; im1facc = 1; jfacc   = 1
           jm1facc = 1; ibottom = 1
           kt_upper_fac(i,k,j) = 0
           if(k.ge.kmt_big(i+1,j)) ifacc   = 0
           if(k.ge.kmt_big(i-1,j)) im1facc = 0
           if(k.ge.kmt_big(i,j+1)) jfacc   = 0
           if(k.ge.kmt_big(i,j-1)) jm1facc = 0
           if(k.ge.kmt_big(i,j)) then
              ifacc   = 0; im1facc = 0
              jfacc   = 0; jm1facc = 0; ibottom = 0
           endif

           if(merkupper(i).eq.0) then
            ijall = ifacc*im1facc*jfacc*jm1facc
            if(ijall.eq.0) then
             if( 0.eq.kmt_big(i+1,j  ).or.
     &           0.eq.kmt_big(i-1,j  ).or.
     &           0.eq.kmt_big(i  ,j-1).or.
     &           0.eq.kmt_big(i  ,j+1) ) then
              merkupper(i) = 1
             endif
            endif
           endif

           if(merkupper(i).eq.0) then

            ijall = ifacc*im1facc*jfacc*jm1facc
            if(ijall.eq.0) then
             if( k.eq.kmt_big(i+1,j  ).or.
     &           k.eq.kmt_big(i-1,j  ).or.
     &           k.eq.kmt_big(i  ,j-1).or.
     &           k.eq.kmt_big(i  ,j+1) ) then
              kt_upper_fac(i,k,j) = 1
              merkupper(i) = 1
             endif

             if( k.le.kmt_big(i+1,j )-1.and.
     &           k.le.kmt_big(i-1,j )-1.and.
     &           k.le.kmt_big(i  ,j-1)-1.and.
     &           k.le.kmt_big(i  ,j+1)-1 ) then
              kt_upper_fac(i,k,j) = 1
              merkupper(i) = 1
             endif
            endif

            if( k.eq.kmt_big(i  ,j  )-1.and.
     &          k.le.kmt_big(i+1,j  )-1.and.
     &          k.le.kmt_big(i-1,j  )-1.and.
     &          k.le.kmt_big(i  ,j-1)-1.and.
     &          k.le.kmt_big(i  ,j+1)-1 ) then
              kt_upper_fac(i,k,j) = 1
              merkupper(i) = 1
            endif

            if( k.eq.kmt_big(i+1,j  )-1.and.
     &          k.eq.kmt_big(i-1,j  )-1.and.
     &          k.eq.kmt_big(i  ,j-1)-1.and.
     &          k.eq.kmt_big(i  ,j+1)-1 ) then
              kt_upper_fac(i,k,j) = 1
              merkupper(i) = 1
            endif

           endif
         
          enddo
         enddo
       enddo

      endif

      if (my_pe == 0) print*,' done'

      end subroutine bbl_init


      subroutine adv_vel_bbl
c
c=======================================================================
c     calculate advection velocities for tracer equations     
c     compute advective terms in BBL formulation
c=======================================================================
c
      use spflame_module
      use bbl_module
      implicit none
      integer :: is,ie,js,je
      integer :: i,j,k,kn,klo,kk

#ifdef vector_host
      integer :: kvec(is_pe-1:ie_pe)
      real    :: tvec(is_pe-1:ie_pe)
      real    :: svec(is_pe-1:ie_pe)
      real    :: rvec(is_pe-1:ie_pe)
      real    :: tpvec(is_pe-1:ie_pe)
      real    :: spvec(is_pe-1:ie_pe)
      real    :: rpvec(is_pe-1:ie_pe)
#ifdef partial_cell
      real    :: ztpvec(is_pe-1:ie_pe)
#endif
#endif

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

c     compute cross-slope density gradients for a sigma diffusion condition

#ifdef vector_host
c
c     vector version of the code
c
      kvec=1;tvec=0.;svec=0.;spvec=0.;tpvec=0.;rvec=0.;rpvec=0.

      do j=js-1,je

       do i=is-1,ie
        k=kmt_big(i,j)
        if ( k /= 0 ) then
         kn=kmt_big(i+1,j)
         if ( kn > 1 ) then
          kvec(i) = max(k,kn)
#ifdef partial_cell
          ztpvec(i) = ztp(i,kvec(i),j)
#endif
          tvec(i) = t(i  ,k ,j,1,tau)
          svec(i) = t(i  ,k ,j,2,tau)
          tpvec(i)= t(i+1,kn,j,1,tau)
          spvec(i)= t(i+1,kn,j,2,tau)
         endif
        endif 
       enddo
       call model_dens_kvec(tvec ,svec ,rvec ,kvec,ie-is+2
#ifdef partial_cell
     &                      ,ztpvec
#endif
     &    )
       call model_dens_kvec(tpvec,spvec,rpvec,kvec,ie-is+2
#ifdef partial_cell
     &                      ,ztpvec
#endif
     &    )
       do i=is-1,ie
        rhodiffxp(i,j)=rvec(i)-rpvec(i)
       enddo

       do i=is-1,ie
        k=kmt_big(i,j)
        if ( k /= 0 ) then
         kn=kmt_big(i,j+1)
         if ( kn > 1 ) then
          kvec(i) = max(k,kn)
#ifdef partial_cell
          ztpvec(i) = ztp(i,kvec(i),j)
#endif
          tvec(i) = t(i,k ,j  ,1,tau)
          svec(i) = t(i,k ,j  ,2,tau)
          tpvec(i)= t(i,kn,j+1,1,tau)
          spvec(i)= t(i,kn,j+1,2,tau)
         endif
        endif 
       enddo
       call model_dens_kvec(tvec ,svec ,rvec ,kvec,ie-is+2 
#ifdef partial_cell
     &                      ,ztpvec
#endif
     &    )
       call model_dens_kvec(tpvec,spvec,rpvec,kvec,ie-is+2 
#ifdef partial_cell
     &                      ,ztpvec
#endif
     &    )
       do i=is-1,ie
        rhodiffyp(i,j)=rvec(i)-rpvec(i)
       enddo

      enddo

#else
c
c     scalar version of the code
c
      do j=js-1,je
          do i=is-1,ie
            k=kmt_big(i,j)
            if ( k /= 0 ) then

              kn=kmt_big(i+1,j)
              if ( kn > 1 ) then
                klo=max(k,kn)
                rhodiffxp(i,j)=
     &    model_dens_scalar(t(i  ,k ,j,1,tau),t(i  ,k ,j,2,tau),klo
#ifdef partial_cell
     &                       ,ztp(i,klo,j)
#endif
     &                       )
     &   -model_dens_scalar(t(i+1,kn,j,1,tau),t(i+1,kn,j,2,tau),klo
#ifdef partial_cell
     &                       ,ztp(i+1,klo,j)
#endif
     &                       )
              endif

              kn=kmt_big(i,j+1)
              if ( kn > 1 ) then
                klo=max(k,kn)
                rhodiffyp(i,j)=
     &    model_dens_scalar(t(i,k ,j  ,1,tau),t(i,k ,j  ,2,tau),klo
#ifdef partial_cell
     &                       ,ztp(i,klo,j)
#endif
     &                       )
     &   -model_dens_scalar(t(i,kn,j+1,1,tau),t(i,kn,j+1,2,tau),klo
#ifdef partial_cell
     &                       ,ztp(i,klo,j+1)
#endif
     &                       )
              endif

            endif 
          enddo
      enddo
#endif


       if (enable_bbl_only_north) then
        do j=js-1,je
c          if (j .lt. 828) then
          if ( abs(yt(j)) < 45.1 ) then
            rhodiffyp(i,j)=0.
            rhodiffxp(i,j)=0.
          endif
        enddo
       endif


c     the advective velocity adv_vnt is located at the northern
c     interface of the ij-box. Consequently, advective velocities 
c     in the bbl (adv_vnt_s) are also infered from the northern 
c     interface.
c     Viewed from the ij-position, adv_vnt is multiplied by 
c     dzt(upper)/dzt(lower) only if the adjacent northern neigbor 
c     bottombox is located above the level of the ij-bottombox.
c     Consequential treatment in ADV_Ty_s (in fdift.h) leads
c     to mass conservation.
c     In short: 
c     v1 *  interface_area_upper =
c     v1 * (interface_area_upper/interface_area_lower) 
c        * interface_area_lower

      if (enable_bbl_advection) then
c
c     sigma_advection velocity on northern face of "T" cells.

       do j=js-1,je
        do i=is-1,ie
          adv_vnt_s(i,j) = 0.
          if (kmt_big(i,j+1)>0 .and. kmt_big(i,j)>0) then
            adv_vnt_s(i,j)  = adv_vnt(i,kmt_big(i,j),j)
            if (kmt_big(i,j) > kmt_big(i,j+1)) then
              adv_vnt_s(i,j) = adv_vnt(i,kmt_big(i,j+1),j)
     &            *dzt(kmt_big(i,j+1))/dzt(kmt_big(i,j))
            endif
          endif
        enddo
       end do

       call set_cyclic(adv_vnt_s,1,1)
       if (my_blk_i == 1 .and. enable_obc_west) 
     &              adv_vnt_s(1,:) = adv_vnt_s(2,:)
       if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &              adv_vnt_s(imt,:) = adv_vnt_s(imt-1,:)
c
c     sigma_advection velocity on eastern face of "T" cells.
c
       do j=js-1,je
        do i=is-1,ie
          adv_vet_s(i,j) = 0.
          if(kmt_big(i+1,j)>0.and.kmt_big(i,j)>0) then
            adv_vet_s(i,j) = adv_vet(i,kmt_big(i,j),j)
            if(kmt_big(i,j) > kmt_big(i+1,j)) then
              adv_vet_s(i,j) = adv_vet(i,kmt_big(i+1,j),j) 
     &           * dzt(kmt_big(i+1,j)) / dzt(kmt_big(i,j)) 
            endif
          endif
        enddo
       enddo

       call set_cyclic(adv_vet_s,1,1)
       if (my_blk_i == 1 .and. enable_obc_west) 
     &              adv_vet_s(1,:) = adv_vet_s(2,:)
       if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &              adv_vet_s(imt,:) = adv_vet_s(imt-1,:)
c
c     sigma_advection vertical velocity on upper face of bottom cells.
c      construct divergence of advection velocity * level thickness
c
       do j=js-1,je
          do i=is-1,ie

            jpdown(i,j)=1
            ipdown(i,j)=1
            k=kmt_big(i,j)
            if (k > kmt_big(i,j+1)) then         ! flag downhill flow
              jpdown(i,j) = max(0.,-1.*sign(1.,adv_vnt_s(i,j)))
            else if (k<kmt_big(i,j+1)) then
              jpdown(i,j) = max(0.,+1.*sign(1.,adv_vnt_s(i,j)))
            endif
            if (jpdown(i,j) == 1) then   ! cross-check with density difference
              if ((k>kmt_big(i,j+1).and.rhodiffyp(i,j).lt.0.).or.
     &            (k<kmt_big(i,j+1).and.rhodiffyp(i,j).gt.0.)) then
                continue                        ! denser water upstream
              else
                jpdown(i,j)=0
              endif
            endif

            if (k .gt. kmt_big(i+1,j)) then
              ipdown(i,j) = max(0.,-1.*sign(1.,adv_vet_s(i,j)))
            else if (k.lt.kmt_big(i+1,j)) then
              ipdown(i,j) = max(0.,+1.*sign(1.,adv_vet_s(i,j)))
            endif
            if (ipdown(i,j) .eq. 1) then
              if ((k.gt.kmt_big(i+1,j).and.rhodiffxp(i,j).lt.0.).or.
     &            (k.lt.kmt_big(i+1,j).and.rhodiffxp(i,j).gt.0.)) then
                continue
              else
                ipdown(i,j)=0
              endif
            endif

          enddo
       enddo

       call set_cyclic_int(ipdown,1,1)
       call set_cyclic_int(jpdown,1,1)
       if (my_blk_i == 1 .and. enable_obc_west) then
             ipdown(1,:) = ipdown(2,:)
             jpdown(1,:) = jpdown(2,:)
       endif
       if (my_blk_i == n_pes_i .and. enable_obc_east)  then
             ipdown(imt,:) = ipdown(imt-1,:)
             jpdown(imt,:) = jpdown(imt-1,:)
       endif
c
c     construct vertical velocity on the bottom face of outer
c     "T" cells adjacent to the sigma bottom layer
c
       do j=js,je
c
c       determine masks for blending out "redirected" downslope flow:
c
        do k=1,km                ! initialize whole water column with 1
          do i=is,ie
            ifac(i,k,j)   = 1; im1fac(i,k,j) = 1
            jfac(i,k,j)   = 1; jm1fac(i,k,j) = 1
          enddo
        enddo

#ifdef vector_host
        do k=1,km
         do i=is,ie          ! set sides of topo-steps and bathymetry to 0
          if (k>=kmt_big(i+1,j)+1)    ifac(i,k,j) = 0
          if (k>=kmt_big(i-1,j)+1)  im1fac(i,k,j) = 0
          if (k>=kmt_big(i,j+1)+1)    jfac(i,k,j) = 0
          if (k>=kmt_big(i,j-1)+1)  jm1fac(i,k,j) = 0
         enddo
        enddo
#else
        do i=is,ie          ! set sides of topo-steps and bathymetry to 0
          do k=min(kmt_big(i+1,j)+1,km),km
            ifac(i,k,j)   = 0
          enddo
          do k=min(kmt_big(i-1,j)+1,km),km
            im1fac(i,k,j) = 0
          enddo
          do k=min(kmt_big(i,j+1)+1,km),km
            jfac(i,k,j)   = 0
          enddo
          do k=min(kmt_big(i,j-1)+1,km),km
            jm1fac(i,k,j) = 0
          enddo
        enddo
#endif

        do i=is,ie         ! determine sign at k-level of topo step
          k=kmt_big(i+1,j)
          if (k.gt.0) then
!           import (0) or export (1) ?
            ifac(i,k,j)  = max(0.,+1.*sign(1.,adv_vet(i,k,j)))
!           if import, set mask to 0 for dense downslope flow:
!              (adv_vet=0 -> ifac=1, which is correct, 
!                                    strange as it may seem j.d.)
!              (note: mask stays 1 if density diff. across step =0!)
            if (ifac(i,k,j).eq.0) then
              kk=kmt_big(i,j)
              if ((kk.gt.k .and. rhodiffxp(i,j).lt.0.) .or.
     &            (kk.lt.k .and. rhodiffxp(i,j).gt.0.)) then
                continue
              else
                ifac(i,k,j)=1
              endif
            endif
          endif
        enddo

        do i=is,ie
          k=kmt_big(i-1,j)
          if (k.gt.0) then
            im1fac(i,k,j)= max(0.,-1.*sign(1.,adv_vet(i-1,k,j)))
            if (im1fac(i,k,j).eq.0) then
              kk=kmt_big(i,j)
              if ((kk.gt.k .and. -1.0*rhodiffxp(i-1,j).lt.0.) .or.
     &            (kk.lt.k .and. -1.0*rhodiffxp(i-1,j).gt.0.)) then
                continue
              else
                im1fac(i,k,j)=1
              endif
            endif
          endif
        enddo

        do i=is,ie
          k=kmt_big(i,j+1)
          if (k.gt.0) then
            jfac(i,k,j)  = max(0.,+1.*sign(1.,adv_vnt(i,k,j)))
            if (jfac(i,k,j).eq.0) then
              kk=kmt_big(i,j)
              if ((kk.gt.k .and. rhodiffyp(i,j).lt.0.) .or.
     &            (kk.lt.k .and. rhodiffyp(i,j).gt.0.)) then
                continue
              else
                jfac(i,k,j)=1
             endif
            endif
          endif
        enddo

        do i=is,ie
          k=kmt_big(i,j-1)
          if (k.gt.0) then
            jm1fac(i,k,j)= max(0.,-1.*sign(1.,adv_vnt(i,k,j-1)))
            if (jm1fac(i,k,j).eq.0) then
              kk=kmt_big(i,j)
              if ((kk.gt.k .and. -1.0*rhodiffyp(i,j-1).lt.0.) .or.
     &            (kk.lt.k .and. -1.0*rhodiffyp(i,j-1).gt.0.)) then
                continue
              else
                jm1fac(i,k,j)=1
              endif
            endif
          endif
        enddo


#ifdef vector_host
        do k=1,km
         do i=is,ie                 ! set local bathymetry points to 0
          if (k>=kmt_big(i,j)) then
            ifac(i,k,j)   = 0; im1fac(i,k,j) = 0
            jfac(i,k,j)   = 0; jm1fac(i,k,j) = 0
          endif
         enddo
        enddo
#else
        do i=is,ie                 ! set local bathymetry points to 0
          do k=max(1,kmt_big(i,j)),km
            ifac(i,k,j)   = 0; im1fac(i,k,j) = 0
            jfac(i,k,j)   = 0; jm1fac(i,k,j) = 0
          enddo
        enddo
#endif

       enddo


       call set_cyclic_int(ifac,km,1)
       call set_cyclic_int(jfac,km,1)
       call set_cyclic_int(im1fac,km,1)
       call set_cyclic_int(jm1fac,km,1)

       if (my_blk_i == 1 .and. enable_obc_west)  then
              ifac(1,:,:) = ifac(2,:,:)
              jfac(1,:,:) = jfac(2,:,:)
              im1fac(1,:,:) = im1fac(2,:,:)
              jm1fac(1,:,:) = jm1fac(2,:,:)
       endif
       if (my_blk_i == n_pes_i .and. enable_obc_east)  then
              ifac(imt,:,:) = ifac(imt-1,:,:)
              jfac(imt,:,:) = jfac(imt-1,:,:)
              im1fac(imt,:,:) = im1fac(imt-1,:,:)
              jm1fac(imt,:,:) = jm1fac(imt-1,:,:)
       endif

       do j=js,je

c       construct divergence of advection velocity * level thickness

c       set "adv_vbt_su" at surface to 0.0 (rigid-lid) or dh/dt (free surf)

        do i=is,ie
          adv_vbt_su(i,0,j)   = 0.
        enddo

        do k=1,km
          do i=is,ie
            adv_vbt_su(i,k,j) =
     &      ((adv_vet(i,k,j)*ifac(i,k,j) - 
     &                           adv_vet(i-1,k,j)*im1fac(i,k,j))*dxtr(i)
     &      +(adv_vnt(i,k,j)*jfac(i,k,j) - 
     &                        adv_vnt(i,k,j-1)*jm1fac(i,k,j))*dytr(j)
     &                   )*cstr(j)*dzt(k)
          enddo
        enddo

c
c       integrate downward to define "adv_vbt" at the bottom of levels
c
        do k=1,km
          do i=is,ie
            adv_vbt_su(i,k,j) = adv_vbt_su(i,k  ,j)
     &                        + adv_vbt_su(i,k-1,j)
          enddo
        enddo

        do i=is,ie        ! integrate and delete afterwards?
          k=kmt_big(i,j)
          if(k.ge.2) then
            adv_vbt_su(i,k,j) = 0.
          endif
        enddo

       enddo


       call set_cyclic(adv_vbt_su,km+1,1)
       if (my_blk_i == 1 .and. enable_obc_west) 
     &              adv_vbt_su(1,:,:) = adv_vbt_su(2,:,:)
       if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &              adv_vbt_su(imt,:,:) = adv_vbt_su(imt-1,:,:)

      endif

      end subroutine adv_vel_bbl
      




      subroutine add_bbl(n,is,ie,js,je,adv_fe,adv_fb)
c=======================================================================
c     add what necessary for BBL to tracer tendency
c=======================================================================
      use spflame_module
      use bbl_module
      implicit none
      integer, intent(in) :: n,is,ie,js,je
      real, intent(in) ::  adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real, intent(in) ::  adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)

      integer i,j,k,kn,ibottom,kb
      integer kmt_jp1,kmt_jm1
      integer kmt_ip1,kmt_im1
      real factxp,factxm
      real factyp,factym
c     the following 2 arrays are needed for sigma_bottom_advection only
      real adv_fe_s(is_pe-1:ie_pe,js_pe:je_pe)  
      real adv_fb_su(is_pe:ie_pe,0:km,js_pe:je_pe) 
c
c     redefined operators
c
      real ADV_Ty_s
      ADV_Ty_s(i,j) = (
     &                ( 1.0 * 
     &         adv_vnt_s(i,j) * jpdown(i,j) *  
     &                         h_sigma(i,j) *
     &                        (t(i,kmt_jp1,j+1,n,tau) +
     &                         t(i,kmt(i,j ),j  ,n,tau)) 
     &       + adv_vnt(i,k,j) * (1-jpdown(i,j)) * 
     &                         h_sigma(i,j) *
     &                        (t(i,k,j+1,n,tau) +
     &                         t(i,k,j  ,n,tau)) )
     &              - (  1.0 * 
     &         adv_vnt_s(i,j-1) * jpdown(i,j-1) * 
     &                         h_sigma(i,j-1) *
     &                        (t(i,kmt(i,j  ),j  ,n,tau) +
     &                         t(i,kmt_jm1,j-1,n,tau))
     &     + adv_vnt(i,k,j-1) * (1-jpdown(i,j-1)) * 
     &                         h_sigma(i,j) *
     &                        (t(i,k,j  ,n,tau) +
     &                         t(i,k,j-1,n,tau)) )
     &                 ) * cstdyt2r(j) / h_sigma(i,j)

      real ADV_Tx_s
      ADV_Tx_s(i,j) = (
     &                ( h_sigma(i  ,j) * ipdown(i,j) *
     &         adv_fe_s(i,j) )
     &              + ( h_sigma(i  ,j) * (1-ipdown(i,j)) *
     &         adv_fe(i,k,j) )
     &              - ( h_sigma(i-1,j) * ipdown(i-1,j) *
     &         adv_fe_s(i-1,j) )
     &              - ( h_sigma(i  ,j) * (1-ipdown(i-1,j)) *
     &         adv_fe(i-1,k,j) )
     &                ) * cstdxt2r(i,j) / h_sigma(i,j)


      real ADV_Tz_s
      ADV_Tz_s(i,j) =  (adv_fb_su(i,k-1,j)  )
     &                                       *dzt2r(k)

      real ADV_Tx_su
      ADV_Tx_su(i,k,j) = 
     &  (adv_fe(i,k,j)*ifac(i,k,j) - adv_fe(i-1,k,j)*im1fac(i,k,j))
     &                   *cstdxt2r(i,j)

      real ADV_Ty_su
      ADV_Ty_su(i,k,j) = 
     &   (jfac(i,k,j)*adv_vnt(i,k,j)*(t(i,k,j,n,tau)+t(i,k,j+1,n,tau))
     &  - jm1fac(i,k,j)*adv_vnt(i,k,j-1)*
     &                               (t(i,k,j-1,n,tau)+t(i,k,j,n,tau))
     &                )*cstdyt2r(j)

      real ADV_Tz_su
      ADV_Tz_su(i,k,j) =ibottom*(
     &    adv_fb_su(i,k-1,j) * (1-kt_upper_fac(i,k,j)) 
     &  + adv_fb(i,k-1,j)    *    kt_upper_fac(i,k,j) 
     &  - adv_fb_su(i,k,j) )
     &               *dzt2r(k)


c     this array is needed for sigma_bottom_diffusion only
      real diff_fe_s(is_pe-1:ie_pe,js_pe:je_pe) 

      real DIFF_Tx_s
      DIFF_Tx_s(i,j) = (
     &   ( 0.5 * (h_sigma(i+1,j)+h_sigma(i,j)) *
     &   factxp*diff_fe_s(i  ,j)*tmask(i+1,kmt_ip1,j) )
     &     - ( 0.5 * (h_sigma(i-1,j)+h_sigma(i,j)) *
     &   factxm*diff_fe_s(i-1,j)*tmask(i-1,kmt_im1,j) )
     &   ) * cstdxtr(i,j) / h_sigma(i,j)

      real DIFF_Ty_s
      DIFF_Ty_s(i,j) =
     &   (
     &   (
     &   factyp*ahc_north_sigma(j) * tmask(i,kmt_jp1,j+1) *
     &   ( 0.5 * (h_sigma(i,j+1)+h_sigma(i,j)) ) *
     &   (t(i,kmt_jp1,j+1,n,taum1) - t(i,kmt(i,j),j,n,taum1))
     &   )
     & - (
     &   factym*ahc_south_sigma(j) * tmask(i,kmt_jm1,j-1) *
     &   ( 0.5 * (h_sigma(i,j)+h_sigma(i,j-1)) ) *
     &   (t(i,kmt(i,j),j,n,taum1) - t(i,kmt_jm1,j-1,n,taum1))
     &   )
     &   ) / h_sigma(i,j)



      if (enable_bbl_advection) then
c
c      calculate 2*advective flux in sigma bottom layer across
c      eastern face of "T" cells. north face component is built
c      into ADV_Ty_s
c
       do j=js,je
        do i=is-1,ie
            if(kmt_big(i+1,j).ne.0.and.kmt_big(i,j).ne.0) then
              adv_fe_s(i,j) = 
     &        adv_vet_s(i,j) *
     &       ( t(i,  kmt_big(i  ,j),j,n,tau) +
     &         t(i+1,kmt_big(i+1,j),j,n,tau) )
            else
              adv_fe_s(i,j) = 0.0
            endif

        enddo
       enddo
c
c      calculate 2*advective flux (for speed) across bottom face of
c      "T" cells. 
c
       do j=js,je
        do k=1,km-1
         do i=is,ie
           adv_fb_su(i,k,j) =
     &         adv_vbt_su(i,k,j) *(t(i,k,  j,n,tau  ) +
     &                             t(i,k+1,j,n,tau  ))
         enddo
        enddo
#ifndef vector_host
        do i=is,ie
         kb=kmt_big(i,j)
         adv_fb_su(i,kb,j) = 0.
         adv_fb_su(i,0 ,j) = 0.
         do k=kb,km
          adv_fb_su(i,k,j) = 0.
         enddo
        enddo
#endif
#ifdef vector_host
        adv_fb_su(:,0,j) = 0.
        do k=1,km
         do i=is,ie
          if (k>=kmt_big(i,j)) adv_fb_su(i,k,j) = 0.
         enddo
        enddo
#endif
       enddo
c
c      Time integration of sigma advection within the
c      inner bottom layer (big Ts, bottom boxes)
c
       do j=js,je
        do i=is,ie
          k=kmt_big(i,j)
          if(k.ne.0) then   
           kmt_im1 = max(1,kmt_big(i-1,j))
           kmt_ip1 = max(1,kmt_big(i+1,j))
           kmt_jp1 = max(1,kmt_big(i,j+1))
           kmt_jm1 = max(1,kmt_big(i,j-1))
           t(i,k,j,n,taup1) = t(i,k,j,n,taup1)+c2dt*
     &                  ( -1.0 * ADV_Tx_s(i,j)*sigma_weight(i,k,j)
     &                    -1.0 * ADV_Ty_s(i,j)*sigma_weight(i,k,j)
     &                    -1.0 * ADV_Tz_s(i,j)*sigma_weight(i,k,j)
     &                  ) * tmask(i,k,j)
          endif       
        enddo
       enddo
c
c      advection in the outer sigma influenced near bottom boxes
c      (small ts, no bottom boxes)
c
       do j=js,je
         do k=1,km
          do i=is,ie
           if(sigma_weight(i,k,j).ne.0.0) then
            ibottom = 1; if(k.ge.kmt_big(i,j)) ibottom = 0
            t(i,k,j,n,taup1) = t(i,k,j,n,taup1)
     &                 + sigma_weight(i,k,j)*c2dt*
     &                  ( -1.0 * ADV_Tx_su(i,k,j)
     &                    -1.0 * ADV_Ty_su(i,k,j)
     &                    -1.0 * ADV_Tz_su(i,k,j)
     &                  )* tmask(i,k,j)
           endif
          enddo
         enddo
       enddo


      endif

      if (enable_bbl_diffusion) then

c       compute diffusive flux across eastern face of "T" cells,
c       diffusive flux across northern face of "T" cells is
c       computed in DIFF_Ty_s 

         do j=js,je
          do i=is-1,ie
            if(kmt_big(i+1,j).ne.0.and.kmt_big(i,j).ne.0) then
              diff_fe_s(i,j) = ah_sigma*cstr(j)*dxur(i)*
     &            ( t(i+1,kmt_big(i+1,j),j,n,taum1) -
     &              t(i  ,kmt_big(i  ,j),j,n,taum1) )
            else
              diff_fe_s(i,j) = 0.0
            endif
          enddo
         enddo
c
c      Time integration of sigma diffusion within the
c      inner bottom layer (big Ts, bottom boxes)
c
c      compute factors for a sigma diffusion condition
c
       do j=js,je
        do i=is,ie
         k=kmt_big(i,j)
         if(k.ne.0) then       ! start k.ne.0

            factxm=0.0; factxp=0.0; factym=0.0; factyp=0.0
            kn=kmt_big(i-1,j)
            factxm=ah_sigma_min/ah_sigma
            if(kn.gt.1) then
              if((k.gt.kn.and.-1.0*rhodiffxp(i-1,j).lt.0.0).or.
     &           (k.lt.kn.and.-1.0*rhodiffxp(i-1,j).gt.0.0)) then
                factxm=1.
              else
                factxm=ah_sigma_min/ah_sigma
              endif
            endif

            kn=kmt_big(i+1,j)
            factxp=ah_sigma_min/ah_sigma
            if(kn.gt.1) then
              if((k.gt.kn.and.rhodiffxp(i,j).lt.0.0).or.
     &           (k.lt.kn.and.rhodiffxp(i,j).gt.0.0)) then
                factxp=1.
              else
                factxp=ah_sigma_min/ah_sigma
              endif
            endif

            kn=kmt_big(i,j-1)
            factym=ah_sigma_min/ah_sigma
            if(kn.gt.1) then
              if((k.gt.kn.and.-1.0*rhodiffyp(i,j-1).lt.0.0).or.
     &           (k.lt.kn.and.-1.0*rhodiffyp(i,j-1).gt.0.0)) then
                factym=1.
              else
                factym=ah_sigma_min/ah_sigma
              endif
            endif

            kn=kmt_big(i,j+1)
            factyp=ah_sigma_min/ah_sigma
            if(kn.gt.1) then
              if((k.gt.kn.and.rhodiffyp(i,j).lt.0.0).or.
     &           (k.lt.kn.and.rhodiffyp(i,j).gt.0.0)) then
                factyp=1.
              else
                factyp=ah_sigma_min/ah_sigma
              endif
            endif

            kmt_im1 = max(1,kmt_big(i-1,j))
            kmt_ip1 = max(1,kmt_big(i+1,j))
            kmt_jp1 = max(1,kmt_big(i,j+1))
            kmt_jm1 = max(1,kmt_big(i,j-1))

            t(i,k,j,n,taup1) = t(i,k,j,n,taup1)+c2dt*
     &          (DIFF_Tx_s(i,j) + DIFF_Ty_s(i,j))*tmask(i,k,j)

         endif       ! end k.ne.0
        enddo
       enddo
      endif



      end subroutine add_bbl



