C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_diag_balance.F,v 1.3 2012/09/18 00:07:36 jmc Exp $
C $Name:  $

#include "OBCS_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: OBCS_DIAG_BALANCE

C     !INTERFACE:
      SUBROUTINE OBCS_DIAG_BALANCE(
     I                              k,
     U                              uTrans, vTrans,
     I                              myTime, myIter, myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE OBCS_DIAG_BALANCE
C     | o For diagnostics purpose, modify transport normal to OB
C     |   to ensure no net inflow
C     | o use same scheme and same parameters as OBCS_BALANCE_FLOW
C     *==========================================================*

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     k       :: current level index
C     uTrans  :: horizontal transport to balance [m^3/s]
C     vTrans  :: horizontal transport to balance [m^3/s]
C     myTime  :: current time of simulation (s)
C     myIter  :: current iteration number
C     myThid  :: my Thread Id number
      INTEGER k
      _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_OBCS
#ifdef ALLOW_DIAGNOSTICS

C     !FUNCTIONS:

C     !LOCAL VARIABLES:
C     bi, bj       :: tile indices
C     i,j,k        :: loop indices
C     iB, jB       :: local index of open boundary
C     msgBuf       :: Informational/error message buffer
      INTEGER bi, bj
      INTEGER i, j, iB, jB
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      _RL areaOB, areaE, areaW, areaN, areaS, tmpA
      _RL inFlow, flowE, flowW, flowN, flowS
      _RL tileArea(nSx,nSy)
      _RL tileFlow(nSx,nSy)
      _RL tileAreaOB(nSx,nSy)
      _RL tileInFlow(nSx,nSy)
      LOGICAL flag

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Old method (OBCS_balanceFac < 0): balance each OB separately
C-- New method applies to all OB with BCS_balanceFac >=0 :
C      ensure that the net inflow through all OB is balanced by correcting
C      each OB normal flow with a uniform velocity, using the corresponding
C      weight factor OBCS_balanceFac.
C  e.g., OBCS_balanceFac_E,W,N,S= 1, -1, 2, 0 :
C      => correct Western OBWu (by substracting a uniform velocity) to ensure
C         zero net transport through Western OB
C      => correct Eastern and Northern normal flow (twice larger Northern
C         velocity correction than Eastern correction) to ensure that
C         the total inflow through E+N+S OB is balanced

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

C--   Integrate the transport through each OB
      flag = .FALSE.
      areaOB = 0. _d 0
      inFlow = 0. _d 0
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        tileAreaOB(bi,bj) = 0.
        tileInFlow(bi,bj) = 0.
       ENDDO
      ENDDO

#ifdef ALLOW_OBCS_EAST
      areaE = 0. _d 0
      flowE = 0. _d 0
      flag = flag .OR. ( OBCS_balanceFacE.GT.0. )
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        tileArea(bi,bj) = 0.
        tileFlow(bi,bj) = 0.
        IF ( tileHasOBE(bi,bj) ) THEN
         DO j=1,sNy
           iB = OB_Ie(j,bi,bj)
           IF ( iB.NE.OB_indexNone .AND. iB.GT.1 ) THEN
            tmpA = drF(k)*hFacW(iB,j,k,bi,bj)*dyG(iB,j,bi,bj)
     &                   *maskInW(iB,j,bi,bj)
            tileArea(bi,bj) = tileArea(bi,bj) + tmpA
            tileFlow(bi,bj) = tileFlow(bi,bj) + uTrans(iB,j,bi,bj)
           ENDIF
         ENDDO
         IF ( OBCS_balanceFacE.GE.0. ) THEN
           tileInFlow(bi,bj) = tileInFlow(bi,bj) - tileFlow(bi,bj)
           tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
     &                       + tileArea(bi,bj)*OBCS_balanceFacE
         ENDIF
         areaE = areaE + tileArea(bi,bj)
         flowE = flowE + tileFlow(bi,bj)
        ENDIF
       ENDDO
      ENDDO
c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
c    &   'OBCS_balance it,areaE,flowE=', myIter, areaE, flowE
      IF ( OBCS_balanceFacE.LT.0. ) THEN
        CALL GLOBAL_SUM_TILE_RL( tileArea, areaE, myThid )
        IF ( areaE.GT.0. ) THEN
          CALL GLOBAL_SUM_TILE_RL( tileFlow, flowE, myThid )
          IF ( debugLevel.GE.debLevC ) THEN
            WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_DIAG_balance (it=',
     &            myIter, ' ) correct OBEu:', flowE, -flowE/areaE
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
          ENDIF
          flowE = -flowE/areaE
        ENDIF
      ENDIF
#endif /* ALLOW_OBCS_EAST */

#ifdef ALLOW_OBCS_WEST
      areaW = 0. _d 0
      flowW = 0. _d 0
      flag = flag .OR. ( OBCS_balanceFacW.GT.0. )
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        tileArea(bi,bj) = 0.
        tileFlow(bi,bj) = 0.
        IF ( tileHasOBW(bi,bj) ) THEN
         DO j=1,sNy
           iB = OB_Iw(j,bi,bj)
           IF ( iB.NE.OB_indexNone .AND. iB.LT.sNx ) THEN
            iB = 1+iB
            tmpA = drF(k)*hFacW(iB,j,k,bi,bj)*dyG(iB,j,bi,bj)
     &                   *maskInW(iB,j,bi,bj)
            tileArea(bi,bj) = tileArea(bi,bj) + tmpA
            tileFlow(bi,bj) = tileFlow(bi,bj) + uTrans(iB,j,bi,bj)
           ENDIF
         ENDDO
         IF ( OBCS_balanceFacW.GE.0. ) THEN
           tileInFlow(bi,bj) = tileInFlow(bi,bj) + tileFlow(bi,bj)
           tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
     &                       + tileArea(bi,bj)*OBCS_balanceFacW
         ENDIF
         areaW = areaW + tileArea(bi,bj)
         flowW = flowW + tileFlow(bi,bj)
        ENDIF
       ENDDO
      ENDDO
c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
c    &   'OBCS_balance it,areaW,flowW=', myIter, areaW, flowW
      IF ( OBCS_balanceFacW.LT.0. ) THEN
        CALL GLOBAL_SUM_TILE_RL( tileArea, areaW, myThid )
        IF ( areaW.GT.0. ) THEN
          CALL GLOBAL_SUM_TILE_RL( tileFlow, flowW, myThid )
          IF ( debugLevel.GE.debLevC ) THEN
            WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_DIAG_balance (it=',
     &            myIter, ' ) correct OBWu:', flowW, -flowW/areaW
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
          ENDIF
          flowW = -flowW/areaW
        ENDIF
      ENDIF
#endif /* ALLOW_OBCS_WEST */

#ifdef ALLOW_OBCS_NORTH
      areaN = 0. _d 0
      flowN = 0. _d 0
      flag = flag .OR. ( OBCS_balanceFacN.GT.0. )
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        tileArea(bi,bj) = 0.
        tileFlow(bi,bj) = 0.
        IF ( tileHasOBN(bi,bj) ) THEN
         DO i=1,sNx
           jB = OB_Jn(i,bi,bj)
           IF ( jB.NE.OB_indexNone .AND. jB.GT.1 ) THEN
            tmpA = drF(k)*hFacS(i,jB,k,bi,bj)*dxG(i,jB,bi,bj)
     &                   *maskInS(i,jB,bi,bj)
            tileArea(bi,bj) = tileArea(bi,bj) + tmpA
            tileFlow(bi,bj) = tileFlow(bi,bj) + vTrans(i,jB,bi,bj)
           ENDIF
         ENDDO
         IF ( OBCS_balanceFacN.GE.0. ) THEN
           tileInFlow(bi,bj) = tileInFlow(bi,bj) - tileFlow(bi,bj)
           tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
     &                       + tileArea(bi,bj)*OBCS_balanceFacN
         ENDIF
         areaN = areaN + tileArea(bi,bj)
         flowN = flowN + tileFlow(bi,bj)
        ENDIF
       ENDDO
      ENDDO
c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
c    &   'OBCS_balance it,areaN,flowN=', myIter, areaN, flowN
      IF ( OBCS_balanceFacN.LT.0. ) THEN
        CALL GLOBAL_SUM_TILE_RL( tileArea, areaN, myThid )
        IF ( areaN.GT.0. ) THEN
          CALL GLOBAL_SUM_TILE_RL( tileFlow, flowN, myThid )
          IF ( debugLevel.GE.debLevC ) THEN
            WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_DIAG_balance (it=',
     &            myIter, ' ) correct OBNv:', flowN, -flowN/areaN
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
          ENDIF
          flowN = -flowN/areaN
        ENDIF
      ENDIF
#endif /* ALLOW_OBCS_NORTH */

#ifdef ALLOW_OBCS_SOUTH
      areaS = 0. _d 0
      flowS = 0. _d 0
      flag = flag .OR. ( OBCS_balanceFacS.GT.0. )
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        tileArea(bi,bj) = 0.
        tileFlow(bi,bj) = 0.
        IF ( tileHasOBS(bi,bj) ) THEN
         DO i=1,sNx
           jB = OB_Js(i,bi,bj)
           IF ( jB.NE.OB_indexNone .AND. jB.LT.sNy ) THEN
            jB = 1+jB
            tmpA = drF(k)*hFacS(i,jB,k,bi,bj)*dxG(i,jB,bi,bj)
     &                   *maskInS(i,jB,bi,bj)
            tileArea(bi,bj) = tileArea(bi,bj) + tmpA
            tileFlow(bi,bj) = tileFlow(bi,bj) + vTrans(i,jB,bi,bj)
           ENDIF
         ENDDO
         IF ( OBCS_balanceFacS.GE.0. ) THEN
           tileInFlow(bi,bj) = tileInFlow(bi,bj) + tileFlow(bi,bj)
           tileAreaOB(bi,bj) = tileAreaOB(bi,bj)
     &                       + tileArea(bi,bj)*OBCS_balanceFacS
         ENDIF
         areaS = areaS + tileArea(bi,bj)
         flowS = flowS + tileFlow(bi,bj)
        ENDIF
       ENDDO
      ENDDO
c     WRITE(standardMessageUnit,'(A,I9,1P2E16.8)')
c    &   'OBCS_balance it,areaS,flowS=', myIter, areaS, flowS
      IF ( OBCS_balanceFacS.LT.0. ) THEN
        CALL GLOBAL_SUM_TILE_RL( tileArea, areaS, myThid )
        IF ( areaS.GT.0. ) THEN
          CALL GLOBAL_SUM_TILE_RL( tileFlow, flowS, myThid )
          IF ( debugLevel.GE.debLevC ) THEN
            WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_DIAG_balance (it=',
     &            myIter, ' ) correct OBSv:', flowS, -flowS/areaS
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
          ENDIF
          flowS = -flowS/areaS
        ENDIF
      ENDIF
#endif /* ALLOW_OBCS_SOUTH */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   Calculate a unique velocity correction for all (OBCS_balanceFac>0) OB
C     and correct each OB using corresponding weight factor OBCS_balanceFac

      IF ( flag ) CALL GLOBAL_SUM_TILE_RL( tileAreaOB, areaOB, myThid )
      IF ( areaOB.GT.0. ) THEN
        CALL GLOBAL_SUM_TILE_RL( tileInFlow, inFlow, myThid )
        IF ( debugLevel.GE.debLevC ) THEN
          WRITE(msgBuf,'(A,I9,A,1P2E16.8)') 'OBCS_DIAG_balance (it=',
     &       myIter, ' ) correct for inFlow:', inFlow, inFlow/areaOB
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
        ENDIF
        inFlow = inFlow / areaOB
      ENDIF
      IF ( OBCS_balanceFacE.GE.0. ) flowE =  inFlow*OBCS_balanceFacE
      IF ( OBCS_balanceFacW.GE.0. ) flowW = -inFlow*OBCS_balanceFacW
      IF ( OBCS_balanceFacN.GE.0. ) flowN =  inFlow*OBCS_balanceFacN
      IF ( OBCS_balanceFacS.GE.0. ) flowS = -inFlow*OBCS_balanceFacS

      IF ( debugLevel.GE.debLevC  .AND. areaOB.GT.0. ) THEN
          WRITE(msgBuf,'(A,1P2E16.8)')
     &     'OBCS_DIAG_balance correction to uTrans:', flowE, flowW
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
          WRITE(msgBuf,'(A,1P2E16.8)')
     &     'OBCS_DIAG_balance correction to vTrans:', flowN, flowS
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   Add correction:

#ifdef ALLOW_OBCS_EAST
      IF ( OBCS_balanceFacE.NE.0. ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         IF ( tileHasOBE(bi,bj) ) THEN
           DO j=1-OLy,sNy+OLy
            IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
             iB = OB_Ie(j,bi,bj)
             uTrans(iB,j,bi,bj) = uTrans(iB,j,bi,bj)
     &                    + flowE*drF(k)*hFacW(iB,j,k,bi,bj)
     &                           *dyG(iB,j,bi,bj)*maskInW(iB,j,bi,bj)
            ENDIF
           ENDDO
         ENDIF
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_OBCS_EAST */

#ifdef ALLOW_OBCS_WEST
      IF ( OBCS_balanceFacW.NE.0. ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         IF ( tileHasOBW(bi,bj) ) THEN
           DO j=1-OLy,sNy+OLy
            IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
             iB = 1+OB_Iw(j,bi,bj)
             uTrans(iB,j,bi,bj) = uTrans(iB,j,bi,bj)
     &                    + flowW*drF(k)*hFacW(iB,j,k,bi,bj)
     &                           *dyG(iB,j,bi,bj)*maskInW(iB,j,bi,bj)
            ENDIF
           ENDDO
         ENDIF
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_OBCS_WEST */

#ifdef ALLOW_OBCS_NORTH
      IF ( OBCS_balanceFacW.NE.0. ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         IF ( tileHasOBN(bi,bj) ) THEN
           DO i=1-OLx,sNx+OLx
            IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
             jB = OB_Jn(i,bi,bj)
             vTrans(i,jB,bi,bj) = vTrans(i,jB,bi,bj)
     &                    + flowN*drF(k)*hFacS(i,jB,k,bi,bj)
     &                           *dxG(i,jB,bi,bj)*maskInS(i,jB,bi,bj)
            ENDIF
           ENDDO
         ENDIF
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_OBCS_NORTH */

#ifdef ALLOW_OBCS_SOUTH
      IF ( OBCS_balanceFacS.NE.0. ) THEN
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         IF ( tileHasOBS(bi,bj) ) THEN
           DO i=1-OLx,sNx+OLx
            IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
             jB = 1+OB_Js(i,bi,bj)
             vTrans(i,jB,bi,bj) = vTrans(i,jB,bi,bj)
     &                    + flowS*drF(k)*hFacS(i,jB,k,bi,bj)
     &                           *dxG(i,jB,bi,bj)*maskInS(i,jB,bi,bj)
            ENDIF
           ENDDO
         ENDIF
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_OBCS_SOUTH */

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

#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_OBCS */

      RETURN
      END