#include "BBL_OPTIONS.h"

CBOP
C     !ROUTINE: BBL_CALC_RHS

C     !INTERFACE:
      SUBROUTINE BBL_CALC_RHS(
     I        myTime, myIter, myThid )

C     !DESCRIPTION:
C     Calculate tendency of tracers due to bottom boundary layer.

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "BBL.h"

C     !INPUT PARAMETERS:
C     myTime    :: Current time in simulation
C     myIter    :: Current time-step number
C     myThid    :: my Thread Id number
      _RL     myTime
      INTEGER myIter, myThid

C     !OUTPUT PARAMETERS:

C     !LOCAL VARIABLES:
C     bi,bj     :: Tile indices
C     i,j       :: Loop indices
C     d,r       :: Donnor/Receiver indices
C     kBot      :: k index of bottommost wet grid box
C     kLowC1    :: k index of bottommost (i,j) cell
C     kLowC2    :: k index of bottommost (i+1,j) or (i,j+1) cell
C     kl        :: k index at which to compare 2 cells
C     thk_d     :: thickness of donnor bottommost wet grid cell
C     thk_r     :: thickness of receiver bottommost wet grid cell
C     bblEta_d  :: bbl_eta of donnor cell
C     bblEta_r  :: bbl_eta of receiver cell
C     resThk_r  :: thk_r - bblEta_r
C     Theta_r   :: Theta of receiver cell
C     bblTheta_d:: Theta of donnor bbl
C     bblTheta_r:: Theta of receiver bbl
C     resTheta_r:: Theta of resThk_r
C     Salt_r    :: Salt of receiver cell
C     bblSalt_d :: Salt of donnor bbl
C     bblSalt_r :: Salt of receiver bbl
C     resSalt_r :: Salt of resThk_r
C     deltaRho  :: density change
C     deltaDpt  :: depth change
C     dVol      :: horizontal volume transport
C     bbl_tend  :: temporary variable for tendency terms
C     sloc      :: salinity of bottommost wet grid box
C     tloc      :: temperature of bottommost wet grid box
C     rholoc    :: in situ density of bottommost wet grid box
C     rhoBBL    :: in situ density of bottom boundary layer
C     bbl_rho1  :: local (i,j) density
C     bbl_rho2  :: local (i+1, j) or (i,j+1) density
      INTEGER bi, bj
      INTEGER i, j, d, r, kBot, kLowC1, kLowC2, kl
      _RL     thk_d, thk_r, bblEta_d, bblEta_r, resThk_r
      _RL     Theta_r, bblTheta_d, bblTheta_r, resTheta_r
      _RL     Salt_r,  bblSalt_d,  bblSalt_r,  resSalt_r
      _RL     deltaRho, deltaDpt, dVol, bbl_tend
      _RL     sloc   ( 0:sNx+1, 0:sNy+1 )
      _RL     tloc   ( 0:sNx+1, 0:sNy+1 )
      _RL     rholoc ( 0:sNx+1, 0:sNy+1 )
      _RL     rhoBBL ( 0:sNx+1, 0:sNy+1 )
      _RL     bbl_rho1, bbl_rho2
CEOP

C--   Loops on tile indices bi,bj
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

C     Initialize tendency terms, make local copy of
C     bottomost temperature, salinity, in-situ density
C     and in-situ BBL density.
        DO j=0,sNy+1
         DO i=0,sNx+1
          bbl_TendTheta(i,j,bi,bj) = 0. _d 0
          bbl_TendSalt (i,j,bi,bj) = 0. _d 0
          kBot        = max(1,kLowC(i,j,bi,bj))
          tLoc(i,j)   = theta(i,j,kBot,bi,bj)
          sLoc(i,j)   = salt (i,j,kBot,bi,bj)
          rholoc(i,j) = rhoInSitu(i,j,kBot,bi,bj)
          IF ( kBot .EQ. Nr ) THEN
           rhoBBL(i,j) = bbl_rho_nr(i,j,bi,bj)
          ELSE
           rhoBBL(i,j) = rhoInSitu(i,j,kBot+1,bi,bj)
          ENDIF
         ENDDO
        ENDDO
 
