pkg/obcs/obcs_balance_flow.F

20c0bcbffa Jean*0001 #include "OBCS_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: OBCS_BALANCE_FLOW
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE OBCS_BALANCE_FLOW( myTime, myIter, myThid )
                0009
                0010 C     !DESCRIPTION:
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE OBCS_BALANCE_FLOW
                0013 C     | o Modify OB normal flow to ensure no net inflow
                0014 C     *==========================================================*
                0015
                0016 C     !USES:
                0017       IMPLICIT NONE
                0018
                0019 C     === Global variables ===
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
9b4f2a04e2 Jean*0024 #include "OBCS_PARAMS.h"
                0025 #include "OBCS_GRID.h"
                0026 #include "OBCS_FIELDS.h"
abfe198bce Mart*0027 #include "FFIELDS.h"
                0028
20c0bcbffa Jean*0029 C     !INPUT/OUTPUT PARAMETERS:
                0030       _RL myTime
                0031       INTEGER myIter
                0032       INTEGER myThid
                0033 CEOP
                0034
af61e5eb16 Mart*0035 #if ( defined ALLOW_OBCS && defined ALLOW_OBCS_BALANCE )
20c0bcbffa Jean*0036
                0037 C     !FUNCTIONS:
                0038
                0039 C     !LOCAL VARIABLES:
                0040 C     bi, bj       :: tile indices
                0041 C     i,j,k        :: loop indices
                0042 C     iB, jB       :: local index of open boundary
                0043 C     msgBuf       :: Informational/error message buffer
                0044       INTEGER bi, bj
005af54e38 Jean*0045       INTEGER i, j, k
                0046 #if (defined ALLOW_OBCS_EAST ) || (defined ALLOW_OBCS_WEST )
                0047       INTEGER iB
                0048 #endif
                0049 #if (defined ALLOW_OBCS_NORTH) || (defined ALLOW_OBCS_SOUTH)
                0050       INTEGER jB
                0051 #endif
20c0bcbffa Jean*0052       CHARACTER*(MAX_LEN_MBUF) msgBuf
005af54e38 Jean*0053       _RL areaOB, tmpA
20c0bcbffa Jean*0054       _RL inFlow, flowE, flowW, flowN, flowS
                0055       _RL tileArea(nSx,nSy)
                0056       _RL tileFlow(nSx,nSy)
                0057       _RL tileAreaOB(nSx,nSy)
                0058       _RL tileInFlow(nSx,nSy)
                0059       LOGICAL flag
abfe198bce Mart*0060       _RL netFreshWaterFlux
                0061       _RL shelfIceNetMassFlux
005af54e38 Jean*0062 #ifdef ALLOW_OBCS_EAST
                0063       _RL areaE
                0064 #endif
                0065 #ifdef ALLOW_OBCS_WEST
                0066       _RL areaW
                0067 #endif
                0068 #ifdef ALLOW_OBCS_NORTH
                0069       _RL areaN
                0070 #endif
                0071 #ifdef ALLOW_OBCS_SOUTH
                0072       _RL areaS
                0073 #endif
20c0bcbffa Jean*0074
                0075 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0076 C-- Old method (OBCS_balanceFac < 0): balance each OB separately
                0077 C-- New method applies to all OB with BCS_balanceFac >=0 :
                0078 C      ensure that the net inflow through all OB is balanced by correcting
                0079 C      each OB normal flow with a uniform velocity, using the corresponding
                0080 C      weight factor OBCS_balanceFac.
