!#################################################################
!
! this module holds all specific routines needed to run MPI
! WP january 2003
! 
!### macro's #####################################################
!
#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_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

module mpi_comm

  use GO       , only : gol, goErr, goPr
  use mpi_const

  implicit none

  ! --- in/out -------------------------------

  private

  !public :: startmpi
  public :: abortmpi
  !public :: stopmpi
  public :: initialize_domains
  public :: barrier_t, barrier_k, barrier
  public :: scatter_after_read_k_3d, scatter_after_read_k
  public :: scatter_after_read_t
  public :: gather_tracer_k_3d, gather_tracer_k
  public :: gather_tracer_t
  public :: tracer_to_k, k_to_tracer
  public :: escape_tm_mpi
  !public :: summasr
  public :: dump_field3d, dump_field4d

  ! --- const ----------------------------------

  character(len=*), parameter  ::  mname = 'mpi_comm'


contains


  ! ================================================================


  subroutine tracer_to_k( array_k, array_tracer, im, jm, lm, &
       ih, jh, lh, ntracer, ntracerloc, tracerar )
    !
    ! re-distribute an array which is distributed over the k-index
    ! to an array which is distributed over the tracers
    ! uses MPI_ALLTOALLV and a help array
    ! we have to use MPI_ALLTOALLV: the number of data received from the PEs
    ! can be different from PE to PE
    !
    use omp_partools, only : my_omp_domain

    implicit none

    ! input/output
    integer,intent(in) :: im, jm, lm
    integer,intent(in) :: ih, jh, lh
    integer,intent(in) :: ntracer, ntracerloc
    integer,dimension(0:npes-1),intent(in)  :: tracerar  
    real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracerloc) :: array_tracer
    real,dimension(1-ih:im+ih,1-jh:jm+jh,1   :lmloc,ntracer) :: array_k

    ! local
    real,dimension((im+2*ih)*(jm+2*jh)*(lm+2*lh)*ntracerloc) :: array_tracer_help
    integer,dimension(0:npes-1) :: sdispls   
    integer,dimension(0:npes-1) :: rdispls   
    integer,dimension(0:npes-1) ::  sendcounts
    integer,dimension(0:npes-1) ::  recvcounts
    integer :: recvcount
    integer :: my_vector
    integer :: slab2d
    integer :: itracer, i, j, k, ii, jj, iproc, kadd, kloop
    logical :: okdebug=.false.
    ! Openmp parameters
    integer :: ns, ne

    ! start

    slab2d = (im+2*ih)*(jm+2*jh)
    !
    ! fill necessary arrays
    !
    do i=0,npes-1
       recvcounts(i) = lmloc*slab2d*tracerar(i)
    end do

    rdispls(0)       = 0
    !
    do i=1,npes-1
       rdispls   (i) = rdispls(i-1)+recvcounts(i-1)
    end do
    !
    do j=0,npes-1
       sendcounts(j) = lmar(j)*slab2d*ntracerloc
    end do
    !
    sdispls(0)       = 0
    !
    do j=1,npes-1
       sdispls(j)    = sdispls(j-1)+sendcounts(j-1)
    end do
    !
    ! data received from npes PE's, numbered 0,... npes-1
    !
    !$OMP PARALLEL &
    !$OMP  default (none) &
    !$OMP  shared  (npes, im, ih, jm, jh, ntracerloc, lmar, array_tracer, &
    !$OMP           array_tracer_help, slab2d) &
    !$OMP  private (i, j, k, ns, ne, ii, jj, kadd, kloop)
    call my_omp_domain(1,ntracerloc,ns,ne)
    jj = 0
    do iproc  =0,npes-1
       !
       ! data received from each processor for ntracerloc tracers
       !
       !
       ii = jj + (ns - 1) * lmar(iproc) * slab2d
       kadd = 0
       do kloop = 0, iproc - 1
          kadd = kadd + lmar(kloop)
       end do
       do itracer=ns,ne
          do k=kadd+1, kadd+lmar(iproc)
             do j=1-jh,jm+jh
                do i=1-ih,im+ih
                   ii = ii + 1
                   array_tracer_help(ii) = array_tracer(i,j,k,itracer) 
                end do
             end do
          end do
       end do
       jj = jj + ntracerloc * lmar(iproc) * slab2d
    end do
    !$OMP END PARALLEL
    !
    !
    !
    ! fill array_tracer_help
    !
    !if(okdebug)then
    !     write(6,*)'COUNTCHECK t_to_k',myid,ntracerloc,lmloc
    !     write(6,*)'COUNTCHECK t_to_k',myid,tracerar
    !     write(6,*)'COUNTCHECK t_to_k',myid,sdispls
    !     write(6,*)'COUNTCHECK t_to_k',myid,sendcounts
    !     write(6,*)'COUNTCHECK t_to_k',myid,recvcounts
    !     write(6,*)'COUNTCHECK t_to_k',myid,'to send ',sum(array_tracer_help)
    !end if

    call MPI_ALLTOALLV( array_tracer_help, sendcounts, sdispls, MY_REAL, &
         array_k          , recvcounts, rdispls, MY_REAL, &
         localComm, ierr )

    !
    ! array array_tracer_help contains all the data, which are now 
    ! distributed properly over the local tracer array
    !
    ! the structure of array array_tracer_help is now as follows:
    !
    !
    !     PE0            itrace1      lmar(0)          2D slabs  ---
    !     PE0            itrace2      lmar(0)          2D slabs    |
    ! itracer=1,ntracerloc from PE0:
    !     PE0            itrace3      lmar(0)          2D slabs    |
    !         ...                                                  |
    !     PE0            ntracerloc   lmar(0)          2D slabs  ---
    !
    !     PE1            itrace1      lmar(1)          2D slabs  ---
    !     PE1            itrace2      lmar(1)          2D slabs    |
    ! itracer=1,ntracerloc from PE1:
    !     PE1            itrace3      lmar(1)          2D slabs    |
    !         ...                                                  |
    !     PE1            ntracerloc   lmar(0)          2D slabs  ---
    !
    !     ....
    !     ....
    !     ....
    !     ....
    !
    !     PEntracerloc   itrace1      lmar(ntracerloc) 2D slabs  ---
    !     PEntracerloc   itrace2      lmar(ntracerloc) 2D slabs    |
    ! itracer=1,ntracerloc from PEnpes-1
    !     PEntracerloc   itrace3      lmar(ntracerloc) 2D slabs    |
    !         ...                                                  |
    !     PEntracerloc   ntracerloc   lmar(ntracerloc) 2D slabs  ---
    !
  end subroutine tracer_to_k



  subroutine k_to_tracer( array_k, array_tracer, im, jm, lm, &
       ih, jh, lh, ntracer, ntracerloc, tracerar )
    !
    ! re-distribute an array which is distributed over the k-index
    ! to an array which is distributed over the tracers
    ! uses MPI_ALLTOALLV and a help array
    ! we have to use MPI_ALLTOALLV: the number of data received from the PEs
    ! can be different from PE to PE

    use omp_partools, only : my_omp_domain

    implicit none

    ! input/output
    integer,intent(in) :: im,jm,lm,ih,jh,lh,ntracer,ntracerloc
    integer,dimension(0:npes-1),intent(in)        :: tracerar  
    real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracerloc), intent(out) :: array_tracer
    real,dimension(1-ih:im+ih,1-jh:jm+jh,1   :lmloc,ntracer), intent(inout) :: array_k

    ! local
    real,dimension((im+2*ih)*(jm+2*jh)*(lm+2*lh)*ntracerloc) :: array_tracer_help
    integer sdispls   (0:npes-1)
    integer rdispls   (0:npes-1)
    integer sendcounts(0:npes-1)
    integer recvcounts(0:npes-1)
    integer recvcount
    integer my_vector
    integer slab2d
    integer itracer, i, j, k, ii, jj, iproc, kadd, kloop
    ! Openmp parameters
    integer ns, ne

    ! start

    slab2d = (im+2*ih)*(jm+2*jh)
    !
    ! fill necessary arrays
    !
    do i=0,npes-1
       sendcounts(i) = lmloc*slab2d*tracerar(i)
    end do

    sdispls(0)       = 0
    !
    do i=1,npes-1
       sdispls   (i) = sdispls(i-1)+sendcounts(i-1)
    end do
    !
    do j=0,npes-1
       recvcounts(j) = lmar(j)*slab2d*ntracerloc
    end do
    !
    rdispls(0)       = 0
    !
    do j=1,npes-1
       rdispls(j)    = rdispls(j-1)+recvcounts(j-1)
    end do

    !
    ! fill array_tracer_help
    !
    call MPI_ALLTOALLV( array_k          , sendcounts, sdispls, MY_REAL, &
         array_tracer_help, recvcounts, rdispls, MY_REAL, &
         localComm, ierr )
    !
    ! array array_tracer_help contains all the data, which are now distributed properly over
    ! the local tracer array
    !
    ! the structure of array array_tracer_help is now as follows:
    !
    !
    !     PE0            itrace1      lmar(0)          2D slabs  ---
    !     PE0            itrace2      lmar(0)          2D slabs    |
    ! itracer=1,ntracerloc from PE0
    !     PE0            itrace3      lmar(0)          2D slabs    |
    !         ...                                                  |
    !     PE0            ntracerloc   lmar(0)          2D slabs  ---
    !     PE1            itrace1      lmar(1)          2D slabs  ---
    !     PE1            itrace2      lmar(1)          2D slabs    |
    ! itracer=1,ntracerloc from PE1
    !     PE1            itrace3      lmar(1)          2D slabs    |
    !         ...                                                  |
    !     PE1            ntracerloc   lmar(0)          2D slabs  ---
    !     ....
    !     ....
    !     ....
    !     ....
    !     PEntracerloc   itrace1      lmar(ntracerloc) 2D slabs  ---
    !     PEntracerloc   itrace2      lmar(ntracerloc) 2D slabs    |
    ! itracer=1,ntracerloc from PEnpes-1
    !     PEntracerloc   itrace3      lmar(ntracerloc) 2D slabs    |
    !         ...                                                  |
    !     PEntracerloc   ntracerloc   lmar(ntracerloc) 2D slabs  ---
    !
    !
    ! data received from npes PE's, numbered 0,... npes-1
    !
    !$OMP PARALLEL &
    !$OMP  default (none) &
    !$OMP  shared  (npes, im, ih, jm, jh, ntracerloc, lmar, array_tracer, &
    !$OMP           array_tracer_help, slab2d) &
    !$OMP  private (i, j, k, ns, ne, ii, jj, kadd, kloop)
    call my_omp_domain (1,ntracerloc,ns,ne)
    jj = 0
    do iproc  =0,npes-1
       !
       ! data received from each processor for ntracerloc tracers
       !
       ii = jj + (ns - 1) * lmar(iproc) * slab2d
       kadd = 0
       do kloop = 0, iproc - 1
          kadd = kadd + lmar(kloop)
       end do
       do itracer=ns,ne
          do k=kadd+1,kadd+lmar(iproc)
             do j=1-jh,jm+jh
                do i=1-ih,im+ih
                   ii = ii + 1
                   array_tracer(i,j,k,itracer) = array_tracer_help(ii)
                end do
             end do
          end do   !  kloop
       end do  ! itracer
       jj = jj + ntracerloc * lmar(iproc) * slab2d
    end do  ! iproc
    !$OMP END PARALLEL
    !
    !
  end subroutine k_to_tracer



  !  ! subroutine initilizes mpi, 
  !  ! sets up localComm, 
  !  ! checks the number of processors, and
  !  ! proceeds to set the paralel mode of the model to 'tracers'
  !
  !  subroutine startmpi( status )
  !
  !    ! --- in/out ------------------------------
  !    
  !    integer, intent(out)  ::  status
  !    
  !    ! --- const ------------------------------
  !
  !    character(len=*), parameter ::  rname = mname//'/Par_Done'
  !    
  !    ! --- begin -------------------------------
  !
  !    !AJS: moved to par_tools >>>
  !    !call MPI_INIT(ierr)
  !    !call MPI_COMM_RANK( localComm, myid, ierr )
  !    !call MPI_COMM_SIZE( localComm, npes, ierr )
  !    !<<<
  !    
  !    !WP! my_real comes from mpif.h, other declarations may be needed sometimes
  !    my_real = MPI_DOUBLE_PRECISION 
  !
  !    ! allocate arrays
  !    allocate( lmar      (0:npes-1) )
  !    allocate( ntracet_ar(0:npes-1) )
  !
  !    call initialize_domains( status )
  !    if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if
  !
  !    which_par='tracer'
  !    previous_par(:)='tracer'
  !
  !    write (gol,'("startmpi: MPI initialized, starting")'); call goPr
  !    write (gol,'("startmpi: data is now paralel over ",a)') trim(which_par); call goPr
  !
  !    call barrier
  !    
  !    ! ok
  !    status = 0
  !
  !  end subroutine startmpi


  !AJS:>>> now in partools
  !  subroutine stopmpi
  !    implicit none
  !
  !    ! deallocate arrays
  !    deallocate( lmar       )
  !    deallocate( ntracet_ar )
  !
  !    call barrier
  !    call mpi_finalize(ierr)
  !    STOP 'stopmpi: EINDE'
  !
  !  end subroutine stopmpi
  !<<<


  subroutine abortmpi
    implicit none

    ! deallocate arrays

    call mpi_abort(localComm, 999, ierr)
    STOP 'abortmpi: FORCED quit'

  end subroutine abortmpi




  subroutine escape_tm_mpi(msg)
    !--------------------------------------------------------------
    !
    ! abnormal termination of a model run
    ! 
    ! msg       character string to be printed on unit kmain
    !
    !-------------------------------------------------------------
    use dims,        only : okdebug, kmain, itau
    use datetime,    only : tstamp
    implicit none
    !      
    character(*),intent(in) :: msg
    !
    write (gol,'(///4(1x,39("--")/))'); call goPr
    call tstamp(kmain,itau,msg)
    write (gol,'(///4(1x,39("--")/))'); call goPr

    call abortmpi

  end subroutine escape_tm_mpi
  
  !-------------------------------------------------------------------------
  !                    TM5                                                 !
  !-------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  initialize_domains
  !
  ! !DESCRIPTION: <Add description here !!>
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine initialize_domains( status )
    !
    ! !USES:
    !
    use dims      , only : lm
    use chem_param, only : ntracet
