Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:42:41 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
e4a300ec91 Jean*0001 #include "OBCS_OPTIONS.h"
                0002 
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: OBCS_SET_CONNECT
                0006 
                0007 C     !INTERFACE:
                0008       SUBROUTINE OBCS_SET_CONNECT( myThid )
                0009 
                0010 C     !DESCRIPTION:
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE OBCS_SET_CONNECT
                0013 C     | o Set OB connected piece Id for each level
                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"
                0024 #include "OBCS_PARAMS.h"
                0025 #include "OBCS_GRID.h"
                0026 
                0027 C     !INPUT/OUTPUT PARAMETERS:
                0028 C     myThid   :: my Thread Id. number
                0029       INTEGER myThid
                0030 CEOP
                0031 
                0032 #ifdef ALLOW_OBCS
                0033 C     !LOCAL VARIABLES:
                0034 C     msgBuf   :: Informational/error message buffer
                0035 C     bi,bj    :: tile indices
                0036 C     i, j, k  :: Loop counters
                0037       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0038       INTEGER bi, bj
                0039       INTEGER i, j, k
                0040       INTEGER idN, idS, idE, idW
                0041       INTEGER fp, prtMsg
                0042       INTEGER n, newConnect, maxConnect
                0043       INTEGER numConnect, listConnect(OBCS_maxConnect)
                0044       INTEGER numLocal(nSx,nSy), listLocal(OBCS_maxConnect,nSx,nSy)
                0045       INTEGER tmpConnect(sNx+sNx+sNy+sNy)
                0046       _RS tmpXZ(1-OLx:sNx+OLx,Nr,nSx,nSy)
                0047       _RS tmpYZ(1-OLy:sNy+OLy,Nr,nSx,nSy)
                0048       _RL tmpRL
                0049 
                0050 #ifdef ALLOW_DEBUG
                0051       IF (debugMode) CALL DEBUG_ENTER('OBCS_SET_CONNECT',myThid)
                0052 #endif
                0053 
                0054 C--   Initialise domain connected-piece Id for OB grid points:
                0055       DO bj = myByLo(myThid), myByHi(myThid)
                0056        DO bi = myBxLo(myThid), myBxHi(myThid)
                0057         idN = 0
                0058         idS = 0
                0059         idE = 0
                0060         idW = 0
                0061         IF ( tileHasOBN(bi,bj) ) idN = 1
                0062         IF ( tileHasOBS(bi,bj) ) idS = 1
                0063         IF ( tileHasOBE(bi,bj) ) idE = 1
                0064         IF ( tileHasOBW(bi,bj) ) idW = 1
                0065         DO k=1,Nr
                0066          DO i=1,sNx
                0067           OBN_connect(i,k,bi,bj) = idN
                0068           OBS_connect(i,k,bi,bj) = idS
                0069          ENDDO
                0070          DO j=1,sNy
                0071           OBE_connect(j,k,bi,bj) = idE
                0072           OBW_connect(j,k,bi,bj) = idW
                0073          ENDDO
                0074         ENDDO
                0075        ENDDO
                0076       ENDDO
                0077 
                0078 C--   Read from files domain connected-piece Id for OB grid points:
b5bd2d619e Jean*0079 C-    Note: current Section READ routines (MDS_READ_SEC_XZ,_YZ) are
                0080 C     single (Master) thread (output only available to Master Thread)
e4a300ec91 Jean*0081       prtMsg = 0
                0082       fp = readBinaryPrec
b5bd2d619e Jean*0083       _BARRIER
e4a300ec91 Jean*0084       IF ( OBNconnectFile.NE.' ' ) THEN
                0085         CALL READ_REC_XZ_RS( OBNconnectFile, fp,Nr, tmpXZ, 1,0,myThid )
b5bd2d619e Jean*0086         _BEGIN_MASTER(myThid)
                0087         DO bj = 1,nSy
                0088          DO bi = 1,nSx
e4a300ec91 Jean*0089           IF ( tileHasOBN(bi,bj) ) THEN
                0090            DO k=1,Nr
                0091             DO i=1,sNx
                0092              OBN_connect(i,k,bi,bj) = NINT( tmpXZ(i,k,bi,bj) )
                0093             ENDDO
                0094            ENDDO
                0095           ENDIF
                0096          ENDDO
                0097         ENDDO
