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