ecearth3/sources/tm5mp/base/tm5_partools.F90

!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MPI(action) if (ierr/=MPI_SUCCESS) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      PARTOOLS
!
! !DESCRIPTION: MPI general tools. This module:
!
! - holds MPI constants, specifically defining some if we are not running with MPI
! 
! - initializes, finalizes and aborts MPI (wrappers, so it works if not runnning MPI):
!
!      TM5_MPI_Init, TM5_MPI_Init2 (init row/col communicators)
!      TM5_MPI_Done
!      TM5_MPI_Abort
!
! - provides wrappers around some of MPI broadcast/barrier/reduce calls :
! 
!      PAR_BARRIER          : 
!      PAR_STOPMPI          : 
!      PAR_BROADCAST_STATUS : 
!      PAR_BROADCAST        : 
!      PAR_REDUCE           : reduction is done across the array. Result size = 1
!      PAR_REDUCE_ELEMENT   : reduction is done for each element. Result size = size of original array
!
! These are general tools, that deal with MPI, regardless of the domain decomposition used.
!
! - Some of these call can be limited to the 'row' or 'col' communicators, but
!   this is partially implemented (although trivial to extent if needed)
! 
!\\
!\\
! !INTERFACE: 
!
MODULE PARTOOLS
  !
  ! !USES:
  !
  use GO,   only : gol, goPr, goErr
  use dims, only : nregions_all

#ifdef MPI
  use mpi
#endif   
  IMPLICIT NONE
  !
  ! !PUBLIC DATA MEMBERS:
  !
  integer, parameter :: PAR_OPER_SUM = 100   ! mpi reduce flag 
  integer            :: localComm  ! global communicator (equal to MPI_COMM_WORLD if not coupled model)
  integer            :: row_comm   ! row    communicator
  integer            :: col_comm   ! column communicator
  integer            :: my_real    ! platform dependent reference to real values for MPI 
  integer            :: myid       ! rank in localComm
  integer            :: npes       ! total number of PE's
  integer            :: npe_lat    ! #pes in lat direction
  integer            :: npe_lon    ! #pes in long direction
  integer            :: ierr       ! return status of MPI routine calls
  integer            :: root       ! myid of root in localComm
  integer            :: RowRootId  ! myid of row_root of my row (in localComm)
#ifndef MPI
  integer, parameter :: MPI_SUCCESS   = 0
  integer, parameter :: MPI_INFO_NULL = 0
  integer, parameter :: MPI_CHARACTER = 0
  integer, parameter :: MPI_INTEGER   = 0
  integer, parameter :: MPI_PROC_NULL = -999
#endif
  logical            :: isRoot, isRowRoot     ! is process global root, row root?
  character(len=6)   :: procname   ! character keys for each processor
  !
  ! !PRIVATE DATA MEMBERS:
  !
  integer,                     private :: irow, jcol       ! rank in column and row communicators
  character(len=*), parameter, private :: mname='ParTools'
  !
  ! !INTERFACE:
  !
  INTERFACE Par_Broadcast
     MODULE PROCEDURE Par_Broadcast_i0
     MODULE PROCEDURE Par_Broadcast_i1
     MODULE PROCEDURE Par_Broadcast_s
     MODULE PROCEDURE Par_Broadcast_l0
     MODULE PROCEDURE Par_Broadcast_l1
     MODULE PROCEDURE Par_Broadcast_r0
     MODULE PROCEDURE Par_Broadcast_r1
     MODULE PROCEDURE Par_Broadcast_r2
     MODULE PROCEDURE Par_Broadcast_r3
  END INTERFACE

  INTERFACE Par_Reduce
     MODULE PROCEDURE Par_Reduce_i0
     MODULE PROCEDURE Par_Reduce_r0
     MODULE PROCEDURE Par_Reduce_r1    ! r1, r2, ... are useful only to get results that do not depend