C==== Compute and apply vertical exchange between BBL and
C     residual volume of botommost wet grid box.
C     This operation does not change total tracer quantity
C     in botommost wet grid box.

        DO j=0,sNy+1
         DO i=0,sNx+1
c        DO j=-oly,sNy+oly
c         DO i=-olx,sNx+olx
          kBot = kLowC(i,j,bi,bj)
          IF ( kBot .GT. 0 ) THEN
C     If model density is lower than BBL, slowly diffuse upward.
           IF ( rhoLoc(i,j) .LT. rhoBBL(i,j) )
     &            bbl_eta(i,j,bi,bj) = max ( 0. _d 0 ,
     &            bbl_eta(i,j,bi,bj) - bbl_wvel * dTtracerLev(kBot) )
C     If model density is higher than BBL then mix instantly.
           IF ( rhoLoc(i,j) .GE. rhoBBL(i,j) .OR.
     &          bbl_eta(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
            bbl_theta(i,j,bi,bj) = tLoc(i,j)
            bbl_salt (i,j,bi,bj) = sLoc(i,j)
            bbl_eta  (i,j,bi,bj) = 0. _d 0
           ENDIF
          ENDIF
         ENDDO
        ENDDO

C==== Compute meridional bbl exchange at northern edge.
        j=sNy
        DO i=0,sNx+1
         kLowC1 = kLowC(i,j  ,bi,bj)
         kLowC2 = kLowC(i,j+1,bi,bj)
         IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN

C     Compare the bbl densities at the higher pressure
C     (highest possible density of given t,s)
C     bbl in situ density is stored in k > kLowC indices
          kl = MAX(kLowC1,kLowC2) + 1
          deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
     &               R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
          IF ( deltaDpt .GT. 0. ) THEN
           IF ( kl .GT. Nr ) THEN
            bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
           ELSE
            bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
           ENDIF
           bbl_rho2 = rhoInSitu(i,j+1,kLowC2,bi,bj)
          ELSE
           bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
           IF ( kl .GT. Nr ) THEN
            bbl_rho2 = bbl_rho_nr(i,j+1,bi,bj)
           ELSE
            bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)
           ENDIF
          ENDIF
          deltaRho = bbl_rho2 - bbl_rho1
          IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
C     If heavy BBL water is higher than light BBL water,
C     exchange properties laterally.

C     Determine donnor and receiver cells.
           IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
            d = j
            r = j + 1
           ELSE
            d = j + 1
            r = j
           ENDIF

C     Replenish thickness of donor cell, if needed.
           thk_d = drF(kLowC(i,d,bi,bj)) *
     &          hFacC(i,d,kLowC(i,d,bi,bj),bi,bj)
           IF ( ( bbl_theta(i,d,bi,bj) .EQ. tloc(i,d) ) .AND.
     &          ( bbl_salt (i,d,bi,bj) .EQ. sloc(i,d) ) .AND.
     &          ( bbl_eta  (i,d,bi,bj) .LT. bbl_initEta ) )
     &          bbl_eta(i,d,bi,bj) = min ( bbl_initEta, thk_d )

C     Compute some donnor and receiver cell properties.
           thk_r      = drF(kLowC(i,r,bi,bj)) *
     &                  hFacC(i,r,kLowC(i,r,bi,bj),bi,bj)
           Theta_r    = tLoc(i,r)
           Salt_r     = sLoc(i,r)
           bblTheta_d = bbl_theta(i,d,bi,bj)
           bblTheta_r = bbl_theta(i,r,bi,bj)
           bblSalt_d  = bbl_salt (i,d,bi,bj)
           bblSalt_r  = bbl_salt (i,r,bi,bj)
           bblEta_d   = bbl_eta  (i,d,bi,bj)
           bblEta_r   = bbl_eta  (i,r,bi,bj)
           resThk_r   = thk_r - bblEta_r
           resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
           resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r

C     Compute volume transport from donnor to receiver.
           dVol = min ( bblEta_d * rA(i,d,bi,bj) / 2. _d 0,
     &          resThk_r * rA(i,r,bi,bj) / 2. _d 0,
     &          dxG(i,j+1,bi,bj) * bblEta_d * bbl_hvel * deltaT )

C     Compute temperature tracer tendencies for donor and receiver cell.
           bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
           bbl_TendTheta(i,d,bi,bj) = bbl_TendTheta(i,d,bi,bj) -
     &                                bbl_tend / rA(i,d,bi,bj) / thk_d
           bbl_TendTheta(i,r,bi,bj) = bbl_TendTheta(i,r,bi,bj) +
     &                                bbl_tend / rA(i,r,bi,bj) / thk_r

C     Compute salinity tracer tendencies for donor and receiver cell.
           bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
           bbl_TendSalt(i,d,bi,bj) = bbl_TendSalt(i,d,bi,bj) -
     &                               bbl_tend / rA(i,d,bi,bj) / thk_d
           bbl_TendSalt(i,r,bi,bj) = bbl_TendSalt(i,r,bi,bj) +
     &                               bbl_tend / rA(i,r,bi,bj) / thk_r

C     Adjust pbl thickness and tracer properties.
           bbl_eta(i,d,bi,bj) = bblEta_d - dVol / rA(i,d,bi,bj)
           IF ( bbl_eta(i,d,bi,bj) .LT. 0.0001 ) THEN
            bbl_eta(i,d,bi,bj) = 0. _d 0
            bbl_theta(i,d,bi,bj) = tLoc(i,d)
            bbl_salt (i,d,bi,bj) = sLoc(i,d)
           ENDIF
           bbl_eta(i,r,bi,bj) = bblEta_r + dVol / rA(i,r,bi,bj)
           bbl_theta(i,r,bi,bj) = ( dVol * bblTheta_d +
     &          bblEta_r * rA(i,r,bi,bj) * bblTheta_r ) /
     &          ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
           bbl_salt(i,r,bi,bj)  = ( dVol * bblSalt_d  +
     &          bblEta_r * rA(i,r,bi,bj) * bblSalt_r  ) /
     &          ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
          ENDIF
         ENDIF
        ENDDO

C==== Compute meridional bbl exchange inside tile.
        DO j=0,sNy-1
         DO i=0,sNx+1
          kLowC1 = kLowC(i,j  ,bi,bj)
          kLowC2 = kLowC(i,j+1,bi,bj)
          IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN

C     Compare the bbl densities at the higher pressure
C     (highest possible density of given t,s)
C     bbl in situ density is stored in k > kLowC indices
           kl = MAX(kLowC1,kLowC2) + 1
           deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
     &                R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
           IF ( deltaDpt .GT. 0. ) THEN
            IF ( kl .GT. Nr ) THEN
             bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
            ELSE
             bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
            ENDIF
            bbl_rho2 = rhoInSitu(i,j+1,kLowC2,bi,bj)
           ELSE
            bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
            IF ( kl .GT. Nr ) THEN
             bbl_rho2 = bbl_rho_nr(i,j+1,bi,bj)
            ELSE
             bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)
            ENDIF
           ENDIF
           deltaRho = bbl_rho2 - bbl_rho1
           IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
