C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_set_connect.F,v 1.2 2014/11/25 15:26:10 jmc Exp $
C $Name:  $

#include "OBCS_OPTIONS.h"

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

C     !INTERFACE:
      SUBROUTINE OBCS_SET_CONNECT( myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | SUBROUTINE OBCS_SET_CONNECT
C     | o Set OB connected piece Id for each level
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"

C     !INPUT/OUTPUT PARAMETERS:
C     myThid   :: my Thread Id. number
      INTEGER myThid
CEOP

#ifdef ALLOW_OBCS
C     !LOCAL VARIABLES:
C     msgBuf   :: Informational/error message buffer
C     bi,bj    :: tile indices
C     i, j, k  :: Loop counters
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER bi, bj
      INTEGER i, j, k
      INTEGER idN, idS, idE, idW
      INTEGER fp, prtMsg
      INTEGER n, newConnect, maxConnect
      INTEGER numConnect, listConnect(OBCS_maxConnect)
      INTEGER numLocal(nSx,nSy), listLocal(OBCS_maxConnect,nSx,nSy)
      INTEGER tmpConnect(sNx+sNx+sNy+sNy)
      _RS tmpXZ(1-OLx:sNx+OLx,Nr,nSx,nSy)
      _RS tmpYZ(1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL tmpRL

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

C--   Initialise domain connected-piece Id for OB grid points:
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        idN = 0
        idS = 0
        idE = 0
        idW = 0
        IF ( tileHasOBN(bi,bj) ) idN = 1
        IF ( tileHasOBS(bi,bj) ) idS = 1
        IF ( tileHasOBE(bi,bj) ) idE = 1
        IF ( tileHasOBW(bi,bj) ) idW = 1
        DO k=1,Nr
         DO i=1,sNx
          OBN_connect(i,k,bi,bj) = idN
          OBS_connect(i,k,bi,bj) = idS
         ENDDO
         DO j=1,sNy
          OBE_connect(j,k,bi,bj) = idE
          OBW_connect(j,k,bi,bj) = idW
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Read from files domain connected-piece Id for OB grid points:
C-    Note: current Section READ routines (MDS_READ_SEC_XZ,_YZ) are
C     single (Master) thread (output only available to Master Thread)
      prtMsg = 0
      fp = readBinaryPrec
      _BARRIER
      IF ( OBNconnectFile.NE.' ' ) THEN
        CALL READ_REC_XZ_RS( OBNconnectFile, fp,Nr, tmpXZ, 1,0,myThid )
        _BEGIN_MASTER(myThid)
        DO bj = 1,nSy
         DO bi = 1,nSx
          IF ( tileHasOBN(bi,bj) ) THEN
           DO k=1,Nr
            DO i=1,sNx
             OBN_connect(i,k,bi,bj) = NINT( tmpXZ(i,k,bi,bj) )
            ENDDO
           ENDDO
          ENDIF
         ENDDO
        ENDDO
        _END_MASTER(myThid)
        prtMsg = 1
      ENDIF
      IF ( OBSconnectFile.NE.' ' ) THEN
        CALL READ_REC_XZ_RS( OBSconnectFile, fp,Nr, tmpXZ, 1,0,myThid )
        _BEGIN_MASTER(myThid)
        DO bj = 1,nSy
         DO bi = 1,nSx
          IF ( tileHasOBS(bi,bj) ) THEN
           DO k=1,Nr
            DO i=1,sNx
             OBS_connect(i,k,bi,bj) = NINT( tmpXZ(i,k,bi,bj) )
            ENDDO
           ENDDO
          ENDIF
         ENDDO
        ENDDO
        _END_MASTER(myThid)
        prtMsg = 1
      ENDIF
      IF ( OBEconnectFile.NE.' ' ) THEN
        CALL READ_REC_YZ_RS( OBEconnectFile, fp,Nr, tmpYZ, 1,0,myThid )
        _BEGIN_MASTER(myThid)
        DO bj = 1,nSy
         DO bi = 1,nSx
          IF ( tileHasOBE(bi,bj) ) THEN
           DO k=1,Nr
            DO j=1,sNy
             OBE_connect(j,k,bi,bj) = NINT( tmpYZ(j,k,bi,bj) )
            ENDDO
           ENDDO
          ENDIF
         ENDDO
        ENDDO
        _END_MASTER(myThid)
        prtMsg = 1
      ENDIF
      IF ( OBWconnectFile.NE.' ' ) THEN
        CALL READ_REC_YZ_RS( OBWconnectFile, fp,Nr, tmpYZ, 1,0,myThid )
        _BEGIN_MASTER(myThid)
        DO bj = 1,nSy
         DO bi = 1,nSx
          IF ( tileHasOBW(bi,bj) ) THEN
           DO k=1,Nr
            DO j=1,sNy
             OBW_connect(j,k,bi,bj) = NINT( tmpYZ(j,k,bi,bj) )
            ENDDO
           ENDDO
          ENDIF
         ENDDO
        ENDDO
        _END_MASTER(myThid)
        prtMsg = 1
      ENDIF
      _BARRIER

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
         DO i=1,sNx
          IF (OB_Jn(i,bi,bj).EQ.OB_indexNone) OBN_connect(i,k,bi,bj)=0
          IF (OB_Js(i,bi,bj).EQ.OB_indexNone) OBS_connect(i,k,bi,bj)=0
         ENDDO
         DO j=1,sNy
          IF (OB_Ie(j,bi,bj).EQ.OB_indexNone) OBE_connect(j,k,bi,bj)=0
          IF (OB_Iw(j,bi,bj).EQ.OB_indexNone) OBW_connect(j,k,bi,bj)=0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Count how many connected parts there are  for each level:
      prtMsg = prtMsg*debugLevel
      DO k=1,Nr

        maxConnect = 0
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)

C-    make a local copy
          DO i=1,sNx
            tmpConnect(i) = OBN_connect(i,k,bi,bj)
            tmpConnect(sNx+i) = OBS_connect(i,k,bi,bj)
          ENDDO
          DO j=1,sNy
            tmpConnect(sNx*2+j) = OBW_connect(j,k,bi,bj)
            tmpConnect(sNx*2+sNy+j) = OBE_connect(j,k,bi,bj)
          ENDDO

C-   make a list for each tile
          numLocal(bi,bj) = 0
          DO n=1,OBCS_maxConnect
            listLocal(n,bi,bj) = 0
          ENDDO
          newConnect = 1
          DO WHILE ( newConnect.NE. 0 )
           newConnect = 0
           DO i=1,(sNx+sNy)*2
            IF ( tmpConnect(i).GE.1 ) THEN
              IF ( newConnect.EQ.0 ) THEN
               newConnect = tmpConnect(i)
               numLocal(bi,bj) = numLocal(bi,bj) + 1
               IF ( numLocal(bi,bj).LE.OBCS_maxConnect )
     &           listLocal(numLocal(bi,bj),bi,bj) = newConnect
              ENDIF
              IF ( tmpConnect(i).EQ.newConnect )
     &          tmpConnect(i) = 0
            ENDIF
           ENDDO
          ENDDO
          IF ( numLocal(bi,bj).GT.OBCS_maxConnect ) THEN
            WRITE(msgBuf,'(A,3(A,I4),2(A,I10))') 'OBCS_SET_CONNECT: ',
     &       'k=', k, ' numLocal(', bi,',',bj,')=', numLocal(bi,bj),
     &       ' exceeds OBCS_maxConnect=', OBCS_maxConnect
            CALL PRINT_ERROR( msgBuf, myThid )
            STOP 'ABNORMAL END: S/R OBCS_SET_CONNECT'
          ENDIF
          IF ( prtMsg.GE.debLevC ) THEN
           IF ( numLocal(bi,bj).EQ.0 ) THEN
            WRITE(msgBuf,'(A,2I4,A,I8)') 'OBCS_SET_CONNECT: bi,bj=',
     &       bi, bj, ' , numLocal=', numLocal(bi,bj)
           ELSE
            WRITE(msgBuf,'(A,2I4,2(A,I8))') 'OBCS_SET_CONNECT: bi,bj=',
     &       bi, bj, ' , numLocal=', numLocal(bi,bj),
     &           ' , listLocal:', listLocal(1,bi,bj)
           ENDIF
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
            DO j=2,numLocal(bi,bj),15
              n = MIN(numLocal(bi,bj),j+14)
              WRITE(msgBuf,'(A,15I8)')
     &            ' ... ', (listLocal(i,bi,bj),i=j,n)
              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                            SQUEEZE_RIGHT, myThid )
            ENDDO
          ENDIF
          DO n=1,numLocal(bi,bj)
            maxConnect = MAX( maxConnect, listLocal(n,bi,bj) )
          ENDDO

         ENDDO
        ENDDO

        tmpRL = maxConnect
        _GLOBAL_MAX_RL( tmpRL, myThid )
        maxConnect = NINT( tmpRL )