b5bd2d619e Jean*0098         _END_MASTER(myThid)
e4a300ec91 Jean*0099         prtMsg = 1
                0100       ENDIF
                0101       IF ( OBSconnectFile.NE.' ' ) THEN
                0102         CALL READ_REC_XZ_RS( OBSconnectFile, fp,Nr, tmpXZ, 1,0,myThid )
b5bd2d619e Jean*0103         _BEGIN_MASTER(myThid)
                0104         DO bj = 1,nSy
                0105          DO bi = 1,nSx
e4a300ec91 Jean*0106           IF ( tileHasOBS(bi,bj) ) THEN
                0107            DO k=1,Nr
                0108             DO i=1,sNx
                0109              OBS_connect(i,k,bi,bj) = NINT( tmpXZ(i,k,bi,bj) )
                0110             ENDDO
                0111            ENDDO
                0112           ENDIF
                0113          ENDDO
                0114         ENDDO
b5bd2d619e Jean*0115         _END_MASTER(myThid)
e4a300ec91 Jean*0116         prtMsg = 1
                0117       ENDIF
                0118       IF ( OBEconnectFile.NE.' ' ) THEN
                0119         CALL READ_REC_YZ_RS( OBEconnectFile, fp,Nr, tmpYZ, 1,0,myThid )
b5bd2d619e Jean*0120         _BEGIN_MASTER(myThid)
                0121         DO bj = 1,nSy
                0122          DO bi = 1,nSx
e4a300ec91 Jean*0123           IF ( tileHasOBE(bi,bj) ) THEN
                0124            DO k=1,Nr
                0125             DO j=1,sNy
                0126              OBE_connect(j,k,bi,bj) = NINT( tmpYZ(j,k,bi,bj) )
                0127             ENDDO
                0128            ENDDO
                0129           ENDIF
                0130          ENDDO
                0131         ENDDO
b5bd2d619e Jean*0132         _END_MASTER(myThid)
e4a300ec91 Jean*0133         prtMsg = 1
                0134       ENDIF
                0135       IF ( OBWconnectFile.NE.' ' ) THEN
                0136         CALL READ_REC_YZ_RS( OBWconnectFile, fp,Nr, tmpYZ, 1,0,myThid )
b5bd2d619e Jean*0137         _BEGIN_MASTER(myThid)
                0138         DO bj = 1,nSy
                0139          DO bi = 1,nSx
e4a300ec91 Jean*0140           IF ( tileHasOBW(bi,bj) ) THEN
                0141            DO k=1,Nr
                0142             DO j=1,sNy
                0143              OBW_connect(j,k,bi,bj) = NINT( tmpYZ(j,k,bi,bj) )
                0144             ENDDO
                0145            ENDDO
                0146           ENDIF
                0147          ENDDO
                0148         ENDDO
b5bd2d619e Jean*0149         _END_MASTER(myThid)
e4a300ec91 Jean*0150         prtMsg = 1
                0151       ENDIF