abfe198bce Mart*0081 C  e.g., OBCS_balanceFacE,W,N,S= 1, -1, 2, 0 :
20c0bcbffa Jean*0082 C      => correct Western OBWu (by substracting a uniform velocity) to ensure
                0083 C         zero net transport through Western OB
                0084 C      => correct Eastern and Northern normal flow (twice larger Northern
                0085 C         velocity correction than Eastern correction) to ensure that
                0086 C         the total inflow through E+N+S OB is balanced
                0087
                0088 #ifdef ALLOW_DEBUG
                0089       IF (debugMode) CALL DEBUG_ENTER('OBCS_BALANCE_FLOW',myThid)
                0090 #endif
                0091
                0092 C--   Integrate the transport through each OB
                0093       flag = .FALSE.
                0094       areaOB = 0. _d 0
                0095       inFlow = 0. _d 0
                0096       DO bj=myByLo(myThid),myByHi(myThid)
                0097        DO bi=myBxLo(myThid),myBxHi(myThid)
                0098         tileAreaOB(bi,bj) = 0.
                0099         tileInFlow(bi,bj) = 0.
                0100        ENDDO
                0101       ENDDO
                0102
                0103 #ifdef ALLOW_OBCS_EAST
                0104       areaE = 0. _d 0
                0105       flowE = 0. _d 0
                0106       flag = flag .OR. ( OBCS_balanceFacE.GT.0. )
                0107       DO bj=myByLo(myThid),myByHi(myThid)
                0108        DO bi=myBxLo(myThid),myBxHi(myThid)
                0109         tileArea(bi,bj) = 0.
                0110         tileFlow(bi,bj) = 0.
                0111         IF ( tileHasOBE(bi,bj) ) THEN
                0112          DO k=1,Nr
                0113           DO j=1,sNy
                0114            iB = OB_Ie(j,bi,bj)
838c7f9401 Jean*0115 C-    If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
                0116 C     communicates with tile interior (sNx+1) rather than with halo region (i=1)
                0117            IF ( iB.NE.OB_indexNone .AND. iB.GT.1 ) THEN
f441182300 Jean*0118             tmpA = drF(k)*hFacW(iB,j,k,bi,bj)*dyG(iB,j,bi,bj)
                0119      &                   *maskInW(iB,j,bi,bj)
47f36df0c2 Jean*0120             tileArea(bi,bj) = tileArea(bi,bj) + tmpA
                0121             tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBEu(j,k,bi,bj)
20c0bcbffa Jean*0122            ENDIF
                0123           ENDDO
                0124          ENDDO
                0125          IF ( OBCS_balanceFacE.GE.0. ) THEN
                0126            tileInFlow(bi,bj) = tileInFlow(bi,bj) - tileFlow(bi,bj)
                0127            tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
                0128      &                       + tileArea(bi,bj)*OBCS_balanceFacE
                0129          ENDIF
                0130          areaE = areaE + tileArea(bi,bj)
                0131          flowE = flowE + tileFlow(bi,bj)
                0132         ENDIF
                0133        ENDDO
                0134       ENDDO
                0135 c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
                0136 c    &   'OBCS_balance it,areaE,flowE=', myIter, areaE, flowE
                0137       IF ( OBCS_balanceFacE.LT.0. ) THEN
                0138         CALL GLOBAL_SUM_TILE_RL( tileArea, areaE, myThid )
                0139         IF ( areaE.GT.0. ) THEN
                0140           CALL GLOBAL_SUM_TILE_RL( tileFlow, flowE, myThid )
8830b8f970 Jean*0141           IF ( debugLevel.GE.debLevC ) THEN
20c0bcbffa Jean*0142             WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
                0143      &            myIter, ' ) correct OBEu:', flowE, -flowE/areaE
                0144             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0145      &                          SQUEEZE_RIGHT, myThid )
                0146           ENDIF
                0147           flowE = -flowE/areaE
                0148         ENDIF
                0149       ENDIF
                0150 #endif /* ALLOW_OBCS_EAST */
                0151
                0152 #ifdef ALLOW_OBCS_WEST
                0153       areaW = 0. _d 0
                0154       flowW = 0. _d 0
                0155       flag = flag .OR. ( OBCS_balanceFacW.GT.0. )
                0156       DO bj=myByLo(myThid),myByHi(myThid)
                0157        DO bi=myBxLo(myThid),myBxHi(myThid)
                0158         tileArea(bi,bj) = 0.
                0159         tileFlow(bi,bj) = 0.
                0160         IF ( tileHasOBW(bi,bj) ) THEN
                0161          DO k=1,Nr
                0162           DO j=1,sNy
                0163            iB = OB_Iw(j,bi,bj)
838c7f9401 Jean*0164 C-    If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
                0165 C     communicates with tile interior (i=0) rather than with halo region (i=sNx)
                0166            IF ( iB.NE.OB_indexNone .AND. iB.LT.sNx ) THEN
f441182300 Jean*0167             tmpA = drF(k)*hFacW(1+iB,j,k,bi,bj)*dyG(1+iB,j,bi,bj)
                0168      &                   *maskInW(1+iB,j,bi,bj)
