C $Header: /u/gcmpack/MITgcm/pkg/exch2/w2_set_map_cumsum.F,v 1.2 2011/09/21 16:27:42 jmc Exp $
C $Name:  $

#include "CPP_EEOPTIONS.h"
#include "W2_OPTIONS.h"

CBOP 0
C !ROUTINE: W2_SET_MAP_CUMSUM

C !INTERFACE:
      SUBROUTINE W2_SET_MAP_CUMSUM( myThid )

C     !DESCRIPTION:
C     Set-up mapping for global cumulated sum.

C     !USES:
      IMPLICIT NONE

C      Tile topology settings data structures
#include "SIZE.h"
#include "EEPARAMS.h"
#include "W2_EXCH2_SIZE.h"
#include "W2_EXCH2_PARAMS.h"
#include "W2_EXCH2_TOPOLOGY.h"

C     !INPUT PARAMETERS:
C     myThid  :: my Thread Id number
C               (Note: not relevant since threading has not yet started)
      INTEGER myThid

C     !LOCAL VARIABLES:
C     === Local variables ===
C     facetXYSum   :: Tile to Facet Matrix for facet-Increment in X & Y dir.
C     facet_CSum   :: Tile to Facet Matrix for CumSum @ facet-origin
C     msgBuf       :: Informational/error message buffer
      INTEGER fNx, fNy, nbTx, nbTy
      INTEGER nActiveFacets
      INTEGER fCnt, prev_fCnt
      INTEGER npass, nType
      LOGICAL prtFlag
      LOGICAL fIsSet(0:W2_maxNbFacets)
      INTEGER tN, i, j, k, ii, jj
#ifdef W2_CUMSUM_USE_MATRIX
      INTEGER tS, bi, bj, l, is, ie
      INTEGER facetXYSum(2,W2_maxNbTiles,W2_maxNbFacets)
      INTEGER facet_CSum(2,W2_maxNbTiles,W2_maxNbFacets)
#endif
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
c      DATA edge / 'N' , 'S' , 'E' , 'W' /

      WRITE(msgBuf,'(2A)') 'W2_SET_MAP_CUMSUM: ',
     &       'setting Facet Matrix for CUMUL-SUM'
      CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
      prtFlag = ABS(W2_printMsg).GE.2
     &       .OR. ( W2_printMsg .NE.0 .AND. myProcId.EQ.0 )

C--   Initialise Common-block
      W2_tMC1 = 0
      W2_tMC2 = 0
      DO j=1,W2_maxNbFacets
       DO i=1,W2_maxNbFacets
         W2_cumSum_facet(1,i,j) = 0
         W2_cumSum_facet(2,i,j) = 0
       ENDDO
      ENDDO
#ifdef W2_CUMSUM_USE_MATRIX
      DO j=1,W2_maxNbTiles
       DO i=1,W2_maxNbTiles
         W2_cumSum_tiles(1,i,j) = 0
         W2_cumSum_tiles(2,i,j) = 0
       ENDDO
      ENDDO
#endif /* W2_CUMSUM_USE_MATRIX */

C--   Start setting cumul-sum Face mapping
      fCnt = 0
      fIsSet(0) = .TRUE.
      DO j=1,nFacets
        fIsSet(j) = .FALSE.
      ENDDO

C--   Start with first non-empty face:
      nActiveFacets = 0
      DO j=1,nFacets
        IF ( facet_dims(2*j-1)*facet_dims(2*j).GE.1 ) THEN
          nActiveFacets = nActiveFacets + 1
          IF ( fCnt.EQ.0 ) THEN
            fIsSet(j) = .TRUE.
            fCnt = 1
            IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(A,I4)')
     &        ' CumSum starts @ SW.corner of facet #', j
          ENDIF
        ENDIF
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Go through list of connections:
      prev_fCnt = 0
      npass = 0
      DO WHILE ( fCnt.GT.prev_fCnt )
       npass = npass + 1
       prev_fCnt = fCnt
       DO nType=1,3
        IF ( fCnt.EQ.prev_fCnt ) THEN
         DO j=1,nFacets
          IF ( fIsSet(j) ) THEN
           DO i=1,4
C-    connected to:
            jj = INT(facet_link(i,j))
            ii = MOD( NINT(facet_link(i,j)*10.), 10 )
C--   1) with same orientation ( N <-> S or E <-> W ), forward progression
            IF ( .NOT.fIsSet(jj) .AND. nType.EQ.1 ) THEN
C-    N <-- S ;
             IF ( i.EQ.W2_NORTH .AND. ii.EQ.W2_SOUTH ) THEN
              DO k=1,nFacets
               W2_cumSum_facet(1,k,jj) = W2_cumSum_facet(1,k,j)
               W2_cumSum_facet(2,k,jj) = W2_cumSum_facet(2,k,j)
              ENDDO
              W2_cumSum_facet(2,j,jj) = W2_cumSum_facet(2,j,jj) + 1
              fCnt = fCnt + 1
              fIsSet(jj) = .TRUE.
              IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(5(A,I4))')
     &         ' CumSum SW.corner of facet #', jj,' set from facet',j,
     &         ' (pass,type=', npass,',', nType,')'
             ENDIF
C-    E <-- W ;
             IF ( i.EQ.W2_EAST .AND. ii.EQ.W2_WEST ) THEN
              DO k=1,nFacets
               W2_cumSum_facet(1,k,jj) = W2_cumSum_facet(1,k,j)
               W2_cumSum_facet(2,k,jj) = W2_cumSum_facet(2,k,j)
              ENDDO
              W2_cumSum_facet(1,j,jj) = W2_cumSum_facet(1,j,jj) + 1
              fCnt = fCnt + 1
              fIsSet(jj) = .TRUE.
              IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(5(A,I4))')
     &         ' CumSum SW.corner of facet #', jj,' set from facet',j,
     &         ' (pass,type=', npass,',', nType,')'
             ENDIF
            ENDIF
C--   2) with same orientation ( N <-> S or E <-> W ), backward progression
            IF ( .NOT.fIsSet(jj) .AND. nType.EQ.2 ) THEN
C-    S <-- N ;
             IF ( i.EQ.W2_SOUTH .AND. ii.EQ.W2_NORTH ) THEN
              DO k=1,nFacets
               W2_cumSum_facet(1,k,jj) = W2_cumSum_facet(1,k,j)
               W2_cumSum_facet(2,k,jj) = W2_cumSum_facet(2,k,j)
              ENDDO
              W2_cumSum_facet(2,jj,jj) = W2_cumSum_facet(2,jj,jj) - 1
              fCnt = fCnt + 1
              fIsSet(jj) = .TRUE.
              IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(5(A,I4))')
     &         ' CumSum SW.corner of facet #', jj,' set from facet',j,
     &         ' (pass,type=', npass,',', nType,')'
             ENDIF
C-    W <-- E ;
             IF ( i.EQ.W2_WEST .AND. ii.EQ.W2_EAST ) THEN
              DO k=1,nFacets
               W2_cumSum_facet(1,k,jj) = W2_cumSum_facet(1,k,j)
               W2_cumSum_facet(2,k,jj) = W2_cumSum_facet(2,k,j)
              ENDDO
              W2_cumSum_facet(1,jj,jj) = W2_cumSum_facet(1,jj,jj) - 1
              fCnt = fCnt + 1
              fIsSet(jj) = .TRUE.
              IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(5(A,I4))')
     &         ' CumSum SW.corner of facet #', jj,' set from facet',j,
     &         ' (pass,type=', npass,',', nType,')'
             ENDIF
            ENDIF
C--   3) with different orientation ( N <-> W or S <-> E )
C- Note: cannot rely on these connections for Cumul-Sum @ grid-cell center
            IF ( .NOT.fIsSet(jj) .AND. nType.EQ.3
     &                           .AND. fCnt.EQ.prev_fCnt ) THEN
C-    N <-- W or W <-- N ;
             IF (  ( i.EQ.W2_NORTH .AND. ii.EQ.W2_WEST  ) .OR.
     &             ( i.EQ.W2_WEST  .AND. ii.EQ.W2_NORTH ) ) THEN
              DO k=1,nFacets
               W2_cumSum_facet(1,k,jj) = W2_cumSum_facet(1,k,j)
               W2_cumSum_facet(2,k,jj) = W2_cumSum_facet(2,k,j)
              ENDDO
              W2_cumSum_facet(2,j, jj) = W2_cumSum_facet(2,j, jj) + 1
              W2_cumSum_facet(2,jj,jj) = W2_cumSum_facet(2,jj,jj) - 1
              fCnt = fCnt + 1
              fIsSet(jj) = .TRUE.
              IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(5(A,I4))')
     &         ' CumSum SW.corner of facet #', jj,' set from facet',j,
     &         ' (pass,type=', npass,',', nType,')'
             ENDIF
C-    E <-- S or S <-- E ;
             IF (  ( i.EQ.W2_EAST  .AND. ii.EQ.W2_SOUTH ) .OR.
     &             ( i.EQ.W2_SOUTH .AND. ii.EQ.W2_EAST  ) ) THEN
              DO k=1,nFacets
               W2_cumSum_facet(1,k,jj) = W2_cumSum_facet(1,k,j)
               W2_cumSum_facet(2,k,jj) = W2_cumSum_facet(2,k,j)
              ENDDO
              W2_cumSum_facet(1,j, jj) = W2_cumSum_facet(1,j, jj) + 1
              W2_cumSum_facet(1,jj,jj) = W2_cumSum_facet(1,jj,jj) - 1
              fCnt = fCnt + 1
              fIsSet(jj) = .TRUE.
              IF ( ABS(W2_printMsg).GE.2 ) WRITE(W2_oUnit,'(5(A,I4))')
     &         ' CumSum SW.corner of facet #', jj,' set from facet',j,
     &         ' (pass,type=', npass,',', nType,')'
             ENDIF
            ENDIF
C          end i loop
           ENDDO
          ENDIF
C        end facets loop
         ENDDO
         IF ( fCnt.GT.prev_fCnt ) THEN
          WRITE(msgBuf,'(2A,3(I4,A),I2,A)') 'W2_SET_MAP_CUMSUM: ',
     &     'set ', fCnt - prev_fCnt, ' /', nActiveFacets,
     &     ' active facets (pass,type=', npass, ',', nType, ')'
          CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
         ENDIF
        ENDIF
C      end nType loop
       ENDDO
C-    end do while (npass count)
      ENDDO
      IF ( fCnt.LT.nActiveFacets ) THEN
        WRITE(msgBuf,'(2A,2(I4,A))') 'W2_SET_MAP_CUMSUM: ',
     &  'Only get', fCnt, ' /', nActiveFacets,' active facets done'
        CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
        WRITE(msgBuf,'(2A)') '** WARNING ** W2_SET_MAP_CUMSUM: ',
     &    ' missing connections in Cumulated Sum'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
      ENDIF
      IF ( myProcId.EQ.0 ) THEN
        WRITE(W2_oUnit,'(2A,2(I4,A))')' Facet Matrix for CUMUL-SUM (',
     &         'nFacets=',nFacets, ', nActive=', nActiveFacets, ' ):'
        DO j=1,nFacets
          WRITE(W2_oUnit,'(A,I3,A,30(2I3,A))') '- facet', j, ' :',
     &         (W2_cumSum_facet(1,i,j),W2_cumSum_facet(2,i,j),' ,',
     &                                 i=1,nFacets)
        ENDDO
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   record "missing corner" tile

      IF ( useCubedSphereExchange ) THEN
       DO j=1,nFacets
        fNx = facet_dims(2*j-1)
        fNy = facet_dims( 2*j )
        IF ( fNx*fNy .GE. 1 ) THEN
         nbTx = fNx/sNx
         nbTy = fNy/sNy
         tN = facet_owns(1,j) - 1 + nbTx
         IF ( W2_tMC2.EQ.0 .AND. MOD(j,2).EQ.0
     &                     .AND. exch2_myFace(tN).NE.0 ) W2_tMC2 = tN
         tN = facet_owns(1,j) + (nbTy-1)*nbTx
         IF ( W2_tMC1.EQ.0 .AND. MOD(j,2).EQ.1
     &                     .AND. exch2_myFace(tN).NE.0 ) W2_tMC1 = tN
        ENDIF
       ENDDO
       IF ( myProcId.EQ.0 ) WRITE(W2_oUnit,'(3(A,I6))')
     &   ' missing-corner Tile for CUMUL-SUM (nTiles=', exch2_nTiles,
     &   ' ): W2_tMC1=', W2_tMC1, ' , W2_tMC2=', W2_tMC2
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--   Now Set cumul-sum Tile mapping

#ifdef W2_CUMSUM_USE_MATRIX
      DO j=1,W2_maxNbFacets
       DO i=1,W2_maxNbTiles
        facetXYSum(1,i,j) = 0
        facetXYSum(2,i,j) = 0
        facet_CSum(1,i,j) = 0
        facet_CSum(2,i,j) = 0
       ENDDO
      ENDDO

C-    First within each face:
      DO j=1,nFacets
       fNx = facet_dims(2*j-1)
       fNy = facet_dims( 2*j )
       IF ( fNx*fNy .GE. 1 ) THEN
        nbTx = fNx/sNx
        nbTy = fNy/sNy

        DO bj=1,nbTy
         DO bi=1,nbTx-1
           tS = facet_owns(1,j) - 1 + bi
           tN = tS + 1 + (bj-1)*nbTx
           DO k=facet_owns(1,j),tS
            W2_cumSum_tiles(1,k,tN) = 1
           ENDDO
         ENDDO
        ENDDO
        tN = facet_owns(1,j) - 1 + nbTx
        facetXYSum(1,tN,j) = 1
        DO k=facet_owns(1,j),tN-1
          facetXYSum(1,k,j) = W2_cumSum_tiles(1,k,tN)
        ENDDO

        DO bj=1,nbTy-1
         DO bi=1,nbTx
          tS = facet_owns(1,j) - 1 + bi
          tN = tS + bj*nbTx
          DO k=1,bj
           l = tS + (k-1)*nbTx
           W2_cumSum_tiles(2,l,tN) = 1
          ENDDO
         ENDDO
        ENDDO
        tN = facet_owns(1,j) + (nbTy-1)*nbTx
        facetXYSum(2,tN,j) = 1
        DO k=facet_owns(1,j),tN-1
          facetXYSum(2,k,j) = W2_cumSum_tiles(2,k,tN)
        ENDDO

       ENDIF
      ENDDO

C-    Then across facet:
      DO j=1,nFacets
       DO k=1,exch2_nTiles
        DO i=1,nFacets
        facet_CSum(1,k,j) = facet_CSum(1,k,j)
     &                    + W2_cumSum_facet(1,i,j)*facetXYSum(1,k,i)
        facet_CSum(2,k,j) = facet_CSum(2,k,j)
     &                    + W2_cumSum_facet(2,i,j)*facetXYSum(2,k,i)
        ENDDO
       ENDDO
      ENDDO

C-    Finally, account for cumulated sum at facet origin:
      DO j=1,nFacets
       DO tN=facet_owns(1,j),facet_owns(2,j)
        DO k=1,exch2_nTiles
         W2_cumSum_tiles(1,k,tN) = W2_cumSum_tiles(1,k,tN)
     &                           + facet_CSum(1,k,j)
         W2_cumSum_tiles(2,k,tN) = W2_cumSum_tiles(2,k,tN)
     &                           + facet_CSum(2,k,j)
        ENDDO
       ENDDO
      ENDDO

      IF (prtFlag) THEN
       WRITE(W2_oUnit,'(A,I6,A)')
     &    ' Tile Matrix for CUMUL-SUM (nTiles=', exch2_nTiles, ' ):'
       DO j=1,exch2_nTiles
        DO is=1,exch2_nTiles,10
         ie = MIN(is+9,exch2_nTiles)
         IF ( is.EQ.1 ) THEN
          WRITE(W2_oUnit,'(3(I6,A),10(2I3,A))') j,' ,',is,' ->',ie,' :',
     &     (W2_cumSum_tiles(1,i,j),W2_cumSum_tiles(2,i,j),' ,',i=is,ie)
         ELSE
          WRITE(W2_oUnit,'(8X,2(I6,A),10(2I3,A))') is,' ->',ie,' :',
     &     (W2_cumSum_tiles(1,i,j),W2_cumSum_tiles(2,i,j),' ,',i=is,ie)
         ENDIF
        ENDDO
       ENDDO
      ENDIF

      WRITE(msgBuf,'(2A)') 'W2_SET_MAP_CUMSUM: ',
     &  'setting Tile Matrix for CUMUL-SUM : done'
      CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )

#else /* W2_CUMSUM_USE_MATRIX */
      WRITE(msgBuf,'(2A)') 'W2_SET_MAP_CUMSUM: ',
     &  'done (skip Tile Matrix setting)'
      CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
#endif /* W2_CUMSUM_USE_MATRIX */

      RETURN
      END