C-   combine set of list (1 per tile) into 1 global list
        numConnect = 0
        DO j=1,maxConnect
          tmpRL = zeroRL
          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO n=1,numLocal(bi,bj)
             IF ( listLocal(n,bi,bj).EQ.j ) tmpRL = oneRL
            ENDDO
           ENDDO
          ENDDO
          _GLOBAL_MAX_RL( tmpRL, myThid )
          IF ( tmpRL.EQ.oneRL ) THEN
            numConnect = numConnect + 1
            IF ( numConnect.LE.OBCS_maxConnect )
     &        listConnect(numConnect) = j
          ENDIF
        ENDDO
        IF ( numConnect.GT.OBCS_maxConnect ) THEN
          WRITE(msgBuf,'(A,I4,2(A,I10))') 'OBCS_SET_CONNECT: @ k=', k,
     &     ' numConnect=', numConnect,
     &     ' exceeds OBCS_maxConnect=', OBCS_maxConnect
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R OBCS_SET_CONNECT'
        ENDIF
        IF ( prtMsg.GE.debLevA ) THEN
          _BEGIN_MASTER(myThid)
          WRITE(msgBuf,'(A,I4,2(A,I10))') 'OBCS_SET_CONNECT: @ k=', k,
     &     ', maxConnect=', maxConnect, ' , numConnect=', numConnect
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
          DO j=1,numConnect,15
            n = MIN(numConnect,j+14)
            WRITE(msgBuf,'(A,15I8)')
     &          ' listConnect:', (listConnect(i),i=j,n)
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
          ENDDO
          _END_MASTER(myThid)
        ENDIF