b5bd2d619e Jean*0152       _BARRIER
e4a300ec91 Jean*0153 
                0154       DO bj = myByLo(myThid), myByHi(myThid)
                0155        DO bi = myBxLo(myThid), myBxHi(myThid)
                0156         DO k=1,Nr
                0157          DO i=1,sNx
                0158           IF (OB_Jn(i,bi,bj).EQ.OB_indexNone) OBN_connect(i,k,bi,bj)=0
                0159           IF (OB_Js(i,bi,bj).EQ.OB_indexNone) OBS_connect(i,k,bi,bj)=0
                0160          ENDDO
                0161          DO j=1,sNy
                0162           IF (OB_Ie(j,bi,bj).EQ.OB_indexNone) OBE_connect(j,k,bi,bj)=0
                0163           IF (OB_Iw(j,bi,bj).EQ.OB_indexNone) OBW_connect(j,k,bi,bj)=0
                0164          ENDDO
                0165         ENDDO
                0166        ENDDO
                0167       ENDDO
                0168 
                0169 C--   Count how many connected parts there are  for each level:
                0170       prtMsg = prtMsg*debugLevel
                0171       DO k=1,Nr
                0172 
                0173         maxConnect = 0
                0174         DO bj = myByLo(myThid), myByHi(myThid)
                0175          DO bi = myBxLo(myThid), myBxHi(myThid)
                0176 
                0177 C-    make a local copy
                0178           DO i=1,sNx
                0179             tmpConnect(i) = OBN_connect(i,k,bi,bj)
                0180             tmpConnect(sNx+i) = OBS_connect(i,k,bi,bj)
                0181           ENDDO
                0182           DO j=1,sNy
                0183             tmpConnect(sNx*2+j) = OBW_connect(j,k,bi,bj)
                0184             tmpConnect(sNx*2+sNy+j) = OBE_connect(j,k,bi,bj)
                0185           ENDDO
                0186 
                0187 C-   make a list for each tile
                0188           numLocal(bi,bj) = 0
                0189           DO n=1,OBCS_maxConnect
                0190             listLocal(n,bi,bj) = 0
                0191           ENDDO
                0192           newConnect = 1
                0193           DO WHILE ( newConnect.NE. 0 )
                0194            newConnect = 0
                0195            DO i=1,(sNx+sNy)*2
                0196             IF ( tmpConnect(i).GE.1 ) THEN
                0197               IF ( newConnect.EQ.0 ) THEN
                0198                newConnect = tmpConnect(i)
                0199                numLocal(bi,bj) = numLocal(bi,bj) + 1
                0200                IF ( numLocal(bi,bj).LE.OBCS_maxConnect )
                0201      &           listLocal(numLocal(bi,bj),bi,bj) = newConnect
                0202               ENDIF
                0203               IF ( tmpConnect(i).EQ.newConnect )
                0204      &          tmpConnect(i) = 0
                0205             ENDIF
                0206            ENDDO
                0207           ENDDO
                0208           IF ( numLocal(bi,bj).GT.OBCS_maxConnect ) THEN
                0209             WRITE(msgBuf,'(A,3(A,I4),2(A,I10))') 'OBCS_SET_CONNECT: ',
                0210      &       'k=', k, ' numLocal(', bi,',',bj,')=', numLocal(bi,bj),
                0211      &       ' exceeds OBCS_maxConnect=', OBCS_maxConnect
                0212             CALL PRINT_ERROR( msgBuf, myThid )
                0213             STOP 'ABNORMAL END: S/R OBCS_SET_CONNECT'
                0214           ENDIF
                0215           IF ( prtMsg.GE.debLevC ) THEN
                0216            IF ( numLocal(bi,bj).EQ.0 ) THEN
                0217             WRITE(msgBuf,'(A,2I4,A,I8)') 'OBCS_SET_CONNECT: bi,bj=',
                0218      &       bi, bj, ' , numLocal=', numLocal(bi,bj)
                0219            ELSE
                0220             WRITE(msgBuf,'(A,2I4,2(A,I8))') 'OBCS_SET_CONNECT: bi,bj=',
                0221      &       bi, bj, ' , numLocal=', numLocal(bi,bj),
                0222      &           ' , listLocal:', listLocal(1,bi,bj)
                0223            ENDIF
                0224             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0225      &                          SQUEEZE_RIGHT, myThid )
                0226             DO j=2,numLocal(bi,bj),15
                0227               n = MIN(numLocal(bi,bj),j+14)
                0228               WRITE(msgBuf,'(A,15I8)')
                0229      &            ' ... ', (listLocal(i,bi,bj),i=j,n)
                0230               CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0231      &                            SQUEEZE_RIGHT, myThid )
                0232             ENDDO
                0233           ENDIF
                0234           DO n=1,numLocal(bi,bj)
                0235             maxConnect = MAX( maxConnect, listLocal(n,bi,bj) )
                0236           ENDDO
                0237 
                0238          ENDDO
                0239         ENDDO
                0240 
                0241         tmpRL = maxConnect
                0242         _GLOBAL_MAX_RL( tmpRL, myThid )
                0243         maxConnect = NINT( tmpRL )
                0244 
                0245 C-   combine set of list (1 per tile) into 1 global list
                0246         numConnect = 0
                0247         DO j=1,maxConnect
                0248           tmpRL = zeroRL
                0249           DO bj = myByLo(myThid), myByHi(myThid)
                0250            DO bi = myBxLo(myThid), myBxHi(myThid)
                0251             DO n=1,numLocal(bi,bj)
                0252              IF ( listLocal(n,bi,bj).EQ.j ) tmpRL = oneRL
                0253             ENDDO
                0254            ENDDO
                0255           ENDDO
                0256           _GLOBAL_MAX_RL( tmpRL, myThid )
                0257           IF ( tmpRL.EQ.oneRL ) THEN
                0258             numConnect = numConnect + 1
                0259             IF ( numConnect.LE.OBCS_maxConnect )
                0260      &        listConnect(numConnect) = j
                0261           ENDIF
                0262         ENDDO
                0263         IF ( numConnect.GT.OBCS_maxConnect ) THEN
                0264           WRITE(msgBuf,'(A,I4,2(A,I10))') 'OBCS_SET_CONNECT: @ k=', k,
                0265      &     ' numConnect=', numConnect,
                0266      &     ' exceeds OBCS_maxConnect=', OBCS_maxConnect
                0267           CALL PRINT_ERROR( msgBuf, myThid )
                0268           STOP 'ABNORMAL END: S/R OBCS_SET_CONNECT'
                0269         ENDIF
                0270         IF ( prtMsg.GE.debLevA ) THEN
