mpi_mod.f90

! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later

module mpi_mod

!*****************************************************************************
!                                                                            *
!  DESCRIPTION                                                               *
!    This module contains subroutines and common variables used for the      *
!    MPI parallelization of FLEXPART.                                        *
!                                                                            *
!  NOTE                                                                      *
!    Depending on the MPI library installed on your system (e.g. mpich2,     *
!    OpenMPI) you may need to choose below in this file between              *
!      use mpi                                                               *
!    (if the MPI library comes with the file 'mpi.mod'); or                  *
!      include 'mpif.h'                                                      *
!                                                                            *
!                                                                            *
!*****************************************************************************
!                                                                            *
! Variables:                                                                 *
!                                                                            *
! mp_ierr                 MPI error code                                     *
! mp_np                   Number of MPI processes                            *
! mp_pid                  Process ID of each MPI process                     *
! mp_seed                 Parameter for random number seed                   *
! read_grp_min            Minimum number of processes at which one will be   *
!                         used as reader                                     *
! numpart_mpi,            Number of particles per node                       *
! maxpart_mpi                                                                *
! mp_partid               MPI process ID for particle calculation            *
! mp_partgroup_           Refers to the subset of processors performing      *
!                         loops over particles. Will be all processes        *
!                         unless a dedicated process runs getfields/readwind *
! lmp_sync                If .false., use asynchronous MPI                   *
! mp_cp                   Real precision to use for deposition fields        *
!                                                                            *
!                                                                            *
!                                                                            *
!                                                                            *
!*****************************************************************************
!                                                                            *
! Modification by DJM, 2017-05-09 - added #ifdef USE_MPIINPLACE cpp          *
! directive to mpif_tm_reduce_grid() to insure that MPI_IN_PLACE is          *
! used in the MPI_Reduce() only if specifically compiled with that           *
! directive.                                                                 *
!                                                                            *
!*****************************************************************************

  use mpi 
  use par_mod, only: dp,sp
  use com_mod, only: lroot

  implicit none

!  include 'mpif.h'

  public

! Set aside a process for reading windfields if using at least these many processes
!==================================================
  integer, parameter, private :: read_grp_min=4
!==================================================

! Variables for each MPI process in the world group
  integer :: mp_ierr, mp_np, mp_pid, mp_partid
  integer, private :: world_group_id

! Variables for MPI processes in the 'particle' group
  integer, allocatable, dimension(:) :: mp_partgroup_rank
  integer, allocatable, dimension(:) :: npart_per_process
  integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np

  integer :: mp_seed=0
  integer, parameter :: mp_sp=MPI_REAL4, mp_dp=MPI_REAL8
  integer :: mp_cp
  integer, parameter :: id_root=0 ! master process

! MPI tags/requests for send/receive operation
  integer :: tm1
  integer, parameter :: nvar_async=26
!integer, dimension(:), allocatable :: tags
  integer, dimension(:), allocatable :: reqs

! Status array used for certain MPI operations (MPI_RECV)
  integer, dimension(MPI_STATUS_SIZE) :: mp_status


  integer :: id_read   ! readwind/getfield process
  integer :: numpart_mpi,maxpart_mpi ! number of particles per node
  integer :: tot_numpart=0
  integer :: mp_comm_used ! global or subgroup communicator

  logical :: lmpreader=.false. ! is set to true for reading process(es) only.
  logical :: lmp_use_reader=.false. ! true if separate readwind process is used

! .true. if only using synchronous MPI send/recv (default)
! If setting this to .false., numwfmem must be set to 3
!===============================================================================
  logical :: lmp_sync=.true. 
!===============================================================================

! mp_dbg_mode       Used for debugging MPI.
! mp_dev_mode       various checks related to debugging the parallel code
! mp_dbg_out        write some arrays to data file for debugging
! mp_measure_time   Measure cpu/wall time, write out at end of run
! mp_time_barrier   Measure MPI barrier time
! mp_exact_numpart  Use an extra MPI communication to give the exact number of particles
!                   to standard output (this does not otherwise affect the simulation)
  logical, parameter :: mp_dbg_mode = .false.
  logical, parameter :: mp_dev_mode = .false.
  logical, parameter :: mp_dbg_out = .false.
  logical, parameter :: mp_time_barrier=.true.
  logical, parameter :: mp_measure_time=.false.
  logical, parameter :: mp_exact_numpart=.true.

