Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:38:02 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
53b68b7823 Dimi*0001 #include "BBL_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: BBL_CALC_RHS
                0005 
                0006 C     !INTERFACE:
                0007       SUBROUTINE BBL_CALC_RHS(
                0008      I        myTime, myIter, myThid )
                0009 
                0010 C     !DESCRIPTION:
                0011 C     Calculate tendency of tracers due to bottom boundary layer.
                0012 
                0013 C     !USES:
                0014       IMPLICIT NONE
                0015 #include "SIZE.h"
                0016 #include "EEPARAMS.h"
                0017 #include "PARAMS.h"
                0018 #include "GRID.h"
                0019 #include "DYNVARS.h"
                0020 #include "BBL.h"
                0021 
                0022 C     !INPUT PARAMETERS:
                0023 C     myTime    :: Current time in simulation
                0024 C     myIter    :: Current time-step number
                0025 C     myThid    :: my Thread Id number
                0026       _RL     myTime
                0027       INTEGER myIter, myThid
                0028 
                0029 C     !OUTPUT PARAMETERS:
                0030 
                0031 C     !LOCAL VARIABLES:
                0032 C     bi,bj     :: Tile indices
                0033 C     i,j       :: Loop indices
5f5cc705f4 Dimi*0034 C     d,r       :: Donnor/Receiver indices
53b68b7823 Dimi*0035 C     kBot      :: k index of bottommost wet grid box
                0036 C     kLowC1    :: k index of bottommost (i,j) cell
                0037 C     kLowC2    :: k index of bottommost (i+1,j) or (i,j+1) cell
                0038 C     kl        :: k index at which to compare 2 cells
5f5cc705f4 Dimi*0039 C     thk_d     :: thickness of donnor bottommost wet grid cell
                0040 C     thk_r     :: thickness of receiver bottommost wet grid cell
                0041 C     bblEta_d  :: bbl_eta of donnor cell
                0042 C     bblEta_r  :: bbl_eta of receiver cell
                0043 C     resThk_r  :: thk_r - bblEta_r
                0044 C     Theta_r   :: Theta of receiver cell
                0045 C     bblTheta_d:: Theta of donnor bbl
                0046 C     bblTheta_r:: Theta of receiver bbl
                0047 C     resTheta_r:: Theta of resThk_r
                0048 C     Salt_r    :: Salt of receiver cell
                0049 C     bblSalt_d :: Salt of donnor bbl
                0050 C     bblSalt_r :: Salt of receiver bbl
                0051 C     resSalt_r :: Salt of resThk_r
53b68b7823 Dimi*0052 C     deltaRho  :: density change
                0053 C     deltaDpt  :: depth change
5f5cc705f4 Dimi*0054 C     dVol      :: horizontal volume transport
53b68b7823 Dimi*0055 C     bbl_tend  :: temporary variable for tendency terms
                0056 C     sloc      :: salinity of bottommost wet grid box
                0057 C     tloc      :: temperature of bottommost wet grid box
                0058 C     rholoc    :: in situ density of bottommost wet grid box
                0059 C     rhoBBL    :: in situ density of bottom boundary layer
                0060 C     bbl_rho1  :: local (i,j) density
ddd44d55cd Jean*0061 C     bbl_rho2  :: local (i+1, j) or (i,j+1) density
53b68b7823 Dimi*0062       INTEGER bi, bj
5f5cc705f4 Dimi*0063       INTEGER i, j, d, r, kBot, kLowC1, kLowC2, kl
                0064       _RL     thk_d, thk_r, bblEta_d, bblEta_r, resThk_r
                0065       _RL     Theta_r, bblTheta_d, bblTheta_r, resTheta_r
                0066       _RL     Salt_r,  bblSalt_d,  bblSalt_r,  resSalt_r
                0067       _RL     deltaRho, deltaDpt, dVol, bbl_tend
                0068       _RL     sloc   ( 0:sNx+1, 0:sNy+1 )
                0069       _RL     tloc   ( 0:sNx+1, 0:sNy+1 )
                0070       _RL     rholoc ( 0:sNx+1, 0:sNy+1 )
                0071       _RL     rhoBBL ( 0:sNx+1, 0:sNy+1 )
