pkg/obcs/obcs_calc_stevens.F

b2cb1ccb9a Mart*0001 #include "OBCS_OPTIONS.h"
4e4ad91a39 Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
eb8aff1784 Jean*0005 #undef CHECK_BALANCE
b2cb1ccb9a Mart*0006
4e4ad91a39 Jean*0007 C--   File obcs_calc_stevens.F:
0daa63ec35 Mart*0008 C--    Contents
                0009 C--    o OBCS_CALC_STEVENS
                0010 C--    o OBCS_STEVENS_CALC_TRACER_EAST
                0011 C--    o OBCS_STEVENS_CALC_TRACER_WEST
                0012 C--    o OBCS_STEVENS_CALC_TRACER_NORTH
                0013 C--    o OBCS_STEVENS_CALC_TRACER_SOUTH
                0014 C--    o OBCS_STEVENS_SAVE_TRACER
                0015
b2cb1ccb9a Mart*0016 CBOP
                0017 C     !ROUTINE: OBCS_CALC_STEVENS
                0018 C     !INTERFACE:
eb8aff1784 Jean*0019       SUBROUTINE OBCS_CALC_STEVENS(
b2cb1ccb9a Mart*0020      I     futureTime, futureIter,
                0021      I     myThid )
2c64c37c71 Mart*0022 C     !DESCRIPTION: \bv
eb8aff1784 Jean*0023 C     *==========================================================*
b2cb1ccb9a Mart*0024 C     | SUBROUTINE OBCS_CALC_STEVENS
                0025 C     | o Calculate future boundary data at open boundaries
eb8aff1784 Jean*0026 C     |   at time = futureTime
                0027 C     |   from input data following Stevens(1990), and some
b2cb1ccb9a Mart*0028 C     |   MOM3 legacy code
                0029 C     |
                0030 C     | o the code works like this
eb8aff1784 Jean*0031 C     |  - the "barotropic" (= vertically averaged) velocity
                0032 C     |    normal to the boundary is assumed to be in
b2cb1ccb9a Mart*0033 C     |    OBE/W/N/Su/v (normal) when this routine is entered
eb8aff1784 Jean*0034 C     |  - the vertically averaged velocity is corrected
                0035 C     |    by the "baroclinic" (= deviation from vertically
                0036 C     |    averaged velocity) velocity to give a new OB?u/v;
2c64c37c71 Mart*0037 C     |    the "barolinic" velocity is estimated from the previous
                0038 C     |    time step which makes this boundary condition depend on
                0039 C     |    a restart file. If OBCS_STEVENS_USE_INTERIOR_VELOCITY
                0040 C     |    is defined the velocity is simply copied from the model
                0041 C     |    interior to the boundary, thereby avoiding a restart
                0042 C     |    file or complicated reconstruction, but this solution
                0043 C     |    can give unexpected results.
b2cb1ccb9a Mart*0044 C     |    (Note: in this context the terms barotropic and baroclinic
                0045 C     |    are MOM jargon and --- to my mind ---- should not be used)
                0046 C     |  - a wave phase speed is estimated from temporal and
                0047 C     |    horizontal variations of the tracer fields for each
b1353a35da Mart*0048 C     |    tracer individually, this similar to Orlanski BCs,
74019f026d Jean*0049 C     |    but for simplicity the fields of the previous time step
                0050 C     |    are used, and the time derivative is estimated
b1353a35da Mart*0051 C     |    independently of the time stepping procedure by simple
                0052 C     |    differencing
b2cb1ccb9a Mart*0053 C     |  - velocity tangential to the boundary is always zero
                0054 C     |    (although this could be changed)
0daa63ec35 Mart*0055 C     |  - a new tracer is computed from a local advection equation
b2cb1ccb9a Mart*0056 C     |    with an upwind scheme: tracer from the interior is
                0057 C     |    advected out of the domain, and tracer from the boundary
                0058 C     |    is "advected" into the domain by a restoring mechanism
eb8aff1784 Jean*0059 C     |  - for the advection equation only values from the
b2cb1ccb9a Mart*0060 C     |    the current (not the updated) time level are used
eb8aff1784 Jean*0061 C     |
                0062 C     *==========================================================*
b2cb1ccb9a Mart*0063 C     | Feb, 2009: started by Martin Losch (Martin.Losch@awi.de)
eb8aff1784 Jean*0064 C     *==========================================================*
2c64c37c71 Mart*0065 C     \ev
b2cb1ccb9a Mart*0066 C     !USES:
eb8aff1784 Jean*0067       IMPLICIT NONE
b2cb1ccb9a Mart*0068 C     === Global variables ===
                0069 #include "SIZE.h"
                0070 #include "EEPARAMS.h"
                0071 #include "PARAMS.h"
eb8aff1784 Jean*0072 #include "GRID.h"
9b4f2a04e2 Jean*0073 #include "OBCS_PARAMS.h"
                0074 #include "OBCS_GRID.h"
                0075 #include "OBCS_FIELDS.h"
b2cb1ccb9a Mart*0076 #include "DYNVARS.h"
                0077 #ifdef ALLOW_PTRACERS
                0078 #include "PTRACERS_SIZE.h"
                0079 #include "PTRACERS_PARAMS.h"
                0080 #include "PTRACERS_FIELDS.h"
                0081 #include "OBCS_PTRACERS.h"
                0082 #endif /* ALLOW_PTRACERS */
38b71af82a Patr*0083 #ifdef ALLOW_AUTODIFF_TAMC
4e4ad91a39 Jean*0084 # include "tamc.h"
38b71af82a Patr*0085 #endif /* ALLOW_AUTODIFF_TAMC */
b2cb1ccb9a Mart*0086
                0087 C     !INPUT/OUTPUT PARAMETERS:
                0088 C     == Routine arguments ==
                0089       _RL futureTime
eb8aff1784 Jean*0090       INTEGER futureIter
b2cb1ccb9a Mart*0091       INTEGER myThid
                0092
                0093 #ifdef ALLOW_OBCS_STEVENS
                0094
                0095 C     !LOCAL VARIABLES:
                0096 C     == Local variables ==
eb8aff1784 Jean*0097 C     I,J,K        :: loop indices
                0098 C     msgBuf       :: Informational/error message buffer
                0099 C     uMer/vZonBar :: vertically averaged velocity at open boundary
                0100 C     drFBar       :: local depth for vertical average
                0101 C     uMer/vZonPri :: velocity anomalies applied to the open boundaries
                0102 C     gammat/s     :: restoring parameters (1./(T/SrelaxStevens - time scale))
b2cb1ccb9a Mart*0103 C     auxillary variables
b1353a35da Mart*0104 C     cflMer/Zon   :: ratio of grid spacing and time step
eb8aff1784 Jean*0105 C     aFac         :: switch (0 or 1) that turns on advective contribution
2c64c37c71 Mart*0106 C     gFacM/Z      :: switch (0 or 1) that turns on restoring boundary condition
eb8aff1784 Jean*0107 C     pFac         :: switch that turns on/off phase velocity contribution
b2cb1ccb9a Mart*0108       INTEGER bi, bj
4e4ad91a39 Jean*0109       INTEGER i, j, k
eb8aff1784 Jean*0110 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
2c64c37c71 Mart*0111       _RL cflMer (1-OLy:sNy+OLy,Nr)
                0112       _RL gFacM  (1-OLy:sNy+OLy,Nr)
e37ac28e89 Mart*0113       _RL uMerPri(Nr)
                0114       _RL uMerBar
16a91a6f85 Mart*0115       _RL cflZon (1-OLx:sNx+OLx,Nr)
3cf13754ce Mart*0116       _RL gFacZ  (1-OLx:sNx+OLx,Nr)
e37ac28e89 Mart*0117       _RL vZonPri(Nr)
                0118       _RL vZonBar
74019f026d Jean*0119       _RL drFBar
2c64c37c71 Mart*0120       _RL gammat, gammas, pFac, aFac
b2cb1ccb9a Mart*0121 #ifdef ALLOW_PTRACERS
eb8aff1784 Jean*0122 c     INTEGER iTracer
b2cb1ccb9a Mart*0123 #endif /* ALLOW_PTRACERS */
4e4ad91a39 Jean*0124 #ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0125 C     tkey :: tape key (depends on tiles)
                0126       INTEGER tkey
4e4ad91a39 Jean*0127 #endif /* ALLOW_AUTODIFF_TAMC */
eb8aff1784 Jean*0128 #ifdef CHECK_BALANCE
                0129       _RL uVelLoc, vVelLoc
                0130       _RL vPhase
                0131 #endif