! for measuring CPU/Wall time
  real(sp),private :: mp_comm_time_beg, mp_comm_time_end, mp_comm_time_total=0.
  real(dp),private :: mp_comm_wtime_beg, mp_comm_wtime_end, mp_comm_wtime_total=0.
  real(sp),private :: mp_root_time_beg, mp_root_time_end, mp_root_time_total=0.
  real(dp),private :: mp_root_wtime_beg, mp_root_wtime_end, mp_root_wtime_total=0.
  real(sp),private :: mp_barrier_time_beg, mp_barrier_time_end, mp_barrier_time_total=0.
  real(dp),private :: mp_barrier_wtime_beg, mp_barrier_wtime_end, mp_barrier_wtime_total=0.
  real(sp),private :: tm_nploop_beg, tm_nploop_end, tm_nploop_total=0.
  real(sp),private :: tm_tot_beg, tm_tot_end, tm_tot_total=0.
  real(dp),private :: mp_getfields_wtime_beg, mp_getfields_wtime_end, mp_getfields_wtime_total=0.
  real(sp),private :: mp_getfields_time_beg, mp_getfields_time_end, mp_getfields_time_total=0.
  real(dp),private :: mp_readwind_wtime_beg, mp_readwind_wtime_end, mp_readwind_wtime_total=0.
  real(sp),private :: mp_readwind_time_beg, mp_readwind_time_end, mp_readwind_time_total=0.
  real(dp),private :: mp_io_wtime_beg, mp_io_wtime_end, mp_io_wtime_total=0.
  real(sp),private :: mp_io_time_beg, mp_io_time_end, mp_io_time_total=0.
  real(dp),private :: mp_wetdepo_wtime_beg, mp_wetdepo_wtime_end, mp_wetdepo_wtime_total=0.
  real(sp),private :: mp_wetdepo_time_beg, mp_wetdepo_time_end, mp_wetdepo_time_total=0.
  real(dp),private :: mp_advance_wtime_beg, mp_advance_wtime_end, mp_advance_wtime_total=0.
  real(dp),private :: mp_conccalc_time_beg, mp_conccalc_time_end, mp_conccalc_time_total=0.
  real(dp),private :: mp_total_wtime_beg, mp_total_wtime_end, mp_total_wtime_total=0.
  real(dp),private :: mp_vt_wtime_beg, mp_vt_wtime_end, mp_vt_wtime_total
  real(sp),private :: mp_vt_time_beg, mp_vt_time_end, mp_vt_time_total

! dat_lun           logical unit number for i/o
  integer, private :: dat_lun 

! Temporary arrays for particles (allocated and deallocated as needed)
  integer, allocatable, dimension(:) :: nclass_tmp, npoint_tmp, itra1_tmp, idt_tmp, &
       & itramem_tmp, itrasplit_tmp
  real(kind=dp), allocatable, dimension(:) :: xtra1_tmp, ytra1_tmp
  real, allocatable, dimension(:) :: ztra1_tmp 
  real, allocatable, dimension(:,:) :: xmass1_tmp

! mp_redist_fract        Exchange particles between processes if relative difference
!                        is larger. A good value is between 0.0 and 0.5
! mp_maxpart_factor      Allocate more memory per process than strictly needed
!                        (safety factor). Recommended value between 1.5 and 2.5
! mp_min_redist          Do not redistribute particles if below this limit
  real, parameter :: mp_redist_fract=0.2, mp_maxpart_factor=1.5
  integer,parameter :: mp_min_redist=100000


contains

  subroutine mpif_init
!***********************************************************************
! DESCRIPTION
!   Initialize MPI.
!   
!   Create the global communicator MPI_COMM_WORLD
!   If using a separate MPI process for getfields/readwind, a subgroup
!   is created for the other processes.
!
! VARIABLES
!   mpi_mode    default 0, set to 2/3 if running MPI version
!   mp_np       number of running processes, decided at run-time
!***********************************************************************
    use par_mod, only: maxpart, numwfmem, dep_prec, maxreceptor, maxspec
    use com_mod, only: mpi_mode, verbosity, creceptor0

    implicit none

    integer :: i,j,s,addmaxpart=0

! Each process gets an ID (mp_pid) in the range 0,..,mp_np-1
    call MPI_INIT(mp_ierr)
    if (mp_ierr /= 0) goto 100
    call MPI_COMM_RANK(MPI_COMM_WORLD, mp_pid, mp_ierr)
    if (mp_ierr /= 0) goto 100
    call MPI_COMM_SIZE(MPI_COMM_WORLD, mp_np, mp_ierr)
    if (mp_ierr /= 0) goto 100


! Variable mpi_mode is used to handle subroutines common to parallel/serial version
    if (lmp_sync) then
      mpi_mode=2 ! hold 2 windfields in memory
    else
      mpi_mode=3 ! hold 3 windfields in memory
    end if

    if (mp_pid.ne.0) then
      lroot = .false.
    end if

! Set MPI precision to use for transferring deposition fields
!************************************************************
    if (dep_prec==dp) then
      mp_cp = MPI_REAL8
! TODO: write info message for serial version as well 
      if (lroot.and.verbosity>0) write(*,*) 'Using double precision for deposition fields'
    else if (dep_prec==sp) then
      mp_cp = MPI_REAL4
      if (lroot.and.verbosity>0) write(*,*) 'Using single precision for deposition fields'
    else
      write(*,*) 'ERROR: something went wrong setting MPI real precision'
      stop
    end if

! Check for sensible combination of parameters
!*********************************************

    if (.not.lmp_sync.and.numwfmem.ne.3) then
      if (lroot) then
        write(*,FMT='(80("#"))')
        write(*,*) '#### mpi_mod::mpif_init> ERROR: ', &
             & 'numwfmem must be set to 3 for asyncronous reading ####'
        write(*,FMT='(80("#"))')
      end if
      call MPI_FINALIZE(mp_ierr)
      stop
    else if (lmp_sync.and.numwfmem.ne.2.and.lroot) then
      write(*,FMT='(80("#"))')
      write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
           & 'numwfmem should be set to 2 for syncronous'
      write(*,*) ' reading. Results will still be valid, but unneccesary memory &
           &is allocated.'
      write(*,FMT='(80("#"))')