47f36df0c2 Jean*0169             tileArea(bi,bj) = tileArea(bi,bj) + tmpA
                0170             tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBWu(j,k,bi,bj)
20c0bcbffa Jean*0171            ENDIF
                0172           ENDDO
                0173          ENDDO
                0174          IF ( OBCS_balanceFacW.GE.0. ) THEN
                0175            tileInFlow(bi,bj) = tileInFlow(bi,bj) + tileFlow(bi,bj)
                0176            tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
                0177      &                       + tileArea(bi,bj)*OBCS_balanceFacW
                0178          ENDIF
                0179          areaW = areaW + tileArea(bi,bj)
                0180          flowW = flowW + tileFlow(bi,bj)
                0181         ENDIF
                0182        ENDDO
                0183       ENDDO
                0184 c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
                0185 c    &   'OBCS_balance it,areaW,flowW=', myIter, areaW, flowW
                0186       IF ( OBCS_balanceFacW.LT.0. ) THEN
                0187         CALL GLOBAL_SUM_TILE_RL( tileArea, areaW, myThid )
                0188         IF ( areaW.GT.0. ) THEN
                0189           CALL GLOBAL_SUM_TILE_RL( tileFlow, flowW, myThid )
8830b8f970 Jean*0190           IF ( debugLevel.GE.debLevC ) THEN
20c0bcbffa Jean*0191             WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
                0192      &            myIter, ' ) correct OBWu:', flowW, -flowW/areaW
                0193             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0194      &                          SQUEEZE_RIGHT, myThid )
                0195           ENDIF
                0196           flowW = -flowW/areaW
                0197         ENDIF
                0198       ENDIF
                0199 #endif /* ALLOW_OBCS_WEST */
                0200
                0201 #ifdef ALLOW_OBCS_NORTH
                0202       areaN = 0. _d 0
                0203       flowN = 0. _d 0
                0204       flag = flag .OR. ( OBCS_balanceFacN.GT.0. )
                0205       DO bj=myByLo(myThid),myByHi(myThid)
                0206        DO bi=myBxLo(myThid),myBxHi(myThid)
                0207         tileArea(bi,bj) = 0.
                0208         tileFlow(bi,bj) = 0.
                0209         IF ( tileHasOBN(bi,bj) ) THEN
                0210          DO k=1,Nr
                0211           DO i=1,sNx
                0212            jB = OB_Jn(i,bi,bj)
838c7f9401 Jean*0213 C-    If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
                0214 C     communicates with tile interior (sNy+1) rather than with halo region (j=1)
                0215            IF ( jB.NE.OB_indexNone .AND. jB.GT.1 ) THEN
f441182300 Jean*0216             tmpA = drF(k)*hFacS(i,jB,k,bi,bj)*dxG(i,jB,bi,bj)
                0217      &                   *maskInS(i,jB,bi,bj)
47f36df0c2 Jean*0218             tileArea(bi,bj) = tileArea(bi,bj) + tmpA
                0219             tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBNv(i,k,bi,bj)
20c0bcbffa Jean*0220            ENDIF
                0221           ENDDO
                0222          ENDDO
                0223          IF ( OBCS_balanceFacN.GE.0. ) THEN
                0224            tileInFlow(bi,bj) = tileInFlow(bi,bj) - tileFlow(bi,bj)
                0225            tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
                0226      &                       + tileArea(bi,bj)*OBCS_balanceFacN
                0227          ENDIF
                0228          areaN = areaN + tileArea(bi,bj)
                0229          flowN = flowN + tileFlow(bi,bj)
                0230         ENDIF
                0231        ENDDO
                0232       ENDDO
                0233 c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
                0234 c    &   'OBCS_balance it,areaN,flowN=', myIter, areaN, flowN
                0235       IF ( OBCS_balanceFacN.LT.0. ) THEN
                0236         CALL GLOBAL_SUM_TILE_RL( tileArea, areaN, myThid )
                0237         IF ( areaN.GT.0. ) THEN
                0238           CALL GLOBAL_SUM_TILE_RL( tileFlow, flowN, myThid )