C     If heavy BBL water is higher than light BBL water,
C     exchange properties laterally.

C     Determine donnor and receiver cells.
            IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
             d = j
             r = j + 1
            ELSE
             d = j + 1
             r = j
            ENDIF

C     Replenish thickness of donor cell, if needed.
            thk_d = drF(kLowC(i,d,bi,bj)) *
     &              hFacC(i,d,kLowC(i,d,bi,bj),bi,bj)
            IF ( ( bbl_theta(i,d,bi,bj) .EQ. tloc(i,d) ) .AND.
     &           ( bbl_salt (i,d,bi,bj) .EQ. sloc(i,d) ) .AND.
     &           ( bbl_eta  (i,d,bi,bj) .LT. bbl_initEta ) )
     &           bbl_eta(i,d,bi,bj) = min ( bbl_initEta, thk_d )

C     Compute some donnor and receiver cell properties.
            thk_r      = drF(kLowC(i,r,bi,bj)) *
     &                   hFacC(i,r,kLowC(i,r,bi,bj),bi,bj)
            Theta_r    = tLoc(i,r)
            Salt_r     = sLoc(i,r)
            bblTheta_d = bbl_theta(i,d,bi,bj)
            bblTheta_r = bbl_theta(i,r,bi,bj)
            bblSalt_d  = bbl_salt (i,d,bi,bj)
            bblSalt_r  = bbl_salt (i,r,bi,bj)
            bblEta_d   = bbl_eta  (i,d,bi,bj)
            bblEta_r   = bbl_eta  (i,r,bi,bj)
            resThk_r   = thk_r - bblEta_r
            resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
            resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r

