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