8830b8f970 Jean*0239           IF ( debugLevel.GE.debLevC ) THEN
20c0bcbffa Jean*0240             WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
                0241      &            myIter, ' ) correct OBNv:', flowN, -flowN/areaN
                0242             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0243      &                          SQUEEZE_RIGHT, myThid )
                0244           ENDIF
                0245           flowN = -flowN/areaN
                0246         ENDIF
                0247       ENDIF
                0248 #endif /* ALLOW_OBCS_NORTH */
                0249
                0250 #ifdef ALLOW_OBCS_SOUTH
                0251       areaS = 0. _d 0
                0252       flowS = 0. _d 0
                0253       flag = flag .OR. ( OBCS_balanceFacS.GT.0. )
                0254       DO bj=myByLo(myThid),myByHi(myThid)
                0255        DO bi=myBxLo(myThid),myBxHi(myThid)
                0256         tileArea(bi,bj) = 0.
                0257         tileFlow(bi,bj) = 0.
                0258         IF ( tileHasOBS(bi,bj) ) THEN
                0259          DO k=1,Nr
                0260           DO i=1,sNx
                0261            jB = OB_Js(i,bi,bj)
838c7f9401 Jean*0262 C-    If 1 OB location is on 2 tiles (@ edge of 2 tiles), select the one which
                0263 C     communicates with tile interior (j=0) rather than with halo region (j=sNy)
                0264            IF ( jB.NE.OB_indexNone .AND. jB.LT.sNy ) THEN
f441182300 Jean*0265             tmpA = drF(k)*hFacS(i,1+jB,k,bi,bj)*dxG(i,1+jB,bi,bj)
                0266      &                   *maskInS(i,1+jB,bi,bj)
47f36df0c2 Jean*0267             tileArea(bi,bj) = tileArea(bi,bj) + tmpA
                0268             tileFlow(bi,bj) = tileFlow(bi,bj) + tmpA*OBSv(i,k,bi,bj)
20c0bcbffa Jean*0269            ENDIF
                0270           ENDDO
                0271          ENDDO
                0272          IF ( OBCS_balanceFacS.GE.0. ) THEN
                0273            tileInFlow(bi,bj) = tileInFlow(bi,bj) + tileFlow(bi,bj)
                0274            tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
                0275      &                       + tileArea(bi,bj)*OBCS_balanceFacS
                0276          ENDIF
                0277          areaS = areaS + tileArea(bi,bj)
                0278          flowS = flowS + tileFlow(bi,bj)
                0279         ENDIF
                0280        ENDDO
                0281       ENDDO
                0282 c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
                0283 c    &   'OBCS_balance it,areaS,flowS=', myIter, areaS, flowS
                0284       IF ( OBCS_balanceFacS.LT.0. ) THEN
                0285         CALL GLOBAL_SUM_TILE_RL( tileArea, areaS, myThid )
                0286         IF ( areaS.GT.0. ) THEN
                0287           CALL GLOBAL_SUM_TILE_RL( tileFlow, flowS, myThid )
8830b8f970 Jean*0288           IF ( debugLevel.GE.debLevC ) THEN
20c0bcbffa Jean*0289             WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
                0290      &            myIter, ' ) correct OBSv:', flowS, -flowS/areaS
                0291             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0292      &                          SQUEEZE_RIGHT, myThid )
                0293           ENDIF
                0294           flowS = -flowS/areaS
                0295         ENDIF
                0296       ENDIF
                0297 #endif /* ALLOW_OBCS_SOUTH */
                0298