C     Compute volume transport from donnor to receiver.
            dVol = min ( bblEta_d * rA(i,d,bi,bj) / 2. _d 0,
     &           resThk_r * rA(i,r,bi,bj) / 2. _d 0,
     &           dxG(i,j+1,bi,bj) * bblEta_d * bbl_hvel * deltaT )

C     Compute temperature tracer tendencies for donor and receiver cell.
            bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
            bbl_TendTheta(i,d,bi,bj) = bbl_TendTheta(i,d,bi,bj) -
     &                                 bbl_tend / rA(i,d,bi,bj) / thk_d
            bbl_TendTheta(i,r,bi,bj) = bbl_TendTheta(i,r,bi,bj) +
     &                                 bbl_tend / rA(i,r,bi,bj) / thk_r

C     Compute salinity tracer tendencies for donor and receiver cell.
            bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
            bbl_TendSalt(i,d,bi,bj) = bbl_TendSalt(i,d,bi,bj) -
     &                                bbl_tend / rA(i,d,bi,bj) / thk_d
            bbl_TendSalt(i,r,bi,bj) = bbl_TendSalt(i,r,bi,bj) +
     &                                bbl_tend / rA(i,r,bi,bj) / thk_r

C     Adjust pbl thickness and tracer properties.
            bbl_eta(i,d,bi,bj) = bblEta_d - dVol / rA(i,d,bi,bj)
            IF ( bbl_eta(i,d,bi,bj) .LT. 0.0001 ) THEN
             bbl_eta(i,d,bi,bj) = 0. _d 0
             bbl_theta(i,d,bi,bj) = tLoc(i,d)
             bbl_salt (i,d,bi,bj) = sLoc(i,d)
            ENDIF
            bbl_eta(i,r,bi,bj) = bblEta_r + dVol / rA(i,r,bi,bj)
            bbl_theta(i,r,bi,bj) = ( dVol * bblTheta_d +
     &           bblEta_r * rA(i,r,bi,bj) * bblTheta_r ) /
     &           ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
            bbl_salt(i,r,bi,bj)  = ( dVol * bblSalt_d  +
     &           bblEta_r * rA(i,r,bi,bj) * bblSalt_r  ) /
     &           ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
           ENDIF
          ENDIF
         ENDDO
        ENDDO

C==== Compute zonal bbl exchange at Eastern edge.
        i=sNx
        DO j=1,sNy
         kLowC1 = kLowC(i  ,j,bi,bj)
         kLowC2 = kLowC(i+1,j,bi,bj)
         IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN

C     Compare the bbl densities at the higher pressure
C     (highest possible density of given t,s)
C     bbl in situ density is stored in k > kLowC indices
          kl = MAX(kLowC1,kLowC2) + 1
          deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
     &               R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
          IF ( deltaDpt .GT. 0. ) THEN
           IF ( kl .GT. Nr ) THEN
            bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
           ELSE
            bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
           ENDIF
           bbl_rho2 = rhoInSitu(i+1,j,kLowC2,bi,bj)
          ELSE
           bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
           IF ( kl .GT. Nr ) THEN
            bbl_rho2 = bbl_rho_nr(i+1,j,bi,bj)
           ELSE
            bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)
           ENDIF
          ENDIF
          deltaRho = bbl_rho2 - bbl_rho1
          IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
C     If heavy BBL water is higher than light BBL water,
C     exchange properties laterally.

C     Determine donnor and receiver cells.
           IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
            d = i
            r = i + 1
           ELSE
            d = i + 1
            r = i
           ENDIF

C     Replenish thickness of donor cell, if needed.
           thk_d = drF(kLowC(d,j,bi,bj)) *
     &             hFacC(d,j,kLowC(d,j,bi,bj),bi,bj)
           IF ( ( bbl_theta(d,j,bi,bj) .EQ. tloc(d,j) ) .AND.
     &          ( bbl_salt (d,j,bi,bj) .EQ. sloc(d,j) ) .AND.
     &          ( bbl_eta  (d,j,bi,bj) .LT. bbl_initEta ) )
     &          bbl_eta(d,j,bi,bj) = min ( bbl_initEta, thk_d )

C     Compute some donnor and receiver cell properties.
           thk_r      = drF(kLowC(r,j,bi,bj)) *
     &                  hFacC(r,j,kLowC(r,j,bi,bj),bi,bj)
           Theta_r    = tLoc(r,j)
           Salt_r     = sLoc(r,j)
           bblTheta_d = bbl_theta(d,j,bi,bj)
           bblTheta_r = bbl_theta(r,j,bi,bj)
           bblSalt_d  = bbl_salt (d,j,bi,bj)
           bblSalt_r  = bbl_salt (r,j,bi,bj)
           bblEta_d   = bbl_eta  (d,j,bi,bj)
           bblEta_r   = bbl_eta  (r,j,bi,bj)
           resThk_r   = thk_r - bblEta_r
           resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
           resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r

C     Compute volume transport from donnor to receiver.
           dVol = min ( bblEta_d * rA(d,j,bi,bj) / 2. _d 0,
     &          resThk_r * rA(r,j,bi,bj) / 2. _d 0,
     &          dxG(i+1,j,bi,bj) * bblEta_d * bbl_hvel * deltaT )

C     Compute temperature tracer tendencies for donor and receiver cell.
           bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
           bbl_TendTheta(d,j,bi,bj) = bbl_TendTheta(d,j,bi,bj) -
     &                                bbl_tend / rA(d,j,bi,bj) / thk_d
           bbl_TendTheta(r,j,bi,bj) = bbl_TendTheta(r,j,bi,bj) +
     &                                bbl_tend / rA(r,j,bi,bj) / thk_r

C     Compute salinity tracer tendencies for donor and receiver cell.
           bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
           bbl_TendSalt(d,j,bi,bj) = bbl_TendSalt(d,j,bi,bj) -
     &                               bbl_tend / rA(d,j,bi,bj) / thk_d
           bbl_TendSalt(r,j,bi,bj) = bbl_TendSalt(r,j,bi,bj) +
     &                               bbl_tend / rA(r,j,bi,bj) / thk_r

C     Adjust pbl thickness and tracer properties.
           bbl_eta(d,j,bi,bj) = bblEta_d - dVol / rA(d,j,bi,bj)
           IF ( bbl_eta(d,j,bi,bj) .LT. 0.0001 ) THEN
            bbl_eta(d,j,bi,bj) = 0. _d 0
            bbl_theta(d,j,bi,bj) = tLoc(d,j)
            bbl_salt (d,j,bi,bj) = sLoc(d,j)
           ENDIF
           bbl_eta(r,j,bi,bj) = bblEta_r + dVol / rA(r,j,bi,bj)
           bbl_theta(r,j,bi,bj) = ( dVol * bblTheta_d +
     &          bblEta_r * rA(r,j,bi,bj) * bblTheta_r ) /
     &          ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
           bbl_salt(r,j,bi,bj)  = ( dVol * bblSalt_d  +
     &          bblEta_r * rA(r,j,bi,bj) * bblSalt_r  ) /
     &          ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
          ENDIF
         ENDIF
        ENDDO

C==== Compute zonal bbl exchange inside tile.
        DO j=1,sNy
         DO i=0,sNx-1
          kLowC1 = kLowC(i  ,j,bi,bj)
          kLowC2 = kLowC(i+1,j,bi,bj)
          IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN

C     Compare the bbl densities at the higher pressure
C     (highest possible density of given t,s)
C     bbl in situ density is stored in k > kLowC indices
           kl = MAX(kLowC1,kLowC2) + 1
           deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
     &                R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
           IF ( deltaDpt .GT. 0. ) THEN
            IF ( kl .GT. Nr ) THEN
             bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
            ELSE
             bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
            ENDIF
            bbl_rho2 = rhoInSitu(i+1,j,kLowC2,bi,bj)
           ELSE
            bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
            IF ( kl .GT. Nr ) THEN
             bbl_rho2 = bbl_rho_nr(i+1,j,bi,bj)
            ELSE
             bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)
            ENDIF
           ENDIF
           deltaRho = bbl_rho2 - bbl_rho1
           IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
C     If heavy BBL water is higher than light BBL water,
C     exchange properties laterally.

C     Determine donnor and receiver cells.
            IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
             d = i
             r = i + 1
            ELSE
             d = i + 1
             r = i
            ENDIF

C     Replenish thickness of donor cell, if needed.
            thk_d = drF(kLowC(d,j,bi,bj)) *
     &              hFacC(d,j,kLowC(d,j,bi,bj),bi,bj)
            IF ( ( bbl_theta(d,j,bi,bj) .EQ. tloc(d,j) ) .AND.
     &           ( bbl_salt (d,j,bi,bj) .EQ. sloc(d,j) ) .AND.
     &           ( bbl_eta  (d,j,bi,bj) .LT. bbl_initEta ) )
     &           bbl_eta(d,j,bi,bj) = min ( bbl_initEta, thk_d )

C     Compute some donnor and receiver cell properties.
            thk_r      = drF(kLowC(r,j,bi,bj)) *
     &                   hFacC(r,j,kLowC(r,j,bi,bj),bi,bj)
            Theta_r    = tLoc(r,j)
            Salt_r     = sLoc(r,j)
            bblTheta_d = bbl_theta(d,j,bi,bj)
            bblTheta_r = bbl_theta(r,j,bi,bj)
            bblSalt_d  = bbl_salt (d,j,bi,bj)
            bblSalt_r  = bbl_salt (r,j,bi,bj)
            bblEta_d   = bbl_eta  (d,j,bi,bj)
            bblEta_r   = bbl_eta  (r,j,bi,bj)
            resThk_r   = thk_r - bblEta_r
            resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
            resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r

C     Compute volume transport from donnor to receiver.
            dVol = min ( bblEta_d * rA(d,j,bi,bj) / 2. _d 0,
     &           resThk_r * rA(r,j,bi,bj) / 2. _d 0,
     &           dxG(i+1,j,bi,bj) * bblEta_d * bbl_hvel * deltaT )

C     Compute temperature tracer tendencies for donor and receiver cell.
            bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
            bbl_TendTheta(d,j,bi,bj) = bbl_TendTheta(d,j,bi,bj) -
     &                                 bbl_tend / rA(d,j,bi,bj) / thk_d
            bbl_TendTheta(r,j,bi,bj) = bbl_TendTheta(r,j,bi,bj) +
     &                                 bbl_tend / rA(r,j,bi,bj) / thk_r

C     Compute salinity tracer tendencies for donor and receiver cell.
            bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
            bbl_TendSalt(d,j,bi,bj) = bbl_TendSalt(d,j,bi,bj) -
     &                                bbl_tend / rA(d,j,bi,bj) / thk_d
            bbl_TendSalt(r,j,bi,bj) = bbl_TendSalt(r,j,bi,bj) +
     &                                bbl_tend / rA(r,j,bi,bj) / thk_r

C     Adjust pbl thickness and tracer properties.
            bbl_eta(d,j,bi,bj) = bblEta_d - dVol / rA(d,j,bi,bj)
            IF ( bbl_eta(d,j,bi,bj) .LT. 0.0001 ) THEN
             bbl_eta(d,j,bi,bj) = 0. _d 0
             bbl_theta(d,j,bi,bj) = tLoc(d,j)
             bbl_salt (d,j,bi,bj) = sLoc(d,j)
            ENDIF
            bbl_eta(r,j,bi,bj) = bblEta_r + dVol / rA(r,j,bi,bj)
            bbl_theta(r,j,bi,bj) = ( dVol * bblTheta_d +
     &           bblEta_r * rA(r,j,bi,bj) * bblTheta_r ) /
     &           ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
            bbl_salt(r,j,bi,bj)  = ( dVol * bblSalt_d  +
     &           bblEta_r * rA(r,j,bi,bj) * bblSalt_r  ) /
     &           ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
           ENDIF
          ENDIF
         ENDDO
        ENDDO

C--   end bi,bj loops.
       ENDDO
      ENDDO

      CALL EXCH_XY_RL( bbl_eta      , myThid )
      CALL EXCH_XY_RL( bbl_theta    , myThid )
      CALL EXCH_XY_RL( bbl_salt     , myThid )
      CALL EXCH_XY_RL( bbl_TendTheta, myThid )
      CALL EXCH_XY_RL( bbl_TendSalt , myThid )

      RETURN
      END
