C $Header: /u/gcmpack/MITgcm/pkg/exch2/w2_set_map_tiles.F,v 1.6 2012/07/16 20:25:10 jmc Exp $
C $Name:  $

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

C--  File w2_set_map_tiles.F:
C--   Contents
C--   o W2_SET_MAP_TILES :: Set tiles and IO mapping
C--   o FIND_GCD_N       :: Returns the Greatest Common Divisor

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

C !INTERFACE:
      SUBROUTINE W2_SET_MAP_TILES( myThid )

C     !DESCRIPTION:
C     Set-up tiles mapping and IO global mapping

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     !FUNCTIONS:
      INTEGER  FIND_GCD_N
      EXTERNAL 

C     !LOCAL VARIABLES:
C     === Local variables ===
C     msgBuf     :: Informational/error message buffer
      CHARACTER*(MAX_LEN_MBUF) msgBuf
      INTEGER tNx, tNy, fNx, fNy, nbPts, fBaseX, fBaseY
      INTEGER nbTx, nbTy
      INTEGER j, ii, k, tId, tx, ty
      INTEGER divide, nnx(W2_maxNbFacets)
      INTEGER errCnt, tCnt
      LOGICAL tileIsActive, prtFlag
CEOP

C     Set-up tiles mapping and IO global mapping
      WRITE(msgBuf,'(2A)') 'W2_SET_MAP_TILES:',
     &       ' tile mapping within facet and global Map:'
      CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
      prtFlag = ABS(W2_printMsg).GE.2
     &       .OR. ( W2_printMsg .NE.0 .AND. myProcId.EQ.0 )

      tNx = sNx
      tNy = sNy
C--   Check that tile dims divide facet dims
      errCnt = 0
      tCnt = 0
      nbPts = 0
      DO j=1,nFacets
        fNx = facet_dims(2*j-1)
        fNy = facet_dims( 2*j )
        nbTx = fNx/tNx
        nbTy = fNy/tNy
        IF ( nbTx*tNx .NE. fNx ) THEN
          WRITE(msgBuf,'(A,I3,2(A,I7))') 'Facet',j,
     &      ' : X-size=', fNx, ' not multiple of sNx=', tNx
          CALL PRINT_ERROR( msgBuf, myThid )
          errCnt = errCnt + 1
        ENDIF
        IF ( nbTy*tNy .NE. fNy ) THEN
          WRITE(msgBuf,'(A,I3,2(A,I7))') 'Facet',j,
     &      ' : Y-size=', fNy, ' not multiple of sNy=', tNy
          CALL PRINT_ERROR( msgBuf, myThid )
          errCnt = errCnt + 1
        ENDIF
        facet_owns(1,j) = tCnt+1
        tCnt = tCnt + nbTx*nbTy
        facet_owns(2,j) = tCnt
        nbPts = nbPts + fNx*fNy
      ENDDO
      IF ( errCnt.GT.0 ) THEN
        WRITE(msgBuf,'(A,I3,A)')
     &   ' W2_SET_MAP_TILES: found', errCnt, ' Fatal errors'
        CALL PRINT_ERROR( msgBuf, myThid )
        STOP 'ABNORMAL END: S/R W2_SET_MAP_TILES'
      ENDIF
C--   Check that domain size and (SIZE.h + blankList) match:
      IF ( tCnt.NE.exch2_nTiles ) THEN
        WRITE(msgBuf,'(A,I6,A)')
     &   'W2_SET_MAP_TILES: Domain Total # of tiles =', tCnt, ' does'
        CALL PRINT_ERROR( msgBuf, myThid )
        WRITE(msgBuf,'(A,I6)')
     &   'W2_SET_MAP_TILES:  not match (SIZE.h+blankList)=',exch2_nTiles
        CALL PRINT_ERROR( msgBuf, myThid )
        STOP 'ABNORMAL END: S/R W2_SET_MAP_TILES'
      ENDIF

      IF ( W2_mapIO.EQ.1 ) THEN