#ifdef with_tracerorder
    Use chem_param, only : distribution
    Use chem_param, only : tracerorder_procs
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REVISION HISTORY: 
    !   30 Mar 2010 - P. Le Sager - included the with_tracerorder from JadB
    !   30 Jun 2010 - Arjo Segers - added check on number of processors for
    !                                 tracerorder 
    !    4 Feb 2011 - WP - allow for npes > nlevels
    !
    ! !REMARKS:
    !
    !EOP
    !-----------------------------------------------------------------------
    !BOC

    ! --- const ----------------------------------

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

    ! --- begin ----------------------------------

!#ifndef without_chemistry
!    
!    if ( npes > lm(1) .or. npes > ntracet ) then
!       write (gol,'("npes should not exceed ntracet or lm:")'); call goErr
!       write (gol,'("  npes    : ",i4)') npes; call goErr
!       write (gol,'("  lm      : ",i4)') lm(1); call goErr
!       write (gol,'("  ntracet : ",i4)') ntracet; call goErr
!       write (gol,'("in ",a)') rname; call goErr; status=1; return
!    end if
!
!#else
!    ! allow for more processors than #levels
    if ( npes > ntracet ) then
       write (gol,'("npes should not exceed ntracet:")'); call goErr
       write (gol,'("  npes    : ",i4)') npes; call goErr
       write (gol,'("  ntracet : ",i4)') ntracet; call goErr
       write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