abfe198bce Mart*0299 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0300 C--   Calculate a flux throught the surface from the global change of etaH
                0301 C     to take into account both freshwater fluxes like precipitation and
                0302 C     displacement fluxes from advancing ice sheets or imported sea ice
                0303
                0304       netFreshWaterFlux = 0. _d 0
                0305       shelfIceNetMassFlux = 0. _d 0
                0306       IF ( OBCSbalanceSurf ) THEN
                0307        DO bj=myByLo(myThid),myByHi(myThid)
                0308         DO bi=myBxLo(myThid),myBxHi(myThid)
                0309          tileFlow(bi,bj) = 0.
                0310          DO j=1,sNy
                0311           DO i=1,sNx
                0312            tileFlow(bi,bj) = tileFlow(bi,bj)
                0313      &          - EmPmR(i,j,bi,bj)
                0314      &          * _rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
                0315           ENDDO
                0316          ENDDO
                0317         ENDDO
                0318        ENDDO
                0319        CALL GLOBAL_SUM_TILE_RL( tileFlow, netFreshWaterFlux, myThid )
                0320        IF ( debugLevel.GE.debLevC ) THEN
                0321         WRITE(msgBuf,'(A,I9,A,1P1E16.8)') 'OBCS_balance (it=',
                0322      &       myIter, ' ) correct for netFreshWaterFlux:',
                0323      &       netFreshWaterFlux
                0324         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0325      &       SQUEEZE_RIGHT, myThid )
                0326        ENDIF
                0327 #ifdef ALLOW_SHELFICE
                0328        IF ( useShelfIce ) THEN
                0329         CALL SHELFICE_NETMASSFLUX_SURF(
                0330      O                        shelfIceNetMassFlux,
                0331      I                        myTime, myIter, myThid )
                0332         IF ( debugLevel.GE.debLevC ) THEN
                0333          WRITE(msgBuf,'(A,I9,A,1P1E16.8)') 'OBCS_balance (it=',
                0334      &        myIter, ' ) correct for shelfIceNetMassFlux:',
                0335      &        shelfIceNetMassFlux
                0336          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0337      &        SQUEEZE_RIGHT, myThid )
                0338         ENDIF
                0339        ENDIF
                0340 #endif /* ALLOW_SHELFICE */
                0341       ENDIF
                0342
20c0bcbffa Jean*0343 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0344 C--   Calculate a unique velocity correction for all (OBCS_balanceFac>0) OB
                0345 C     and correct each OB using corresponding weight factor OBCS_balanceFac
                0346
                0347       IF ( flag ) CALL GLOBAL_SUM_TILE_RL( tileAreaOB, areaOB, myThid )
                0348       IF ( areaOB.GT.0. ) THEN
                0349         CALL GLOBAL_SUM_TILE_RL( tileInFlow, inFlow, myThid )
8830b8f970 Jean*0350         IF ( debugLevel.GE.debLevC ) THEN
20c0bcbffa Jean*0351           WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_balance (it=',
                0352      &       myIter, ' ) correct for inFlow:', inFlow, inFlow/areaOB
                0353           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0354      &                        SQUEEZE_RIGHT, myThid )
abfe198bce Mart*0355         ENDIF
                0356         IF ( OBCSbalanceSurf ) THEN
                0357          inFlow = inFlow + netFreshWaterFlux*mass2rUnit
                0358 #ifdef ALLOW_SHELFICE
                0359          IF ( useShelfIce )
                0360      &        inFlow = inFlow + shelfIceNetMassFlux*mass2rUnit
                0361 #endif
20c0bcbffa Jean*0362         ENDIF
                0363         inFlow = inFlow / areaOB
                0364       ENDIF
                0365       IF ( OBCS_balanceFacE.GE.0. ) flowE =  inFlow*OBCS_balanceFacE
                0366       IF ( OBCS_balanceFacW.GE.0. ) flowW = -inFlow*OBCS_balanceFacW
                0367       IF ( OBCS_balanceFacN.GE.0. ) flowN =  inFlow*OBCS_balanceFacN
                0368       IF ( OBCS_balanceFacS.GE.0. ) flowS = -inFlow*OBCS_balanceFacS
                0369
8830b8f970 Jean*0370       IF ( debugLevel.GE.debLevC  .AND. areaOB.GT.0. ) THEN
20c0bcbffa Jean*0371           WRITE(msgBuf,'(A,1P2E16.8)')
                0372      &     'OBCS_balance correction to OBEu,OBWu:', flowE, flowW
                0373           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0374      &                        SQUEEZE_RIGHT, myThid )
                0375           WRITE(msgBuf,'(A,1P2E16.8)')
                0376      &     'OBCS_balance correction to OBNv,OBSv:', flowN, flowS
                0377           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0378      &                        SQUEEZE_RIGHT, myThid )
                0379       ENDIF
                0380
                0381 c     IF ( .NOT.useOBCSbalance ) RETURN
                0382
                0383 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0384 C--   Add correction:
                0385
                0386 #ifdef ALLOW_OBCS_EAST
                0387       IF ( OBCS_balanceFacE.NE.0. ) THEN
                0388        DO bj=myByLo(myThid),myByHi(myThid)
                0389         DO bi=myBxLo(myThid),myBxHi(myThid)
                0390          IF ( tileHasOBE(bi,bj) ) THEN
                0391           DO k=1,Nr