C--   Compact IO map (mostly in Y dir): search for Greatest Common Divisor
C     of all x-size (faster to apply GCD to Nb of Tiles in X):
        k = 0
        nnx(1) = 0
        DO j=1,nFacets
C     skip empty facet
          IF ( facet_dims(2*j-1).GT.0 ) THEN
            k = k + 1
            nnx(k) = facet_dims(2*j-1)/tNx
          ENDIF
        ENDDO
        divide = FIND_GCD_N( nnx, k )
        W2_mapIO = divide*tNx
        WRITE(msgBuf,'(A,2(I5,A))') ' W2_mapIO =', W2_mapIO,
     &                              ' (=', divide, '*sNx)'
        CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
      ENDIF

C--   Global Map size:
C     facets stacked in x direction
      exch2_xStack_Nx = 0
      exch2_xStack_Ny = 0
      DO j=1,nFacets
        exch2_xStack_Nx =   exch2_xStack_Nx + facet_dims(2*j-1)
        exch2_xStack_Ny = MAX( exch2_xStack_Ny, facet_dims(2*j) )
      ENDDO
C     facets stacked in y direction
      exch2_yStack_Nx = 0
      exch2_yStack_Ny = 0
      DO j=1,nFacets
        exch2_yStack_Nx = MAX( exch2_yStack_Nx, facet_dims(2*j-1) )
        exch2_yStack_Ny =   exch2_yStack_Ny + facet_dims(2*j)
      ENDDO
      IF ( W2_mapIO.EQ.-1 ) THEN
        exch2_global_Nx = exch2_xStack_Nx
        exch2_global_Ny = exch2_xStack_Ny
      ELSEIF ( W2_mapIO.EQ.0 ) THEN
        exch2_global_Nx = nbPts
        exch2_global_Ny = 1
      ELSE
        exch2_global_Nx = W2_mapIO
        exch2_global_Ny = nbPts/W2_mapIO
      ENDIF
      WRITE(msgBuf,'(A,2(A,I8))') ' Global Map (IO):',
     &  ' X-size=', exch2_global_Nx, ' , Y-size=', exch2_global_Ny
      CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )

C--   Set tiles mapping within facet (sub-domain) and within Global Map
      WRITE(msgBuf,'(2A)') 'W2_SET_MAP_TILES:',
     &       ' tile offset within facet and global Map:'
      CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
      tId = 0
      nbPts = 0
      fBaseX = 0
      fBaseY = 0
      DO j=1,nFacets
        fNx = facet_dims(2*j-1)
        fNy = facet_dims( 2*j )
        nbTx = fNx/tNx
        nbTy = fNy/tNy
        WRITE(W2_oUnit,'(A,I3,2(A,I6),A,I5,2(A,I4),A)')
     &    '- facet', j, ' : X-size=', fNx, ' , Y-size=', fNy,
     &    ' ;', nbTx*nbTy, ' tiles (Tx,Ty=', nbTx,',',nbTy,')'
c       CALL PRINT_MESSAGE( msgBuf, W2_oUnit, SQUEEZE_RIGHT, myThid )
        DO ty=1,nbTy
         DO tx=1,nbTx
          tId = tId + 1
C--   Tags blank tile by removing facet # (exch2_myFace) but keeps its location
          tileIsActive = .TRUE.
          DO k=1,nBlankTiles
           IF ( blankList(k).EQ.tId ) tileIsActive = .FALSE.
          ENDDO
          IF ( tileIsActive ) exch2_myFace(tId) = j
          exch2_mydNx ( tId ) = fNx
          exch2_mydNy ( tId ) = fNy
          exch2_tNx   ( tId ) = tNx
          exch2_tNy   ( tId ) = tNy
          exch2_tBasex( tId ) = (tx-1)*tNx
          exch2_tBasey( tId ) = (ty-1)*tNy