C-   reset connected Id in order to use all number from 1 to numConnect
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)

C-    make a local copy
          DO i=1,sNx
            tmpConnect(i) = OBN_connect(i,k,bi,bj)
            tmpConnect(sNx+i) = OBS_connect(i,k,bi,bj)
          ENDDO
          DO j=1,sNy
            tmpConnect(sNx*2+j) = OBW_connect(j,k,bi,bj)
            tmpConnect(sNx*2+sNy+j) = OBE_connect(j,k,bi,bj)
          ENDDO
          DO n=1,numConnect
C-    change Id value: listConnect(n) to n
           IF ( listConnect(n).NE.n ) THEN
            DO i=1,(sNx+sNy)*2
             IF ( tmpConnect(i).EQ.listConnect(n) ) tmpConnect(i) = n
            ENDDO
           ENDIF
          ENDDO
C-    copy back into OB[N,S,E,W]_connect arrays
          DO i=1,sNx
            OBN_connect(i,k,bi,bj) = tmpConnect(i)
            OBS_connect(i,k,bi,bj) = tmpConnect(sNx+i)
          ENDDO
          DO j=1,sNy
            OBW_connect(j,k,bi,bj) = tmpConnect(sNx*2+j)
            OBE_connect(j,k,bi,bj) = tmpConnect(sNx*2+sNy+j)
          ENDDO

         ENDDO
        ENDDO

C-    store numConnect in common block
        _BEGIN_MASTER(myThid)
        OB_connectNumber(k) = numConnect
        _END_MASTER(myThid)

C--   end k loop
      ENDDO

      _BARRIER

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

#endif /* ALLOW_OBCS */
      RETURN
      END