838c7f9401 Jean*0392            DO j=1-OLy,sNy+OLy
20c0bcbffa Jean*0393             iB = OB_Ie(j,bi,bj)
838c7f9401 Jean*0394             IF ( iB.NE.OB_indexNone ) THEN
20c0bcbffa Jean*0395              OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
                0396      &                       + flowE*maskW(iB,j,k,bi,bj)
                0397             ENDIF
                0398            ENDDO
                0399           ENDDO
                0400          ENDIF
                0401         ENDDO
                0402        ENDDO
                0403       ENDIF
                0404 #endif /* ALLOW_OBCS_EAST */
                0405
                0406 #ifdef ALLOW_OBCS_WEST
                0407       IF ( OBCS_balanceFacW.NE.0. ) THEN
                0408        DO bj=myByLo(myThid),myByHi(myThid)
                0409         DO bi=myBxLo(myThid),myBxHi(myThid)
                0410          IF ( tileHasOBW(bi,bj) ) THEN
                0411           DO k=1,Nr
838c7f9401 Jean*0412            DO j=1-OLy,sNy+OLy
20c0bcbffa Jean*0413             iB = OB_Iw(j,bi,bj)
838c7f9401 Jean*0414             IF ( iB.NE.OB_indexNone ) THEN
20c0bcbffa Jean*0415              OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
                0416      &                       + flowW*maskW(1+iB,j,k,bi,bj)
                0417             ENDIF
                0418            ENDDO
                0419           ENDDO
                0420          ENDIF
                0421         ENDDO
                0422        ENDDO
                0423       ENDIF
                0424 #endif /* ALLOW_OBCS_WEST */
                0425
                0426 #ifdef ALLOW_OBCS_NORTH
f2bfd4c10a Jean*0427       IF ( OBCS_balanceFacN.NE.0. ) THEN
20c0bcbffa Jean*0428        DO bj=myByLo(myThid),myByHi(myThid)
                0429         DO bi=myBxLo(myThid),myBxHi(myThid)
                0430          IF ( tileHasOBN(bi,bj) ) THEN
                0431           DO k=1,Nr
838c7f9401 Jean*0432            DO i=1-OLx,sNx+OLx
20c0bcbffa Jean*0433             jB = OB_Jn(i,bi,bj)
838c7f9401 Jean*0434             IF ( jB.NE.OB_indexNone ) THEN
20c0bcbffa Jean*0435              OBNv(i,k,bi,bj) = OBNv(i,k,bi,bj)
                0436      &                       + flowN*maskS(i,jB,k,bi,bj)
                0437             ENDIF
                0438            ENDDO
                0439           ENDDO
                0440          ENDIF
                0441         ENDDO
                0442        ENDDO
                0443       ENDIF
                0444 #endif /* ALLOW_OBCS_NORTH */
                0445
                0446 #ifdef ALLOW_OBCS_SOUTH
                0447       IF ( OBCS_balanceFacS.NE.0. ) THEN
                0448        DO bj=myByLo(myThid),myByHi(myThid)
                0449         DO bi=myBxLo(myThid),myBxHi(myThid)
                0450          IF ( tileHasOBS(bi,bj) ) THEN
                0451           DO k=1,Nr
838c7f9401 Jean*0452            DO i=1-OLx,sNx+OLx
20c0bcbffa Jean*0453             jB = OB_Js(i,bi,bj)
838c7f9401 Jean*0454             IF ( jB.NE.OB_indexNone ) THEN
20c0bcbffa Jean*0455              OBSv(i,k,bi,bj) = OBSv(i,k,bi,bj)
                0456      &                       + flowS*maskS(i,1+jB,k,bi,bj)
                0457             ENDIF
                0458            ENDDO
                0459           ENDDO
                0460          ENDIF
                0461         ENDDO
                0462        ENDDO
                0463       ENDIF
                0464 #endif /* ALLOW_OBCS_SOUTH */
                0465
                0466 #ifdef ALLOW_DEBUG
                0467       IF (debugMode) CALL DEBUG_LEAVE('OBCS_BALANCE_FLOW',myThid)
                0468 #endif
                0469
af61e5eb16 Mart*0470 #endif /* ALLOW_OBCS and ALLOW_OBCS_BALANCE */
20c0bcbffa Jean*0471
                0472       RETURN
                0473       END