!#endif
    
    ! lmar    = lm(1) ! Jadb: This is silly because lmar will be redefined in the function.
    call determine_lmar(lmar,npes,lm(1))

#ifdef with_tracerorder
    ! JadB:
    ! I pre-determined the distribution. In normal cases, the result will be the same
    ! But if you have many PEs, it could be different.
    ! This because one interdependent block of six tracers may
    ! be much bigger than what an average PE has in total.
    ntracet_ar = distribution
    ! check ...
    if ( npes /= tracerorder_procs ) then
       write (gol,'("number of processors (",i2,") should be equal to ")') npes; call goErr
       write (gol,'("  number for which tracer order was generated (",i2,")")') tracerorder_procs; call goErr
       write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if
#else
    call determine_lmar(ntracet_ar,npes,ntracet)
#endif

    lmloc   = lmar(myid)
    ntracetloc   = ntracet_ar(myid)

    ! JadB: commented, since it simply sets first arg to zero.
    ! call determine_first(PE_FIRST_TRACER,ntracet_ar)
    ! call determine_first(PE_FIRST_L     ,lmar   )
    PE_FIRST_TRACER = 0
    PE_FIRST_L = 0

    ! JadB: At least not write down level distribution after the word 'Tracer'
    ! JadB: Now it writes down its ID, then zero and then the distribution.
    write (gol,*) 'initialize_domains: myid, PE_FIRST_TRACER', myid, PE_FIRST_TRACER, ntracet_ar; call goPr
    write (gol,*) 'initialize_domains: myid, PE_FIRST_L     ', myid, PE_FIRST_L,      lmar      ; call goPr  
    call set_domain('tracer')
    call set_domain('levels')
    if ( root_t /= 0 .or. root_k /= 0 .or. root /= 0 ) then
       write (gol,'("initialize_domains: ERROR  root_t /= root_k")'); call goPr
       write (gol,'("initialize_domains: problems for budget calculations")'); call goPr
       write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    call determine_tracer_active

    ! ok
    status = 0

  end subroutine initialize_domains
  !EOC

  !-------------------------------------------------------------------------
  !                    TM5                                                 !
  !-------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  determine_lmar
  !
  ! !DESCRIPTION: <Add description here !!>
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine determine_lmar(xlmar,npesmax,xlm)
    !
    ! !USES:
    !
    implicit none
    !
    ! !REVISION HISTORY: 
    !    4 Feb 2011 - WP - allow for npes > nlevels
    !
    ! !REMARKS:
    !
    !EOP
    !----------------------------------------------------------------------
    !BOC

    integer resdiv, remainder, i, count
    integer npesmax
    integer xlm
    integer xlmar(0:npesmax-1)
    !
    ! this subroutine fills elements 0,npes-1 of array lmar
    ! with a number equal to lm div npes = resdiv (integer division)
    ! if resdiv*npes not equal to lm, there is a remainder
    ! the remainder is divided over the last elements 
    ! (npes-1, npes-remainder)
    !
    ! JadB:
    ! I would like to add interdependent groups of tracers. The point
    ! is that some tracers need one another during the sedimentation
    ! step. An example of these are the aerosol numbers and compound
    ! masses. Those tracers must be handled by the same processor. The
    ! independent tracers will be distributed, resulting in a total
    ! distribution as even as possible. This is done somewhere else.
    ! Therefore this routine is not used for tracers.
    !
    ! It looks as if it is logical that the tracers handled by a single
    ! processor have adjacent identifiers. This would mean that the order
    ! should be adapted to interdependencies in chem_param. In this
    ! subroutine, the distribution could be adapted if it does not fit.
    ! The problem of this method is that there is communication between
    ! chem_param and this method. And hacking one will cause problems later
    ! on.
    !
    ! It could be that the being adjacent is not that big of a deal. Then,
    ! the distribution can be done here and the chem_param can be left alone.

    if(npes.gt.xlm) then

       xlmar       = 0
       do i=0,xlm-1  ! first xlm pe's get a value, others get nothing
          xlmar(i)=1
       enddo

    else 
       resdiv = xlm/npes
       remainder  = xlm - resdiv*npes

       xlmar       = 0
       do i=0,npes-1
          xlmar(i) = resdiv
       end do
       !
       count = remainder
       ! CMK changed!  PE0 should have less work (surface = more iterations)
       i = npes-1
       do
          if(count == 0)exit 
          xlmar(i) = xlmar(i) + 1
          count   = count   - 1
          i       = i       - 1   ! CMK changed!
       end do

    endif

    do i=0,npes-1
       write (gol,'("distribution over PE s: i, array(i) ",2i4)') i, xlmar(i); call goPr
    end do

  end subroutine determine_lmar
  !EOC


!=============================================================================
!  
!  PLS-JadB : commented for now, since just set pe_first=0
!  
!=============================================================================
!  subroutine determine_first(pe_first,xlmar)
!    !
!    ! we assume a distribution over PEs with:
!    !
!    ! lmar(0)    lmar(1)    lmar(2)     ....   lmar(npes-1)
!    !
!    !
!    implicit none
!    integer :: pe_first, xlmar(0:npes-1)
!    integer :: iproc
!
!    pe_first=0
!
!    !if(xlmar(npes-1) == 0)then
!    !   call barrier
!    !   call mpi_finalize(ierr)
!    !   STOP 'ERROR: NO DATA on LAST PE'
!    !end if
!    !
!    !iproc = npes-1
!    !
!    !do while (iproc.ge.0)
!    !   if(xlmar(iproc) == 0)exit
!    !   iproc = iproc - 1
!    !end do
!    !pe_first = iproc + 1
!
!  end subroutine determine_first


  ! ***


  subroutine determine_tracer_active
    integer :: n,offsetn,offsetnp,ipos
    offsetn = sum(ntracet_ar(0:myid-1)) + 1
    offsetnp = offsetn + ntracetloc - 1 

    do n=1,ntracet
       if (n >= offsetn .and. n <= offsetnp) then
          tracer_active(n) = .true.
          tracer_loc(n)    = n-offsetn+1    ! the counter in the 
       else
          tracer_active(n) = .false.
          tracer_loc(n)    = -9           !  for checks....
       end if
    end do

    ipos = 1
    do n=0,npes-1
       ! proc_tracer - which processor handles tracer i ? 
       proc_tracer(ipos:ipos+ntracet_ar(n)-1) = n
       ipos = ipos + ntracet_ar(n)
    end do
    !write (gol,'("determine_tracer_active: Tracers active for processor ",i4,60l2)') &
    !                  myid, tracer_active; call goPr
    !write (gol,'("determine_tracer_active: Local tracer number on proc   "i4,60i2)') &
    !                  myid, tracer_loc   ; call goPr
    !write (gol,'("determine_tracer_active: Tracer numbers processed by  ",i4,60i2)') &
    !                  myid, proc_tracer  ; call goPr

  end subroutine determine_tracer_active


  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  set_domain
  !
  ! !DESCRIPTION: <Add description here !!>
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine set_domain(switch)
    !
    ! !USES:
    !
    use dims
    implicit none
    !
    ! !INPUT PARAMETERS:
    !
    character(len=6),intent(in) :: switch
    !
    ! !REVISION HISTORY: 
    !    4 Feb 2011 - WP - allow for npes > nlevels
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    integer :: orig_group2,orig_group,com_trac_group,com_lev_group
    integer :: a,b,i
    integer :: new_rank,lev_rank
    integer,allocatable,dimension(:) :: inc_ranks

    myid_t=999 !WP! only processors in the communicator receive a myid
    myid_k=999 !WP! only processors in the communicator receive a myid

    !WP! set up communication world for tracers
    if ( switch == 'tracer' ) then
       if ( okdebug ) then
          write (gol,'("set_domain: Setting up communicator world for tracer domain",i4)') myid; call goPr
       end if
       call MPI_COMM_GROUP(localComm,orig_group,ierr)
       b=min(ntracet,npes)
       allocate(inc_ranks(b)) !WP! to hold ranks of PE's
       do i=1,b
          inc_ranks(i)=PE_FIRST_TRACER+i-1 !WP! fill ranks array
       end do
       call MPI_GROUP_INCL(orig_group,b,inc_ranks,com_trac_group,ierr)
       deallocate(inc_ranks)
       call barrier
       call mpi_comm_create(localComm,com_trac_group,com_trac,ierr) 
       call barrier
       !WP! ranks in new communicator
       if ( ntracetloc /= 0 ) call MPI_COMM_RANK(com_trac,myid_t,ierr)
       root_t=PE_FIRST_TRACER

       write (gol,*) 'set_domain: In tracer domain myid: ',myid, ' has number ',myid_t,' and uses as root: ',root_t; call goPr
       !WP! first free old communciator group
       call MPI_GROUP_FREE(com_trac_group,ierr)
       !WP!
       !WP! com_trac is setup now, proceed with com_lev
       !WP!

    end if
    if ( switch == 'levels' ) then

       if ( okdebug ) then
          write (gol,*) 'set_domain: Setting up communicator world for levels domain ' ; call goPr
       end if
       call barrier
       call MPI_COMM_GROUP(localComm,orig_group2,ierr)
       b=min(npes,lm(1))
       allocate(inc_ranks(b))
       do i=1,b
          inc_ranks(i)=PE_FIRST_L+i-1
       end do
       call MPI_GROUP_INCL(orig_group2,b,inc_ranks,com_lev_group,ierr)
       deallocate(inc_ranks)
       call barrier
       call mpi_comm_create(localComm,com_lev_group,com_lev,ierr)
       call barrier
       !WP! ranks in new communicator
       if ( lmloc /= 0) call MPI_COMM_RANK(com_lev,myid_k,ierr)
       call barrier
       root_k=PE_FIRST_L
       write (gol,*) 'set_domain: In levels domain myid: ',myid, ' has number ',myid_k,' and uses as root: ',root_k; call goPr
       !WP! first free old communciator group
       call MPI_GROUP_FREE(com_lev_group,ierr)
       call barrier
    end if

    !call MPI_COMM_FREE(com_trac,ierr)  !WP! first free old communciator group
    !call MPI_COMM_FREE(com_lev,ierr)  !WP! first free old communciator group

    call MPI_BARRIER(localComm,ierr)

  end subroutine set_domain



  subroutine barrier
    implicit none