53b68b7823 Dimi*0072       _RL     bbl_rho1, bbl_rho2
                0073 CEOP
                0074 
                0075 C--   Loops on tile indices bi,bj
                0076       DO bj=myByLo(myThid),myByHi(myThid)
                0077        DO bi=myBxLo(myThid),myBxHi(myThid)
                0078 
97c15a0afc Dimi*0079 C     Initialize tendency terms, make local copy of
                0080 C     bottomost temperature, salinity, in-situ density
5f5cc705f4 Dimi*0081 C     and in-situ BBL density.
                0082         DO j=0,sNy+1
                0083          DO i=0,sNx+1
53b68b7823 Dimi*0084           bbl_TendTheta(i,j,bi,bj) = 0. _d 0
                0085           bbl_TendSalt (i,j,bi,bj) = 0. _d 0
15338fa568 Dimi*0086           kBot        = max(1,kLowC(i,j,bi,bj))
53b68b7823 Dimi*0087           tLoc(i,j)   = theta(i,j,kBot,bi,bj)
                0088           sLoc(i,j)   = salt (i,j,kBot,bi,bj)
                0089           rholoc(i,j) = rhoInSitu(i,j,kBot,bi,bj)
15338fa568 Dimi*0090           IF ( kBot .EQ. Nr ) THEN
                0091            rhoBBL(i,j) = bbl_rho_nr(i,j,bi,bj)
                0092           ELSE
5f5cc705f4 Dimi*0093            rhoBBL(i,j) = rhoInSitu(i,j,kBot+1,bi,bj)
15338fa568 Dimi*0094           ENDIF
53b68b7823 Dimi*0095          ENDDO
                0096         ENDDO
5f5cc705f4 Dimi*0097  
53b68b7823 Dimi*0098 C==== Compute and apply vertical exchange between BBL and
                0099 C     residual volume of botommost wet grid box.
                0100 C     This operation does not change total tracer quantity
                0101 C     in botommost wet grid box.
                0102 
5f5cc705f4 Dimi*0103         DO j=0,sNy+1
                0104          DO i=0,sNx+1
                0105 c        DO j=-oly,sNy+oly
                0106 c         DO i=-olx,sNx+olx
53b68b7823 Dimi*0107           kBot = kLowC(i,j,bi,bj)
                0108           IF ( kBot .GT. 0 ) THEN
5f5cc705f4 Dimi*0109 C     If model density is lower than BBL, slowly diffuse upward.
                0110            IF ( rhoLoc(i,j) .LT. rhoBBL(i,j) )
                0111      &            bbl_eta(i,j,bi,bj) = max ( 0. _d 0 ,
                0112      &            bbl_eta(i,j,bi,bj) - bbl_wvel * dTtracerLev(kBot) )
                0113 C     If model density is higher than BBL then mix instantly.
                0114            IF ( rhoLoc(i,j) .GE. rhoBBL(i,j) .OR.
                0115      &          bbl_eta(i,j,bi,bj) .EQ. 0. _d 0 ) THEN
53b68b7823 Dimi*0116             bbl_theta(i,j,bi,bj) = tLoc(i,j)
                0117             bbl_salt (i,j,bi,bj) = sLoc(i,j)
5f5cc705f4 Dimi*0118             bbl_eta  (i,j,bi,bj) = 0. _d 0
53b68b7823 Dimi*0119            ENDIF
                0120           ENDIF
                0121          ENDDO
                0122         ENDDO
                0123 
5f5cc705f4 Dimi*0124 C==== Compute meridional bbl exchange at northern edge.
                0125         j=sNy
                0126         DO i=0,sNx+1
                0127          kLowC1 = kLowC(i,j  ,bi,bj)
                0128          kLowC2 = kLowC(i,j+1,bi,bj)
                0129          IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
                0130 
53b68b7823 Dimi*0131 C     Compare the bbl densities at the higher pressure
                0132 C     (highest possible density of given t,s)
5f5cc705f4 Dimi*0133 C     bbl in situ density is stored in k > kLowC indices
                0134           kl = MAX(kLowC1,kLowC2) + 1
                0135           deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
                0136      &               R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
                0137           IF ( deltaDpt .GT. 0. ) THEN