b2cb1ccb9a Mart*0132 CEOP
                0133
                0134 #ifdef ALLOW_DEBUG
eb8aff1784 Jean*0135       IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC_STEVENS',myThid)
b2cb1ccb9a Mart*0136 #endif
                0137
                0138       aFac   = 1. _d 0
                0139       IF (.NOT. useStevensAdvection ) aFac   = 0. _d 0
                0140       pFac   = 1. _d 0
                0141       IF (.NOT. useStevensPhaseVel )  pFac   = 0. _d 0
                0142       gammat = 0. _d 0
                0143       IF (TrelaxStevens .GT. 0. _d 0) gammat = 1./TrelaxStevens
                0144       gammas = 0. _d 0
                0145       IF (SrelaxStevens .GT. 0. _d 0) gammas = 1./SrelaxStevens
                0146
                0147       DO bj=myByLo(myThid),myByHi(myThid)
                0148        DO bi=myBxLo(myThid),myBxHi(myThid)
eb8aff1784 Jean*0149
38b71af82a Patr*0150 #ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0151         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0152 # ifdef ALLOW_OBCS_EAST
                0153 CADJ STORE OBEt(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0154 CADJ STORE OBEs(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0155 CADJ STORE OBEu(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0156 # endif
                0157 # ifdef ALLOW_OBCS_WEST
                0158 CADJ STORE OBWt(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0159 CADJ STORE OBWs(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0160 CADJ STORE OBWu(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0161 # endif
                0162 # ifdef ALLOW_OBCS_NORTH
                0163 CADJ STORE OBNt(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0164 CADJ STORE OBNs(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0165 CADJ STORE OBNv(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0166 # endif
                0167 # ifdef ALLOW_OBCS_SOUTH
                0168 CADJ STORE OBSt(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0169 CADJ STORE OBSs(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0170 CADJ STORE OBSv(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0171 # endif
38b71af82a Patr*0172 #endif /* ALLOW_AUTODIFF_TAMC */
                0173
b2cb1ccb9a Mart*0174 #ifdef ALLOW_OBCS_EAST
                0175         IF ( useStevensEast ) THEN
                0176 C     Eastern OB
                0177 #ifdef ALLOW_DEBUG
eb8aff1784 Jean*0178          IF (debugMode)
b2cb1ccb9a Mart*0179      &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: East',myThid)
                0180 #endif
eb8aff1784 Jean*0181 C     compute vertical average and deviation from vertical
9f68a71aa4 Mart*0182 C     average for I_obe
4e4ad91a39 Jean*0183          DO j=1-OLy,sNy+OLy
                0184           i = OB_Ie(j,bi,bj)
                0185           IF ( i.NE.OB_indexNone ) THEN
e37ac28e89 Mart*0186 C     first initialize some fields
                0187            drFbar  = 0. _d 0
                0188            uMerBar = 0. _d 0
4e4ad91a39 Jean*0189            DO k=1,Nr
                0190             uMerPri(k) = 0. _d 0
e37ac28e89 Mart*0191            ENDDO
4e4ad91a39 Jean*0192            DO k=1,Nr
2c64c37c71 Mart*0193 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0194             uMerBar = uMerBar + uVel(i-1,j,k,bi,bj)
2c64c37c71 Mart*0195 #else
4e4ad91a39 Jean*0196             uMerBar = uMerBar + OBEuStevens(j,k,bi,bj)
2c64c37c71 Mart*0197 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0198      &           *drF(k)* _hFacW(i,j,k,bi,bj)
                0199             drFBar = drFBar + drF(k)* _hFacW(i,j,k,bi,bj)
9f68a71aa4 Mart*0200            ENDDO
e37ac28e89 Mart*0201            IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
4e4ad91a39 Jean*0202            DO k=1,Nr
2c64c37c71 Mart*0203 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0204             uMerPri(k) = (uVel(i-1,j,k,bi,bj)-uMerBar)
74019f026d Jean*0205 #else
4e4ad91a39 Jean*0206             uMerPri(k) = (OBEuStevens(j,k,bi,bj)-uMerBar)
2c64c37c71 Mart*0207 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0208      &           * _maskW(i,j,k,bi,bj)
9f68a71aa4 Mart*0209            ENDDO
eb8aff1784 Jean*0210 C     vertical average of input field
e37ac28e89 Mart*0211            drFbar  = 0. _d 0
                0212            uMerBar = 0. _d 0
4e4ad91a39 Jean*0213            DO k=1,Nr
                0214             uMerBar = uMerBar + OBEu(j,k,bi,bj)
                0215      &           *drF(k)* _hFacW(i,j,k,bi,bj)
                0216             drFBar = drFBar + drF(k)* _hFacW(i,j,k,bi,bj)
9f68a71aa4 Mart*0217            ENDDO
e37ac28e89 Mart*0218            IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
b2cb1ccb9a Mart*0219 C     Now the absolute velocity normal to the boundary is
e37ac28e89 Mart*0220 C     uMerBar + uMerPri(K).
4e4ad91a39 Jean*0221            DO k=1,Nr
                0222             OBEu(j,k,bi,bj) = (uMerBar + uMerPri(k))
                0223      &           * _maskW(i,j,k,bi,bj)
9f68a71aa4 Mart*0224 CML            OBEv(J,K,bi,bj) = 0. _d 0
b2cb1ccb9a Mart*0225 #ifdef ALLOW_NONHYDROSTATIC
4e4ad91a39 Jean*0226             OBEw(j,k,bi,bj)=0.
b2cb1ccb9a Mart*0227 #endif
9f68a71aa4 Mart*0228            ENDDO
                0229           ENDIF
b2cb1ccb9a Mart*0230          ENDDO
2c64c37c71 Mart*0231 #ifdef NONLIN_FRSURF
                0232 C     this is a bit of a hack
                0233          IF ( nonlinFreeSurf.GT.0 ) THEN
4e4ad91a39 Jean*0234           DO j=1-OLy,sNy+OLy
                0235            i = OB_Ie(j,bi,bj)
                0236            IF ( i.NE.OB_indexNone ) THEN
                0237             OBEeta(j,bi,bj) = etaN(i-1,j,bi,bj)
2c64c37c71 Mart*0238            ENDIF
                0239           ENDDO
                0240          ENDIF
                0241 #endif /* NONLIN_FRSURF */
b2cb1ccb9a Mart*0242 C     Next, we compute the phase speed correction, which depends on the
2c64c37c71 Mart*0243 C     tracer!
4e4ad91a39 Jean*0244          DO k=1,Nr
                0245           DO j=1-OLy,sNy+OLy
                0246            i = OB_Ie(j,bi,bj)
                0247            IF ( i.NE.OB_indexNone ) THEN
                0248             cflMer(j,k) = 0.5 _d 0 * _dxC(i-1,j,bi,bj)/dTtracerLev(k)
2c64c37c71 Mart*0249 CML         gFacM(J,K)  = 0. _d 0
                0250 CML         IF ( uVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacM(J,K) = 1. _d 0
4e4ad91a39 Jean*0251             gFacM(j,k)  = ABS(MIN(SIGN(1.D0,uVel(i,j,k,bi,bj)),0.D0))
2c64c37c71 Mart*0252            ELSE
4e4ad91a39 Jean*0253             cflMer(j,k) = 0. _d 0
                0254             gFacM (j,k) = 0. _d 0
2c64c37c71 Mart*0255            ENDIF
                0256           ENDDO
                0257          ENDDO
4e4ad91a39 Jean*0258 # ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0259 CADJ STORE cflMer, gFacM = comlev1_bibj, key=tkey, kind=isbyte
4e4ad91a39 Jean*0260 # endif
0daa63ec35 Mart*0261 C     theta
                0262          CALL OBCS_STEVENS_CALC_TRACER_EAST(
4e4ad91a39 Jean*0263      U        OBEt,
                0264      I        OBEtStevens, theta, gammat,
0daa63ec35 Mart*0265      I        uVel, cflMer, gFacM, pFac, aFac,
                0266      I        OB_Ie, OB_indexNone, bi, bj,
                0267      I        futureTime, futureIter,
                0268      I        myThid )
                0269 C     salinity
                0270          CALL OBCS_STEVENS_CALC_TRACER_EAST(
4e4ad91a39 Jean*0271      U        OBEs,
0daa63ec35 Mart*0272      I        OBEsStevens, salt, gammas,
                0273      I        uVel, cflMer, gFacM, pFac, aFac,
                0274      I        OB_Ie, OB_indexNone, bi, bj,
                0275      I        futureTime, futureIter,
                0276      I        myThid )
                0277 C     Template for passive tracers, requires work
                0278 CML#ifdef ALLOW_PTRACERS
                0279 CMLC     ptracers
                0280 CML         IF ( usePtracers ) THEN
                0281 CML          DO itracer = 1, PTRACERnumInUse
                0282 CML           CALL OBCS_STEVENS_CALC_TRACER_EAST(
4e4ad91a39 Jean*0283 CML     O          OBEptr       (1-OLy,1,1,1,iTracer),
                0284 CML     I          OBEpStevens  (1-OLy,1,1,1,iTracer),
0daa63ec35 Mart*0285 CML     I          pTracer(1-OLx,1-OLy,1,1,1,iTracer), gammas,
                0286 CML     I          uVel, cflMer, gFacM, pFac, aFac,
                0287 CML     I          OB_Ie, OB_indexNone, bi, bj,
                0288 CML     I          futureTime, futureIter,
                0289 CML     I          myThid )
                0290 CML          ENDDO
                0291 CML         ENDIF
                0292 CML#endif /* ALLOW_PTRACERS */
b2cb1ccb9a Mart*0293 C     IF ( useStevensEast ) THEN
eb8aff1784 Jean*0294         ENDIF
b2cb1ccb9a Mart*0295 #endif /* ALLOW_OBCS_EAST */
eb8aff1784 Jean*0296
b2cb1ccb9a Mart*0297 C ------------------------------------------------------------------------------
                0298
                0299 #ifdef ALLOW_OBCS_WEST
eb8aff1784 Jean*0300         IF ( useStevensWest ) THEN
b2cb1ccb9a Mart*0301 C     Western OB
                0302 #ifdef ALLOW_DEBUG
eb8aff1784 Jean*0303          IF (debugMode)
                0304      &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: West',myThid)
b2cb1ccb9a Mart*0305 #endif
eb8aff1784 Jean*0306 C     compute vertical average and deviation from vertical
9f68a71aa4 Mart*0307 C     average for I_obw+1
4e4ad91a39 Jean*0308          DO j=1-OLy,sNy+OLy
                0309           i = OB_Iw(j,bi,bj)
                0310           IF ( i.NE.OB_indexNone ) THEN
e37ac28e89 Mart*0311 C     first initialize some fields
                0312            drFBar  = 0. _d 0
                0313            uMerBar = 0. _d 0
4e4ad91a39 Jean*0314            DO k=1,Nr
                0315             uMerPri(k) = 0. _d 0
e37ac28e89 Mart*0316            ENDDO
4e4ad91a39 Jean*0317            DO k=1,Nr
2c64c37c71 Mart*0318 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0319             uMerBar = uMerBar + uVel(i+2,j,k,bi,bj)
2c64c37c71 Mart*0320 #else
4e4ad91a39 Jean*0321             uMerBar = uMerBar + OBWuStevens(j,k,bi,bj)
2c64c37c71 Mart*0322 #endif /* OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0323      &           *drF(k)* _hFacW(i+1,j,k,bi,bj)
                0324             drFBar = drFBar + drF(k)* _hFacW(i+1,j,k,bi,bj)
9f68a71aa4 Mart*0325            ENDDO
e37ac28e89 Mart*0326            IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
4e4ad91a39 Jean*0327            DO k=1,Nr
2c64c37c71 Mart*0328 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0329             uMerPri(k) = (uVel(i+2,j,k,bi,bj)-uMerBar)
2c64c37c71 Mart*0330 #else
4e4ad91a39 Jean*0331             uMerPri(k) = (OBWuStevens(j,k,bi,bj)-uMerBar)
2c64c37c71 Mart*0332 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0333      &           * _maskW(i+1,j,k,bi,bj)
9f68a71aa4 Mart*0334            ENDDO
eb8aff1784 Jean*0335 C     vertical average of input field
9f68a71aa4 Mart*0336            drFBar  = 0. _d 0
e37ac28e89 Mart*0337            uMerBar = 0. _d 0
4e4ad91a39 Jean*0338            DO k=1,Nr
                0339             uMerBar = uMerBar + OBWu(j,k,bi,bj)
                0340      &           *drF(k)* _hFacW(i+1,j,k,bi,bj)
                0341             drFBar = drFBar + drF(k)* _hFacW(i+1,j,k,bi,bj)
9f68a71aa4 Mart*0342            ENDDO
e37ac28e89 Mart*0343            IF ( drFBar .GT. 0. _d 0 ) uMerBar = uMerBar/drFBar
b2cb1ccb9a Mart*0344 C     Now the absolute velocity normal to the boundary is
e37ac28e89 Mart*0345 C     uMerBar + uMerPri(K).
4e4ad91a39 Jean*0346            DO k=1,Nr
                0347             OBWu(j,k,bi,bj) = (uMerBar + uMerPri(k))
                0348      &          * _maskW(i+1,j,k,bi,bj)
9f68a71aa4 Mart*0349 CML            OBWv(J,K,bi,bj) = 0. _d 0
b2cb1ccb9a Mart*0350 #ifdef ALLOW_NONHYDROSTATIC
4e4ad91a39 Jean*0351             OBWw(j,k,bi,bj)=0.
b2cb1ccb9a Mart*0352 #endif
9f68a71aa4 Mart*0353            ENDDO
                0354           ENDIF
eb8aff1784 Jean*0355          ENDDO
2c64c37c71 Mart*0356 #ifdef NONLIN_FRSURF
                0357 C     this is a bit of a hack
                0358          IF ( nonlinFreeSurf.GT.0 ) THEN
4e4ad91a39 Jean*0359           DO j=1-OLy,sNy+OLy
                0360            i = OB_Iw(j,bi,bj)
                0361            IF ( i.NE.OB_indexNone ) THEN
                0362             OBWeta(j,bi,bj) = etaN(i+1,j,bi,bj)
2c64c37c71 Mart*0363            ENDIF
                0364           ENDDO
                0365          ENDIF
                0366 #endif /* NONLIN_FRSURF */
b2cb1ccb9a Mart*0367 C     Next, we compute the phase speed correction, which depends on the
74019f026d Jean*0368 C     tracer!
4e4ad91a39 Jean*0369          DO k=1,Nr
                0370           DO j=1-OLy,sNy+OLy
                0371            i = OB_Iw(j,bi,bj)
                0372            IF ( i.NE.OB_indexNone ) THEN
                0373             cflMer(j,k) = 0.5 _d 0 * _dxC(i+2,j,bi,bj)/dTtracerLev(k)
2c64c37c71 Mart*0374 CML         gFacM = 0. _d 0
                0375 CML         IF ( uVel(I+1,J,K,bi,bj) .GT. 0. _d 0 ) gFacM = 1. _d 0
4e4ad91a39 Jean*0376             gFacM(j,k)  = ABS(MAX(SIGN(1.D0,uVel(i+1,j,k,bi,bj)),0.D0))
2c64c37c71 Mart*0377            ELSE
4e4ad91a39 Jean*0378             cflMer(j,k) = 0. _d 0
                0379             gFacM (j,k) = 0. _d 0
2c64c37c71 Mart*0380            ENDIF
                0381           ENDDO
                0382          ENDDO
4e4ad91a39 Jean*0383 # ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0384 CADJ STORE cflMer, gFacM = comlev1_bibj, key=tkey, kind=isbyte
4e4ad91a39 Jean*0385 # endif
0daa63ec35 Mart*0386 C     theta
                0387          CALL OBCS_STEVENS_CALC_TRACER_WEST(
4e4ad91a39 Jean*0388      U        OBWt,
                0389      I        OBWtStevens, theta, gammat,
0daa63ec35 Mart*0390      I        uVel, cflMer, gFacM, pFac, aFac,
                0391      I        OB_Iw, OB_indexNone, bi, bj,
                0392      I        futureTime, futureIter,
                0393      I        myThid )
                0394 C     salinity
                0395          CALL OBCS_STEVENS_CALC_TRACER_WEST(
4e4ad91a39 Jean*0396      U        OBWs,
0daa63ec35 Mart*0397      I        OBWsStevens, salt, gammas,
                0398      I        uVel, cflMer, gFacM, pFac, aFac,
                0399      I        OB_Iw, OB_indexNone, bi, bj,
                0400      I        futureTime, futureIter,
                0401      I        myThid )
                0402 C     ptracers
                0403 C     IF ( useStevensWest ) THEN
eb8aff1784 Jean*0404         ENDIF
b2cb1ccb9a Mart*0405 #endif /* ALLOW_OBCS_WEST */
                0406
                0407 C ------------------------------------------------------------------------------
                0408
                0409 #ifdef ALLOW_OBCS_NORTH
eb8aff1784 Jean*0410         IF ( useStevensNorth ) THEN
b2cb1ccb9a Mart*0411 C         Northern OB
                0412 #ifdef ALLOW_DEBUG
eb8aff1784 Jean*0413          IF (debugMode)
                0414      &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: North',myThid)
b2cb1ccb9a Mart*0415 #endif
9f68a71aa4 Mart*0416 C     compute vertical average and deviation from vertical
                0417 C     average for J_obn
4e4ad91a39 Jean*0418          DO i=1-OLx,sNx+OLx
                0419           j = OB_Jn(i,bi,bj)
                0420           IF ( j.NE.OB_indexNone ) THEN
e37ac28e89 Mart*0421 C     first initialize some fields
9f68a71aa4 Mart*0422            drFBar  = 0. _d 0
e37ac28e89 Mart*0423            vZonBar = 0. _d 0
4e4ad91a39 Jean*0424            DO k=1,Nr
                0425             vZonPri(k) = 0. _d 0
e37ac28e89 Mart*0426            ENDDO
4e4ad91a39 Jean*0427            DO k=1,Nr
9f68a71aa4 Mart*0428 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0429             vZonBar = vZonBar + vVel(i,j-1,k,bi,bj)
9f68a71aa4 Mart*0430 #else
4e4ad91a39 Jean*0431             vZonBar = vZonBar + OBNvStevens(i,k,bi,bj)
9f68a71aa4 Mart*0432 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0433      &           *drF(k)* _hFacS(i,j,k,bi,bj)
                0434             drFBar = drFBar + drF(k)* _hFacS(i,j,k,bi,bj)
9f68a71aa4 Mart*0435            ENDDO
e37ac28e89 Mart*0436            IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
4e4ad91a39 Jean*0437            DO k=1,Nr
9f68a71aa4 Mart*0438 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0439             vZonPri(k) = (vVel(i,j-1,k,bi,bj)-vZonBar)
9f68a71aa4 Mart*0440 #else
4e4ad91a39 Jean*0441             vZonPri(k) = (OBNvStevens(i,k,bi,bj)-vZonBar)
9f68a71aa4 Mart*0442 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0443      &           * _maskS(i,j,k,bi,bj)
9f68a71aa4 Mart*0444            ENDDO
                0445 C     vertical average of input field
                0446            drFBar  = 0. _d 0
e37ac28e89 Mart*0447            vZonBar = 0. _d 0
4e4ad91a39 Jean*0448            DO k=1,Nr
                0449             vZonBar = vZonBar + OBNv(i,k,bi,bj)
                0450      &           *drF(k)* _hFacS(i,j,k,bi,bj)
                0451             drFBar = drFBar + drF(k)* _hFacS(i,j,k,bi,bj)
9f68a71aa4 Mart*0452            ENDDO
e37ac28e89 Mart*0453            IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
9f68a71aa4 Mart*0454 C     Now the absolute velocity normal to the boundary is
e37ac28e89 Mart*0455 C     vZonBar + vZonPri(K).
4e4ad91a39 Jean*0456            DO k=1,Nr
                0457             OBNv(i,k,bi,bj) = (vZonBar + vZonPri(k))
                0458      &           * _maskS(i,j,k,bi,bj)
9f68a71aa4 Mart*0459 CML            OBNu(I,K,bi,bj) = 0. _d 0
b2cb1ccb9a Mart*0460 #ifdef ALLOW_NONHYDROSTATIC
4e4ad91a39 Jean*0461             OBNw(i,k,bi,bj)=0.
b2cb1ccb9a Mart*0462 #endif
9f68a71aa4 Mart*0463            ENDDO
                0464           ENDIF
                0465          ENDDO
b2cb1ccb9a Mart*0466 #ifdef NONLIN_FRSURF
9f68a71aa4 Mart*0467 C     this is a bit of a hack
                0468          IF ( nonlinFreeSurf.GT.0 ) THEN
4e4ad91a39 Jean*0469           DO i=1-OLx,sNx+OLx
                0470            j = OB_Jn(i,bi,bj)
                0471            IF ( j.NE.OB_indexNone ) THEN
                0472             OBNeta(i,bi,bj) = etaN(i,j-1,bi,bj)
9f68a71aa4 Mart*0473            ENDIF
                0474           ENDDO
                0475          ENDIF
                0476 #endif /* NONLIN_FRSURF */
                0477 C     Next, we compute the phase speed correction, which depends on the
74019f026d Jean*0478 C     tracer!
4e4ad91a39 Jean*0479          DO k=1,Nr
                0480           DO i=1-OLx,sNx+OLx
                0481            j = OB_Jn(i,bi,bj)
                0482            IF ( j.NE.OB_indexNone ) THEN
                0483             cflZon(i,k) = 0.5 _d 0 * _dyC(i,j-1,bi,bj)/dTtracerLev(k)
9f68a71aa4 Mart*0484 CML         gFacZ(I,K) = 0. _d 0
                0485 CML         IF ( vVel(I,J,K,bi,bj) .LT. 0. _d 0 ) gFacZ(I,K) = 1. _d 0
4e4ad91a39 Jean*0486             gFacZ(i,k)  = ABS(MIN(SIGN(1.D0,vVel(i,j,k,bi,bj)),0.D0))
9f68a71aa4 Mart*0487            ELSE
4e4ad91a39 Jean*0488             cflZon(i,k) = 0. _d 0
                0489             gFacZ (i,k) = 0. _d 0
9f68a71aa4 Mart*0490            ENDIF
                0491           ENDDO
                0492          ENDDO
4e4ad91a39 Jean*0493 # ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0494 CADJ STORE cflZon, gFacZ = comlev1_bibj, key=tkey, kind=isbyte
4e4ad91a39 Jean*0495 # endif
0daa63ec35 Mart*0496 C     theta
                0497          CALL OBCS_STEVENS_CALC_TRACER_NORTH(
4e4ad91a39 Jean*0498      U        OBNt,
0daa63ec35 Mart*0499      I        OBNtStevens, theta, gammat,
4e4ad91a39 Jean*0500      I        vVel, cflZon, gFacZ, pFac, aFac,
0daa63ec35 Mart*0501      I        OB_Jn, OB_indexNone, bi, bj,
                0502      I        futureTime, futureIter,
                0503      I        myThid )
                0504 C     salinity
                0505          CALL OBCS_STEVENS_CALC_TRACER_NORTH(
4e4ad91a39 Jean*0506      U        OBNs,
0daa63ec35 Mart*0507      I        OBNsStevens, salt, gammas,
4e4ad91a39 Jean*0508      I        vVel, cflZon, gFacZ, pFac, aFac,
0daa63ec35 Mart*0509      I        OB_Jn, OB_indexNone, bi, bj,
                0510      I        futureTime, futureIter,
                0511      I        myThid )
                0512 C     ptracers
                0513 C     IF ( useStevensNorth ) THEN
eb8aff1784 Jean*0514         ENDIF
b2cb1ccb9a Mart*0515 #endif /* ALLOW_OBCS_NORTH */
                0516
                0517 C ------------------------------------------------------------------------------
                0518
                0519 #ifdef ALLOW_OBCS_SOUTH
eb8aff1784 Jean*0520         IF ( useStevensSouth ) THEN
b2cb1ccb9a Mart*0521 C         Southern OB
                0522 #ifdef ALLOW_DEBUG
eb8aff1784 Jean*0523          IF (debugMode)
                0524      &        CALL DEBUG_MSG('OBCS_CALC_STEVENS: South',myThid)
b2cb1ccb9a Mart*0525 #endif
9f68a71aa4 Mart*0526 C     compute vertical average and deviation from vertical
                0527 C     average for J_obs+1
4e4ad91a39 Jean*0528          DO i=1-OLx,sNx+OLx
                0529           j = OB_Js(i,bi,bj)
                0530           IF ( j.NE.OB_indexNone ) THEN
e37ac28e89 Mart*0531 C     first initialize some fields
9f68a71aa4 Mart*0532            drFBar  = 0. _d 0
e37ac28e89 Mart*0533            vZonBar = 0. _d 0
4e4ad91a39 Jean*0534            DO k=1,Nr
                0535             vZonPri(k) = 0. _d 0
e37ac28e89 Mart*0536            ENDDO
4e4ad91a39 Jean*0537            DO k=1,Nr
9f68a71aa4 Mart*0538 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0539             vZonBar = vZonBar + vVel(i,j+2,k,bi,bj)
9f68a71aa4 Mart*0540 #else
4e4ad91a39 Jean*0541             vZonBar = vZonBar + OBSvStevens(i,k,bi,bj)
9f68a71aa4 Mart*0542 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0543      &           *drF(k)* _hFacS(i,j+1,k,bi,bj)
                0544             drFBar = drFBar + drF(k)* _hFacS(i,j+1,k,bi,bj)
9f68a71aa4 Mart*0545            ENDDO
e37ac28e89 Mart*0546            IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
4e4ad91a39 Jean*0547            DO k=1,Nr
9f68a71aa4 Mart*0548 #ifdef OBCS_STEVENS_USE_INTERIOR_VELOCITY
4e4ad91a39 Jean*0549             vZonPri(k) = (vVel(i,j+2,k,bi,bj)-vZonBar)
9f68a71aa4 Mart*0550 #else
4e4ad91a39 Jean*0551             vZonPri(k) = (OBSvStevens(i,k,bi,bj)-vZonBar)
9f68a71aa4 Mart*0552 #endif /*  OBCS_STEVENS_USE_INTERIOR_VELOCITY */
4e4ad91a39 Jean*0553      &           * _maskS(i,j+1,k,bi,bj)
9f68a71aa4 Mart*0554            ENDDO
                0555 C     vertical average of input field
                0556            drFBar  = 0. _d 0
e37ac28e89 Mart*0557            vZonBar = 0. _d 0
4e4ad91a39 Jean*0558            DO k=1,Nr
                0559             vZonBar = vZonBar + OBSv(i,k,bi,bj)
                0560      &           *drF(k)* _hFacS(i,j+1,k,bi,bj)
                0561             drFBar = drFBar + drF(k)* _hFacS(i,j+1,k,bi,bj)
9f68a71aa4 Mart*0562            ENDDO
e37ac28e89 Mart*0563            IF ( drFBar .GT. 0. _d 0 ) vZonBar = vZonBar/drFBar
9f68a71aa4 Mart*0564 C     Now the absolute velocity normal to the boundary is
e37ac28e89 Mart*0565 C     vZonBar + vZonPri(K).
4e4ad91a39 Jean*0566            DO k=1,Nr
                0567             OBSv(i,k,bi,bj) = (vZonBar + vZonPri(k))
                0568      &          * _maskS(i,j+1,k,bi,bj)
9f68a71aa4 Mart*0569 CML            OBSu(I,K,bi,bj) = 0. _d 0
b2cb1ccb9a Mart*0570 #ifdef ALLOW_NONHYDROSTATIC
4e4ad91a39 Jean*0571             OBSw(i,k,bi,bj)=0.
b2cb1ccb9a Mart*0572 #endif
9f68a71aa4 Mart*0573            ENDDO
                0574           ENDIF
                0575          ENDDO
b2cb1ccb9a Mart*0576 #ifdef NONLIN_FRSURF
9f68a71aa4 Mart*0577 C     this is a bit of a hack
                0578          IF ( nonlinFreeSurf.GT.0 ) THEN
4e4ad91a39 Jean*0579           DO i=1-OLx,sNx+OLx
                0580            j = OB_Js(i,bi,bj)
                0581            IF ( j.NE.OB_indexNone ) THEN
                0582             OBSeta(i,bi,bj) = etaN(i,j+1,bi,bj)
9f68a71aa4 Mart*0583            ENDIF
                0584           ENDDO
                0585          ENDIF
                0586 #endif /* NONLIN_FRSURF */
                0587 C     Next, we compute the phase speed correction, which depends on the
74019f026d Jean*0588 C     tracer!
4e4ad91a39 Jean*0589          DO k=1,Nr
                0590           DO i=1-OLx,sNx+OLx
                0591            j = OB_Js(i,bi,bj)
                0592            IF ( j.NE.OB_indexNone ) THEN
                0593             cflZon(i,k) = 0.5 _d 0 * _dyC(i,j+2,bi,bj)/dTtracerLev(k)
9f68a71aa4 Mart*0594 CML         gFacZ = 0. _d 0
                0595 CML         IF ( vVel(I,J+1,K,bi,bj) .GT. 0. _d 0 ) gFacZ = 1. _d 0
4e4ad91a39 Jean*0596             gFacZ(i,k)  = ABS(MAX(SIGN(1.D0,vVel(i,j+1,k,bi,bj)),0.D0))
9f68a71aa4 Mart*0597            ELSE
4e4ad91a39 Jean*0598             cflZon(i,k) = 0. _d 0
                0599             gFacZ (i,k) = 0. _d 0
9f68a71aa4 Mart*0600            ENDIF
                0601           ENDDO
                0602          ENDDO
4e4ad91a39 Jean*0603 # ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0604 CADJ STORE cflZon, gFacZ = comlev1_bibj, key=tkey, kind=isbyte
4e4ad91a39 Jean*0605 # endif
0daa63ec35 Mart*0606 C     theta
                0607          CALL OBCS_STEVENS_CALC_TRACER_SOUTH(
4e4ad91a39 Jean*0608      U        OBSt,
0daa63ec35 Mart*0609      I        OBStStevens, theta, gammat,
4e4ad91a39 Jean*0610      I        vVel, cflZon, gFacZ, pFac, aFac,
0daa63ec35 Mart*0611      I        OB_Js, OB_indexNone, bi, bj,
                0612      I        futureTime, futureIter,
                0613      I        myThid )
                0614 C     salinity
                0615          CALL OBCS_STEVENS_CALC_TRACER_SOUTH(
4e4ad91a39 Jean*0616      U        OBSs,
0daa63ec35 Mart*0617      I        OBSsStevens, salt, gammas,
4e4ad91a39 Jean*0618      I        vVel, cflZon, gFacZ, pFac, aFac,
0daa63ec35 Mart*0619      I        OB_Js, OB_indexNone, bi, bj,
                0620      I        futureTime, futureIter,
                0621      I        myThid )
                0622 C     ptracers
9f68a71aa4 Mart*0623 C     IF ( useStevensSouth ) THEN
eb8aff1784 Jean*0624         ENDIF
b2cb1ccb9a Mart*0625 #endif /* ALLOW_OBCS_SOUTH */
                0626
                0627 C     end bi/bj-loops
                0628        ENDDO
                0629       ENDDO
                0630
b1353a35da Mart*0631 C     save the tracer fields of the previous time step for the next time step
                0632       CALL OBCS_STEVENS_SAVE_TRACERS(
                0633      I     futureTime, futureIter,
                0634      I     myThid )
b2cb1ccb9a Mart*0635 C ------------------------------------------------------------------------------
                0636
                0637 #ifdef CHECK_BALANCE
                0638       DO bj=myByLo(myThid),myByHi(myThid)
                0639        DO bi=myBxLo(myThid),myBxHi(myThid)
                0640         uPhase=0.
                0641         vPhase=0.
eb8aff1784 Jean*0642         uVelLoc = 0.
4e4ad91a39 Jean*0643         DO j=1-OLy,sNy+OLy
e37ac28e89 Mart*0644          uMerBar=0. _d 0
4e4ad91a39 Jean*0645          DO k=1,Nr
                0646           i = OB_Ie(j,bi,bj)
                0647           IF ( i.EQ.OB_indexNone ) i = 1
                0648           uPhase = uPhase + OBEu(j,k,bi,bj)
                0649      &         *drF(k)* _hFacW(i,j,k,bi,bj)*dyG(i,j,bi,bj)
                0650           i = OB_Iw(j,bi,bj)
                0651           IF ( i.EQ.OB_indexNone ) i = 1
                0652           vPhase = vPhase + OBWu(j,k,bi,bj)
                0653      &         *drF(k)* _hFacW(i+1,j,k,bi,bj)*dyG(i+1,j,bi,bj)
e37ac28e89 Mart*0654 CML          uVelLoc = uVelLoc + uMerPri(J,K)
                0655 CML     &         *drF(k)* _hFacW(I+1,J,K,bi,bj)*dyG(I+1,J,bi,bj)
                0656 CML          uMerBar(J)=uMerBar(J) + uMerPri(J,K)
                0657 CML     &         *drF(k)* _hFacW(I+1,J,K,bi,bj)
b2cb1ccb9a Mart*0658          ENDDO
e37ac28e89 Mart*0659 CML         print *, 'ml-obcs: uBar = ', j,uMerBar(J)
eb8aff1784 Jean*0660         ENDDO
b2cb1ccb9a Mart*0661 C     end bi/bj-loops
                0662        ENDDO
                0663       ENDDO
eb8aff1784 Jean*0664       _GLOBAL_SUM_RL( uPhase, myThid )
                0665       _GLOBAL_SUM_RL( vPhase, myThid )
e37ac28e89 Mart*0666 CML      _GLOBAL_SUM_RL( uVelLoc, myThid )
eb8aff1784 Jean*0667       print *, 'ml-obcs: OBE  = ',  uPhase*1 _d -6, ' Sv'
                0668       print *, 'ml-obcs: OBW  = ',  vPhase*1 _d -6, ' Sv'
e37ac28e89 Mart*0669 CML      print *, 'ml-obcs: OBWp = ', uVelLoc*1 _d -6, ' Sv'
eb8aff1784 Jean*0670 #endif /* CHECK_BALANCE */
b2cb1ccb9a Mart*0671
                0672 #ifdef ALLOW_DEBUG
eb8aff1784 Jean*0673          IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC_STEVENS',myThid)
b2cb1ccb9a Mart*0674 #endif
eb8aff1784 Jean*0675
b2cb1ccb9a Mart*0676 #endif /* ALLOW_OBCS_STEVENS */
                0677       RETURN
                0678       END
b1353a35da Mart*0679
                0680 CBOP
0daa63ec35 Mart*0681 C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_EAST
                0682 C     !INTERFACE:
                0683       SUBROUTINE OBCS_STEVENS_CALC_TRACER_EAST(
4e4ad91a39 Jean*0684      U     OBEf,
0daa63ec35 Mart*0685      I     OBE_Stevens, tracer, gammaf,
4e4ad91a39 Jean*0686      I     uVel, cflMer, gFacM, pFac, aFac,
0daa63ec35 Mart*0687      I     OB_I, OB_indexNone, bi, bj,
                0688      I     futureTime, futureIter,
                0689      I     myThid )
                0690
                0691 C     !DESCRIPTION: \bv
                0692 C     *==========================================================*
                0693 C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_EAST
                0694 C     | Calculate tracer value at the eastern OB location
                0695 C     *==========================================================*
                0696 C     \ev
                0697 C     !USES:
                0698       IMPLICIT NONE
                0699 C     == Global variables ==
                0700 #include "SIZE.h"
                0701 #include "GRID.h"
                0702
                0703 C     !INPUT/OUTPUT PARAMETERS:
                0704 C     == Routine arguments ==
                0705 C    myThid    :: my Thread Id number
                0706 C    bi, bj    :: indices of current tile
                0707       _RL futureTime
                0708       INTEGER futureIter
                0709       INTEGER myThid
                0710       INTEGER bi, bj
                0711       INTEGER OB_indexNone
                0712       INTEGER OB_I             (1-OLy:sNy+OLy,nSx,nSy)
                0713       _RL cflMer               (1-OLy:sNy+OLy,Nr)
                0714       _RL gFacM                (1-OLy:sNy+OLy,Nr)
                0715       _RL gammaf, pFac, aFac
4e4ad91a39 Jean*0716       _RL OBEf                 (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0717       _RL OBE_Stevens          (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0718       _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0719       _RL uVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0daa63ec35 Mart*0720
                0721 #ifdef ALLOW_OBCS_STEVENS
                0722 C     !LOCAL VARIABLES:
                0723 C     i,j,k        :: loop indices
                0724 C     uPhase       :: estimate of phase velocity for radiation condition
                0725 C     dtracSpace   :: horizontal difference of tracer
                0726 C     dtracTime    :: temporal difference of tracer
                0727       INTEGER i,j,k
                0728       _RL uPhase
                0729       _RL dtracSpace
                0730       _RL dTracTime
                0731 CEOP
4e4ad91a39 Jean*0732       DO k=1,Nr
                0733        DO j=1-OLy,sNy+OLy
                0734         i = OB_I(j,bi,bj)
                0735         IF ( i.NE.OB_indexNone ) THEN
                0736          dTracSpace = (tracer(i-1,j,k,bi,bj)-tracer(i-2,j,k,bi,bj))
                0737      &        * _maskW(i-1,j,k,bi,bj)
                0738          dTracTime  = (tracer(i-1,j,k,bi,bj)-OBE_Stevens(j,k,bi,bj))
                0739          uPhase = cflMer(j,k) * pFac
0daa63ec35 Mart*0740          IF ( dTracSpace .NE. 0. _d 0 ) THEN
4e4ad91a39 Jean*0741           uPhase = MIN( cflMer(j,k),
                0742      &         MAX( 0.D0, -cflMer(j,k)*dTracTime/dTracSpace )
0daa63ec35 Mart*0743      &         ) * pFac
                0744          ENDIF
                0745 C     Compute the tracer tendency here, the tracer will be updated
                0746 C     with a simple Euler forward step in S/R obcs_apply_ts
4e4ad91a39 Jean*0747          OBEf(j,k,bi,bj) = _maskW(i,j,k,bi,bj) * (
                0748      &        - ( aFac*MAX(0.D0,uVel(i,j,k,bi,bj)) + uPhase )
                0749      &        *(tracer(i,j,k,bi,bj)-tracer(i-1,j,k,bi,bj))
                0750      &        * _recip_dxC(i,j,bi,bj)
                0751      &        - gFacM(j,k) * gammaf
                0752      &        * (tracer(i,j,k,bi,bj)-OBEf(j,k,bi,bj)) )
0daa63ec35 Mart*0753         ENDIF
                0754        ENDDO
                0755       ENDDO
                0756
                0757 #endif /* ALLOW_OBCS_STEVENS */
                0758       RETURN
                0759       END
                0760
                0761 CBOP
                0762 C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_WEST
                0763 C     !INTERFACE:
                0764       SUBROUTINE OBCS_STEVENS_CALC_TRACER_WEST(
4e4ad91a39 Jean*0765      U     OBWf,
0daa63ec35 Mart*0766      I     OBW_Stevens, tracer, gammaf,
4e4ad91a39 Jean*0767      I     uVel, cflMer, gFacM, pFac, aFac,
0daa63ec35 Mart*0768      I     OB_I, OB_indexNone, bi, bj,
                0769      I     futureTime, futureIter,
                0770      I     myThid )
                0771
                0772 C     !DESCRIPTION: \bv
                0773 C     *==========================================================*
                0774 C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_WEST
                0775 C     | Calculate tracer value at the western OB location
                0776 C     *==========================================================*
                0777 C     \ev
                0778 C     !USES:
                0779       IMPLICIT NONE
                0780 C     == Global variables ==
                0781 #include "SIZE.h"
                0782 #include "GRID.h"
                0783
                0784 C     !INPUT/OUTPUT PARAMETERS:
                0785 C     == Routine arguments ==
                0786 C    myThid    :: my Thread Id number
                0787 C    bi, bj    :: indices of current tile
                0788       _RL futureTime
                0789       INTEGER futureIter
                0790       INTEGER myThid
                0791       INTEGER bi, bj
                0792       INTEGER OB_indexNone
                0793       INTEGER OB_I             (1-OLy:sNy+OLy,nSx,nSy)
                0794       _RL cflMer               (1-OLy:sNy+OLy,Nr)
                0795       _RL gFacM                (1-OLy:sNy+OLy,Nr)
                0796       _RL gammaf, pFac, aFac
4e4ad91a39 Jean*0797       _RL OBWf                 (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0798       _RL OBW_Stevens          (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0799       _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0800       _RL uVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0daa63ec35 Mart*0801
                0802 #ifdef ALLOW_OBCS_STEVENS
                0803 C     !LOCAL VARIABLES:
                0804 C     i,j,k        :: loop indices
                0805 C     uPhase       :: estimate of phase velocity for radiation condition
                0806 C     dtracSpace   :: horizontal difference of tracer
                0807 C     dtracTime    :: temporal difference of tracer
                0808       INTEGER i,j,k
                0809       _RL uPhase
                0810       _RL dtracSpace
                0811       _RL dTracTime
                0812 CEOP
                0813
4e4ad91a39 Jean*0814       DO k=1,Nr
                0815        DO j=1-OLy,sNy+OLy
                0816         i = OB_I(j,bi,bj)
                0817         IF ( i.NE.OB_indexNone ) THEN
                0818          dTracSpace = (tracer(i+2,j,k,bi,bj)-tracer(i+1,j,k,bi,bj))
                0819      &        * _maskW(i+2,j,k,bi,bj)
                0820          dTracTime  = (tracer(i+1,j,k,bi,bj)-OBW_Stevens(j,k,bi,bj))
                0821          uPhase = -cflMer(j,k) * pFac
0daa63ec35 Mart*0822          IF ( dTracSpace .NE. 0. _d 0 ) THEN
4e4ad91a39 Jean*0823           uPhase = MAX( -cflMer(j,k),
                0824      &         MIN( 0.D0, -cflMer(j,k)*dTracTime/dTracSpace )
0daa63ec35 Mart*0825      &         ) * pFac
                0826          ENDIF
                0827 C     Compute the tracer tendency here, the tracer will be updated
                0828 C     with a simple Euler forward step in S/R obcs_apply_ts
4e4ad91a39 Jean*0829          OBWf(j,k,bi,bj) = _maskW(i+1,j,k,bi,bj) * (
                0830      &        - ( aFac*MIN(0.D0,uVel(i+1,j,k,bi,bj)) + uPhase )
                0831      &        *(tracer(i+1,j,k,bi,bj)-tracer(i,j,k,bi,bj))
                0832      &        * _recip_dxC(i+1,j,bi,bj)
                0833      &        - gFacM(j,k) * gammaf
                0834      &        * (tracer(i,j,k,bi,bj)-OBWf(j,k,bi,bj)) )
0daa63ec35 Mart*0835         ENDIF
                0836        ENDDO
                0837       ENDDO
                0838
                0839 #endif /* ALLOW_OBCS_STEVENS */
                0840       RETURN
                0841       END
                0842
                0843 CBOP
                0844 C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_NORTH
                0845 C     !INTERFACE:
                0846       SUBROUTINE OBCS_STEVENS_CALC_TRACER_NORTH(
4e4ad91a39 Jean*0847      U     OBNf,
0daa63ec35 Mart*0848      I     OBN_Stevens, tracer, gammaf,
4e4ad91a39 Jean*0849      I     vVel, cflZon, gFacZ, pFac, aFac,
0daa63ec35 Mart*0850      I     OB_J, OB_indexNone, bi, bj,
                0851      I     futureTime, futureIter,
                0852      I     myThid )
                0853
                0854 C     !DESCRIPTION: \bv
                0855 C     *==========================================================*
                0856 C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_NORTH
                0857 C     | Calculate tracer value at the northern OB location
                0858 C     *==========================================================*
                0859 C     \ev
                0860 C     !USES:
                0861       IMPLICIT NONE
                0862 C     == Global variables ==
                0863 #include "SIZE.h"
                0864 #include "GRID.h"
                0865
                0866 C     !INPUT/OUTPUT PARAMETERS:
                0867 C     == Routine arguments ==
                0868 C    myThid    :: my Thread Id number
                0869 C    bi, bj    :: indices of current tile
                0870       _RL futureTime
                0871       INTEGER futureIter
                0872       INTEGER myThid
                0873       INTEGER bi, bj
                0874       INTEGER OB_indexNone
                0875       INTEGER OB_J             (1-OLx:sNx+OLx,nSx,nSy)
                0876       _RL cflZon               (1-OLx:sNx+OLx,Nr)
                0877       _RL gFacZ                (1-OLx:sNx+OLx,Nr)
                0878       _RL gammaf, pFac, aFac
4e4ad91a39 Jean*0879       _RL OBNf                 (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0880       _RL OBN_Stevens          (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0881       _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0882       _RL vVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0daa63ec35 Mart*0883
                0884 #ifdef ALLOW_OBCS_STEVENS
                0885 C     !LOCAL VARIABLES:
                0886 C     i,j,k        :: loop indices
                0887 C     vPhase       :: estimate of phase velocity for radiation condition
                0888 C     dtracSpace   :: horizontal difference of tracer
                0889 C     dtracTime    :: temporal difference of tracer
                0890       INTEGER i,j,k
                0891       _RL vPhase
                0892       _RL dtracSpace
                0893       _RL dTracTime
                0894 CEOP
4e4ad91a39 Jean*0895       DO k=1,Nr
                0896        DO i=1-OLx,sNx+OLx
                0897         j = OB_J(i,bi,bj)
                0898         IF ( j.NE.OB_indexNone ) THEN
0daa63ec35 Mart*0899 C     Theta first:
4e4ad91a39 Jean*0900          dTracSpace = (tracer(i,j-1,k,bi,bj)-tracer(i,j-2,k,bi,bj))
                0901      &        * _maskS(i,j-1,k,bi,bj)
                0902          dTracTime  = (tracer(i,j-1,k,bi,bj)-OBN_Stevens(i,k,bi,bj))
                0903          vPhase = cflZon(i,k) * pFac
0daa63ec35 Mart*0904          IF ( dTracSpace .NE. 0. _d 0 ) THEN
4e4ad91a39 Jean*0905           vPhase = MIN( cflZon(i,k),
                0906      &         MAX( 0.D0, -cflZon(i,k)*dTracTime/dTracSpace )
0daa63ec35 Mart*0907      &         ) * pFac
                0908          ENDIF
                0909 C     Compute the tracer tendency here, the tracer will be updated
                0910 C     with a simple Euler forward step in S/R obcs_apply_ts
4e4ad91a39 Jean*0911          OBNf(i,k,bi,bj) = _maskS(i,j,k,bi,bj) * (
                0912      &        - ( aFac*MAX(0.D0,vVel(i,j,k,bi,bj)) + vPhase )
                0913      &        *(tracer(i,j,k,bi,bj)-tracer(i,j-1,k,bi,bj))
                0914      &        * _recip_dyC(i,j,bi,bj)
                0915      &        - gFacZ(i,k) * gammaf
                0916      &        * (tracer(i,j,k,bi,bj)-OBNf(i,k,bi,bj)) )
0daa63ec35 Mart*0917         ENDIF
                0918        ENDDO
                0919       ENDDO
                0920
                0921 #endif /* ALLOW_OBCS_STEVENS */
                0922       RETURN
                0923       END
                0924
                0925 CBOP
                0926 C     !ROUTINE: OBCS_STEVENS_CALC_TRACER_SOUTH
                0927 C     !INTERFACE:
                0928       SUBROUTINE OBCS_STEVENS_CALC_TRACER_SOUTH(
4e4ad91a39 Jean*0929      U     OBSf,
0daa63ec35 Mart*0930      I     OBS_Stevens, tracer, gammaf,
4e4ad91a39 Jean*0931      I     vVel, cflZon, gFacZ, pFac, aFac,
0daa63ec35 Mart*0932      I     OB_J, OB_indexNone, bi, bj,
                0933      I     futureTime, futureIter,
                0934      I     myThid )
                0935
                0936 C     !DESCRIPTION: \bv
                0937 C     *==========================================================*
                0938 C     | SUBROUTINE OBCS_STEVENS_CALC_TRACER_SOUTH
                0939 C     | Calculate tracer value at the southern OB location
                0940 C     *==========================================================*
                0941 C     \ev
                0942 C     !USES:
                0943       IMPLICIT NONE
                0944 C     == Global variables ==
                0945 #include "SIZE.h"
                0946 #include "GRID.h"
                0947
                0948 C     !INPUT/OUTPUT PARAMETERS:
                0949 C     == Routine arguments ==
                0950 C    myThid    :: my Thread Id number
                0951 C    bi, bj    :: indices of current tile
                0952       _RL futureTime
                0953       INTEGER futureIter
                0954       INTEGER myThid
                0955       INTEGER bi, bj
                0956       INTEGER OB_indexNone
                0957       INTEGER OB_J             (1-OLx:sNx+OLx,nSx,nSy)
                0958       _RL cflZon               (1-OLx:sNx+OLx,Nr)
                0959       _RL gFacZ                (1-OLx:sNx+OLx,Nr)
                0960       _RL gammaf, pFac, aFac
4e4ad91a39 Jean*0961       _RL OBSf                 (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0962       _RL OBS_Stevens          (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0963       _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0964       _RL vVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
0daa63ec35 Mart*0965
                0966 #ifdef ALLOW_OBCS_STEVENS
                0967 C     !LOCAL VARIABLES:
                0968 C     i,j,k        :: loop indices
                0969 C     vPhase       :: estimate of phase velocity for radiation condition
                0970 C     dtracSpace   :: horizontal difference of tracer
                0971 C     dtracTime    :: temporal difference of tracer
                0972       INTEGER i,j,k
                0973       _RL vPhase
                0974       _RL dtracSpace
                0975       _RL dTracTime
                0976 CEOP
4e4ad91a39 Jean*0977       DO k=1,Nr
                0978        DO i=1-OLx,sNx+OLx
                0979         j = OB_J(i,bi,bj)
                0980         IF ( j.NE.OB_indexNone ) THEN
                0981          dTracSpace = (tracer(i,j+2,k,bi,bj)-tracer(i,j+1,k,bi,bj))
                0982      &        * _maskS(i,j+2,k,bi,bj)
                0983          dTracTime  = (tracer(i,j+1,k,bi,bj)-OBS_Stevens(i,k,bi,bj))
                0984          vPhase = -cflZon(i,k) * pFac
0daa63ec35 Mart*0985          IF ( dTracSpace .NE. 0. _d 0 ) THEN
4e4ad91a39 Jean*0986           vPhase = MAX( -cflZon(i,k),
                0987      &         MIN( 0.D0, -cflZon(i,k)*dTracTime/dTracSpace )
0daa63ec35 Mart*0988      &         ) * pFac
                0989          ENDIF
                0990 C     Compute the tracer tendency here, the tracer will be updated
                0991 C     with a simple Euler forward step in S/R obcs_apply_ts
4e4ad91a39 Jean*0992          OBSf(i,k,bi,bj) = _maskS(i,j+1,k,bi,bj) * (
                0993      &        - ( aFac*MIN(0.D0,vVel(i,j+1,k,bi,bj)) + vPhase )
                0994      &        *(tracer(i,j+1,k,bi,bj)-tracer(i,j,k,bi,bj))
                0995      &        * _recip_dyC(i,j+1,bi,bj)
                0996      &        - gFacZ(i,k) * gammaf
                0997      &        * (tracer(i,j,k,bi,bj)-OBSf(i,k,bi,bj)) )
0daa63ec35 Mart*0998         ENDIF
                0999        ENDDO
                1000       ENDDO
                1001
                1002 #endif /* ALLOW_OBCS_STEVENS */
                1003       RETURN
                1004       END
                1005
                1006 CBOP
b1353a35da Mart*1007 C     !ROUTINE: OBCS_STEVENS_SAVE_TRACERS
                1008 C     !INTERFACE:
                1009       SUBROUTINE OBCS_STEVENS_SAVE_TRACERS(
                1010      I     futureTime, futureIter,
                1011      I     myThid )
                1012 C     !DESCRIPTION: \bv
                1013 C     *==========================================================*
                1014 C     | SUBROUTINE OBCS_STEVENS_SAVE_TRACERS
74019f026d Jean*1015 C     | Save tracers (of previous time step) at the OB location
                1016 C     | to be used in the next time step for Stevens boundary
b1353a35da Mart*1017 C     | conditions
                1018 C     *==========================================================*
                1019 C     \ev
                1020 C     !USES:
                1021       IMPLICIT NONE
                1022 C     == Global variables ==
                1023 #include "SIZE.h"
                1024 #include "EEPARAMS.h"
                1025 #include "GRID.h"
                1026 #include "OBCS_PARAMS.h"
                1027 #include "OBCS_GRID.h"
                1028 #include "OBCS_FIELDS.h"
                1029 #include "DYNVARS.h"
                1030 CML#ifdef ALLOW_PTRACERS
                1031 CML#include "PTRACERS_SIZE.h"
                1032 CML#include "PTRACERS_PARAMS.h"
                1033 CML#include "PTRACERS_FIELDS.h"
                1034 CML#include "OBCS_PTRACERS.h"
                1035 CML#endif /* ALLOW_PTRACERS */
                1036
                1037 C     !INPUT/OUTPUT PARAMETERS:
                1038 C     == Routine arguments ==
                1039 C    myThid    :: my Thread Id number
                1040       _RL futureTime
                1041       INTEGER futureIter
                1042       INTEGER myThid
                1043
397c34a218 Mart*1044 #ifdef ALLOW_OBCS_STEVENS
b1353a35da Mart*1045 C     !LOCAL VARIABLES:
4e4ad91a39 Jean*1046 C     bi, bj     :: indices of current tile
                1047 C     i,j,k      :: loop indices
                1048 C     Iobc, Jobc :: position-index of open boundary
b1353a35da Mart*1049       INTEGER bi,bj,i,j,k,Iobc,Jobc
                1050 CEOP
                1051
                1052       DO bj=myByLo(myThid),myByHi(myThid)
                1053        DO bi=myBxLo(myThid),myBxHi(myThid)
                1054 #ifdef ALLOW_OBCS_NORTH
                1055         IF ( tileHasOBN(bi,bj) .AND. useStevensNorth ) THEN
                1056 C Northern boundary
74019f026d Jean*1057          DO i=1-OLx,sNx+OLx
b1353a35da Mart*1058           Jobc = OB_Jn(i,bi,bj)
74019f026d Jean*1059           IF ( Jobc.NE.OB_indexNone ) THEN
b1353a35da Mart*1060            DO k = 1,Nr
                1061             OBNtStevens(i,k,bi,bj) = theta(i,Jobc-1,k,bi,bj)
                1062      &           *maskC(i,Jobc+1,k,bi,bj)
                1063             OBNsStevens(i,k,bi,bj) = salt(i,Jobc-1,k,bi,bj)
                1064      &           *maskC(i,Jobc+1,k,bi,bj)
                1065            ENDDO
                1066           ENDIF
                1067          ENDDO
                1068         ENDIF
                1069 #endif /* ALLOW_OBCS_NORTH */
                1070 #ifdef ALLOW_OBCS_SOUTH
                1071         IF ( tileHasOBS(bi,bj) .AND. useStevensSouth ) THEN
                1072 C Southern boundary
74019f026d Jean*1073          DO i=1-OLx,sNx+OLx
b1353a35da Mart*1074           Jobc = OB_Js(i,bi,bj)
74019f026d Jean*1075           IF ( Jobc.NE.OB_indexNone ) THEN
b1353a35da Mart*1076            DO k = 1,Nr
                1077             OBStStevens(i,k,bi,bj) = theta(i,Jobc+1,k,bi,bj)
                1078      &           *maskC(i,Jobc+1,k,bi,bj)
                1079             OBSsStevens(i,k,bi,bj) = salt(i,Jobc+1,k,bi,bj)
                1080      &           *maskC(i,Jobc+1,k,bi,bj)
                1081            ENDDO
                1082           ENDIF
                1083          ENDDO
                1084         ENDIF
                1085 #endif /* ALLOW_OBCS_SOUTH */
                1086 #ifdef ALLOW_OBCS_EAST
                1087         IF ( tileHasOBE(bi,bj) .AND. useStevensEast ) THEN
                1088 C Eastern boundary
74019f026d Jean*1089          DO j=1-OLy,sNy+OLy
b1353a35da Mart*1090           Iobc = OB_Ie(j,bi,bj)
74019f026d Jean*1091           IF ( Iobc.NE.OB_indexNone ) THEN
b1353a35da Mart*1092            DO k = 1,Nr
                1093             OBEtStevens(j,k,bi,bj) = theta(Iobc-1,j,k,bi,bj)
                1094      &           *maskC(Iobc-1,j,k,bi,bj)
                1095             OBEsStevens(j,k,bi,bj) = salt(Iobc-1,j,k,bi,bj)
                1096      &           *maskC(Iobc-1,j,k,bi,bj)
                1097            ENDDO
                1098           ENDIF
                1099          ENDDO
                1100         ENDIF
                1101 #endif /* ALLOW_OBCS_EAST */
                1102 #ifdef ALLOW_OBCS_WEST
                1103         IF ( tileHasOBW(bi,bj) .AND. useStevensWest ) THEN
                1104 C Western boundary
74019f026d Jean*1105          DO j=1-OLy,sNy+OLy
b1353a35da Mart*1106           Iobc = OB_Iw(j,bi,bj)
74019f026d Jean*1107           IF ( Iobc.NE.OB_indexNone ) THEN
b1353a35da Mart*1108            DO k = 1,Nr
                1109             OBWtStevens(j,k,bi,bj) = theta(Iobc+1,j,k,bi,bj)
                1110      &           *maskC(Iobc+1,j,k,bi,bj)
                1111             OBWsStevens(j,k,bi,bj) = salt(Iobc+1,j,k,bi,bj)
                1112      &           *maskC(Iobc+1,j,k,bi,bj)
                1113            ENDDO
                1114           ENDIF
                1115          ENDDO
                1116         ENDIF
                1117 #endif /* ALLOW_OBCS_WEST */
                1118        ENDDO
                1119       ENDDO
                1120 #endif /* ALLOW_OBCS_STEVENS */
                1121       RETURN
                1122       END