#ifdef with_barrier
    call MPI_BARRIER(localComm,ierr)
#endif

  end subroutine barrier



  subroutine barrier_t
    implicit none

#ifdef with_barrier
    call MPI_BARRIER(com_trac,ierr)
#endif

  end subroutine barrier_t



  subroutine barrier_k
    implicit none

    call MPI_BARRIER(com_lev,ierr)

  end subroutine barrier_k


  !
  ! gather an array  
  !    tracer(1-ih,im+ih,1-jh,jm+jh,1-lh,lm+lh,ntracer) 
  ! on root_t
  ! from local arrays  
  !    tracerloc(1-ih,im+ih,1-jh,jm+jh,1-lh:lm+lh,ntracerloc) 
  ! on all PE's where  
  ! the array is distributed over the n-index
  !

  subroutine gather_tracer_t(tracer,im,jm,lm,ih,jh,lh,ntracer,tracerloc,broadcast)

    use dims, only : CheckShape

    ! input/output
    logical,intent(in)                    :: broadcast
    integer,intent(in)                    :: im,jm,lm
    integer,intent(in)                    :: ih,jh,lh
    integer,intent(in)                    :: ntracer

    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracer), &
    !     intent(inout)   :: tracer    
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracetloc), &
    !     intent(in)      :: tracerloc    
    real,dimension(:,:,:,:), intent(inout)   ::  tracer    
    real,dimension(:,:,:,:), intent(in)      ::  tracerloc    

    ! local
    integer :: slab3d
    integer,dimension(0:npes-1) :: displ
    integer,dimension(0:npes-1) :: recvcounts
    integer :: sendcount, nsend
    integer :: itracer, i , ierr
    logical :: okdebug=.false.

    ! start

    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh,ntracer   /), shape(tracer   ) )
    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh,ntracetloc/), shape(tracerloc) )

    if ( okdebug ) print*,'Gather_tracer_t: total on pe ',myid, ' = ',sum(tracerloc)
    call barrier
    slab3d    = (im+2*ih)*(jm+2*jh)*(lm+2*lh)
    displ(:)=0  ! offset

    do i=1,npes-1
       displ(i) = displ(i-1) + ntracet_ar(i-1)*slab3d ! calculate offset slabs
    end do
    !
    do i=0,npes-1
       recvcounts(i) = ntracet_ar(i)*slab3d ! receive nr of slabs from other PE's
    end do
    !
    sendcount = ntracetloc*slab3d  ! send the local nr of slabs
    !
    !
    call barrier
    ! gather on PE root_t
    !call MPI_GATHERV( &
    !       tracerloc(1-ih,1-jh,1-lh,1), sendcount        , my_real, &
    !       tracer   (1-ih,1-jh,1-lh,1), recvcounts, displ, my_real, &
    !       root_t, localComm, ierr )
    call MPI_GATHERV( &
         tracerloc(:,:,:,:), sendcount        , my_real, &
         tracer   (:,:,:,:), recvcounts, displ, my_real, &
         root_t, localComm, ierr )

    if ( ierr /= 0 ) call escape_tm_mpi('gather_tracer_t: ERROR')
    call barrier
    !
    if ( myid==root_t .and. okdebug ) then 
       print*,'Gather_tracer_t: total on root ',myid,' = ',sum(tracer)
    end if

    if ( broadcast ) then
       nsend=(im+2*ih)*(jm+2*jh)*(lm+2*lh)*ntracer
       call mpi_bcast( tracer, nsend, my_real, root_t, localComm, ierr )
    end if

  end subroutine gather_tracer_t



  !
  ! gather   an array  
  !    tracer    (1-ih,im+ih,1-jh,jm+jh,1-lh,lm+lh,ntracer)
  ! on root_k
  ! from local arrays  
  !    tracerloc(1-ih,im+ih,1-jh,jm+jh,1-lh,lmloc+lh,ntracer) 
  ! on all PE's where  
  ! the array is distributed over the l-index
  !

  subroutine gather_tracer_k_3d( tracer, im,jm,lm, ih,jh,lh, tracerloc, broadcast )

    use dims, only : CheckShape

    ! input/output
    logical,intent(in)                        :: broadcast
    integer,intent(in)                        :: im,jm,lm
    integer,intent(in)                        :: ih,jh,lh

    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh), intent(inout)   :: tracer    
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1   :lmloc), intent(in)      :: tracerloc
    real, dimension(:,:,:), intent(inout)   ::  tracer    
    real, dimension(:,:,:), intent(in)      ::  tracerloc

    ! local
    integer :: slab2d
    integer,dimension(0:npes-1) :: displ
    integer,dimension(0:npes-1) :: recvcounts
    integer :: sendcount, nsend
    integer :: i , ierr
    logical :: okdebug=.false.

    ! start

    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh/), shape(tracer   ) )
    call CheckShape( (/im+2*ih,jm+2*jh,lmloc  /), shape(tracerloc) )

    if ( okdebug ) &
         print*,'Gather_tracer_k: total on pe ',myid,' = ',sum(tracerloc)

    slab2d    = (im+2*ih)*(jm+2*jh)
    displ(:)=0  ! offset
    displ(0) = lh*slab2d  ! offset with lh as these do not exist on tracerloc

    do i=1,npes-1
       displ(i) = displ(i-1) + lmar(i-1)*slab2d ! calculate offset slabs
    end do
    !
    do i=0,npes-1
       recvcounts(i) = lmar(i)*slab2d ! receive nr of slabs from other PE's
    end do
    !
    sendcount = lmloc*slab2d  ! send the local nr of slabs
    !
    !
    call barrier
    ! gather on PE root_k
    !call MPI_GATHERV( tracerloc(1-ih,1-jh,   1), sendcount, my_real, &
    !                  tracer   (1-ih,1-jh,1-lh), recvcounts, displ, my_real, &
    !                  root_k, localComm, ierr )
    call MPI_GATHERV( tracerloc(:,:,:), sendcount, my_real, &
         tracer   (:,:,:), recvcounts, displ, my_real, &
         root_k, localComm, ierr )
    if(ierr/=0) call escape_tm_mpi('gather_tracer_k: error ')
    call barrier
    !
    if ( myid==root_k .and. okdebug ) then 
       print*,'Gather_tracer_k: total on root ',myid,' = ',sum(tracer)
    end if

    if ( broadcast ) then
       nsend=(im+2*ih)*(jm+2*jh)*(lm+2*lh)
       call mpi_bcast(tracer,nsend,my_real,root_k,localComm,ierr)
    end if

  end subroutine gather_tracer_k_3d


  ! ***


  subroutine gather_tracer_k(tracer,im,jm,lm,ih,jh,lh,ntracer,tracerloc,broadcast)

    use dims, only : CheckShape

    ! input/output
    logical,intent(in)                        :: broadcast
    integer,intent(in)                        :: im,jm,lm
    integer,intent(in)                        :: ih,jh,lh
    integer,intent(in)                        :: ntracer

    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracer), &
    !     intent(inout)   :: tracer    
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1   :lmloc,ntracer), &
    !     intent(in)      :: tracerloc
    real, dimension(:,:,:,:), intent(inout)   ::  tracer    
    real, dimension(:,:,:,:), intent(in)      ::  tracerloc


    ! local
    integer :: slab2d
    integer,dimension(0:npes-1) :: displ
    integer,dimension(0:npes-1) :: recvcounts
    integer :: sendcount, nsend
    integer :: itracer, i , ierr
    logical :: okdebug=.false.

    ! start

    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh,ntracer/), shape(tracer   ) )
    call CheckShape( (/im+2*ih,jm+2*jh,  lmloc,ntracer/), shape(tracerloc) )

    if ( okdebug ) &
         print*,'Gather_tracer_k: total on pe ',myid,' = ',sum(tracerloc)

    slab2d    = (im+2*ih)*(jm+2*jh)
    displ(:)=0  ! offset
    displ(0) = lh*slab2d  ! offset with lh as these do not exist on tracerloc

    do i=1,npes-1
       displ(i) = displ(i-1) + lmar(i-1)*slab2d ! calculate offset slabs
    end do
    !
    do i=0,npes-1
       recvcounts(i) = lmar(i)*slab2d ! receive nr of slabs from other PE's
    end do
    !
    sendcount = lmloc*slab2d  ! send the local nr of slabs
    !
    !
    call barrier
    do itracer = 1, ntracer
       ! gather on PE root_k
       !call MPI_GATHERV( tracerloc(1-ih,1-jh,   1,itracer), sendcount        , my_real, &
       !                  tracer   (1-ih,1-jh,1-lh,itracer), recvcounts, displ, my_real, &
       !                  root_k,localComm, ierr)
       call MPI_GATHERV( tracerloc(:,:,:,itracer), sendcount        , my_real, &
            tracer   (:,:,:,itracer), recvcounts, displ, my_real, &
            root_k,localComm, ierr)
       if(ierr/=0) call escape_tm_mpi('gather_tracer_k: error ')
    end do
    call barrier
    !
    if ( myid==root_k .and. okdebug ) then 
       print*,'Gather_tracer_k: total on root ',myid,' = ',sum(tracer)
    end if

    if ( broadcast ) then
       nsend=(im+2*ih)*(jm+2*jh)*(lm+2*lh)*ntracer
       call mpi_bcast(tracer,nsend,my_real,root_k,localComm,ierr)
    end if

  end subroutine gather_tracer_k


  ! *********************************************


  !
  ! scatter an array  
  !    tracer   (1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracer)
  ! on root (pe=0)
  ! over local arrays  
  !    tracerloc(1-ih:im+ih,1-jh:jm+jh,1:lmloc   ,ntracer)
  ! on all PE's so that 
  ! the array is distributed over the l-index
  !

  subroutine scatter_after_read_k_3d( tracer_read, im,jm,lm, ih,jh,lh, &
       tracerloc, send_id )

    use mpi_const,only : lmar,myid,lmloc
    use dims, only : CheckShape

    ! input/output
    integer, intent(in)              :: im,jm,lm,ih,jh,lh
    integer, intent(in)              :: send_id 
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh) :: tracer_read
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1   :lmloc) :: tracerloc   
    real, dimension(:,:,:), intent(in)    :: tracer_read
    real, dimension(:,:,:), intent(inout) :: tracerloc   

    ! local
    real, dimension(1-ih:im+ih,1-jh:jm+jh) :: sum_one,sum_all
    integer                         :: slab2d
    integer, dimension(0:npes-1)    :: displ     
    integer, dimension(0:npes-1)    :: sendcounts
    integer                         :: recvcount
    integer                         :: nsum
    integer                         :: i 

    ! start

    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh/), shape(tracer_read) )
    call CheckShape( (/im+2*ih,jm+2*jh,lmloc  /), shape(tracerloc  ) )

    slab2d    = (im+2*ih)*(jm+2*jh)
    !
    ! offset for gather-scatter
    !
    displ(:) = 0
    displ(0) = lh*slab2d
    do i=1,npes-1
       displ(i) = displ(i-1) + lmar(i-1)*slab2d
    end do
    !
    do i=0,npes-1
       sendcounts(i) = lmar(i)*slab2d
    end do
    !
    recvcount = lmloc*slab2d
    !
    ! call MPI_SCATTERV
    !
    !call MPI_SCATTERV( tracer_read(1-ih,1-jh,1-lh), sendcounts, displ, MY_REAL, & 
    !                   tracerloc  (1-ih,1-jh,1   ), recvcount , MY_REAL, &
    !                   send_id, localComm, ierr )
    call MPI_SCATTERV( tracer_read(:,:,:), sendcounts, displ, MY_REAL, & 
         tracerloc  (:,:,:), recvcount        , MY_REAL, &
         send_id, localComm, ierr )

  end subroutine scatter_after_read_k_3d


  ! ***


  subroutine scatter_after_read_k( tracer_read, im,jm,lm, ih,jh,lh, ntracer, &
       tracerloc, send_id )

    use mpi_const,only : lmar,myid,lmloc
    use dims, only : CheckShape

    ! input/output
    integer, intent(in)              :: im,jm,lm,ih,jh,lh,ntracer
    integer, intent(in)              :: send_id 
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracer) :: tracer_read
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1   :lmloc,ntracer) :: tracerloc   
    real, dimension(:,:,:,:), intent(in)    :: tracer_read
    real, dimension(:,:,:,:), intent(inout) :: tracerloc   

    ! local
    real,dimension(1-ih:im+ih,1-jh:jm+jh,ntracer) :: sum_one,sum_all
    integer                         :: slab2d
    integer,dimension(0:npes-1)   :: displ     
    integer,dimension(0:npes-1)   :: sendcounts
    integer                         :: recvcount
    integer                         ::nsum
    integer                         :: itracer, i 

    ! start

    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh,ntracer/), shape(tracer_read) )
    call CheckShape( (/im+2*ih,jm+2*jh,lmloc  ,ntracer/), shape(tracerloc  ) )

    slab2d    = (im+2*ih)*(jm+2*jh)
    !
    ! offset for gather-scatter
    !
    displ(:) = 0
    displ(0) = lh*slab2d
    do i=1,npes-1
       displ(i) = displ(i-1) + lmar(i-1)*slab2d
    end do
    !
    do i=0,npes-1
       sendcounts(i) = lmar(i)*slab2d
    end do
    !
    recvcount = lmloc*slab2d
    !
    ! call MPI_SCATTERV
    !
    do itracer = 1, ntracer
       !call MPI_SCATTERV( tracer_read(1-ih,1-jh,1-lh,itracer), sendcounts, displ, MY_REAL, & 
       !                   tracerloc  (1-ih,1-jh,1   ,itracer), recvcount , MY_REAL, &
       !                   send_id, localComm, ierr )
       call MPI_SCATTERV( tracer_read(:,:,:,itracer), sendcounts, displ, MY_REAL, & 
            tracerloc  (:,:,:,itracer), recvcount        , MY_REAL, &
            send_id, localComm, ierr )
    end do

  end subroutine scatter_after_read_k


  ! **********************************


  subroutine scatter_after_read_t(tracer_read,im,jm,lm,ih,jh,lh,tracerloc,send_id)

    use mpi_const,only : lmar,myid,lmloc,ntracetloc,ntracet_ar
    use dims, only : CheckShape

    !
    ! scatter an array
    !    tracer    (1-ih,im+ih,1-jh,jm+jh,1-lh,lm+lh,ntracer)
    ! on root (pe=0)
    ! over local arrays
    !    tracerloc(1-ih,im+ih,1-jh,jm+jh,1,  ,lmloc,ntracer) 
    ! on all PE's so that 
    ! the array is distributed over the l-index
    !
    ! input/output
    integer,intent(in)              :: im,jm,lm,ih,jh,lh
    integer,intent(in)              :: send_id 
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracet) :: tracer_read
    !real,dimension(1-ih:im+ih,1-jh:jm+jh,1-lh:lm+lh,ntracetloc) :: tracerloc   
    real, dimension(:,:,:,:)        :: tracer_read
    real, dimension(:,:,:,:)        :: tracerloc   

    ! local
    integer                         :: slab3d
    integer,dimension(0:npes-1)     :: displ     
    integer,dimension(0:npes-1)     :: sendcounts
    integer                         :: recvcount,i
    logical                         :: okdebug=.false.

    ! start

    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh,ntracet   /), shape(tracer_read) )
    call CheckShape( (/im+2*ih,jm+2*jh,lm+2*lh,ntracetloc/), shape(tracerloc  ) )

    slab3d    = (im+2*ih)*(jm+2*jh)*(lm+2*lh)
    !
    ! offset for gather-scatter
    !
    displ(:) = 0
    do i=1,npes-1
       displ(i) = displ(i-1) + ntracet_ar(i-1)*slab3d ! calculate offset slabs
    end do
    !
    do i=0,npes-1
       sendcounts(i) = ntracet_ar(i)*slab3d
    end do
    !
    recvcount = ntracetloc*slab3d
    !
    ! call MPI_SCATTERV
    !
    call barrier
    !call MPI_SCATTERV( tracer_read(1-ih,1-jh,1-lh,1), sendcounts, displ, MY_REAL, &
    !                   tracerloc  (1-ih,1-jh,1-lh,1), recvcount        , MY_REAL, &
    !                   send_id, localComm, ierr )
    call MPI_SCATTERV( tracer_read(:,:,:,:), sendcounts, displ, MY_REAL, &
         tracerloc  (:,:,:,:), recvcount        , MY_REAL, &
         send_id, localComm, ierr )
    call barrier
    !

  end subroutine scatter_after_read_t



  subroutine dump_field4d(region,fieldloc4d,lmr,ntracel,add,name,filename)
    !
    ! This subroutine is used to collect 3d arrays from all 
    ! the processors and write them to a HDF file. 
    ! The array add contains the number of halo cells in each dimension (3). 
    ! depending on the size of lmr, data is gathered and written 
    ! from root_k, or copied and written from root_t. 
    ! If lmloc==lm(region), then npes=0 and gathering is also not necessary.
    ! WP january 2003
    !
    use io_hdf,      only : io_write, DFACC_CREATE, DFNT_INT32
    use dims,        only : im, jm, lm
    use chem_param,  only : ntracet
    use dims       , only : CheckShape

    ! input/output
    character(len=*), intent(in)   :: name,filename
    integer, intent(in)            :: region,lmr,ntracel
    integer, dimension(3)          :: add  !WP! to expand arrays, only one or zero
    !real,dimension(& 
    !     1-add(1):im(region)+add(1),&
    !     1-add(2):jm(region)+add(2),&
    !     1:lmr,ntracel)   :: fieldloc4d
    real, dimension(:,:,:,:), intent(in)  ::  fieldloc4d

    ! local
    integer :: is,ie,js,je
    integer :: io,istat,sfstart,sfend,sfsnatt,root_id
    real,dimension(& 
         1-add(1):im(region)+add(1),&
         1-add(2):jm(region)+add(2),&
         1-add(3):lm(region)+add(3),ntracet) :: fieldglob4d

    ! start

    call CheckShape( (/im(region)+2*add(1),jm(region)+2*add(2),lmr,ntracel/), shape(fieldloc4d) )

    which_par=previous_par(region)
    fieldglob4d=0.0
    if ( which_par == 'levels' ) then ! parallel over levels
       call  gather_tracer_k(fieldglob4d,im(region),jm(region),lm(region), &
            add(1),add(2),add(3),ntracet, fieldloc4d,.false.)
       ! no need to broadcast result
       root_id=root_k                 
    else if(which_par=='tracer') then  !WP! parallel over tracers
       call  gather_tracer_t(fieldglob4d,im(region),jm(region),lm(region), &
            add(1),add(2),add(3),ntracet, fieldloc4d,.false.)
       ! no need to broadcast result
       root_id=root_t                 
    end if

    if ( myid == root_id ) then
       io=SFSTART(filename, DFACC_CREATE)
       istat = sfsnatt(io,'im',    DFNT_INT32, 1, im(region))
       istat = sfsnatt(io,'jm',    DFNT_INT32, 1, jm(region))
       istat = sfsnatt(io,'lm',    DFNT_INT32, 1, lm(region))
       istat = sfsnatt(io,'ntracet',    DFNT_INT32, 1, ntracet)
       call io_write( io, im(region), 'lon', jm(region), 'lat', &
            lm(region),'hybrid',ntracet,'ntracet',&
            fieldglob4d(&
            1:im(region),&
            1:jm(region),&
            1:lm(region),1:ntracet) ,name)
       print *,'dump_field4d: closing output file',SFend(io)
    end if

  end subroutine dump_field4d



  subroutine dump_field3d(region,fieldloc3d,lmr,add,name,filename)
    !
    ! This subroutine is used to collect 3d arrays from all the 
    ! processors and write them to a HDF file. The array add 
    ! contains the number of halo cells in each dimension (3). 
    ! depending on the size of lmr, data is gathered and written 
    ! from root_k, or copied and written from root_t. 
    ! If lmloc==lm(region), then npes=0 and gathering is also not necessary.
    !    WP january 2003
    !
    use io_hdf,      only : io_write, DFACC_CREATE, DFNT_INT32
    use dims,        only : im, jm, lm
    use dims,        only : CheckShape

    ! input/output
    character(len=*), intent(in)         :: name,filename
    integer, intent(in)                  :: region,lmr
    integer, dimension(3), intent(in)    :: add  !WP! to expand arrays, only one or zero
    !real,dimension(& 
    !     1-add(1):im(region)+add(1),&
    !     1-add(2):jm(region)+add(2),&
    !     1:lmr)   :: fieldloc3d  !WP! watch out for am!
    real,dimension(:,:,:), intent(in)    :: fieldloc3d  !WP! watch out for am!

    ! local
    integer :: is,ie,js,je
    integer :: io,istat,sfstart,sfend,sfsnatt,root_id
    real,dimension(& 
         1-add(1):im(region)+add(1),&
         1-add(2):jm(region)+add(2),&
         1-add(3):lm(region)+add(3)) :: fieldglob3d

    ! start

    call CheckShape( (/im(region)+2*add(1),jm(region)+2*add(2),lmr/), shape(fieldloc3d) )

    fieldglob3d=0.0
    if ( lmr == lmloc ) then
       call gather_tracer_k_3d( &
            fieldglob3d, im(region),jm(region),lm(region), add(1),add(2),add(3), &
            fieldloc3d, .false. )
       ! no need to broadcast result
       root_id=root_k                 
    else
       fieldglob3d=fieldloc3d
       root_id=root_t
    end if

    if ( myid == root_id ) then
       io=sfstart(filename, DFACC_CREATE)
       istat = sfsnatt(io,'im',    DFNT_INT32, 1, im(region))
       istat = sfsnatt(io,'jm',    DFNT_INT32, 1, jm(region))
       istat = sfsnatt(io,'lm',    DFNT_INT32, 1, lm(region))
       call io_write(io,im(region),'lon',jm(region),'lat', &
            lm(region),'hybrid',&
            fieldglob3d(&
            1:im(region),&
            1:jm(region),&
            1:lm(region)) ,name)
       print *,'dump_field3d: closing output file',sfend(io)
    end if

  end subroutine dump_field3d



  ! ***


  !  subroutine summasr( region, m, rm )
  !
  !    use dims, only : im, jm, lm
  !    use chem_param,   only : fscale
  !    
  !    ! --- in/out ---------------------------
  !    
  !    integer, intent(in)                ::  region
  !    real,dimension(:,:,:), pointer     ::  m
  !    real,dimension(:,:,:,:), pointer   ::  rm
  !    
  !    ! --- local ---------------------------
  !
  !    integer i, j, k,n,root_id
  !    real summas,summin,summax
  !    real summas_all,summin_all,summax_all
  !    real,allocatable,dimension(:,:,:) :: mglob
  !    real,allocatable,dimension(:,:,:,:) :: rmglob
  !
  !    ! --- begin ------------------------------
  !
  !    allocate(mglob(-1:im(region)+2,-1:jm(region)+2,1:lm(region))) 
  !    allocate(rmglob(-1:im(region)+2,-1:jm(region)+2,1:lm(region),ntracet)) 
  !
  !    mglob=0.0
  !    rmglob=0.0
  !
  !    if ( previous_par(region) == 'levels' ) then
  !       call gather_tracer_k_3d( &
  !               mglob, im(region),jm(region),lm(region), 2,2,0, &
  !               m, .false. ) ! no need to broadcast result
  !       call  gather_tracer_k( &
  !               rmglob, im(region),jm(region),lm(region), 2,2,0, ntracet, &
  !               rm, .false. ) ! no need to broadcast result
  !       root_id = root_k                 
  !    else
  !       mglob=m
  !       call gather_tracer_t( &
  !               rmglob, im(region),jm(region),lm(region), 2,2,0, ntracet, &
  !               rm, .false. ) ! no need to broadcast result
  !       root_id = root_t
  !    end if
  !    summin=1.e20
  !    summax=0.0
  !    if ( myid == root_id ) then
  !       summas=sum(mglob(1:im(region),1:jm(region),1:lm(region)) )
  !       do k=1,lm(region)
  !          do j=1,jm(region)
  !             do i=1,im(region)
  !                summin = min(summin,mglob(i,j,k))
  !                summax = max(summax,mglob(i,j,k))
  !             end do
  !          end do
  !       end do
  !       print *, ' Total mass            ',summas
  !       print *, ' Maximum value  m      ',summax
  !       print *, ' Minimum value  m      ',summin    
  !       summin=1.0
  !       summax=0.0
  !       do k=1,lm(region)
  !          do j=1,jm(region)
  !             do i=1,im(region)
  !                summin = min(summin,fscale(1)*rmglob(i,j,k,1)/mglob(i,j,k))
  !                summax = max(summax,fscale(1)*rmglob(i,j,k,1)/mglob(i,j,k))
  !             end do
  !          end do
  !       end do
  !       print *, ' Maximum value mixing ratio tracer 1 ',summax
  !       print *, ' Minimum value mixing ratio tracer 1 ',summin
  !    end if
  !    deallocate(mglob)
  !    deallocate(rmglob)
  !    call barrier
  !
  !  end subroutine summasr


end module mpi_comm