! Force "syncronized" version if all processes will call getfields
    else if ((.not.lmp_sync.and.mp_np.lt.read_grp_min).or.(mp_np.eq.1)) then
      if (lroot) then
        write(*,FMT='(80("#"))')
        write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
             & 'all procs call getfields. Setting lmp_sync=.true.'
        write(*,FMT='(80("#"))')
      end if
      lmp_sync=.true.
    end if

! TODO: Add more warnings for unimplemented flexpart features

! Set ID of process that calls getfield/readwind. 
! Using the last process in the group ensures statistical identical results
! as running with one process less but not using separate read process
!**********************************************************************

!    id_read = min(mp_np-1, 1)
    id_read = mp_np-1

    if (mp_pid.eq.id_read) lmpreader=.true.

    call MPI_Comm_group (MPI_COMM_WORLD, world_group_id, mp_ierr)

! Create a MPI group/communicator that will calculate trajectories. 
! Skip this step if program is run with only a few processes
!************************************************************************
    allocate(mp_partgroup_rank(0:mp_np-2))

! This allows checking for allocation of arrays when no subgroup is used
    mp_partgroup_pid=0

    if (read_grp_min.lt.2) then
      write(*,*) '#### mpi_mod::mpif_init> ERROR ####', &
           & 'read_grp_min should be at least 2. Exiting'
      stop
    end if

    if (mp_np.ge.read_grp_min) then
      lmp_use_reader = .true.

! Map the subgroup IDs= 0,1,2,...,mp_np-2, skipping reader process
      j=-1 
      do i=0, mp_np-2 !loop over all (subgroup) IDs
        j=j+1
        if (i.eq.id_read) then
          j=j+1
        end if
        mp_partgroup_rank(i) = j
      end do

      call MPI_Group_incl (world_group_id, mp_np-1, mp_partgroup_rank, &
           &mp_partgroup_pid, mp_ierr)
      if (mp_ierr /= 0) goto 100
      call MPI_Comm_create (MPI_COMM_WORLD, mp_partgroup_pid, mp_partgroup_comm, mp_ierr)
      if (mp_ierr /= 0) goto 100

      if (mp_pid.ne.id_read) then
        call MPI_Comm_rank (mp_partgroup_comm, mp_partgroup_pid, mp_ierr)
        if (mp_ierr /= 0) goto 100

! Get size of new group (=mp_np-1)
        call MPI_COMM_SIZE(mp_partgroup_comm, mp_partgroup_np, mp_ierr)
        if (mp_ierr /= 0) goto 100
        if (mp_partgroup_np.ne.mp_np-1) then
          write(*,*)  '#### mpi_mod:: mpif_init> ERROR ####&
               & mp_partgroup_np.ne.mp_np-1'
          stop
        endif

      else
        mp_partgroup_pid = -1
      end if
    end if

    if (lmp_use_reader) then
      mp_comm_used = mp_partgroup_comm
      mp_partid = mp_partgroup_pid
      mp_partgroup_np=mp_np-1
    else
      mp_comm_used = MPI_COMM_WORLD
      mp_partgroup_np = mp_np
      mp_partid = mp_pid
    end if

! Set maxpart per process
! ESO 08/2016: Increase maxpart per process, in case of unbalanced distribution
    maxpart_mpi=int(mp_maxpart_factor*real(maxpart)/real(mp_partgroup_np))
    if (mp_np == 1) maxpart_mpi = maxpart

! Do not allocate particle data arrays for readwind process
    if (lmpreader.and.lmp_use_reader) then
      maxpart_mpi=0
    end if

! Set random seed for each non-root process
    if (mp_pid.gt.0) then
      s = 244
      mp_seed = -abs(mod((s*181)*((mp_pid-83)*359), 104729))
    end if
    if (mp_dbg_mode) then
      mp_seed=0
      if (lroot) write(*,*) 'MPI: setting seed=0'
    end if

! Allocate request array for asynchronous MPI
    if (.not.lmp_sync) then
      allocate(reqs(0:nvar_async*mp_np-1))
      reqs(:)=MPI_REQUEST_NULL
    else
      allocate(reqs(0:1))
      reqs(:)=MPI_REQUEST_NULL
    end if

! Allocate array for number of particles per process    
    allocate(npart_per_process(0:mp_partgroup_np-1))

! Write whether MPI_IN_PLACE is used or not
#ifdef USE_MPIINPLACE
    if (lroot) write(*,*) 'Using MPI_IN_PLACE operations'
#else
    if (lroot) allocate(creceptor0(maxreceptor,maxspec))
    if (lroot) write(*,*) 'Not using MPI_IN_PLACE operations'
#endif
    goto 101

100 write(*,*) '#### mpi_mod::mpif_init> ERROR ####', mp_ierr
    stop

101 end subroutine mpif_init


  subroutine mpif_mpi_barrier
!***********************************************************************
! Set MPI barrier (all processes are synchronized here), with the
! possible exception of a separate 'readwind' process.
! Optionally measure cpu/wall time.
!
!***********************************************************************
    implicit none

    if (mp_time_barrier) then
      call cpu_time(mp_barrier_time_beg)
      mp_barrier_wtime_beg = mpi_wtime()
    endif

    call MPI_BARRIER(mp_comm_used, mp_ierr)
    if (mp_ierr /= 0) goto 100

    if (mp_time_barrier) then
      call cpu_time(mp_barrier_time_end)
      mp_barrier_wtime_end = mpi_wtime()

      mp_barrier_time_total=mp_barrier_time_total+&
           &(mp_barrier_time_end-mp_barrier_time_beg)
      mp_barrier_wtime_total=mp_barrier_wtime_total+&
           &(mp_barrier_wtime_end-mp_barrier_wtime_beg)
    end if

    goto 101