15338fa568 Dimi*0138            IF ( kl .GT. Nr ) THEN
                0139             bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
5f5cc705f4 Dimi*0140            ELSE
15338fa568 Dimi*0141             bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
5f5cc705f4 Dimi*0142            ENDIF
                0143            bbl_rho2 = rhoInSitu(i,j+1,kLowC2,bi,bj)
                0144           ELSE
                0145            bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
                0146            IF ( kl .GT. Nr ) THEN
                0147             bbl_rho2 = bbl_rho_nr(i,j+1,bi,bj)
                0148            ELSE
                0149             bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)
                0150            ENDIF
                0151           ENDIF
                0152           deltaRho = bbl_rho2 - bbl_rho1
                0153           IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
                0154 C     If heavy BBL water is higher than light BBL water,
                0155 C     exchange properties laterally.
                0156 
                0157 C     Determine donnor and receiver cells.
                0158            IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
                0159             d = j
                0160             r = j + 1
                0161            ELSE
                0162             d = j + 1
                0163             r = j
                0164            ENDIF
                0165 
                0166 C     Replenish thickness of donor cell, if needed.
                0167            thk_d = drF(kLowC(i,d,bi,bj)) *
                0168      &          hFacC(i,d,kLowC(i,d,bi,bj),bi,bj)
                0169            IF ( ( bbl_theta(i,d,bi,bj) .EQ. tloc(i,d) ) .AND.
                0170      &          ( bbl_salt (i,d,bi,bj) .EQ. sloc(i,d) ) .AND.
                0171      &          ( bbl_eta  (i,d,bi,bj) .LT. bbl_initEta ) )
                0172      &          bbl_eta(i,d,bi,bj) = min ( bbl_initEta, thk_d )