C--   Global IO Mappings
C       these are for OBCS (vertical slices)
          exch2_txXStackLo( tId ) = 1 + exch2_tBasex(tId) + fBaseX
          exch2_tyXStackLo( tId ) = 1 + exch2_tBasey(tId)
          exch2_txYStackLo( tId ) = 1 + exch2_tBasex(tId)
          exch2_tyYStackLo( tId ) = 1 + exch2_tBasey(tId) + fBaseY
C       and these for global files (3d files/horizontal 2d files)
          IF ( W2_mapIO.EQ.-1 ) THEN
C-        Old format
            exch2_txGlobalo( tId ) = 1 + exch2_tBasex(tId) + fBaseX
            exch2_tyGlobalo( tId ) = 1 + exch2_tBasey(tId)
          ELSEIF ( W2_mapIO.EQ.0 ) THEN
C-        Compact format = 1 long line
            ii = nbPts + exch2_tBasex(tId) + exch2_tBasey(tId)*fNx
            exch2_txGlobalo( tId ) = 1 + ii
            exch2_tyGlobalo( tId ) = 1
          ELSE
C         Compact format: piled in the Y direction
            ii = nbPts + exch2_tBasex(tId) + exch2_tBasey(tId)*fNx
            exch2_txGlobalo( tId ) = 1 + MOD(ii,W2_mapIO)
            exch2_tyGlobalo( tId ) = 1 + ii/W2_mapIO
          ENDIF
          IF ( prtFlag )
     &    WRITE(W2_oUnit,'(A,I5,3(A,I3),2A,2I5,2A,2I8)') '  tile',tId,
     &    ' on facet', exch2_myFace(tId),' (',tx,',',ty,'):',
     &         ' offset=', exch2_tBasex(tId), exch2_tBasey(tId),' ;',
     &    ' on Glob.Map=', exch2_txGlobalo(tId),exch2_tyGlobalo(tId)
         ENDDO
        ENDDO
        fBaseX = fBaseX + fNx
        fBaseY = fBaseY + fNy
        nbPts = nbPts + fNx*fNy
      ENDDO

      RETURN
      END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: FIND_GCD_N C !INTERFACE: INTEGER FUNCTION FIND_GCD_N( fldList, nFld ) C !DESCRIPTION: C *==========================================================* C | FUNCTION FIND_GCD_N C | o Find the Greatest Common Divisor of N integers C *==========================================================* C !USES: IMPLICIT NONE C !INPUT PARAMETERS: C fldList :: list of integers to search for GCD C nFLd :: length of the input integer list. INTEGER nFLd INTEGER fldList(nFld) C !LOCAL VARIABLES: INTEGER mnFld, divide INTEGER j, ii LOGICAL flag LOGICAL localDBg CEOP PARAMETER ( localDBg = .FALSE. ) c PARAMETER ( localDBg = .TRUE. ) mnFld = fldList(1) DO j=1,nFld mnFld = MIN( mnFld, fldList(j) ) ENDDO IF (localDBg) WRITE(0,'(A,I8)') 'FIND_GCD_N: mnFld=',mnFld IF (mnFld.GT.1 ) THEN divide = 1 ii = 2 DO WHILE ( ii.LE.mnFld ) IF (localDBg) WRITE(0,'(A,I8)') ' GCD : try',ii flag = .TRUE. DO j=1,nFld flag = flag.AND.(MOD(fldList(j),ii).EQ.0 ) ENDDO IF ( flag ) THEN divide = divide*ii DO j=1,nFld fldList(j) = fldList(j)/ii ENDDO IF (localDBg) WRITE(0,'(A,I8)') & 'FIND_GCD_N: com.fact=',ii mnFld = mnFld/ii ELSE ii = ii+2 IF (ii.EQ.4) ii=3 ENDIF ENDDO C- Put back the original Nb: IF (localDBg) WRITE(0,'(10I8)') (fldList(j),j=1,nFld) DO j=1,nFld fldList(j) = fldList(j)*divide ENDDO ELSE divide = MAX( 0, mnFld ) ENDIF FIND_GCD_N = divide RETURN END