100 write(*,*) '#### mpi_mod::mpif_mpi_barrier> ERROR ####', mp_ierr
    stop

101 end subroutine mpif_mpi_barrier


  subroutine mpif_com_mod_allocate
!*******************************************************************************    
! Dynamic allocation of arrays (in serial code these have size maxpart)
!
! Not in use anymore, moved to com_mod for interoperability with serial version
! 
! TODO: remove 
!*******************************************************************************
    use com_mod

    implicit none 

    integer :: array_size

    array_size = maxpart_mpi

    allocate(itra1(array_size),npoint(array_size),&
         & nclass(array_size),idt(array_size),&
         & itramem(array_size),itrasplit(array_size),&
         & xtra1(array_size),ytra1(array_size),&
         & ztra1(array_size),xmass1(array_size, maxspec))

    allocate(uap(array_size),ucp(array_size),&
         & uzp(array_size),us(array_size),&
         & vs(array_size),&
         & ws(array_size),cbt(array_size))


  end subroutine mpif_com_mod_allocate


  subroutine mpif_tm_send_dims
!***********************************************************************
! Distribute array dimensions from pid0 to all processes.
! numpart_mpi/numpart is sent to allow dynamic allocation
!
! Currently not in use.
!
! Variables of similar type (integer, real) are placed in an array
! and sent collectively, to avoid the overhead associated with individual
! MPI send/receives
!
!
!***********************************************************************
    use com_mod

    implicit none

    integer :: add_part=0

    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)

! MPI subgroup does the particle-loop 
    if (lmp_use_reader) then
      if (mod(numpart,mp_partgroup_np).ne.0) add_part=1
      numpart_mpi=int(numpart/mp_partgroup_np)+add_part
    else
      if (mod(numpart,mp_np).ne.0) add_part=1
      numpart_mpi=int(numpart/mp_np)+add_part
    end if


! redefine numpart as 'numpart per process' throughout the code
!**************************************************************
    numpart = numpart_mpi

  end subroutine mpif_tm_send_dims


  subroutine mpif_send_part_properties(num_part)
!***********************************************************************
! Distribute particle properties from root to all processes.
!  
! Used by init_domainfill_mpi
!
! Variables:
! 
! num_part        input, number of particles per process (rounded up)
!
!***********************************************************************
    use com_mod

    implicit none

    integer,intent(in) :: num_part
    integer :: i,jj, addone

! Exit if too many particles
    if (num_part.gt.maxpart_mpi) then
      write(*,*) '#####################################################'
      write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED  ####'
      write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE  ####'
      write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
      write(*,*) '#### OR INCREASE MAXPART.                        ####'
      write(*,*) '#####################################################'
!      call MPI_FINALIZE(mp_ierr)
      stop
    end if


! Time for MPI communications
!****************************
    if (mp_measure_time) call mpif_mtime('commtime',0)

