bbl/code/mypackage_calc_rhs.F

C $Header: /u/gcmpack/MITgcm_contrib/bbl/code/mypackage_calc_rhs.F,v 1.2 2010/11/23 07:06:25 dimitri Exp $
C $Name:  $

#include "BBL_OPTIONS.h"

CBOP
C     !ROUTINE: BBL_CALC_RHS

C     !INTERFACE:
      SUBROUTINE MYPACKAGE_CALC_RHS(
     I        bi, bj, myTime, myIter, myThid )

C     !DESCRIPTION:
C     Calculate tendency of tracers due to bottom boundary layer.
C     Scheme is currently coded for dTtracerLev(k) .EQ. deltaT.

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     bi,bj     :: Tile indices
C     myTime    :: Current time in simulation
C     myIter    :: Current time-step number
C     myThid    :: my Thread Id number
      INTEGER bi, bj
      _RL     myTime
      INTEGER myIter, myThid

C     !OUTPUT PARAMETERS:

C     !LOCAL VARIABLES:
C     i,j       :: Loop 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     resThk    :: thickness of bottommost wet grid box minus bbl_eta
C     resTheta  :: temperature of this residual volume
C     resSalt   :: salinity of this residual volume
C     deltaRho  :: density change
C     deltaDpt  :: depth change
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     fZon      :: Zonal flux
C     fMer      :: Meridional flux
C     bbl_rho1  :: local (i,j) density
C     bbl_rho2  :: local (i+1, j) or (i,j+1) density 
      INTEGER i, j, kBot, kLowC1, kLowC2, kl
      _RL     resThk, resTheta, resSalt
      _RL     deltaRho, deltaDpt, bbl_tend
      _RL     bbl_rho1, bbl_rho2
      _RL     sloc    ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
      _RL     tloc    ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
      _RL     rholoc  ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
      _RL     rhoBBL  ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
      _RL     fZon    ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
      _RL     fMer    ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

C     Initialize tendency terms and make local copy of
C     bottomost temperature, salinity, and in-situ desnity
C     and of in-situ BBL density
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        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)
        rhoBBL(i,j) = rhoInSitu(i,j,kBot+1, bi,bj)
       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=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        kBot = kLowC(i,j,bi,bj)
        IF ( kBot .GT. 0 ) THEN
         resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj)

C     If bbl occupies the complete bottom model grid box or
C     if model density is higher than BBL then mix instantly.
         IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.rhoBBL(i,j)) ) THEN
          bbl_theta(i,j,bi,bj) = tLoc(i,j)
          bbl_salt (i,j,bi,bj) = sLoc(i,j)

C     If model density is lower than BBL, slowly diffuse upward.
         ELSE
          resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
     &         (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
          resSalt  = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
     &         (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
          bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
     &         deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR
          bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
     &         deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR
         ENDIF
        ENDIF
       ENDDO
      ENDDO

C==== Compute zonal bbl exchange.
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx-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 (highest possible density of given t,s)
C bbl in situ density is stored in kLowC + 1 index  
         kl = MAX(kLowC1, kLowC2) + 1
         bbl_rho1 = rhoInSitu(i,j,kl,bi,bj) 
         bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)         
         deltaRho = bbl_rho2 - bbl_rho1
         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)

C     If heavy BBL water is higher than light BBL water,
C     exchange properties laterally.
         IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
          bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
     &         ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) /
     &         bbl_RelaxH
          bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) -
     &         ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) *
     &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
     &         ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
          bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
     &         ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) /
     &         bbl_RelaxH
          bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) -
     &         ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) *
     &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
     &         ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
         ENDIF
        ENDIF
       ENDDO
      ENDDO

C==== Compute meridional bbl exchange.
      DO j=1-Oly,sNy+Oly-1
       DO i=1-Olx,sNx+Olx
        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 (highest possible density of given t,s)
C bbl in situ density is stored in kLowC + 1 index  
         kl = MAX(kLowC1, kLowC2) + 1
         bbl_rho1 = rhoInSitu(i,j,kl,bi,bj) 
         bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)         
         deltaRho = bbl_rho2 - bbl_rho1
         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)

C     If heavy BBL water is higher than light BBL water,
C     exchange properties laterally.
         IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
          bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
     &         ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) /
     &         bbl_RelaxH
          bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) -
     &         ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) *
     &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
     &         ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) /
     &         bbl_RelaxH
          bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
     &         ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) /
     &         bbl_RelaxH
          bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) -
     &         ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) *
     &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
     &         ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH )
         ENDIF
        ENDIF
       ENDDO
      ENDDO

C==== Apply lateral BBL exchange then scale tendency term
C     for botommost wet grid box.
      DO j=1-Oly,sNy+Oly-1
       DO i=1-Olx,sNx+Olx-1
        kBot = kLowC(i,j,bi,bj)
        IF ( kBot .GT. 0 ) THEN
         bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
     &        deltaT * bbl_TendTheta(i,j,bi,bj)
         bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
     &        deltaT * bbl_TendSalt (i,j,bi,bj)
         bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) *
     &        bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
         bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) *
     &        bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevB ) THEN
C     Check salinity conservation
       bbl_tend=0
       DO j=1,sNy
        DO i=1,sNx
         kBot = kLowC(i,j,bi,bj)
         IF ( kBot .GT. 0 ) THEN
          bbl_tend = bbl_tend + bbl_TendSalt(i,j,bi,bj) *
     &         hFacC(i,j,kBot,bi,bj) * drF(kBot) *rA(i,j,bi,bj)
         ENDIF
        ENDDO
       ENDDO
       _GLOBAL_SUM_RL( bbl_tend, myThid )
       WRITE(msgBuf,'(A,E10.2)') 'total salt tendency = ', bbl_tend
       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &      SQUEEZE_RIGHT, myThid )
      ENDIF
#endif /* ALLOW_DEBUG */


      CALL EXCH_XY_RL( bbl_theta, myThid )
      CALL EXCH_XY_RL( bbl_salt , myThid )

      RETURN
      END
1	dimitri	1.3	C $Header: /u/gcmpack/MITgcm_contrib/bbl/code/mypackage_calc_rhs.F,v 1.2 2010/11/23 07:06:25 dimitri Exp $
2	dimitri	1.1	C $Name: $
3
4			#include "BBL_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: BBL_CALC_RHS
8
9			C !INTERFACE:
10			SUBROUTINE MYPACKAGE_CALC_RHS(
11			I bi, bj, myTime, myIter, myThid )
12
13			C !DESCRIPTION:
14			C Calculate tendency of tracers due to bottom boundary layer.
15			C Scheme is currently coded for dTtracerLev(k) .EQ. deltaT.
16
17			C !USES:
18			IMPLICIT NONE
19			#include "SIZE.h"
20			#include "EEPARAMS.h"
21			#include "PARAMS.h"
22			#include "GRID.h"
23			#include "DYNVARS.h"
24			#include "BBL.h"
25
26			C !INPUT PARAMETERS:
27			C bi,bj :: Tile indices
28			C myTime :: Current time in simulation
29			C myIter :: Current time-step number
30			C myThid :: my Thread Id number
31			INTEGER bi, bj
32			_RL myTime
33			INTEGER myIter, myThid
34
35			C !OUTPUT PARAMETERS:
36
37			C !LOCAL VARIABLES:
38			C i,j :: Loop indices
39			C kBot :: k index of bottommost wet grid box
40	dimitri	1.3	C kLowC1 :: k index of bottommost (i,j) cell
41			C kLowC2 :: k index of bottommost (i+1,j) or (i,j+1) cell
42			C kl :: k index at which to compare 2 cells
43	dimitri	1.1	C resThk :: thickness of bottommost wet grid box minus bbl_eta
44			C resTheta :: temperature of this residual volume
45			C resSalt :: salinity of this residual volume
46			C deltaRho :: density change
47			C deltaDpt :: depth change
48			C bbl_tend :: temporary variable for tendency terms
49			C sloc :: salinity of bottommost wet grid box
50			C tloc :: temperature of bottommost wet grid box
51	dimitri	1.3	C rholoc :: in situ density of bottommost wet grid box
52			C rhoBBL :: in situ density of bottom boundary layer
53	dimitri	1.1	C fZon :: Zonal flux
54			C fMer :: Meridional flux
55	dimitri	1.3	C bbl_rho1 :: local (i,j) density
56			C bbl_rho2 :: local (i+1, j) or (i,j+1) density
57			INTEGER i, j, kBot, kLowC1, kLowC2, kl
58	dimitri	1.1	_RL resThk, resTheta, resSalt
59			_RL deltaRho, deltaDpt, bbl_tend
60	dimitri	1.3	_RL bbl_rho1, bbl_rho2
61	dimitri	1.1	_RL sloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
62			_RL tloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
63			_RL rholoc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
64	dimitri	1.3	_RL rhoBBL ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
65	dimitri	1.1	_RL fZon ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
66			_RL fMer ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly )
67	dimitri	1.2	CHARACTER*(MAX_LEN_MBUF) msgBuf
68	dimitri	1.1	CEOP
69
70	dimitri	1.3	C Initialize tendency terms and make local copy of
71			C bottomost temperature, salinity, and in-situ desnity
72			C and of in-situ BBL density
73	dimitri	1.1	DO j=1-OLy,sNy+OLy
74			DO i=1-OLx,sNx+OLx
75			bbl_TendTheta(i,j,bi,bj) = 0. _d 0
76			bbl_TendSalt (i,j,bi,bj) = 0. _d 0
77			kBot = max(1,kLowC(i,j,bi,bj))
78	dimitri	1.3	tLoc(i,j) = theta(i,j,kBot,bi,bj)
79			sLoc(i,j) = salt (i,j,kBot,bi,bj)
80			rholoc(i,j) = rhoInSitu(i,j,kBot,bi,bj)
81			rhoBBL(i,j) = rhoInSitu(i,j,kBot+1, bi,bj)
82	dimitri	1.1	ENDDO
83			ENDDO
84
85			C==== Compute and apply vertical exchange between BBL and
86			C residual volume of botommost wet grid box.
87			C This operation does not change total tracer quantity
88			C in botommost wet grid box.
89
90			DO j=1-Oly,sNy+Oly
91			DO i=1-Olx,sNx+Olx
92			kBot = kLowC(i,j,bi,bj)
93			IF ( kBot .GT. 0 ) THEN
94			resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj)
95
96			C If bbl occupies the complete bottom model grid box or
97			C if model density is higher than BBL then mix instantly.
98	dimitri	1.3	IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.rhoBBL(i,j)) ) THEN
99	dimitri	1.1	bbl_theta(i,j,bi,bj) = tLoc(i,j)
100			bbl_salt (i,j,bi,bj) = sLoc(i,j)
101
102			C If model density is lower than BBL, slowly diffuse upward.
103			ELSE
104			resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
105			& (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
106			resSalt = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
107			& (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
108			bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
109			& deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR
110			bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
111			& deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR
112			ENDIF
113			ENDIF
114			ENDDO
115			ENDDO
116
117			C==== Compute zonal bbl exchange.
118			DO j=1-Oly,sNy+Oly
119			DO i=1-Olx,sNx+Olx-1
120	dimitri	1.3	kLowC1 = kLowC(i,j,bi,bj)
121			kLowC2 = kLowC(i+1,j,bi,bj)
122			IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
123			C compare the bbl densities at the higher pressure (highest possible density of given t,s)
124			C bbl in situ density is stored in kLowC + 1 index
125			kl = MAX(kLowC1, kLowC2) + 1
126			bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
127			bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)
128			deltaRho = bbl_rho2 - bbl_rho1
129	dimitri	1.1	deltaDpt = R_low(i ,j,bi,bj) + bbl_eta(i ,j,bi,bj) -
130			& R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
131
132			C If heavy BBL water is higher than light BBL water,
133			C exchange properties laterally.
134			IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
135			bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
136			& ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) /
137			& bbl_RelaxH
138			bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) -
139			& ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) *
140			& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
141			& ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
142			bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
143			& ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) /
144			& bbl_RelaxH
145			bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) -
146			& ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) *
147			& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
148			& ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
149			ENDIF
150			ENDIF
151			ENDDO
152			ENDDO
153
154			C==== Compute meridional bbl exchange.
155			DO j=1-Oly,sNy+Oly-1
156			DO i=1-Olx,sNx+Olx
157	dimitri	1.3	kLowC1 = kLowC(i,j,bi,bj)
158			kLowC2 = kLowC(i,j+1, bi,bj)
159			IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
160			C compare the bbl densities at the higher pressure (highest possible density of given t,s)
161			C bbl in situ density is stored in kLowC + 1 index
162			kl = MAX(kLowC1, kLowC2) + 1
163			bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
164			bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)
165			deltaRho = bbl_rho2 - bbl_rho1
166	dimitri	1.1	deltaDpt = R_low(i,j ,bi,bj) + bbl_eta(i,j ,bi,bj) -
167			& R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
168
169			C If heavy BBL water is higher than light BBL water,
170			C exchange properties laterally.
171			IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
172			bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
173			& ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) /
174			& bbl_RelaxH
175			bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) -
176			& ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) *
177			& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
178			& ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) /
179			& bbl_RelaxH
180			bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
181			& ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) /
182			& bbl_RelaxH
183			bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) -
184			& ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) *
185			& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) /
186			& ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH )
187			ENDIF
188			ENDIF
189			ENDDO
190			ENDDO
191
192			C==== Apply lateral BBL exchange then scale tendency term
193			C for botommost wet grid box.
194			DO j=1-Oly,sNy+Oly-1
195			DO i=1-Olx,sNx+Olx-1
196			kBot = kLowC(i,j,bi,bj)
197			IF ( kBot .GT. 0 ) THEN
198			bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
199			& deltaT * bbl_TendTheta(i,j,bi,bj)
200			bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
201			& deltaT * bbl_TendSalt (i,j,bi,bj)
202			bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) *
203			& bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
204			bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) *
205			& bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
206			ENDIF
207			ENDDO
208			ENDDO
209
210	dimitri	1.2	#ifdef ALLOW_DEBUG
211			IF ( debugLevel .GE. debLevB ) THEN
212			C Check salinity conservation
213			bbl_tend=0
214			DO j=1,sNy
215			DO i=1,sNx
216			kBot = kLowC(i,j,bi,bj)
217			IF ( kBot .GT. 0 ) THEN
218			bbl_tend = bbl_tend + bbl_TendSalt(i,j,bi,bj) *
219			& hFacC(i,j,kBot,bi,bj) * drF(kBot) *rA(i,j,bi,bj)
220			ENDIF
221			ENDDO
222			ENDDO
223			_GLOBAL_SUM_RL( bbl_tend, myThid )
224			WRITE(msgBuf,'(A,E10.2)') 'total salt tendency = ', bbl_tend
225			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
226			& SQUEEZE_RIGHT, myThid )
227			ENDIF
228			#endif /* ALLOW_DEBUG */
229
230
231	dimitri	1.1	CALL EXCH_XY_RL( bbl_theta, myThid )
232			CALL EXCH_XY_RL( bbl_salt , myThid )
233
234			RETURN
235			END