53b68b7823 Dimi*0173 
5f5cc705f4 Dimi*0174 C     Compute some donnor and receiver cell properties.
                0175            thk_r      = drF(kLowC(i,r,bi,bj)) *
                0176      &                  hFacC(i,r,kLowC(i,r,bi,bj),bi,bj)
                0177            Theta_r    = tLoc(i,r)
                0178            Salt_r     = sLoc(i,r)
                0179            bblTheta_d = bbl_theta(i,d,bi,bj)
                0180            bblTheta_r = bbl_theta(i,r,bi,bj)
                0181            bblSalt_d  = bbl_salt (i,d,bi,bj)
                0182            bblSalt_r  = bbl_salt (i,r,bi,bj)
                0183            bblEta_d   = bbl_eta  (i,d,bi,bj)
                0184            bblEta_r   = bbl_eta  (i,r,bi,bj)
                0185            resThk_r   = thk_r - bblEta_r
                0186            resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
                0187            resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r
                0188 
                0189 C     Compute volume transport from donnor to receiver.
                0190            dVol = min ( bblEta_d * rA(i,d,bi,bj) / 2. _d 0,
                0191      &          resThk_r * rA(i,r,bi,bj) / 2. _d 0,
                0192      &          dxG(i,j+1,bi,bj) * bblEta_d * bbl_hvel * deltaT )
                0193 
                0194 C     Compute temperature tracer tendencies for donor and receiver cell.
                0195            bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
                0196            bbl_TendTheta(i,d,bi,bj) = bbl_TendTheta(i,d,bi,bj) -
                0197      &                                bbl_tend / rA(i,d,bi,bj) / thk_d
                0198            bbl_TendTheta(i,r,bi,bj) = bbl_TendTheta(i,r,bi,bj) +
                0199      &                                bbl_tend / rA(i,r,bi,bj) / thk_r
                0200 
                0201 C     Compute salinity tracer tendencies for donor and receiver cell.
                0202            bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
                0203            bbl_TendSalt(i,d,bi,bj) = bbl_TendSalt(i,d,bi,bj) -
                0204      &                               bbl_tend / rA(i,d,bi,bj) / thk_d
                0205            bbl_TendSalt(i,r,bi,bj) = bbl_TendSalt(i,r,bi,bj) +
                0206      &                               bbl_tend / rA(i,r,bi,bj) / thk_r
                0207 
                0208 C     Adjust pbl thickness and tracer properties.
                0209            bbl_eta(i,d,bi,bj) = bblEta_d - dVol / rA(i,d,bi,bj)
                0210            IF ( bbl_eta(i,d,bi,bj) .LT. 0.0001 ) THEN
                0211             bbl_eta(i,d,bi,bj) = 0. _d 0
                0212             bbl_theta(i,d,bi,bj) = tLoc(i,d)
                0213             bbl_salt (i,d,bi,bj) = sLoc(i,d)
                0214            ENDIF
                0215            bbl_eta(i,r,bi,bj) = bblEta_r + dVol / rA(i,r,bi,bj)
                0216            bbl_theta(i,r,bi,bj) = ( dVol * bblTheta_d +
                0217      &          bblEta_r * rA(i,r,bi,bj) * bblTheta_r ) /
                0218      &          ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
                0219            bbl_salt(i,r,bi,bj)  = ( dVol * bblSalt_d  +
                0220      &          bblEta_r * rA(i,r,bi,bj) * bblSalt_r  ) /
                0221      &          ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
                0222           ENDIF
                0223          ENDIF
                0224         ENDDO
                0225 
                0226 C==== Compute meridional bbl exchange inside tile.
                0227         DO j=0,sNy-1
                0228          DO i=0,sNx+1
                0229           kLowC1 = kLowC(i,j  ,bi,bj)
                0230           kLowC2 = kLowC(i,j+1,bi,bj)
                0231           IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
                0232 
                0233 C     Compare the bbl densities at the higher pressure
                0234 C     (highest possible density of given t,s)
                0235 C     bbl in situ density is stored in k > kLowC indices
                0236            kl = MAX(kLowC1,kLowC2) + 1
                0237            deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
                0238      &                R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
                0239            IF ( deltaDpt .GT. 0. ) THEN
                0240             IF ( kl .GT. Nr ) THEN
                0241              bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
                0242             ELSE
                0243              bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
                0244             ENDIF
                0245             bbl_rho2 = rhoInSitu(i,j+1,kLowC2,bi,bj)
                0246            ELSE
                0247             bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
                0248             IF ( kl .GT. Nr ) THEN
                0249              bbl_rho2 = bbl_rho_nr(i,j+1,bi,bj)
                0250             ELSE
                0251              bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)
                0252             ENDIF
                0253            ENDIF
                0254            deltaRho = bbl_rho2 - bbl_rho1
                0255            IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
53b68b7823 Dimi*0256 C     If heavy BBL water is higher than light BBL water,
                0257 C     exchange properties laterally.
5f5cc705f4 Dimi*0258 
                0259 C     Determine donnor and receiver cells.
                0260             IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
                0261              d = j
                0262              r = j + 1
                0263             ELSE
                0264              d = j + 1
                0265              r = j
                0266             ENDIF
                0267 
                0268 C     Replenish thickness of donor cell, if needed.
                0269             thk_d = drF(kLowC(i,d,bi,bj)) *
                0270      &              hFacC(i,d,kLowC(i,d,bi,bj),bi,bj)
                0271             IF ( ( bbl_theta(i,d,bi,bj) .EQ. tloc(i,d) ) .AND.
                0272      &           ( bbl_salt (i,d,bi,bj) .EQ. sloc(i,d) ) .AND.
                0273      &           ( bbl_eta  (i,d,bi,bj) .LT. bbl_initEta ) )
                0274      &           bbl_eta(i,d,bi,bj) = min ( bbl_initEta, thk_d )
                0275 
                0276 C     Compute some donnor and receiver cell properties.
                0277             thk_r      = drF(kLowC(i,r,bi,bj)) *
                0278      &                   hFacC(i,r,kLowC(i,r,bi,bj),bi,bj)
                0279             Theta_r    = tLoc(i,r)
                0280             Salt_r     = sLoc(i,r)
                0281             bblTheta_d = bbl_theta(i,d,bi,bj)
                0282             bblTheta_r = bbl_theta(i,r,bi,bj)
                0283             bblSalt_d  = bbl_salt (i,d,bi,bj)
                0284             bblSalt_r  = bbl_salt (i,r,bi,bj)
                0285             bblEta_d   = bbl_eta  (i,d,bi,bj)
                0286             bblEta_r   = bbl_eta  (i,r,bi,bj)
                0287             resThk_r   = thk_r - bblEta_r
                0288             resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
                0289             resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r
                0290 
                0291 C     Compute volume transport from donnor to receiver.
                0292             dVol = min ( bblEta_d * rA(i,d,bi,bj) / 2. _d 0,
                0293      &           resThk_r * rA(i,r,bi,bj) / 2. _d 0,
                0294      &           dxG(i,j+1,bi,bj) * bblEta_d * bbl_hvel * deltaT )
                0295 
                0296 C     Compute temperature tracer tendencies for donor and receiver cell.
                0297             bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
                0298             bbl_TendTheta(i,d,bi,bj) = bbl_TendTheta(i,d,bi,bj) -
                0299      &                                 bbl_tend / rA(i,d,bi,bj) / thk_d
                0300             bbl_TendTheta(i,r,bi,bj) = bbl_TendTheta(i,r,bi,bj) +
                0301      &                                 bbl_tend / rA(i,r,bi,bj) / thk_r
                0302 
                0303 C     Compute salinity tracer tendencies for donor and receiver cell.
                0304             bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
                0305             bbl_TendSalt(i,d,bi,bj) = bbl_TendSalt(i,d,bi,bj) -
                0306      &                                bbl_tend / rA(i,d,bi,bj) / thk_d
                0307             bbl_TendSalt(i,r,bi,bj) = bbl_TendSalt(i,r,bi,bj) +
                0308      &                                bbl_tend / rA(i,r,bi,bj) / thk_r
                0309 
                0310 C     Adjust pbl thickness and tracer properties.
                0311             bbl_eta(i,d,bi,bj) = bblEta_d - dVol / rA(i,d,bi,bj)
                0312             IF ( bbl_eta(i,d,bi,bj) .LT. 0.0001 ) THEN
                0313              bbl_eta(i,d,bi,bj) = 0. _d 0
                0314              bbl_theta(i,d,bi,bj) = tLoc(i,d)
                0315              bbl_salt (i,d,bi,bj) = sLoc(i,d)
                0316             ENDIF
                0317             bbl_eta(i,r,bi,bj) = bblEta_r + dVol / rA(i,r,bi,bj)
                0318             bbl_theta(i,r,bi,bj) = ( dVol * bblTheta_d +
                0319      &           bblEta_r * rA(i,r,bi,bj) * bblTheta_r ) /
                0320      &           ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
                0321             bbl_salt(i,r,bi,bj)  = ( dVol * bblSalt_d  +
                0322      &           bblEta_r * rA(i,r,bi,bj) * bblSalt_r  ) /
                0323      &           ( bbl_eta(i,r,bi,bj) * rA(i,r,bi,bj) )
53b68b7823 Dimi*0324            ENDIF
                0325           ENDIF
                0326          ENDDO
                0327         ENDDO
                0328 
5f5cc705f4 Dimi*0329 C==== Compute zonal bbl exchange at Eastern edge.
                0330         i=sNx
                0331         DO j=1,sNy
                0332          kLowC1 = kLowC(i  ,j,bi,bj)
                0333          kLowC2 = kLowC(i+1,j,bi,bj)
                0334          IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
                0335 
                0336 C     Compare the bbl densities at the higher pressure
53b68b7823 Dimi*0337 C     (highest possible density of given t,s)
5f5cc705f4 Dimi*0338 C     bbl in situ density is stored in k > kLowC indices
                0339           kl = MAX(kLowC1,kLowC2) + 1
                0340           deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
                0341      &               R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
                0342           IF ( deltaDpt .GT. 0. ) THEN
15338fa568 Dimi*0343            IF ( kl .GT. Nr ) THEN
                0344             bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
