C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_calc_stevens.F,v 1.14 2015/01/12 16:06:38 mlosch Exp $
C $Name:  $

#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