!     MODULE PROCEDURE Par_Reduce_r2   ! on the number of processors (when summing). It is cheaper to
     MODULE PROCEDURE Par_Reduce_r3   ! do the operation on each proc, and then reduce a single value.
  END INTERFACE

  INTERFACE Par_Reduce_Element
     MODULE PROCEDURE Par_Reduce_element_r1
     MODULE PROCEDURE Par_Reduce_element_r2
     MODULE PROCEDURE Par_Reduce_element_r3
  END INTERFACE
  !
  ! !REVISION HISTORY: 
  !   18 Jan 2012 - P. Le Sager - revamped for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    TM5_MPI_Init
  !
  ! !DESCRIPTION: initializes MPI
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE TM5_MPI_Init( status, comm )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)            ::  status
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in), optional   ::  comm
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/TM5_MPI_Init'

#ifdef MPI

    ! communicator provided, for example by prism coupler ?
    if ( present(comm) ) then

       localComm = comm

    else

       ! init mpi here to set MPI_COMM_WORLD etc
       call MPI_INIT( ierr )
       IF_NOTOK_MPI(status=1)

       localComm = MPI_COMM_WORLD

    end if

    ! obtain number of proceses:
    call MPI_COMM_SIZE( localComm, npes, ierr )
    IF_NOTOK_MPI(status=1)

    ! obtain proces number:
    call MPI_COMM_RANK( localComm, myid, ierr )
    IF_NOTOK_MPI(status=1)

    ! set root in localComm to PE 0; 'real' default to double precision
    root = 0
    my_real = MPI_DOUBLE_PRECISION 

#else

    localComm = 0   ! dummy communicator
    npes = 1        ! single processor
    root = 0        ! id for root processor    
    myid = root     ! only one processor, so this is always root

#endif

    ! set processor names: pe0000, pe0001, ...
    write (procname,'("pe",i4.4)') myid

    isRoot = ( myid == root )

    ! ok
    status = 0

  END SUBROUTINE TM5_MPI_Init
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    TM5_MPI_INIT2
  !
  ! !DESCRIPTION: Step #2 of initialization, which requires input from rcfile.
  !               It defines the extra row and column communicators. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE TM5_MPI_Init2( nlon, nlat, status )
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   ::  nlat, nlon  ! # PE on each direction
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REVISION HISTORY: 
    !   21 Feb 2012 - P. Le Sager - 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/TM5_MPI_Init2'
    integer :: irow, jcol
    
    npe_lat=nlat
    npe_lon=nlon

#ifdef MPI

    irow = myid/npe_lon
    jcol = mod(myid,npe_lon)

    ! NOTE: could use myid as key for rank designation if that facilitates some comm
    call MPI_COMM_SPLIT(localComm, irow, jcol, row_comm, ierr)
    IF_NOTOK_MPI(status=1)
    isRowRoot = jcol==0
    RowRootId = (myid/npe_lon)*npe_lon
    
    ! define column communicator [not needed yet]
    call MPI_COMM_SPLIT(localComm, jcol, irow, col_comm, ierr)
    IF_NOTOK_MPI(status=1)
    
#else
    row_comm  = 0   ! dummy communicator
    col_comm  = 0   ! dummy communicator
#endif

    status = 0

  END SUBROUTINE TM5_MPI_Init2
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    TM5_MPI_Done
  !
  ! !DESCRIPTION: finalizes MPI
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE TM5_MPI_Done( status, comm )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)            ::  status
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in), optional   ::  comm
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/TM5_MPI_Done'


#ifdef MPI

    ! shut-down communication, only if communicator is not provided
    if (.not. present(comm) ) then

       !write (*,'("call MPI_Finalize from proces ",i6)') myid
       call MPI_Finalize( ierr )
       IF_NOTOK_MPI(status=1)

    end if

#endif

    ! ok
    status = 0

  END SUBROUTINE TM5_MPI_Done
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  TM5_MPI_Abort
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine TM5_MPI_Abort( errorcode, status )
    !
    ! !USES:
    !
    use GO, only : goExit
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)             ::  errorcode
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)            ::  status
    !
    ! !REVISION HISTORY: 
    !   18 Jan 2012 - P. Le Sager -
    !
    ! !REMARKS:
    ! 
    ! (pls, 8-4-2011) Sometimes the code does not return from MPI_Abort, for
    ! example when a problem reading restart files occurs. From the doc:
    !        -------------------------------------------------
    ! "Before the error value is returned, the current MPI error handler is
    ! called. By default, this error handler aborts the MPI job, except for
    ! I/O function errors."
    !        -------------------------------------------------
    ! so, the only way to nicely abort is to close files when an i/o error
    ! occurs. Done with a new macro for problems with reading/writing restart
    ! (see tm5_restart.F90). Check if your module/routine is prone to
    ! i/o error, and apply a similar patch.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/TM5_MPI_Abort'

#ifdef MPI
    ! emergency break ...
    call MPI_Abort( localComm, errorcode, ierr )
    IF_NOTOK_MPI(status=1)

#else
    ! system exit:
    call goExit( errorcode )
#endif

    ! ok
    status = 0

  end subroutine TM5_MPI_Abort
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Barrier
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_BARRIER( ROW, COL)
    !
    ! !INPUT PARAMETERS:
    !
    logical, intent(in), optional :: row, col  ! to limit call to sub-communicator
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/PAR_BARRIER'

#ifdef MPI

    integer :: l_comm

    l_comm=localComm

    if(present(row)) then
       if(row) l_comm=row_comm
    end if

    if(present(col)) then
       if(col) l_comm=col_comm
    end if

    call mpi_barrier( l_comm, ierr )
    IF_NOTOK_MPI(write (gol,*)"error MPI_BARRIER in PAR_BARRIER";call goErr)
#endif

  END SUBROUTINE PAR_BARRIER
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_StopMPI
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_StopMPI
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter ::  rname = mname//'/Par_StopMPI'

    write (*,'("WARNING - (par_)stopmpi should be avoided; please trace back to main program ...")')

#ifdef MPI
    ! shut down mpi communication:
    call mpi_finalize(ierr)
    IF_NOTOK_MPI(write (gol,*)"error MPI_FINALIZE in Par_StopMPI";call goErr)
#endif

    ! fortran stop ....
    stop 'Fortran STOP in Par_StopMPI ...'

  end subroutine Par_StopMPI
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    Par_Broadcast_Status
  !
  ! !DESCRIPTION: Broadcast integer istat from id PE to all. Same as
  !               Par_Broadcast_i(istat, id, istat)
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_Status( istat, id )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)            ::  istat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)               ::  id
    !
    ! !REVISION HISTORY: 
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_Status'
    integer                      ::  status

#ifdef MPI
    ! send the input status to all other processes:
    call Par_Broadcast( istat, status, id=id )
    IF_NOTOK_RETURN(istat=1)
#endif
    
  end subroutine Par_Broadcast_Status
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    Par_Broadcast_i0
  !
  ! !DESCRIPTION: broadcast one integer (scalar case)
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_i0( i, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)        :: i
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_i'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call MPI_BCast( i, 1, MPI_INTEGER, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)
#endif

    status = 0

  end subroutine Par_Broadcast_i0
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    Par_Broadcast_i1
  !
  ! !DESCRIPTION: Broadcast 1-D array of integers.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_i1( i, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)        :: i(:)
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_i'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call MPI_BCast( i, size(i), MPI_INTEGER, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)
#endif

    status = 0

  end subroutine Par_Broadcast_i1
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_s
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_s( s, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    character(len=*), intent(inout)   :: s
    integer, intent(out)              :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_s'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call MPI_BCast( s, len(s), MPI_CHARACTER, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)
#endif

    status = 0
    
  end subroutine Par_Broadcast_s
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_l0
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_l0( x, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    logical, intent(inout)        :: x
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_l0'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call mpi_bcast( x, 1, MPI_LOGICAL, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)

#endif

    status = 0

  end subroutine Par_Broadcast_l0
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_l1
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_l1( x, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    logical, intent(inout)        :: x(:)
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_l1'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call mpi_bcast( x, size(x), MPI_LOGICAL, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)

#endif

    status = 0

  end subroutine Par_Broadcast_l1
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_r0
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_r0( x, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real,    intent(inout)        :: x
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_r0'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call mpi_bcast( x, 1, my_real, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)

#endif

    status = 0

  end subroutine Par_Broadcast_r0
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_r1
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_r1( x, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real,    intent(inout)        :: x(:)
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_r1'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call mpi_bcast( x, size(x), my_real, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)

#endif

    status = 0

  end subroutine Par_Broadcast_r1
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_r2
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_r2( x, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real,    intent(inout)        :: x(:,:)
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_r2'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call mpi_bcast( x, size(x), my_real, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)
#endif

    status = 0

  end subroutine Par_Broadcast_r2
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  Par_Broadcast_r3
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Par_Broadcast_r3( x, status, ID, ROW, COL )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real,    intent(inout)        :: x(:,:,:)
    integer, intent(out)          :: status

    integer, intent(in), optional :: id    ! broadcaster ID         (default is communicator root)
    logical, intent(in), optional :: row   ! limit to PE on row_comm (default is global localComm)
    logical, intent(in), optional :: col   ! limit to PE on col_comm (default is global localComm)
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2012 - P. Le Sager - added ROW, COL options. Made ID optional.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Par_Broadcast_r3'
    logical :: l_row, l_col
    integer :: l_root, l_comm, l_id

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_col=.false.  ; if(present(col)) l_col=col
    
    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
    else
       l_root = root
       l_comm = localComm
    end if

    l_id = l_root
    if(present(id)) l_id=id

    call mpi_bcast( x, size(x), my_real, l_id, l_comm, ierr )
    IF_NOTOK_MPI(status=1)
    
#endif

    status = 0

  end subroutine Par_Broadcast_r3
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_I0
  !
  ! !DESCRIPTION: Wrapper around MPI_REDUCE or MPI_ALLREDUCE.
  !               
  !               Apply to a SINGLE INTEGER.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_I0( DATA, OP, RESULTAT, STATUS, ALL, ROW )
    !
    ! !INPUT PARAMETERS:
    !
    integer,           intent(in) :: data
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! use MPI_ALLREDUCE instead of MPI_REDUCE
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,    intent(out) :: resultat
    integer,    intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   19 Jun 2013 - Ph. Le Sager - v0 (copy from par_reduce_r0)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_I0'
    logical :: l_row, l_all
    integer :: l_root, l_comm, l_id, l_op

#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_all=.false.  ; if(present(all)) l_all=all

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
       l_id   = jcol
    else
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if


    ! degenerate cases first
    if(l_row.and.(npe_lon==1)) then
       resultat = data
    else

       SELECT CASE( OP )
       case('sum', 'Sum', 'SUM')
          l_op = MPI_SUM
       case('max', 'Max', 'MAX')
          l_op = MPI_MAX
       case('min', 'Min', 'MIN')
          l_op = MPI_MIN
       case default
          write(gol,*) 'UNSUPPORTED OPERATION :', OP; status=1
          IF_NOTOK_RETURN(status=1)
       END SELECT

       if (l_all) then
          call MPI_ALLREDUCE(data, resultat, 1, MPI_INTEGER, l_op, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       else
          call MPI_REDUCE(data, resultat, 1, MPI_INTEGER, l_op, l_root, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       end if

    end if

#else
    resultat = data
#endif

    status=0

  END SUBROUTINE PAR_REDUCE_I0
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_R0
  !
  ! !DESCRIPTION: Wrapper around MPI_REDUCE or MPI_ALLREDUCE.
  !               
  !               Apply to a SINGLE REAL.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_R0( DATA, OP, RESULTAT, STATUS, ALL, ROW, COL )
    !
    ! !INPUT PARAMETERS:
    !
    real,              intent(in) :: data
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! use MPI_ALLREDUCE instead of MPI_REDUCE
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    logical, intent(in), optional :: col   ! limit to PE on col_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: resultat
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   01 Nov 2011 - P. Le Sager - v0
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_R0'
    logical :: l_row, l_all, l_col
    integer :: l_root, l_comm, l_id, l_op
    
#ifdef MPI

    l_col=.false.  ; if(present(col)) l_col=col
    l_row=.false.  ; if(present(row)) l_row=row
    l_all=.false.  ; if(present(all)) l_all=all

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
       l_id   = jcol
    else if(l_col) then
       l_root = 0               ! by our own design
       l_comm = col_comm
       l_id   = irow
    else
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if


    ! degenerate cases first
    if(l_row.and.(npe_lon==1)) then
       resultat = data
    else

       SELECT CASE( OP )
         case('sum', 'Sum', 'SUM')
            l_op = MPI_SUM
         case('max', 'Max', 'MAX')
            l_op = MPI_MAX
         case('min', 'Min', 'MIN')
            l_op = MPI_MIN
         case default
            write(gol,*) 'UNSUPPORTED OPERATION :', OP; status=1
            IF_NOTOK_RETURN(status=1)
       END SELECT

       if (l_all) then
          call MPI_ALLREDUCE(data, resultat, 1, my_real, l_op, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       else
          call MPI_REDUCE(data, resultat, 1, my_real, l_op, l_root, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       end if

    end if
    
#else
    resultat = data
#endif

    status=0

  END SUBROUTINE PAR_REDUCE_R0
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_R1
  !
  ! !DESCRIPTION: Global reduction. Data are gathered on root, where OP is
  !                then done. Apply to a 1D REAL.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_R1( DATA, OP, RESULTAT, STATUS, ALL, ROW )
    !
    ! !INPUT PARAMETERS:
    !
    real,              intent(in) :: data(:)
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! mimic mpi_allreduce instead of mpi_reduce
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: resultat
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   01 Nov 2011 - P. Le Sager - FIXME: NOT TESTED
    !
    ! !REMARKS:
    !
    ! (1) this is a convenient tool to get same results independently of the
    !      number of processors. If this is not what you want, you can do OP
    !      on each processor before calling par_reduce_r0 : it is less expensive.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_R1'

#ifdef MPI
    real, allocatable :: glb(:)
    logical :: l_row
    integer :: sz, l_root, l_comm, l_id
    
    sz = size(data)

    l_row=.false.
    if(present(row)) l_row=row

    if(l_row)then
       allocate( glb(npe_lon*sz) )
       l_root = 0
       l_comm = row_comm
       l_id   = jcol
    else
       allocate( glb(npes*sz) )
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if


    CALL MPI_GATHER( data, sz, my_real, glb, sz, my_real, l_root, l_comm, ierr)
    IF_NOTOK_MPI(status=1)

    if (l_root==l_id) then

       select case( op )

       case('sum', 'Sum', 'SUM')
          resultat = sum(glb)          
       case('max', 'Max', 'MAX')
          resultat = maxval(glb)          
       case('min', 'Min', 'MIN')
          resultat = minval(glb)          
       case default
          write(gol,*) 'UNSUPPORTED OPERATION'; status=1
          IF_NOTOK_RETURN(status=1)

       end select
    end if

    deallocate(glb)
    
    if (present(all)) then
       if (all) call MPI_bcast(resultat, 1, my_real, l_root, l_comm, ierr)
    end if

#else

    select case( op )

    case('sum', 'Sum', 'SUM')
       resultat = sum(data)
    case('max', 'Max', 'MAX')
       resultat = maxval(data)
    case('min', 'Min', 'MIN')
       resultat = minval(data)
    case default
       write(gol,*) 'UNSUPPORTED OPERATION'; status=1
       IF_NOTOK_RETURN(status=1)

    end select
    
#endif

    status=0

  END SUBROUTINE PAR_REDUCE_R1
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_R2
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
!   subroutine Par_Reduce_r2( send, recv, oper, id_recv, status )
!    
  ! !REMARKS:
  
  ! in TM5 v3, this routine was a wrapper around MPI_REDUCE( ..., MPI_SUM, ...), and was called
  ! in project CLIMAQS/user_output_pdump.F90 (an updated retro output) and in
  ! various user_output_retro.F90, to sum budget... 
  !         =====> in tm5-mp, this routine became PAR_REDUCE_ELEMENT_R2 ()



  ! TO DEVELOP IF YOU NEED RESULT THAT ARE INDEPENDENT OF THE NUMBER OF
  ! PROCESSOR YOU ARE USING.


  
!   end subroutine Par_Reduce_r2
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_R3
  !
  ! !DESCRIPTION: Global reduction. Data are gathered on root, where OP is
  !                then done. Apply to a 3D REAL.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_R3( DATA, OP, RESULTAT, STATUS, ALL, ROW )
    !
    ! !INPUT PARAMETERS:
    !
    real,              intent(in) :: data(:,:,:)
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! mimic mpi_allreduce instead of mpi_reduce
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: resultat
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !
    ! !REMARKS:
    !
    ! (1) this is a convenient tool to get same results independently of the
    !      number of processors. If this is not what you want, you can do OP
    !      on each processor before calling par_reduce_r0 : it is less expensive.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_R3'

#ifdef MPI
    real, allocatable :: glb(:)
    logical :: l_row
    integer :: sz, l_root, l_comm, l_id

    sz = size(data)

    l_row=.false.
    if(present(row)) l_row=row

    if(l_row)then
       allocate( glb(npe_lon*sz) )
       l_root = 0
       l_comm = row_comm
       l_id   = jcol
    else
       allocate( glb(npes*sz) )
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if


    CALL MPI_GATHER( data, sz, my_real, glb, sz, my_real, l_root, l_comm, ierr)
    IF_NOTOK_MPI(status=1)

    if (l_root==l_id) then

       select case( op )

       case('sum', 'Sum', 'SUM')
          resultat = sum(glb)          
       case('max', 'Max', 'MAX')
          resultat = maxval(glb)          
       case('min', 'Min', 'MIN')
          resultat = minval(glb)          
       case default
          write(gol,*) 'UNSUPPORTED OPERATION'; status=1
          IF_NOTOK_RETURN(status=1)

       end select
    end if

    deallocate(glb)

    if (present(all)) then
       if (all) call MPI_bcast(resultat, 1, my_real, l_root, l_comm, ierr)
    end if

#else

    select case( op )

    case('sum', 'Sum', 'SUM')
       resultat = sum(data)
    case('max', 'Max', 'MAX')
       resultat = maxval(data)
    case('min', 'Min', 'MIN')
       resultat = minval(data)
    case default
       write(gol,*) 'UNSUPPORTED OPERATION'; status=1
       IF_NOTOK_RETURN(status=1)

    end select

#endif

    status=0

  END SUBROUTINE PAR_REDUCE_R3
  !EOC

  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_ELEMENT_R1
  !
  ! !DESCRIPTION: reduce 1D arrays element-wise
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_ELEMENT_R1 (DATA, OP, RESULTAT, STATUS, ALL, ROW)
    !
    ! !USES:
    !
    use dims, only : CheckShape
    !
    ! !INPUT PARAMETERS:
    !
    real,              intent(in) :: data(:)
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! mimic mpi_allreduce instead of mpi_reduce
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: resultat(:)
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    1 Mar 2012 - P. Le Sager - v0
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_ELEMENT_R2'
    logical :: l_row, l_all
    integer :: l_root, l_comm, l_id, l_op


    call CheckShape( shape(data), shape(resultat), status )
    IF_NOTOK_RETURN(status=1)


#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_all=.false.  ; if(present(all)) l_all=all

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
       l_id   = jcol
    else
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if

    ! degenerate cases first
    if(l_row.and.(npe_lon==1)) then
       resultat = data
    else

       SELECT CASE( OP )
       case('sum', 'Sum', 'SUM')
          l_op = MPI_SUM
       case('max', 'Max', 'MAX')
          l_op = MPI_MAX
       case('min', 'Min', 'MIN')
          l_op = MPI_MIN
       case default
          write(gol,*) 'UNSUPPORTED OPERATION :', OP; status=1
          IF_NOTOK_RETURN(status=1)
       END SELECT

       if (l_all) then
          call MPI_ALLREDUCE(data, resultat, size(data), my_real, l_op, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       else
          call MPI_REDUCE(data, resultat, size(data), my_real, l_op, l_root, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       end if

    end if

#else
    resultat = data
#endif

    status=0

  END SUBROUTINE PAR_REDUCE_ELEMENT_R1
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_ELEMENT_R2
  !
  ! !DESCRIPTION: reduce 2D arrays element-wise
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_ELEMENT_R2 (DATA, OP, RESULTAT, STATUS, ALL, ROW)
    !
    ! !USES:
    !
    use dims, only : CheckShape
    !
    ! !INPUT PARAMETERS:
    !
    real,              intent(in) :: data(:,:)
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! mimic mpi_allreduce instead of mpi_reduce
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: resultat(:,:)
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    1 Mar 2012 - P. Le Sager - v0
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_ELEMENT_R2'
    logical :: l_row, l_all
    integer :: l_root, l_comm, l_id, l_op


    call CheckShape( shape(data), shape(resultat), status )
    IF_NOTOK_RETURN(status=1)


#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_all=.false.  ; if(present(all)) l_all=all

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
       l_id   = jcol
    else
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if

    ! degenerate cases first
    if(l_row.and.(npe_lon==1)) then
       resultat = data
    else

       SELECT CASE( OP )
       case('sum', 'Sum', 'SUM')
          l_op = MPI_SUM
       case('max', 'Max', 'MAX')
          l_op = MPI_MAX
       case('min', 'Min', 'MIN')
          l_op = MPI_MIN
       case default
          write(gol,*) 'UNSUPPORTED OPERATION :', OP; status=1
          IF_NOTOK_RETURN(status=1)
       END SELECT

       if (l_all) then
          call MPI_ALLREDUCE(data, resultat, size(data), my_real, l_op, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       else
          call MPI_REDUCE(data, resultat, size(data), my_real, l_op, l_root, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       end if

    end if

#else
    resultat = data
#endif

    status=0

  END SUBROUTINE PAR_REDUCE_ELEMENT_R2
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_REDUCE_ELEMENT_R3
  !
  ! !DESCRIPTION: reduce 3D arrays element-wise
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_REDUCE_ELEMENT_R3 (DATA, OP, RESULTAT, STATUS, ALL, ROW)
    !
    ! !USES:
    !
    use dims, only : CheckShape
    !
    ! !INPUT PARAMETERS:
    !
    real,              intent(in) :: data(:,:,:)
    character(len=3),  intent(in) :: op    ! 'MAX', 'MIN' or 'SUM'
    logical, intent(in), optional :: all   ! mimic mpi_allreduce instead of mpi_reduce
    logical, intent(in), optional :: row   ! limit to PE on row_comm
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: resultat(:,:,:)
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    1 Mar 2012 - P. Le Sager - v0
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/PAR_REDUCE_ELEMENT_R3'
    logical :: l_row, l_all
    integer :: l_root, l_comm, l_id, l_op


    call CheckShape( shape(data), shape(resultat), status )
    IF_NOTOK_RETURN(status=1)


#ifdef MPI

    l_row=.false.  ; if(present(row)) l_row=row
    l_all=.false.  ; if(present(all)) l_all=all

    if(l_row)then
       l_root = 0               ! by our own design
       l_comm = row_comm
       l_id   = jcol
    else
       l_root = root
       l_comm = localComm
       l_id   = myid
    end if

    ! degenerate cases first
    if(l_row.and.(npe_lon==1)) then
       resultat = data
    else

       SELECT CASE( OP )
       case('sum', 'Sum', 'SUM')
          l_op = MPI_SUM
       case('max', 'Max', 'MAX')
          l_op = MPI_MAX
       case('min', 'Min', 'MIN')
          l_op = MPI_MIN
       case default
          write(gol,*) 'UNSUPPORTED OPERATION :', OP; status=1
          IF_NOTOK_RETURN(status=1)
       END SELECT

       if (l_all) then
          call MPI_ALLREDUCE(data, resultat, size(data), my_real, l_op, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       else
          call MPI_REDUCE(data, resultat, size(data), my_real, l_op, l_root, l_comm, ierr)
          IF_NOTOK_MPI(status=1)
       end if

    end if

#else
    resultat = data
#endif

    status=0

  END SUBROUTINE PAR_REDUCE_ELEMENT_R3
  !EOC


END MODULE ParTools