mpi_mod.f90

!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

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        *
!                                                                            *
!                                                                            *
!                                                                            *
!                                                                            *
!*****************************************************************************

  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 :: 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 !27 !29 :DBG:
!integer, dimension(:), allocatable :: tags
  integer, dimension(:), allocatable :: reqs


  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:
! 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 

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
    use com_mod, only: mpi_mode

    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.and.lroot) then
      write(*,FMT='(80("#"))')
      write(*,*) '#### mpi_mod::mpif_init> ERROR: ', &
           & 'numwfmem must be set to 3 for asyncronous reading ####'
      write(*,FMT='(80("#"))')
      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 '
      write(*,*) 'amount of memory is being 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. ! :DBG: eso fix this...
    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
    if (mp_partid.lt.mod(maxpart,mp_partgroup_np)) addmaxpart=1
    maxpart_mpi=int(maxpart/mp_partgroup_np)+addmaxpart

! 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
!    if (mp_pid.ge.0) then
!      call system_clock(s)
      s = 244
      mp_seed = -abs(mod((s*181)*((mp_pid-83)*359), 104729))
    end if
    if (mp_dev_mode) write(*,*) 'PID, mp_seed: ',mp_pid, mp_seed
    if (mp_dbg_mode) then
! :DBG: For debugging, set all seed to 0 and maxrand to e.g. 20
      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

    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)
!
!*******************************************************************************
    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_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 is 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 were 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 /= 0) goto 600
      call MPI_Bcast(numbwf,1,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_Bcast(nmixz,1,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600 
      call MPI_Bcast(height,nzmax,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600

      first_call=.false.
    endif

    call MPI_Bcast(uu(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(vv(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(uupol(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(vvpol(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(ww(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600
    call MPI_Bcast(tt(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(rho(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(drhodz(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(tth(:,:,:,li:ui),d3s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(qvh(:,:,:,li:ui),d3s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(qv(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600
    call MPI_Bcast(pv(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(clouds(:,:,:,li:ui),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600

! cloud water/ice:
    if (readclouds) then
      ! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      ! if (mp_ierr /= 0) goto 600
      call MPI_Bcast(clw4(:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      if (mp_ierr /= 0) goto 600
      ! call MPI_Bcast(clwc(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      ! if (mp_ierr /= 0) goto 600
      ! call MPI_Bcast(ciwc(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
      ! if (mp_ierr /= 0) goto 600
    end if

! 2D fields
    call MPI_Bcast(cloudsh(:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(vdep(:,:,:,li:ui),d2s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(ps(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600
    call MPI_Bcast(sd(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600
    call MPI_Bcast(tcc(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600
    call MPI_Bcast(tt2(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
    if (mp_ierr /= 0) goto 600 
    call MPI_Bcast(td2(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)