#include "OBCS_OPTIONS.h"
#undef CHECK_BALANCE

C--   File obcs_calc_stevens.F: 
C--    Contents
C--    o OBCS_CALC_STEVENS
C--    o OBCS_STEVENS_CALC_TRACER_EAST
C--    o OBCS_STEVENS_CALC_TRACER_WEST
C--    o OBCS_STEVENS_CALC_TRACER_NORTH
C--    o OBCS_STEVENS_CALC_TRACER_SOUTH
C--    o OBCS_STEVENS_SAVE_TRACER

CBOP
C     !ROUTINE: OBCS_CALC_STEVENS
C     !INTERFACE:
      SUBROUTINE OBCS_CALC_STEVENS(
     I     futureTime, futureIter,
     I     myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE OBCS_CALC_STEVENS
C     | o Calculate future boundary data at open boundaries
C     |   at time = futureTime
C     |   from input data following Stevens(1990), and some
C     |   MOM3 legacy code
C     |
C     | o the code works like this
C     |  - the "barotropic" (= vertically averaged) velocity
C     |    normal to the boundary is assumed to be in
C     |    OBE/W/N/Su/v (normal) when this routine is entered
C     |  - the vertically averaged velocity is corrected
C     |    by the "baroclinic" (= deviation from vertically
C     |    averaged velocity) velocity to give a new OB?u/v;
C     |    the "barolinic" velocity is estimated from the previous
C     |    time step which makes this boundary condition depend on
C     |    a restart file. If OBCS_STEVENS_USE_INTERIOR_VELOCITY
C     |    is defined the velocity is simply copied from the model
C     |    interior to the boundary, thereby avoiding a restart
C     |    file or complicated reconstruction, but this solution
C     |    can give unexpected results.
C     |    (Note: in this context the terms barotropic and baroclinic
C     |    are MOM jargon and --- to my mind ---- should not be used)
C     |  - a wave phase speed is estimated from temporal and
C     |    horizontal variations of the tracer fields for each
C     |    tracer individually, this similar to Orlanski BCs,
C     |    but for simplicity the fields of the previous time step
C     |    are used, and the time derivative is estimated
C     |    independently of the time stepping procedure by simple
C     |    differencing
C     |  - velocity tangential to the boundary is always zero
C     |    (although this could be changed)
C     |  - a new tracer is computed from a local advection equation
C     |    with an upwind scheme: tracer from the interior is
C     |    advected out of the domain, and tracer from the boundary
C     |    is "advected" into the domain by a restoring mechanism
C     |  - for the advection equation only values from the
C     |    the current (not the updated) time level are used
C     |
C     *==========================================================*
C     | Feb, 2009: started by Martin Losch (Martin.Losch@awi.de)
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "DYNVARS.h"
#ifdef ALLOW_PTRACERS
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#include "PTRACERS_FIELDS.h"
#include "OBCS_PTRACERS.h"
#endif /* ALLOW_PTRACERS */
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
      _RL futureTime
      INTEGER futureIter
      INTEGER myThid

#ifdef ALLOW_OBCS_STEVENS

C     !LOCAL VARIABLES:
C     == Local variables ==
C     I,J,K        :: loop indices
C     msgBuf       :: Informational/error message buffer
C     uMer/vZonBar :: vertically averaged velocity at open boundary
C     drFBar       :: local depth for vertical average
C     uMer/vZonPri :: velocity anomalies applied to the open boundaries
C     gammat/s     :: restoring parameters (1./(T/SrelaxStevens - time scale))
C     auxillary variables
C     cflMer/Zon   :: ratio of grid spacing and time step
C     aFac         :: switch (0 or 1) that turns on advective contribution
C     gFacM/Z      :: switch (0 or 1) that turns on restoring boundary condition
C     pFac         :: switch that turns on/off phase velocity contribution
      INTEGER bi, bj
      INTEGER I, J, K
c     CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RL cflMer (1-OLy:sNy+OLy,Nr)
      _RL gFacM  (1-OLy:sNy+OLy,Nr)
      _RL uMerPri(Nr)
      _RL uMerBar
      _RL cflZon (1-OLx:sNx+OLx,Nr)
      _RL gFacZ  (1-OLx:sNx+OLx,Nr)
      _RL vZonPri(Nr)
      _RL vZonBar
      _RL drFBar
      _RL gammat, gammas, pFac, aFac
#ifdef ALLOW_PTRACERS
c     INTEGER iTracer
#endif /* ALLOW_PTRACERS */
#ifdef CHECK_BALANCE
      _RL uVelLoc, vVelLoc
      _RL vPhase
#endif
CEOP

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC_STEVENS',myThid)
#endif

      aFac   = 1. _d 0
      IF (.NOT. useStevensAdvection ) aFac   = 0. _d 0
      pFac   = 1. _d 0
      IF (.NOT. useStevensPhaseVel )  pFac   = 0. _d 0
      gammat = 0. _d 0
      IF (TrelaxStevens .GT. 0. _d 0) gammat = 1./TrelaxStevens
      gammas = 0. _d 0
      IF (SrelaxStevens .GT. 0. _d 0) gammas = 1./SrelaxStevens

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1
          ikey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

#ifdef ALLOW_OBCS_EAST

# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBEt(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBEs(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
# endif
        IF ( useStevensEast ) THEN
C     Eastern OB
#ifdef ALLOW_DEBUG
         IF (debugMode)
     &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: East',myThid)
#endif
C     compute vertical average and deviation from vertical
C     average for I_obe
         DO J=1-OLy,sNy+OLy
          I = OB_Ie(J,bi,bj)
          IF ( I.NE.OB_indexNone ) THEN
C     first initialize some fields
           drFbar  = 0. _d 0
           uMerBar = 0. _d 0
           DO K=1,Nr
            uMerPri(K) = 0. _d 0
           ENDDO
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            uMerBar = uMerBar + uVel(I-1,J,K,bi,bj)
#else
            uMerBar = uMerBar + OBEuStevens(J,K,bi,bj)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           *drF(K)* _hFacW(I,J,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacW(I,J,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            uMerPri(K) = (uVel(I-1,J,K,bi,bj)-uMerBar)
#else
            uMerPri(K) = (OBEuStevens(J,K,bi,bj)-uMerBar)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           * _maskW(I,J,K,bi,bj)
           ENDDO
C     vertical average of input field
           drFbar  = 0. _d 0
           uMerBar = 0. _d 0
           DO K=1,Nr
            uMerBar = uMerBar + OBEu(J,K,bi,bj)
     &           *drF(K)* _hFacW(I,J,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacW(I,J,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
C     Now the absolute velocity normal to the boundary is
C     uMerBar + uMerPri(K).
           DO K=1,Nr
            OBEu(J,K,bi,bj) = (uMerBar + uMerPri(K))
     &           * _maskW(I,J,K,bi,bj)
CML            OBEv(J,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
            OBEw(J,K,bi,bj)=0.
#endif
           ENDDO
          ENDIF
         ENDDO
#ifdef NONLIN_FRSURF
C     this is a bit of a hack
         IF ( nonlinFreeSurf.GT.0 ) THEN
          DO J=1-OLy,sNy+OLy
           I = OB_Ie(J,bi,bj)
           IF ( I.NE.OB_indexNone ) THEN
            OBEeta(J,bi,bj) = etaN(I-1,J,bi,bj)
           ENDIF
          ENDDO
         ENDIF
#endif /* NONLIN_FRSURF */
C     Next, we compute the phase speed correction, which depends on the
C     tracer!
         DO K=1,Nr
          DO J=1-OLy,sNy+OLy
           I = OB_Ie(J,bi,bj)
           IF ( I.NE.OB_indexNone ) THEN
            cflMer(J,K) = 0.5 _d 0 * _dxC(I-1,J,bi,bj)/dTtracerLev(K)
CML         gFacM(J,K)  = 0. _d 0
CML         IF ( uVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacM(J,K) = 1. _d 0
            gFacM(J,K)  = ABS(MIN(SIGN(1.D0,uVel(I,J,K,bi,bj)),0.D0))
           ELSE
            cflMer(J,K) = 0. _d 0
            gFacM (J,K) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
C     theta
         CALL OBCS_STEVENS_CALC_TRACER_EAST(
     U        OBEt, 
     I        OBEtStevens, theta, gammat, 
     I        uVel, cflMer, gFacM, pFac, aFac,
     I        OB_Ie, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     salinity
         CALL OBCS_STEVENS_CALC_TRACER_EAST(
     U        OBEs, 
     I        OBEsStevens, salt, gammas,
     I        uVel, cflMer, gFacM, pFac, aFac,
     I        OB_Ie, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     Template for passive tracers, requires work
CML#ifdef ALLOW_PTRACERS
CMLC     ptracers
CML         IF ( usePtracers ) THEN
CML          DO itracer = 1, PTRACERnumInUse
CML           CALL OBCS_STEVENS_CALC_TRACER_EAST(
CML     O          OBEptr       (1-OLy,1,1,1,iTracer), 
CML     I          OBEpStevens  (1-OLy,1,1,1,iTracer), 
CML     I          pTracer(1-OLx,1-OLy,1,1,1,iTracer), gammas,
CML     I          uVel, cflMer, gFacM, pFac, aFac,
CML     I          OB_Ie, OB_indexNone, bi, bj,
CML     I          futureTime, futureIter,
CML     I          myThid )
CML          ENDDO
CML         ENDIF
CML#endif /* ALLOW_PTRACERS */
C     IF ( useStevensEast ) THEN
        ENDIF
#endif /* ALLOW_OBCS_EAST */

C ------------------------------------------------------------------------------

#ifdef ALLOW_OBCS_WEST

# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBWt(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBWs(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
# endif
        IF ( useStevensWest ) THEN
C     Western OB
#ifdef ALLOW_DEBUG
         IF (debugMode)
     &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: West',myThid)
#endif
C     compute vertical average and deviation from vertical
C     average for I_obw+1
         DO J=1-OLy,sNy+OLy
          I = OB_Iw(J,bi,bj)
          IF ( I.NE.OB_indexNone ) THEN
C     first initialize some fields
           drFBar  = 0. _d 0
           uMerBar = 0. _d 0
           DO K=1,Nr
            uMerPri(K) = 0. _d 0
           ENDDO
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            uMerBar = uMerBar + uVel(I+2,J,K,bi,bj)
#else
            uMerBar = uMerBar + OBWuStevens(J,K,bi,bj)
#endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           *drF(K)* _hFacW(I+1,J,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacW(I+1,J,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            uMerPri(K) = (uVel(I+2,J,K,bi,bj)-uMerBar)
#else
            uMerPri(K) = (OBWuStevens(J,K,bi,bj)-uMerBar)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           * _maskW(I+1,J,K,bi,bj)
           ENDDO
C     vertical average of input field
           drFBar  = 0. _d 0
           uMerBar = 0. _d 0
           DO K=1,Nr
            uMerBar = uMerBar + OBWu(J,K,bi,bj)
     &           *drF(K)* _hFacW(I+1,J,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacW(I+1,J,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
C     Now the absolute velocity normal to the boundary is
C     uMerBar + uMerPri(K).
           DO K=1,Nr
            OBWu(J,K,bi,bj) = (uMerBar + uMerPri(K))
     &          * _maskW(I+1,J,K,bi,bj)
CML            OBWv(J,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
            OBWw(J,K,bi,bj)=0.
#endif
           ENDDO
          ENDIF
         ENDDO
#ifdef NONLIN_FRSURF
C     this is a bit of a hack
         IF ( nonlinFreeSurf.GT.0 ) THEN
          DO J=1-OLy,sNy+OLy
           I = OB_Iw(J,bi,bj)
           IF ( I.NE.OB_indexNone ) THEN
            OBWeta(J,bi,bj) = etaN(I+1,J,bi,bj)
           ENDIF
          ENDDO
         ENDIF
#endif /* NONLIN_FRSURF */
C     Next, we compute the phase speed correction, which depends on the
C     tracer!
         DO K=1,Nr
          DO J=1-OLy,sNy+OLy
           I = OB_Iw(J,bi,bj)
           IF ( I.NE.OB_indexNone ) THEN
            cflMer(J,K) = 0.5 _d 0 * _dxC(I+2,J,bi,bj)/dTtracerLev(K)
CML         gFacM = 0. _d 0
CML         IF ( uVel(I+1,J,K,bi,bj) .GT. 0. _d 0 ) gFacM = 1. _d 0
            gFacM(J,K)  = ABS(MAX(SIGN(1.D0,uVel(I+1,J,K,bi,bj)),0.D0))
           ELSE
            cflMer(J,K) = 0. _d 0
            gFacM (J,K) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
C     theta
         CALL OBCS_STEVENS_CALC_TRACER_WEST(
     U        OBWt, 
     I        OBWtStevens, theta, gammat, 
     I        uVel, cflMer, gFacM, pFac, aFac,
     I        OB_Iw, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     salinity
         CALL OBCS_STEVENS_CALC_TRACER_WEST(
     U        OBWs, 
     I        OBWsStevens, salt, gammas,
     I        uVel, cflMer, gFacM, pFac, aFac,
     I        OB_Iw, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     ptracers
C     IF ( useStevensWest ) THEN
        ENDIF
#endif /* ALLOW_OBCS_WEST */

C ------------------------------------------------------------------------------

#ifdef ALLOW_OBCS_NORTH
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBNt(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBNs(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
# endif
        IF ( useStevensNorth ) THEN
C         Northern OB
#ifdef ALLOW_DEBUG
         IF (debugMode)
     &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: North',myThid)
#endif
C     compute vertical average and deviation from vertical
C     average for J_obn
         DO I=1-OLx,sNx+OLx
          J = OB_Jn(I,bi,bj)
          IF ( J.NE.OB_indexNone ) THEN
C     first initialize some fields
           drFBar  = 0. _d 0
           vZonBar = 0. _d 0
           DO K=1,Nr
            vZonPri(K) = 0. _d 0
           ENDDO
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            vZonBar = vZonBar + vVel(I,J-1,K,bi,bj)
#else
            vZonBar = vZonBar + OBNvStevens(I,K,bi,bj)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           *drF(K)* _hFacS(I,J,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacS(I,J,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            vZonPri(K) = (vVel(I,J-1,K,bi,bj)-vZonBar)
#else
            vZonPri(K) = (OBNvStevens(I,K,bi,bj)-vZonBar)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           * _maskS(I,J,K,bi,bj)
           ENDDO
C     vertical average of input field
           drFBar  = 0. _d 0
           vZonBar = 0. _d 0
           DO K=1,Nr
            vZonBar = vZonBar + OBNv(I,K,bi,bj)
     &           *drF(K)* _hFacS(I,J,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacS(I,J,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
C     Now the absolute velocity normal to the boundary is
C     vZonBar + vZonPri(K).
           DO K=1,Nr
            OBNv(I,K,bi,bj) = (vZonBar + vZonPri(K))
     &           * _maskS(I,J,K,bi,bj)
CML            OBNu(I,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
            OBNw(I,K,bi,bj)=0.
#endif
           ENDDO
          ENDIF
         ENDDO
#ifdef NONLIN_FRSURF
C     this is a bit of a hack
         IF ( nonlinFreeSurf.GT.0 ) THEN
          DO I=1-OLx,sNx+OLx
           J = OB_Jn(I,bi,bj)
           IF ( J.NE.OB_indexNone ) THEN
            OBNeta(I,bi,bj) = etaN(I,J-1,bi,bj)
           ENDIF
          ENDDO
         ENDIF
#endif /* NONLIN_FRSURF */
C     Next, we compute the phase speed correction, which depends on the
C     tracer!
         DO K=1,Nr
          DO I=1-OLx,sNx+OLx
           J = OB_Jn(I,bi,bj)
           IF ( J.NE.OB_indexNone ) THEN
            cflZon(I,K) = 0.5 _d 0 * _dyC(I,J-1,bi,bj)/dTtracerLev(K)
CML         gFacZ(I,K) = 0. _d 0
CML         IF ( vVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacZ(I,K) = 1. _d 0
            gFacZ(I,K)  = ABS(MIN(SIGN(1.D0,vVel(I,J,K,bi,bj)),0.D0))
           ELSE
            cflZon(I,K) = 0. _d 0
            gFacZ (I,K) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
C     theta
         CALL OBCS_STEVENS_CALC_TRACER_NORTH(
     U        OBNt, 
     I        OBNtStevens, theta, gammat,
     I        vVel, cflZon, gFacZ, pFac, aFac, 
     I        OB_Jn, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     salinity
         CALL OBCS_STEVENS_CALC_TRACER_NORTH(
     U        OBNs, 
     I        OBNsStevens, salt, gammas,
     I        vVel, cflZon, gFacZ, pFac, aFac, 
     I        OB_Jn, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     ptracers
C     IF ( useStevensNorth ) THEN
        ENDIF
#endif /* ALLOW_OBCS_NORTH */

C ------------------------------------------------------------------------------

#ifdef ALLOW_OBCS_SOUTH

# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE OBSt(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE OBSs(:,:,bi,bj)  = comlev1_bibj, key=ikey, byte=isbyte
# endif
        IF ( useStevensSouth ) THEN
C         Southern OB
#ifdef ALLOW_DEBUG
         IF (debugMode)
     &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: South',myThid)
#endif
C     compute vertical average and deviation from vertical
C     average for J_obs+1
         DO I=1-OLx,sNx+OLx
          J = OB_Js(I,bi,bj)
          IF ( J.NE.OB_indexNone ) THEN
C     first initialize some fields
           drFBar  = 0. _d 0
           vZonBar = 0. _d 0
           DO K=1,Nr
            vZonPri(K) = 0. _d 0
           ENDDO
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            vZonBar = vZonBar + vVel(I,J+2,K,bi,bj)
#else
            vZonBar = vZonBar + OBSvStevens(I,K,bi,bj)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           *drF(K)* _hFacS(I,J+1,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacS(I,J+1,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
           DO K=1,Nr
#ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
            vZonPri(K) = (vVel(I,J+2,K,bi,bj)-vZonBar)
#else
            vZonPri(K) = (OBSvStevens(I,K,bi,bj)-vZonBar)
#endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
     &           * _maskS(I,J+1,K,bi,bj)
           ENDDO
C     vertical average of input field
           drFBar  = 0. _d 0
           vZonBar = 0. _d 0
           DO K=1,Nr
            vZonBar = vZonBar + OBSv(I,K,bi,bj)
     &           *drF(K)* _hFacS(I,J+1,K,bi,bj)
            drFBar = drFBar + drF(K)* _hFacS(I,J+1,K,bi,bj)
           ENDDO
           IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
C     Now the absolute velocity normal to the boundary is
C     vZonBar + vZonPri(K).
           DO K=1,Nr
            OBSv(I,K,bi,bj) = (vZonBar + vZonPri(K))
     &          * _maskS(I,J+1,K,bi,bj)
CML            OBSu(I,K,bi,bj) = 0. _d 0
#ifdef ALLOW_NONHYDROSTATIC
            OBSw(I,K,bi,bj)=0.
#endif
           ENDDO
          ENDIF
         ENDDO
#ifdef NONLIN_FRSURF
C     this is a bit of a hack
         IF ( nonlinFreeSurf.GT.0 ) THEN
          DO I=1-OLx,sNx+OLx
           J = OB_Js(I,bi,bj)
           IF ( J.NE.OB_indexNone ) THEN
            OBSeta(I,bi,bj) = etaN(I,J+1,bi,bj)
           ENDIF
          ENDDO
         ENDIF
#endif /* NONLIN_FRSURF */
C     Next, we compute the phase speed correction, which depends on the
C     tracer!
         DO K=1,Nr
          DO I=1-OLx,sNx+OLx
           J = OB_Js(I,bi,bj)
           IF ( J.NE.OB_indexNone ) THEN
            cflZon(I,K) = 0.5 _d 0 * _dyC(I,J+2,bi,bj)/dTtracerLev(K)
CML         gFacZ = 0. _d 0
CML         IF ( vVel(I,J+1,K,bi,bj) .GT. 0. _d 0 ) gFacZ = 1. _d 0
            gFacZ(I,K)  = ABS(MAX(SIGN(1.D0,vVel(I,J+1,K,bi,bj)),0.D0))
           ELSE
            cflZon(I,K) = 0. _d 0
            gFacZ (I,K) = 0. _d 0
           ENDIF
          ENDDO
         ENDDO
C     theta
         CALL OBCS_STEVENS_CALC_TRACER_SOUTH(
     U        OBSt, 
     I        OBStStevens, theta, gammat,
     I        vVel, cflZon, gFacZ, pFac, aFac, 
     I        OB_Js, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     salinity
         CALL OBCS_STEVENS_CALC_TRACER_SOUTH(
     U        OBSs, 
     I        OBSsStevens, salt, gammas,
     I        vVel, cflZon, gFacZ, pFac, aFac, 
     I        OB_Js, OB_indexNone, bi, bj,
     I        futureTime, futureIter,
     I        myThid )
C     ptracers
C     IF ( useStevensSouth ) THEN
        ENDIF
#endif /* ALLOW_OBCS_SOUTH */

C     end bi/bj-loops
       ENDDO
      ENDDO

C     save the tracer fields of the previous time step for the next time step
      CALL OBCS_STEVENS_SAVE_TRACERS(
     I     futureTime, futureIter,
     I     myThid )
C ------------------------------------------------------------------------------

#ifdef CHECK_BALANCE
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        uPhase=0.
        vPhase=0.
        uVelLoc = 0.
        DO J=1-OLy,sNy+OLy
         uMerBar=0. _d 0
         DO K=1,Nr
          I = OB_Ie(J,bi,bj)
          IF ( I.EQ.OB_indexNone ) I = 1
          uPhase = uPhase + OBEu(J,K,bi,bj)
     &         *drF(k)* _hFacW(I,J,K,bi,bj)*dyG(I,J,bi,bj)
          I = OB_Iw(J,bi,bj)
          IF ( I.EQ.OB_indexNone ) I = 1
          vPhase = vPhase + OBWu(J,K,bi,bj)
     &         *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
CML          uVelLoc = uVelLoc + uMerPri(J,K)
CML     &         *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
CML          uMerBar(J)=uMerBar(J) + uMerPri(J,K)
CML     &         *drF(k)* _hFacW(I+1,J,K,bi,bj)
         ENDDO
CML         print *, 'ml-obcs: uBar = ', j,uMerBar(J)
        ENDDO
C     end bi/bj-loops
       ENDDO
      ENDDO
      _GLOBAL_SUM_RL( uPhase, myThid )
      _GLOBAL_SUM_RL( vPhase, myThid )
CML      _GLOBAL_SUM_RL( uVelLoc, myThid )
      print *, 'ml-obcs: OBE  = ',  uPhase*1 _d -6, ' Sv'
      print *, 'ml-obcs: OBW  = ',  vPhase*1 _d -6, ' Sv'
CML      print *, 'ml-obcs: OBWp = ', uVelLoc*1 _d -6, ' Sv'
#endif /* CHECK_BALANCE */

#ifdef ALLOW_DEBUG
         IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC_STEVENS',myThid)
#endif

#endif /* ALLOW_OBCS_STEVENS */
      RETURN
      END

CBOP
C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_EAST
C     !INTERFACE:
      SUBROUTINE OBCS_STEVENS_CALC_TRACER_EAST(
     U     OBEf, 
     I     OBE_Stevens, tracer, gammaf,
     I     uVel, cflMer, gFacM, pFac, aFac, 
     I     OB_I, OB_indexNone, bi, bj,
     I     futureTime, futureIter,
     I     myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_EAST
C     | Calculate tracer value at the eastern OB location
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C    myThid    :: my Thread Id number
C    bi, bj    :: indices of current tile
      _RL futureTime
      INTEGER futureIter
      INTEGER myThid
      INTEGER bi, bj
      INTEGER OB_indexNone
      INTEGER OB_I             (1-OLy:sNy+OLy,nSx,nSy)
      _RL cflMer               (1-OLy:sNy+OLy,Nr)
      _RL gFacM                (1-OLy:sNy+OLy,Nr)
      _RL gammaf, pFac, aFac
      _RL OBEf                 (1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL OBE_Stevens          (1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL uVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

#ifdef ALLOW_OBCS_STEVENS
C     !LOCAL VARIABLES:
C     i,j,k        :: loop indices
C     uPhase       :: estimate of phase velocity for radiation condition
C     dtracSpace   :: horizontal difference of tracer
C     dtracTime    :: temporal difference of tracer
      INTEGER i,j,k
      _RL uPhase
      _RL dtracSpace
      _RL dTracTime
CEOP
      DO K=1,Nr
       DO J=1-OLy,sNy+OLy
        I = OB_I(J,bi,bj)
        IF ( I.NE.OB_indexNone ) THEN
         dTracSpace = (tracer(I-1,J,K,bi,bj)-tracer(I-2,J,K,bi,bj))
     &        * _maskW(I-1,J,K,bi,bj)
         dTracTime  = (tracer(I-1,J,K,bi,bj)-OBE_Stevens(J,K,bi,bj))
         uPhase = cflMer(J,K) * pFac
         IF ( dTracSpace .NE. 0. _d 0 ) THEN
          uPhase = MIN( cflMer(J,K),
     &         MAX( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
     &         ) * pFac
         ENDIF
C     Compute the tracer tendency here, the tracer will be updated
C     with a simple Euler forward step in S/R obcs_apply_ts
         OBEf(J,K,bi,bj) = _maskW(I,J,K,bi,bj) * (
     &        - ( aFac*MAX(0.D0,uVel(I,J,K,bi,bj)) + uPhase )
     &        *(tracer(I,J,K,bi,bj)-tracer(I-1,J,K,bi,bj))
     &        * _recip_dxC(I,J,bi,bj)
     &        - gFacM(J,K) * gammaf
     &        * (tracer(I,J,K,bi,bj)-OBEf(J,K,bi,bj)) )
        ENDIF
       ENDDO
      ENDDO

#endif /* ALLOW_OBCS_STEVENS */
      RETURN
      END

CBOP
C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_WEST
C     !INTERFACE:
      SUBROUTINE OBCS_STEVENS_CALC_TRACER_WEST(
     U     OBWf, 
     I     OBW_Stevens, tracer, gammaf,
     I     uVel, cflMer, gFacM, pFac, aFac, 
     I     OB_I, OB_indexNone, bi, bj,
     I     futureTime, futureIter,
     I     myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_WEST
C     | Calculate tracer value at the western OB location
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C    myThid    :: my Thread Id number
C    bi, bj    :: indices of current tile
      _RL futureTime
      INTEGER futureIter
      INTEGER myThid
      INTEGER bi, bj
      INTEGER OB_indexNone
      INTEGER OB_I             (1-OLy:sNy+OLy,nSx,nSy)
      _RL cflMer               (1-OLy:sNy+OLy,Nr)
      _RL gFacM                (1-OLy:sNy+OLy,Nr)
      _RL gammaf, pFac, aFac
      _RL OBWf                 (1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL OBW_Stevens          (1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL uVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

#ifdef ALLOW_OBCS_STEVENS
C     !LOCAL VARIABLES:
C     i,j,k        :: loop indices
C     uPhase       :: estimate of phase velocity for radiation condition
C     dtracSpace   :: horizontal difference of tracer
C     dtracTime    :: temporal difference of tracer
      INTEGER i,j,k
      _RL uPhase
      _RL dtracSpace
      _RL dTracTime
CEOP

      DO K=1,Nr
       DO J=1-OLy,sNy+OLy
        I = OB_I(J,bi,bj)
        IF ( I.NE.OB_indexNone ) THEN
         dTracSpace = (tracer(I+2,J,K,bi,bj)-tracer(I+1,J,K,bi,bj))
     &        * _maskW(I+2,J,K,bi,bj)
         dTracTime  = (tracer(I+1,J,K,bi,bj)-OBW_Stevens(J,K,bi,bj))
         uPhase = -cflMer(J,K) * pFac
         IF ( dTracSpace .NE. 0. _d 0 ) THEN
          uPhase = MAX( -cflMer(J,K),
     &         MIN( 0.D0, -cflMer(J,K)*dTracTime/dTracSpace )
     &         ) * pFac
         ENDIF
C     Compute the tracer tendency here, the tracer will be updated
C     with a simple Euler forward step in S/R obcs_apply_ts
         OBWf(J,K,bi,bj) = _maskW(I+1,J,K,bi,bj) * (
     &        - ( aFac*MIN(0.D0,uVel(I+1,J,K,bi,bj)) + uPhase )
     &        *(tracer(I+1,J,K,bi,bj)-tracer(I,J,K,bi,bj))
     &        * _recip_dxC(I+1,J,bi,bj)
     &        - gFacM(J,K) * gammaf
     &        * (tracer(I,J,K,bi,bj)-OBWf(J,K,bi,bj)) )
        ENDIF
       ENDDO
      ENDDO

#endif /* ALLOW_OBCS_STEVENS */
      RETURN
      END

CBOP
C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_NORTH
C     !INTERFACE:
      SUBROUTINE OBCS_STEVENS_CALC_TRACER_NORTH(
     U     OBNf, 
     I     OBN_Stevens, tracer, gammaf,
     I     vVel, cflZon, gFacZ, pFac, aFac, 
     I     OB_J, OB_indexNone, bi, bj,
     I     futureTime, futureIter,
     I     myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_NORTH
C     | Calculate tracer value at the northern OB location
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C    myThid    :: my Thread Id number
C    bi, bj    :: indices of current tile
      _RL futureTime
      INTEGER futureIter
      INTEGER myThid
      INTEGER bi, bj
      INTEGER OB_indexNone
      INTEGER OB_J             (1-OLx:sNx+OLx,nSx,nSy)
      _RL cflZon               (1-OLx:sNx+OLx,Nr)
      _RL gFacZ                (1-OLx:sNx+OLx,Nr)
      _RL gammaf, pFac, aFac
      _RL OBNf                 (1-Olx:sNx+Olx,Nr,nSx,nSy)
      _RL OBN_Stevens          (1-Olx:sNx+Olx,Nr,nSx,nSy)
      _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL vVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

#ifdef ALLOW_OBCS_STEVENS
C     !LOCAL VARIABLES:
C     i,j,k        :: loop indices
C     vPhase       :: estimate of phase velocity for radiation condition
C     dtracSpace   :: horizontal difference of tracer
C     dtracTime    :: temporal difference of tracer
      INTEGER i,j,k
      _RL vPhase
      _RL dtracSpace
      _RL dTracTime
CEOP
      DO K=1,Nr
       DO I=1-OLx,sNx+OLx
        J = OB_J(I,bi,bj)
        IF ( J.NE.OB_indexNone ) THEN
C     Theta first:
         dTracSpace = (tracer(I,J-1,K,bi,bj)-tracer(I,J-2,K,bi,bj))
     &        * _maskS(I,J-1,K,bi,bj)
         dTracTime  = (tracer(I,J-1,K,bi,bj)-OBN_Stevens(I,K,bi,bj))
         vPhase = cflZon(I,K) * pFac
         IF ( dTracSpace .NE. 0. _d 0 ) THEN
          vPhase = MIN( cflZon(I,K),
     &         MAX( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
     &         ) * pFac
         ENDIF
C     Compute the tracer tendency here, the tracer will be updated
C     with a simple Euler forward step in S/R obcs_apply_ts
         OBNf(I,K,bi,bj) = _maskS(I,J,K,bi,bj) * (
     &        - ( aFac*MAX(0.D0,vVel(I,J,K,bi,bj)) + vPhase )
     &        *(tracer(I,J,K,bi,bj)-tracer(I,J-1,K,bi,bj))
     &        * _recip_dyC(I,J,bi,bj)
     &        - gFacZ(I,K) * gammaf
     &        * (tracer(I,J,K,bi,bj)-OBNf(I,K,bi,bj)) )
        ENDIF
       ENDDO
      ENDDO

#endif /* ALLOW_OBCS_STEVENS */
      RETURN
      END

CBOP
C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_SOUTH
C     !INTERFACE:
      SUBROUTINE OBCS_STEVENS_CALC_TRACER_SOUTH(
     U     OBSf, 
     I     OBS_Stevens, tracer, gammaf,
     I     vVel, cflZon, gFacZ, pFac, aFac, 
     I     OB_J, OB_indexNone, bi, bj,
     I     futureTime, futureIter,
     I     myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_SOUTH
C     | Calculate tracer value at the southern OB location
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C    myThid    :: my Thread Id number
C    bi, bj    :: indices of current tile
      _RL futureTime
      INTEGER futureIter
      INTEGER myThid
      INTEGER bi, bj
      INTEGER OB_indexNone
      INTEGER OB_J             (1-OLx:sNx+OLx,nSx,nSy)
      _RL cflZon               (1-OLx:sNx+OLx,Nr)
      _RL gFacZ                (1-OLx:sNx+OLx,Nr)
      _RL gammaf, pFac, aFac
      _RL OBSf                 (1-Olx:sNx+Olx,Nr,nSx,nSy)
      _RL OBS_Stevens          (1-Olx:sNx+Olx,Nr,nSx,nSy)
      _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL vVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)

#ifdef ALLOW_OBCS_STEVENS
C     !LOCAL VARIABLES:
C     i,j,k        :: loop indices
C     vPhase       :: estimate of phase velocity for radiation condition
C     dtracSpace   :: horizontal difference of tracer
C     dtracTime    :: temporal difference of tracer
      INTEGER i,j,k
      _RL vPhase
      _RL dtracSpace
      _RL dTracTime
CEOP
      DO K=1,Nr
       DO I=1-OLx,sNx+OLx
        J = OB_J(I,bi,bj)
        IF ( J.NE.OB_indexNone ) THEN
         dTracSpace = (tracer(I,J+2,K,bi,bj)-tracer(I,J+1,K,bi,bj))
     &        * _maskS(I,J+2,K,bi,bj)
         dTracTime  = (tracer(I,J+1,K,bi,bj)-OBS_Stevens(I,K,bi,bj))
         vPhase = -cflZon(I,K) * pFac
         IF ( dTracSpace .NE. 0. _d 0 ) THEN
          vPhase = MAX( -cflZon(I,K),
     &         MIN( 0.D0, -cflZon(I,K)*dTracTime/dTracSpace )
     &         ) * pFac
         ENDIF
C     Compute the tracer tendency here, the tracer will be updated
C     with a simple Euler forward step in S/R obcs_apply_ts
         OBSf(I,K,bi,bj) = _maskS(I,J+1,K,bi,bj) * (
     &        - ( aFac*MIN(0.D0,vVel(I,J+1,K,bi,bj)) + vPhase )
     &        *(tracer(I,J+1,K,bi,bj)-tracer(I,J,K,bi,bj))
     &        * _recip_dyC(I,J+1,bi,bj)
     &        - gFacZ(I,K) * gammaf
     &        * (tracer(I,J,K,bi,bj)-OBSf(I,K,bi,bj)) )
        ENDIF
       ENDDO
      ENDDO

#endif /* ALLOW_OBCS_STEVENS */
      RETURN
      END

CBOP
C     !ROUTINE: OBCS_STEVENS_SAVE_TRACERS
C     !INTERFACE:
      SUBROUTINE OBCS_STEVENS_SAVE_TRACERS(
     I     futureTime, futureIter,
     I     myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE OBCS_STEVENS_SAVE_TRACERS
C     | Save tracers (of previous time step) at the OB location
C     | to be used in the next time step for Stevens boundary
C     | conditions
C     *==========================================================*
C     \ev
C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "DYNVARS.h"
CML#ifdef ALLOW_PTRACERS
CML#include "PTRACERS_SIZE.h"
CML#include "PTRACERS_PARAMS.h"
CML#include "PTRACERS_FIELDS.h"
CML#include "OBCS_PTRACERS.h"
CML#endif /* ALLOW_PTRACERS */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C    myThid    :: my Thread Id number
      _RL futureTime
      INTEGER futureIter
      INTEGER myThid

#ifdef ALLOW_OBCS_STEVENS
C     !LOCAL VARIABLES:
C    bi, bj    :: indices of current tile
C    i,j,k     :: loop indices
C    Iobc,Jobc :: position-index of open boundary
      INTEGER bi,bj,i,j,k,Iobc,Jobc
CEOP

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef ALLOW_OBCS_NORTH
        IF ( tileHasOBN(bi,bj) .AND. useStevensNorth ) THEN
C Northern boundary
         DO i=1-OLx,sNx+OLx
          Jobc = OB_Jn(i,bi,bj)
          IF ( Jobc.NE.OB_indexNone ) THEN
           DO k = 1,Nr
            OBNtStevens(i,k,bi,bj) = theta(i,Jobc-1,k,bi,bj)
     &           *maskC(i,Jobc+1,k,bi,bj)
            OBNsStevens(i,k,bi,bj) = salt(i,Jobc-1,k,bi,bj)
     &           *maskC(i,Jobc+1,k,bi,bj)
           ENDDO
          ENDIF
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS_NORTH */
#ifdef ALLOW_OBCS_SOUTH
        IF ( tileHasOBS(bi,bj) .AND. useStevensSouth ) THEN
C Southern boundary
         DO i=1-OLx,sNx+OLx
          Jobc = OB_Js(i,bi,bj)
          IF ( Jobc.NE.OB_indexNone ) THEN
           DO k = 1,Nr
            OBStStevens(i,k,bi,bj) = theta(i,Jobc+1,k,bi,bj)
     &           *maskC(i,Jobc+1,k,bi,bj)
            OBSsStevens(i,k,bi,bj) = salt(i,Jobc+1,k,bi,bj)
     &           *maskC(i,Jobc+1,k,bi,bj)
           ENDDO
          ENDIF
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS_SOUTH */
#ifdef ALLOW_OBCS_EAST
        IF ( tileHasOBE(bi,bj) .AND. useStevensEast ) THEN
C Eastern boundary
         DO j=1-OLy,sNy+OLy
          Iobc = OB_Ie(j,bi,bj)
          IF ( Iobc.NE.OB_indexNone ) THEN
           DO k = 1,Nr
            OBEtStevens(j,k,bi,bj) = theta(Iobc-1,j,k,bi,bj)
     &           *maskC(Iobc-1,j,k,bi,bj)
            OBEsStevens(j,k,bi,bj) = salt(Iobc-1,j,k,bi,bj)
     &           *maskC(Iobc-1,j,k,bi,bj)
           ENDDO
          ENDIF
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS_EAST */
#ifdef ALLOW_OBCS_WEST
        IF ( tileHasOBW(bi,bj) .AND. useStevensWest ) THEN
C Western boundary
         DO j=1-OLy,sNy+OLy
          Iobc = OB_Iw(j,bi,bj)
          IF ( Iobc.NE.OB_indexNone ) THEN
           DO k = 1,Nr
            OBWtStevens(j,k,bi,bj) = theta(Iobc+1,j,k,bi,bj)
     &           *maskC(Iobc+1,j,k,bi,bj)
            OBWsStevens(j,k,bi,bj) = salt(Iobc+1,j,k,bi,bj)
     &           *maskC(Iobc+1,j,k,bi,bj)
           ENDDO
          ENDIF
         ENDDO
        ENDIF
#endif /* ALLOW_OBCS_WEST */
       ENDDO
      ENDDO
#endif /* ALLOW_OBCS_STEVENS */
      RETURN
      END