5f5cc705f4 Dimi*0345            ELSE
15338fa568 Dimi*0346             bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
5f5cc705f4 Dimi*0347            ENDIF
                0348            bbl_rho2 = rhoInSitu(i+1,j,kLowC2,bi,bj)
                0349           ELSE
                0350            bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
                0351            IF ( kl .GT. Nr ) THEN
                0352             bbl_rho2 = bbl_rho_nr(i+1,j,bi,bj)
                0353            ELSE
                0354             bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)
                0355            ENDIF
                0356           ENDIF
                0357           deltaRho = bbl_rho2 - bbl_rho1
                0358           IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
53b68b7823 Dimi*0359 C     If heavy BBL water is higher than light BBL water,
                0360 C     exchange properties laterally.
5f5cc705f4 Dimi*0361 
                0362 C     Determine donnor and receiver cells.
                0363            IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
                0364             d = i
                0365             r = i + 1
                0366            ELSE
                0367             d = i + 1
                0368             r = i
53b68b7823 Dimi*0369            ENDIF
5f5cc705f4 Dimi*0370 
                0371 C     Replenish thickness of donor cell, if needed.
                0372            thk_d = drF(kLowC(d,j,bi,bj)) *
                0373      &             hFacC(d,j,kLowC(d,j,bi,bj),bi,bj)
                0374            IF ( ( bbl_theta(d,j,bi,bj) .EQ. tloc(d,j) ) .AND.
                0375      &          ( bbl_salt (d,j,bi,bj) .EQ. sloc(d,j) ) .AND.
                0376      &          ( bbl_eta  (d,j,bi,bj) .LT. bbl_initEta ) )
                0377      &          bbl_eta(d,j,bi,bj) = min ( bbl_initEta, thk_d )
                0378 
                0379 C     Compute some donnor and receiver cell properties.
                0380            thk_r      = drF(kLowC(r,j,bi,bj)) *
                0381      &                  hFacC(r,j,kLowC(r,j,bi,bj),bi,bj)
                0382            Theta_r    = tLoc(r,j)
                0383            Salt_r     = sLoc(r,j)
                0384            bblTheta_d = bbl_theta(d,j,bi,bj)
                0385            bblTheta_r = bbl_theta(r,j,bi,bj)
                0386            bblSalt_d  = bbl_salt (d,j,bi,bj)
                0387            bblSalt_r  = bbl_salt (r,j,bi,bj)
                0388            bblEta_d   = bbl_eta  (d,j,bi,bj)
                0389            bblEta_r   = bbl_eta  (r,j,bi,bj)
                0390            resThk_r   = thk_r - bblEta_r
                0391            resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
                0392            resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r
                0393 
                0394 C     Compute volume transport from donnor to receiver.
                0395            dVol = min ( bblEta_d * rA(d,j,bi,bj) / 2. _d 0,
                0396      &          resThk_r * rA(r,j,bi,bj) / 2. _d 0,
                0397      &          dxG(i+1,j,bi,bj) * bblEta_d * bbl_hvel * deltaT )
                0398 
                0399 C     Compute temperature tracer tendencies for donor and receiver cell.
                0400            bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
                0401            bbl_TendTheta(d,j,bi,bj) = bbl_TendTheta(d,j,bi,bj) -
                0402      &                                bbl_tend / rA(d,j,bi,bj) / thk_d
                0403            bbl_TendTheta(r,j,bi,bj) = bbl_TendTheta(r,j,bi,bj) +
                0404      &                                bbl_tend / rA(r,j,bi,bj) / thk_r
                0405 
                0406 C     Compute salinity tracer tendencies for donor and receiver cell.
                0407            bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
                0408            bbl_TendSalt(d,j,bi,bj) = bbl_TendSalt(d,j,bi,bj) -
                0409      &                               bbl_tend / rA(d,j,bi,bj) / thk_d
                0410            bbl_TendSalt(r,j,bi,bj) = bbl_TendSalt(r,j,bi,bj) +
                0411      &                               bbl_tend / rA(r,j,bi,bj) / thk_r
                0412 
                0413 C     Adjust pbl thickness and tracer properties.
                0414            bbl_eta(d,j,bi,bj) = bblEta_d - dVol / rA(d,j,bi,bj)
                0415            IF ( bbl_eta(d,j,bi,bj) .LT. 0.0001 ) THEN
                0416             bbl_eta(d,j,bi,bj) = 0. _d 0
                0417             bbl_theta(d,j,bi,bj) = tLoc(d,j)
                0418             bbl_salt (d,j,bi,bj) = sLoc(d,j)
                0419            ENDIF
                0420            bbl_eta(r,j,bi,bj) = bblEta_r + dVol / rA(r,j,bi,bj)
                0421            bbl_theta(r,j,bi,bj) = ( dVol * bblTheta_d +
                0422      &          bblEta_r * rA(r,j,bi,bj) * bblTheta_r ) /
                0423      &          ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
                0424            bbl_salt(r,j,bi,bj)  = ( dVol * bblSalt_d  +
                0425      &          bblEta_r * rA(r,j,bi,bj) * bblSalt_r  ) /
                0426      &          ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
53b68b7823 Dimi*0427           ENDIF
5f5cc705f4 Dimi*0428          ENDIF
53b68b7823 Dimi*0429         ENDDO
                0430 
5f5cc705f4 Dimi*0431 C==== Compute zonal bbl exchange inside tile.
                0432         DO j=1,sNy
                0433          DO i=0,sNx-1
                0434           kLowC1 = kLowC(i  ,j,bi,bj)
                0435           kLowC2 = kLowC(i+1,j,bi,bj)
                0436           IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
                0437 
                0438 C     Compare the bbl densities at the higher pressure
                0439 C     (highest possible density of given t,s)
                0440 C     bbl in situ density is stored in k > kLowC indices
                0441            kl = MAX(kLowC1,kLowC2) + 1
                0442            deltaDpt = R_low(i,j,bi,bj)   + bbl_eta(i,j,bi,bj) -
                0443      &                R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
                0444            IF ( deltaDpt .GT. 0. ) THEN
                0445             IF ( kl .GT. Nr ) THEN
                0446              bbl_rho1 = bbl_rho_nr(i,j,bi,bj)
                0447             ELSE
                0448              bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
                0449             ENDIF
                0450             bbl_rho2 = rhoInSitu(i+1,j,kLowC2,bi,bj)
                0451            ELSE
                0452             bbl_rho1 = rhoInSitu(i,j,kLowC1,bi,bj)
                0453             IF ( kl .GT. Nr ) THEN
                0454              bbl_rho2 = bbl_rho_nr(i+1,j,bi,bj)
                0455             ELSE
                0456              bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)
                0457             ENDIF
                0458            ENDIF
                0459            deltaRho = bbl_rho2 - bbl_rho1
                0460            IF ( (deltaRho*deltaDpt) .LT. 0. ) THEN
                0461 C     If heavy BBL water is higher than light BBL water,
                0462 C     exchange properties laterally.
                0463 
                0464 C     Determine donnor and receiver cells.
                0465             IF ( bbl_rho1 .GT. bbl_rho2 ) THEN
                0466              d = i
                0467              r = i + 1
                0468             ELSE
                0469              d = i + 1
                0470              r = i
                0471             ENDIF
                0472 
                0473 C     Replenish thickness of donor cell, if needed.
                0474             thk_d = drF(kLowC(d,j,bi,bj)) *
                0475      &              hFacC(d,j,kLowC(d,j,bi,bj),bi,bj)
                0476             IF ( ( bbl_theta(d,j,bi,bj) .EQ. tloc(d,j) ) .AND.
                0477      &           ( bbl_salt (d,j,bi,bj) .EQ. sloc(d,j) ) .AND.
                0478      &           ( bbl_eta  (d,j,bi,bj) .LT. bbl_initEta ) )
                0479      &           bbl_eta(d,j,bi,bj) = min ( bbl_initEta, thk_d )
                0480 
                0481 C     Compute some donnor and receiver cell properties.
                0482             thk_r      = drF(kLowC(r,j,bi,bj)) *
                0483      &                   hFacC(r,j,kLowC(r,j,bi,bj),bi,bj)
                0484             Theta_r    = tLoc(r,j)
                0485             Salt_r     = sLoc(r,j)
                0486             bblTheta_d = bbl_theta(d,j,bi,bj)
                0487             bblTheta_r = bbl_theta(r,j,bi,bj)
                0488             bblSalt_d  = bbl_salt (d,j,bi,bj)
                0489             bblSalt_r  = bbl_salt (r,j,bi,bj)
                0490             bblEta_d   = bbl_eta  (d,j,bi,bj)
                0491             bblEta_r   = bbl_eta  (r,j,bi,bj)
                0492             resThk_r   = thk_r - bblEta_r
                0493             resTheta_r = (Theta_r*thk_r-bblTheta_r*bblEta_r)/resThk_r
                0494             resSalt_r  = (Salt_r *thk_r-bblSalt_r *bblEta_r)/resThk_r
                0495 
                0496 C     Compute volume transport from donnor to receiver.
                0497             dVol = min ( bblEta_d * rA(d,j,bi,bj) / 2. _d 0,
                0498      &           resThk_r * rA(r,j,bi,bj) / 2. _d 0,
                0499      &           dxG(i+1,j,bi,bj) * bblEta_d * bbl_hvel * deltaT )
                0500 
                0501 C     Compute temperature tracer tendencies for donor and receiver cell.
                0502             bbl_tend = dVol * (bblTheta_d - resTheta_r) / deltaT
                0503             bbl_TendTheta(d,j,bi,bj) = bbl_TendTheta(d,j,bi,bj) -
                0504      &                                 bbl_tend / rA(d,j,bi,bj) / thk_d
                0505             bbl_TendTheta(r,j,bi,bj) = bbl_TendTheta(r,j,bi,bj) +
                0506      &                                 bbl_tend / rA(r,j,bi,bj) / thk_r
                0507 
                0508 C     Compute salinity tracer tendencies for donor and receiver cell.
                0509             bbl_tend = dVol * (bblSalt_d - resSalt_r) / deltaT
                0510             bbl_TendSalt(d,j,bi,bj) = bbl_TendSalt(d,j,bi,bj) -
                0511      &                                bbl_tend / rA(d,j,bi,bj) / thk_d
                0512             bbl_TendSalt(r,j,bi,bj) = bbl_TendSalt(r,j,bi,bj) +
                0513      &                                bbl_tend / rA(r,j,bi,bj) / thk_r
                0514 
                0515 C     Adjust pbl thickness and tracer properties.
                0516             bbl_eta(d,j,bi,bj) = bblEta_d - dVol / rA(d,j,bi,bj)
                0517             IF ( bbl_eta(d,j,bi,bj) .LT. 0.0001 ) THEN
                0518              bbl_eta(d,j,bi,bj) = 0. _d 0
                0519              bbl_theta(d,j,bi,bj) = tLoc(d,j)
                0520              bbl_salt (d,j,bi,bj) = sLoc(d,j)
                0521             ENDIF
                0522             bbl_eta(r,j,bi,bj) = bblEta_r + dVol / rA(r,j,bi,bj)
                0523             bbl_theta(r,j,bi,bj) = ( dVol * bblTheta_d +
                0524      &           bblEta_r * rA(r,j,bi,bj) * bblTheta_r ) /
                0525      &           ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
                0526             bbl_salt(r,j,bi,bj)  = ( dVol * bblSalt_d  +
                0527      &           bblEta_r * rA(r,j,bi,bj) * bblSalt_r  ) /
                0528      &           ( bbl_eta(r,j,bi,bj) * rA(r,j,bi,bj) )
                0529            ENDIF
53b68b7823 Dimi*0530           ENDIF
                0531          ENDDO
                0532         ENDDO
                0533 
81ffbea13a Jean*0534 C--   end bi,bj loops.
                0535        ENDDO
                0536       ENDDO
                0537 
5f5cc705f4 Dimi*0538       CALL EXCH_XY_RL( bbl_eta      , myThid )
                0539       CALL EXCH_XY_RL( bbl_theta    , myThid )
                0540       CALL EXCH_XY_RL( bbl_salt     , myThid )
                0541       CALL EXCH_XY_RL( bbl_TendTheta, myThid )
                0542       CALL EXCH_XY_RL( bbl_TendSalt , myThid )
15338fa568 Dimi*0543 
53b68b7823 Dimi*0544       RETURN
                0545       END