b5bd2d619e Jean*0271           _BEGIN_MASTER(myThid)
e4a300ec91 Jean*0272           WRITE(msgBuf,'(A,I4,2(A,I10))') 'OBCS_SET_CONNECT: @ k=', k,
                0273      &     ', maxConnect=', maxConnect, ' , numConnect=', numConnect
                0274           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0275      &                        SQUEEZE_RIGHT, myThid )
                0276           DO j=1,numConnect,15
                0277             n = MIN(numConnect,j+14)
                0278             WRITE(msgBuf,'(A,15I8)')
                0279      &          ' listConnect:', (listConnect(i),i=j,n)
                0280             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0281      &                          SQUEEZE_RIGHT, myThid )
                0282           ENDDO
b5bd2d619e Jean*0283           _END_MASTER(myThid)
e4a300ec91 Jean*0284         ENDIF
                0285 
                0286 C-   reset connected Id in order to use all number from 1 to numConnect
                0287         DO bj = myByLo(myThid), myByHi(myThid)
                0288          DO bi = myBxLo(myThid), myBxHi(myThid)
                0289 
                0290 C-    make a local copy
                0291           DO i=1,sNx
                0292             tmpConnect(i) = OBN_connect(i,k,bi,bj)
                0293             tmpConnect(sNx+i) = OBS_connect(i,k,bi,bj)
                0294           ENDDO
                0295           DO j=1,sNy
                0296             tmpConnect(sNx*2+j) = OBW_connect(j,k,bi,bj)
                0297             tmpConnect(sNx*2+sNy+j) = OBE_connect(j,k,bi,bj)
                0298           ENDDO
                0299           DO n=1,numConnect
                0300 C-    change Id value: listConnect(n) to n
                0301            IF ( listConnect(n).NE.n ) THEN
                0302             DO i=1,(sNx+sNy)*2
                0303              IF ( tmpConnect(i).EQ.listConnect(n) ) tmpConnect(i) = n
                0304             ENDDO
                0305            ENDIF
                0306           ENDDO
                0307 C-    copy back into OB[N,S,E,W]_connect arrays
                0308           DO i=1,sNx
                0309             OBN_connect(i,k,bi,bj) = tmpConnect(i)
                0310             OBS_connect(i,k,bi,bj) = tmpConnect(sNx+i)
                0311           ENDDO
                0312           DO j=1,sNy
                0313             OBW_connect(j,k,bi,bj) = tmpConnect(sNx*2+j)
                0314             OBE_connect(j,k,bi,bj) = tmpConnect(sNx*2+sNy+j)
                0315           ENDDO
                0316 
                0317          ENDDO
                0318         ENDDO
                0319 
                0320 C-    store numConnect in common block
                0321         _BEGIN_MASTER(myThid)
                0322         OB_connectNumber(k) = numConnect
                0323         _END_MASTER(myThid)
                0324 
                0325 C--   end k loop
                0326       ENDDO
                0327 
                0328       _BARRIER
                0329 
                0330 #ifdef ALLOW_DEBUG
                0331       IF (debugMode) CALL DEBUG_LEAVE('OBCS_SET_CONNECT',myThid)
                0332 #endif
                0333 
                0334 #endif /* ALLOW_OBCS */
                0335       RETURN
                0336       END