! Distribute variables, send from pid 0 to other processes (including itself)
!****************************************************************************

    call MPI_SCATTER(nclass_tmp,num_part,MPI_INTEGER,nclass,&
         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(npoint_tmp,num_part,MPI_INTEGER,npoint,&
         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(itra1_tmp,num_part,MPI_INTEGER,itra1,&
         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(idt_tmp,num_part,MPI_INTEGER,idt,&
         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(itramem_tmp,num_part,MPI_INTEGER,itramem,&
         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(itrasplit_tmp,num_part,MPI_INTEGER,itrasplit,&
         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(xtra1_tmp,num_part,mp_dp,xtra1,&
         &num_part,mp_dp,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(ytra1_tmp,num_part,mp_dp,ytra1,&
         &num_part,mp_dp,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_SCATTER(ztra1_tmp,num_part,mp_sp,ztra1,&
         &num_part,mp_sp,id_root,mp_comm_used,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    do i=1,nspec
      call MPI_SCATTER(xmass1_tmp(:,i),num_part,mp_sp,xmass1(:,i),&
           &num_part,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
    end do

    if (mp_measure_time) call mpif_mtime('commtime',1)

    goto 601

600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
    stop

! After the transfer of particles it it possible that the value of
! "numpart" is larger than the actual used, so we reduce it if there are
! invalid particles at the end of the arrays

601 do i=numpart, 1, -1
      if (itra1(i).eq.-999999999) then
        numpart=numpart-1
      else
        exit
      end if
    end do


  end subroutine mpif_send_part_properties


  subroutine mpif_calculate_part_redist(itime)
!***********************************************************************
! Calculate number of particles to redistribute between processes. This routine 
! can be called at regular time intervals to keep a level number of
! particles on each process.
!
! First, all processes report their local "numpart" to each other, which is
! stored in array "numpart_mpi(np)". This array is sorted from low to
! high values, along with a corresponding process ID array "idx_arr(np)".
! If the relative difference between the highest and lowest "numpart_mpi" value
! is above a threshold, particles are transferred from process idx_arr(np-1)
! to process idx_arr(0). Repeat for processes idx_arr(np-i) and idx_arr(i)
! for all valid i. 
! Note: If np is an odd number, the process with the median
! number of particles is left unchanged
!
! VARIABLES
! itime        input, current time
!***********************************************************************
    use com_mod

    implicit none
    
    integer, intent(in) :: itime
    real :: pmin,z
    integer :: i,jj,nn, num_part=1,m,imin, num_trans
    logical :: first_iter
    integer,allocatable,dimension(:) :: idx_arr
    real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this

! Exit if running with only 1 process
!************************************************************************
    if (mp_np.eq.1) return

! All processes exchange information on number of particles
!****************************************************************************
    allocate( idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))

    call MPI_Allgather(numpart, 1, MPI_INTEGER, npart_per_process, &
         & 1, MPI_INTEGER, mp_comm_used, mp_ierr)


! Sort from lowest to highest
! Initial guess: correct order
    sorted(:) = npart_per_process(:)
    do i=0, mp_partgroup_np-1
      idx_arr(i) = i
    end do

! Do not rebalance particles for ipout=3    
    if (ipout.eq.3) return

! For each successive element in index array, see if a lower value exists
    do i=0, mp_partgroup_np-2
      pmin=sorted(i)
      imin=idx_arr(i)
      m=i+1
      do jj=m, mp_partgroup_np-1
        if (pmin.le.sorted(jj)) cycle
        z=pmin
        pmin=sorted(jj)
        sorted(jj)=z

        nn=imin
        imin=idx_arr(jj)
        idx_arr(jj)=nn
      end do
      sorted(i)=pmin
      idx_arr(i)=imin
    end do

! For each pair of processes, transfer particles if the difference is above a
! limit, and numpart at sending process large enough

    m=mp_partgroup_np-1 ! index for last sorted process (most particles)
    do i=0,mp_partgroup_np/2-1
      num_trans = npart_per_process(idx_arr(m)) - npart_per_process(idx_arr(i))
      if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then
        if ( npart_per_process(idx_arr(m)).gt.mp_min_redist.and.&
             & real(num_trans)/real(npart_per_process(idx_arr(m))).gt.mp_redist_fract) then
! DBG
          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process', &
          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, npart_per_process
! DBG
          call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
        end if
      end if
      m=m-1
    end do

    deallocate(idx_arr, sorted)

  end subroutine mpif_calculate_part_redist


  subroutine mpif_redist_part(itime, src_proc, dest_proc, num_trans)
!***********************************************************************
! Transfer particle properties between two arbitrary processes.
! 
! VARIABLES
!
! itime           input, current time
! src_proc        input, ID of source process
! dest_proc       input, ID of destination process
! num_trans       input, number of particles to transfer
!
!************************************************************************
    use com_mod !TODO:  ,only: nclass etc

    implicit none 

    integer, intent(in) :: itime, src_proc, dest_proc, num_trans
    integer :: ll, ul ! lower and upper indices in arrays
    integer :: arr_sz ! size of temporary arrays
    integer :: mtag   ! MPI message tag
    integer :: i, j, minpart, ipart, maxnumpart
 
! Check for invalid input arguments
!**********************************
 if (src_proc.eq.dest_proc) then
   write(*,*) '<mpi_mod::mpif_redist_part>: Error: &
        &src_proc.eq.dest_proc' 
   stop
 end if

! Measure time for MPI communications
!************************************
    if (mp_measure_time) call mpif_mtime('commtime',0)

! Send particles
!***************
    if (mp_partid.eq.src_proc) then
      mtag = 2000

! Array indices for the transferred particles
      ll=numpart-num_trans+1
      ul=numpart

      if (mp_dev_mode) then
        write(*,FMT='(72("#"))')
        write(*,*) "Sending ", num_trans, "particles (from/to)", src_proc, dest_proc
        write(*,*) "numpart before: ", numpart
      end if

      call MPI_SEND(nclass(ll:ul),num_trans,&
           & MPI_INTEGER,dest_proc,mtag+1,mp_comm_used,mp_ierr)

      call MPI_SEND(npoint(ll:ul),num_trans,&
           & MPI_INTEGER,dest_proc,mtag+2,mp_comm_used,mp_ierr)

      call MPI_SEND(itra1(ll:ul),num_trans, &
           & MPI_INTEGER,dest_proc,mtag+3,mp_comm_used,mp_ierr)

      call MPI_SEND(idt(ll:ul),num_trans, &
           & MPI_INTEGER,dest_proc,mtag+4,mp_comm_used,mp_ierr)

      call MPI_SEND(itramem(ll:ul),num_trans, &
           & MPI_INTEGER,dest_proc,mtag+5,mp_comm_used,mp_ierr)

      call MPI_SEND(itrasplit(ll:ul),num_trans,&
           & MPI_INTEGER,dest_proc,mtag+6,mp_comm_used,mp_ierr)

      call MPI_SEND(xtra1(ll:ul),num_trans, &
           & mp_dp,dest_proc,mtag+7,mp_comm_used,mp_ierr)

      call MPI_SEND(ytra1(ll:ul),num_trans,&
           & mp_dp,dest_proc,mtag+8,mp_comm_used,mp_ierr)

      call MPI_SEND(ztra1(ll:ul),num_trans,&
           & mp_sp,dest_proc,mtag+9,mp_comm_used,mp_ierr)

      do j=1,nspec
        call MPI_SEND(xmass1(ll:ul,j),num_trans,&
             & mp_sp,dest_proc,mtag+(9+j),mp_comm_used,mp_ierr)
      end do

! Terminate transferred particles, update numpart
      itra1(ll:ul) = -999999999

      numpart = numpart-num_trans

      if (mp_dev_mode) then
        write(*,*) "numpart after: ", numpart
        write(*,FMT='(72("#"))')
      end if

    else if (mp_partid.eq.dest_proc) then

! Receive particles
!******************
      mtag = 2000
      ! Array indices for the transferred particles
      ll=numpart+1
      ul=numpart+num_trans

      if (mp_dev_mode) then
        write(*,FMT='(72("#"))')
        write(*,*) "Receiving ", num_trans, "particles (from/to)", src_proc, dest_proc
        write(*,*) "numpart before: ", numpart
      end if

! We could receive the data directly at nclass(ll:ul) etc., but this leaves
! vacant spaces at lower indices. Using temporary arrays instead.
      arr_sz = ul-ll+1
      allocate(itra1_tmp(arr_sz),npoint_tmp(arr_sz),nclass_tmp(arr_sz),&
           & idt_tmp(arr_sz),itramem_tmp(arr_sz),itrasplit_tmp(arr_sz),&
           & xtra1_tmp(arr_sz),ytra1_tmp(arr_sz),ztra1_tmp(arr_sz),&
           & xmass1_tmp(arr_sz, maxspec))
      
      call MPI_RECV(nclass_tmp,num_trans,&
           & MPI_INTEGER,src_proc,mtag+1,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(npoint_tmp,num_trans,&
           & MPI_INTEGER,src_proc,mtag+2,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(itra1_tmp,num_trans, &
           & MPI_INTEGER,src_proc,mtag+3,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(idt_tmp,num_trans, &
           & MPI_INTEGER,src_proc,mtag+4,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(itramem_tmp,num_trans, &
           & MPI_INTEGER,src_proc,mtag+5,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(itrasplit_tmp,num_trans,&
           & MPI_INTEGER,src_proc,mtag+6,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(xtra1_tmp,num_trans, &
           & mp_dp,src_proc,mtag+7,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(ytra1_tmp,num_trans,&
           & mp_dp,src_proc,mtag+8,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(ztra1_tmp,num_trans,&
           & mp_sp,src_proc,mtag+9,mp_comm_used,mp_status,mp_ierr)

      do j=1,nspec
        call MPI_RECV(xmass1_tmp(:,j),num_trans,&
             & mp_sp,src_proc,mtag+(9+j),mp_comm_used,mp_status,mp_ierr)
      end do

! This is the maximum value numpart can possibly have after the transfer
      maxnumpart=numpart+num_trans

! Search for vacant space and copy from temporary storage
!********************************************************
      minpart=1
      do i=1, num_trans
! Take into acount that we may have transferred invalid particles
        if (itra1_tmp(i).ne.itime) cycle
        do ipart=minpart,maxnumpart
          if (itra1(ipart).ne.itime) then
            itra1(ipart) = itra1_tmp(i)
            npoint(ipart) = npoint_tmp(i)
            nclass(ipart) = nclass_tmp(i)
            idt(ipart) = idt_tmp(i)
            itramem(ipart) = itramem_tmp(i)
            itrasplit(ipart) = itrasplit_tmp(i)
            xtra1(ipart) = xtra1_tmp(i)
            ytra1(ipart) = ytra1_tmp(i)
            ztra1(ipart) = ztra1_tmp(i)
            xmass1(ipart,:) = xmass1_tmp(i,:)
            goto 200 ! Storage space has been found, stop searching
          end if
! :TODO: add check for if too many particles requiried
        end do
200     minpart=ipart+1
      end do
! Increase numpart, if necessary
      numpart=max(numpart,ipart)

      deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,&
           &itrasplit_tmp,xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)

      if (mp_dev_mode) then
        write(*,*) "numpart after: ", numpart
        write(*,FMT='(72("#"))')
      end if

    else
! This routine should only be called by the two participating processes
      write(*,*) "ERROR: wrong process has entered mpi_mod::mpif_redist_part"
      stop
      return
    end if

! Measure time for MPI communications
!************************************
    if (mp_measure_time) call mpif_mtime('commtime',1)

  end subroutine mpif_redist_part


  subroutine mpif_tm_send_vars
!***********************************************************************
! Distribute particle variables from pid0 to all processes.
! Called from timemanager
! *NOT IN USE* at the moment, but can be useful for debugging
!
!***********************************************************************
    use com_mod

    implicit none

    integer :: i

! Time for MPI communications
!****************************
    if (mp_measure_time) call mpif_mtime('commtime',0)

! Distribute variables, send from pid 0 to other processes
!**********************************************************************

! integers
    if (lroot) then
      call MPI_SCATTER(npoint,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(idt,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(itra1,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(nclass,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(itramem,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 

! int2
      call MPI_SCATTER(cbt,numpart_mpi,MPI_INTEGER2,MPI_IN_PLACE,&
           & numpart_mpi,MPI_INTEGER2,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600

! real
      call MPI_SCATTER(uap,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(ucp,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(uzp,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600

      call MPI_SCATTER(us,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(vs,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(ws,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600

      call MPI_SCATTER(xtra1,numpart_mpi,mp_dp,MPI_IN_PLACE,&
           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(ytra1,numpart_mpi,mp_dp,MPI_IN_PLACE,&
           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(ztra1,numpart_mpi,mp_sp,MPI_IN_PLACE,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 

      do i=1,nspec
        call MPI_SCATTER(xmass1(:,i),numpart_mpi,mp_sp,MPI_IN_PLACE,&
             & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
        if (mp_ierr /= 0) goto 600 
      end do

    else ! (mp_pid >= 1)
! integers
      call MPI_SCATTER(npoint,numpart_mpi,MPI_INTEGER,npoint,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(idt,numpart_mpi,MPI_INTEGER,idt,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(itra1,numpart_mpi,MPI_INTEGER,itra1,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(nclass,numpart_mpi,MPI_INTEGER,nclass,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(itramem,numpart_mpi,MPI_INTEGER,itramem,&
           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 

! int2
      call MPI_SCATTER(cbt,numpart_mpi,MPI_INTEGER2,cbt,&
           & numpart_mpi,MPI_INTEGER2,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600

! reals
      call MPI_SCATTER(uap,numpart_mpi,mp_sp,uap,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(ucp,numpart_mpi,mp_sp,ucp,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(uzp,numpart_mpi,mp_sp,uzp,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600

      call MPI_SCATTER(us,numpart_mpi,mp_sp,us,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(vs,numpart_mpi,mp_sp,vs,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_SCATTER(ws,numpart_mpi,mp_sp,ws,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600

      call MPI_SCATTER(xtra1,numpart_mpi,mp_dp,xtra1,&
           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(ytra1,numpart_mpi,mp_dp,ytra1,&
           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_SCATTER(ztra1,numpart_mpi,mp_sp,ztra1,&
           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
      if (mp_ierr /= 0) goto 600 

      do i=1,nspec
        call MPI_SCATTER(xmass1(:,i),numpart_mpi,mp_sp,xmass1(:,i),&
             & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
        if (mp_ierr /= 0) goto 600 
      end do

    end if !(lroot)

    if (mp_measure_time) call mpif_mtime('commtime',1)

    goto 601

600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
    stop
601 end subroutine mpif_tm_send_vars


  subroutine mpif_tm_collect_vars
!*******************************************************************************
! Collect particle data to PID 0 from all processes.                           *
!                                                                              *
! This can be called from timemanager_mpi, before caclulating                  *
! concentrations/writing output grids.                                         *
!                                                                              *
! Currently *not in use*, as each process calculates concentrations            *
! separately, but can be useful for debugging                                  *
!                                                                              *
! To use this routine (together with mpif_tm_send_vars) to transfer data       *
! to the MPI root process (useful for debugging), insert code like this        *
! (in timemanager_mpi.f90)                                                     *
!                                                                              *
!      if (lroot) tot_numpart=numpart ! root stores total numpart              *
!      call mpif_tm_send_dims                                                  *
!      if (numpart>1) then                                                     *
!        call mpif_tm_send_vars                                                *
!      end if                                                                  *
!                                                                              *
!    To collect data from all processes to root process after a parallel       *
!    region:                                                                   *
!                                                                              *
!      if (numpart>0) then                                                     *
!          if (lroot) numpart = tot_numpart                                    *
!       call mpif_tm_collect_vars                                              *
!      end if                                                                  *
!                                                                              *
! Then a section that begins with "if (lroot) ..." will behave like            *
! serial flexpart                                                              *
!                                                                              *
!*******************************************************************************
    use com_mod !, only: numpart, nspec, itra1, npoint, nclass

    implicit none

    integer :: i

    logical :: use_mp_vars = .false.! use mp_... type allocated vars
    logical :: use_in_place = .true.! use MPI_IN_PLACE to collect vars.


! Time for MPI communications
    if (mp_measure_time) call mpif_mtime('commtime',0)

! Distribute variables, send from pid 0 to other processes
!**********************************************************************

    if (.not. use_mp_vars.and..not.use_in_place) then

! Integers:
      call MPI_GATHER(npoint, numpart_mpi, MPI_INTEGER, npoint, &
           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(idt, numpart_mpi, MPI_INTEGER, idt, &
           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(itra1, numpart_mpi, MPI_INTEGER, itra1, &
           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(nclass, numpart_mpi, MPI_INTEGER, nclass, &
           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(itramem, numpart_mpi, MPI_INTEGER, itramem, &
           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 

! int2
      call MPI_GATHER(cbt, numpart_mpi, MPI_INTEGER2, cbt, &
           &numpart_mpi, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 

! Reals:
      call MPI_GATHER(uap, numpart_mpi, mp_sp, uap, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(ucp, numpart_mpi, mp_sp, ucp, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(uzp, numpart_mpi, mp_sp, uzp, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 

      call MPI_GATHER(us, numpart_mpi, mp_sp, us, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600
      call MPI_GATHER(vs, numpart_mpi, mp_sp, vs, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(ws, numpart_mpi, mp_sp, ws, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600


      call MPI_GATHER(xtra1, numpart_mpi, mp_dp, xtra1, &
           &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(ytra1, numpart_mpi, mp_dp, ytra1, &
           &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_GATHER(ztra1, numpart_mpi, mp_sp, ztra1, &
           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
      if (mp_ierr /= 0) goto 600 

      do i=1, nspec
        call MPI_GATHER(xmass1(:,i), numpart_mpi, mp_sp, xmass1(:,i), &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

      end do

! Use variables npoint etc directly for communications
!***********************************************************************
    else if (use_in_place.and..not.use_mp_vars) then
      if (lroot) then
! Integers:
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, npoint, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, idt, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, itra1, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, nclass, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, itramem, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

! int2
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER2, cbt, &
             &numpart_mpi, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

! Reals:
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, uap, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, ucp, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, uzp, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, us, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, vs, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, ws, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600


        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_dp, xtra1, &
             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_dp, ytra1, &
             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, ztra1, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

        do i=1, nspec
          call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, xmass1(:,i), &
               &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
          if (mp_ierr /= 0) goto 600 
        end do

      else ! (for mp_pid >= 1)

! Integers:
        call MPI_GATHER(npoint, numpart_mpi, MPI_INTEGER, npoint, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(idt, numpart_mpi, MPI_INTEGER, idt, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(itra1, numpart_mpi, MPI_INTEGER, itra1, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(nclass, numpart_mpi, MPI_INTEGER, nclass, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(itramem, numpart_mpi, MPI_INTEGER, itramem, &
             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

! int2
        call MPI_GATHER(cbt, numpart_mpi, MPI_INTEGER2, cbt, &
             &numpart_mpi, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

! Reals:
        call MPI_GATHER(uap, numpart_mpi, mp_sp, uap, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(ucp, numpart_mpi, mp_sp, ucp, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(uzp, numpart_mpi, mp_sp, uzp, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

        call MPI_GATHER(us, numpart_mpi, mp_sp, us, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600
        call MPI_GATHER(vs, numpart_mpi, mp_sp, vs, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(ws, numpart_mpi, mp_sp, ws, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600


        call MPI_GATHER(xtra1, numpart_mpi, mp_dp, xtra1, &
             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(ytra1, numpart_mpi, mp_dp, ytra1, &
             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 
        call MPI_GATHER(ztra1, numpart_mpi, mp_sp, ztra1, &
             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
        if (mp_ierr /= 0) goto 600 

        do i=1, nspec
          call MPI_GATHER(xmass1(:,i), numpart_mpi, mp_sp, xmass1(:,i), &
               &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
          if (mp_ierr /= 0) goto 600 
        end do
      end if ! (lroot)
    end if !  (use_in_place)

    if (mp_measure_time) call mpif_mtime('commtime',1)

    goto 601

600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
    stop
601 end subroutine mpif_tm_collect_vars


  subroutine mpif_gf_send_vars(memstat)
!*******************************************************************************
! DESCRIPTION
!   Distribute 'getfield' variables from reader process
! 
!   Called from timemanager
!
! NOTE
!   This subroutine distributes windfields read from the reader process to
!   all other processes. Usually one set of fields are transfered, but at
!   first timestep when there are no fields in memory, both are transfered.
!   MPI_Bcast is used, so implicitly all processes are synchronized at this
!   step
!
!
!*******************************************************************************
    use com_mod
    use par_mod,only: numwfmem

    implicit none

    integer, intent(in) :: memstat

! Common array sizes used for communications
    integer, parameter :: d3_size1 = nxmax*nymax*nzmax
    integer, parameter :: d3_size2 = nxmax*nymax*nuvzmax
    integer, parameter :: d2_size1 = nxmax*nymax
    integer, parameter :: d2_size2 = nxmax*nymax*maxspec
    integer, parameter :: d2_size3 = nxmax*nymax
    integer, parameter :: d1_size1 = maxwf

    integer :: d3s1,d3s2,d2s1,d2s2

! li,ui    lower,upper indices of fields to be transfered (1-3, ui>=li) 
    integer :: li,ui

! First time routine is called the unchangeable fields will be transfered    
    logical :: first_call=.true.

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

! Sizes of arrays transfered
    d3s1=d3_size1
    d3s2=d3_size2
    d2s1=d2_size1 
    d2s2=d2_size2

! Decide which field(s) need to be transfered
    if (memstat.ge.32) then ! distribute 2 fields, to 1st/2nd indices 
      li=1; ui=2
      d3s1=2*d3_size1
      d3s2=2*d3_size2
      d2s1=2*d2_size1 
      d2s2=2*d2_size2
    else if (memstat.eq.17) then ! distribute 1 field, place on 1st index
      li=1; ui=1
    else if (memstat.eq.18) then ! distribute 1 field, place on 2nd index
      li=2; ui=2
    else if (memstat.eq.19) then ! distribute 1 field, place on 3nd index
      li=3; ui=3
    else
      write(*,*) '#### mpi_mod::mpif_gf_send_vars> ERROR: ', &
           & 'wrong value of memstat, exiting ####', memstat
      stop
    end if


! Time for MPI communications
    if (mp_measure_time) call mpif_mtime('commtime',0)

! Send variables from getfield process (id_read) to other processes
!**********************************************************************

! The non-reader processes need to know if cloud water was read.
    call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 

! Static fields/variables sent only at startup
    if (first_call) then
      call MPI_Bcast(oro(:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_Bcast(excessoro(:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_Bcast(lsm(:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_Bcast(xlanduse(:,:,:),d2_size3*numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_Bcast(wftime,d1_size1,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /=