source: flexpart.git/src/mpi_mod.f90 @ 16b61a5

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 16b61a5 was 16b61a5, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Reworked the domain-filling option (MPI). Fixed a slow loop which had errors in loop counter (MPI)

  • Property mode set to 100644
File size: 106.7 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22module mpi_mod
23
24!*****************************************************************************
25!                                                                            *
26!  DESCRIPTION                                                               *
27!    This module contains subroutines and common variables used for the      *
28!    MPI parallelization of FLEXPART.                                        *
29!                                                                            *
30!  NOTE                                                                      *
31!    Depending on the MPI library installed on your system (e.g. mpich2,     *
32!    OpenMPI) you may need to choose below in this file between              *
33!      use mpi                                                               *
34!    (if the MPI library comes with the file 'mpi.mod'); or                  *
35!      include 'mpif.h'                                                      *
36!                                                                            *
37!                                                                            *
38!*****************************************************************************
39!                                                                            *
40! Variables:                                                                 *
41!                                                                            *
42! mp_ierr                 MPI error code                                     *
43! mp_np                   Number of MPI processes                            *
44! mp_pid                  Process ID of each MPI process                     *
45! mp_seed                 Parameter for random number seed                   *
46! read_grp_min            Minimum number of processes at which one will be   *
47!                         used as reader                                     *
48! numpart_mpi,            Number of particles per node                       *
49! maxpart_mpi                                                                *
50! mp_partid               MPI process ID for particle calculation            *
51! mp_partgroup_           Refers to the subset of processors performing      *
52!                         loops over particles. Will be all processes        *
53!                         unless a dedicated process runs getfields/readwind *
54! lmp_sync                If .false., use asynchronous MPI                   *
55! mp_cp                   Real precision to use for deposition fields        *
56!                                                                            *
57!                                                                            *
58!                                                                            *
59!                                                                            *
60!*****************************************************************************
61
62  use mpi
63  use par_mod, only: dp,sp
64  use com_mod, only: lroot
65
66  implicit none
67
68!  include 'mpif.h'
69
70  public
71
72! Set aside a process for reading windfields if using at least these many processes
73!==================================================
74  integer, parameter, private :: read_grp_min=4
75!==================================================
76
77! Variables for each MPI process in the world group
78  integer :: mp_ierr, mp_np, mp_pid, mp_partid
79  integer, private :: world_group_id
80
81! Variables for MPI processes in the 'particle' group
82  integer, allocatable, dimension(:) :: mp_partgroup_rank
83  integer :: mp_partgroup_comm, mp_partgroup_pid, mp_partgroup_np
84
85  integer :: mp_seed=0
86  integer, parameter :: mp_sp=MPI_REAL4, mp_dp=MPI_REAL8
87  integer :: mp_cp
88  integer, parameter :: id_root=0 ! master process
89
90! MPI tags/requests for send/receive operation
91  integer :: tm1
92  integer, parameter :: nvar_async=26
93!integer, dimension(:), allocatable :: tags
94  integer, dimension(:), allocatable :: reqs
95
96! Status array used for certain MPI operations (MPI_RECV)
97  integer, dimension(MPI_STATUS_SIZE) :: mp_status
98
99
100  integer :: id_read   ! readwind/getfield process
101  integer :: numpart_mpi,maxpart_mpi ! number of particles per node
102  integer :: tot_numpart=0
103  integer :: mp_comm_used ! global or subgroup communicator
104
105  logical :: lmpreader=.false. ! is set to true for reading process(es) only.
106  logical :: lmp_use_reader=.false. ! true if separate readwind process is used
107
108! .true. if only using synchronous MPI send/recv (default)
109! If setting this to .false., numwfmem must be set to 3
110!===============================================================================
111  logical :: lmp_sync=.true.
112!===============================================================================
113
114! mp_dbg_mode       Used for debugging MPI.
115! mp_dev_mode       various checks related to debugging the parallel code
116! mp_dbg_out        write some arrays to data file for debugging
117! mp_measure_time   Measure cpu/wall time, write out at end of run
118! mp_time_barrier   Measure MPI barrier time
119! mp_exact_numpart  Use an extra MPI communication to give the exact number of particles
120!                   to standard output (this does *not* otherwise affect the simulation)
121  logical, parameter :: mp_dbg_mode = .false.
122  logical, parameter :: mp_dev_mode = .false.
123  logical, parameter :: mp_dbg_out = .false.
124  logical, parameter :: mp_time_barrier=.true.
125  logical, parameter :: mp_measure_time=.false.
126  logical, parameter :: mp_exact_numpart=.true.
127
128! for measuring CPU/Wall time
129  real(sp),private :: mp_comm_time_beg, mp_comm_time_end, mp_comm_time_total=0.
130  real(dp),private :: mp_comm_wtime_beg, mp_comm_wtime_end, mp_comm_wtime_total=0.
131  real(sp),private :: mp_root_time_beg, mp_root_time_end, mp_root_time_total=0.
132  real(dp),private :: mp_root_wtime_beg, mp_root_wtime_end, mp_root_wtime_total=0.
133  real(sp),private :: mp_barrier_time_beg, mp_barrier_time_end, mp_barrier_time_total=0.
134  real(dp),private :: mp_barrier_wtime_beg, mp_barrier_wtime_end, mp_barrier_wtime_total=0.
135  real(sp),private :: tm_nploop_beg, tm_nploop_end, tm_nploop_total=0.
136  real(sp),private :: tm_tot_beg, tm_tot_end, tm_tot_total=0.
137  real(dp),private :: mp_getfields_wtime_beg, mp_getfields_wtime_end, mp_getfields_wtime_total=0.
138  real(sp),private :: mp_getfields_time_beg, mp_getfields_time_end, mp_getfields_time_total=0.
139  real(dp),private :: mp_readwind_wtime_beg, mp_readwind_wtime_end, mp_readwind_wtime_total=0.
140  real(sp),private :: mp_readwind_time_beg, mp_readwind_time_end, mp_readwind_time_total=0.
141  real(dp),private :: mp_io_wtime_beg, mp_io_wtime_end, mp_io_wtime_total=0.
142  real(sp),private :: mp_io_time_beg, mp_io_time_end, mp_io_time_total=0.
143  real(dp),private :: mp_wetdepo_wtime_beg, mp_wetdepo_wtime_end, mp_wetdepo_wtime_total=0.
144  real(sp),private :: mp_wetdepo_time_beg, mp_wetdepo_time_end, mp_wetdepo_time_total=0.
145  real(dp),private :: mp_advance_wtime_beg, mp_advance_wtime_end, mp_advance_wtime_total=0.
146  real(dp),private :: mp_conccalc_time_beg, mp_conccalc_time_end, mp_conccalc_time_total=0.
147  real(dp),private :: mp_total_wtime_beg, mp_total_wtime_end, mp_total_wtime_total=0.
148  real(dp),private :: mp_vt_wtime_beg, mp_vt_wtime_end, mp_vt_wtime_total
149  real(sp),private :: mp_vt_time_beg, mp_vt_time_end, mp_vt_time_total
150
151! dat_lun           logical unit number for i/o
152  integer, private :: dat_lun
153
154! Temporary arrays for particles (allocated and deallocated as needed)
155  integer, allocatable, dimension(:) :: nclass_tmp, npoint_tmp, itra1_tmp, idt_tmp, &
156       & itramem_tmp, itrasplit_tmp
157  real(kind=dp), allocatable, dimension(:) :: xtra1_tmp, ytra1_tmp
158  real, allocatable, dimension(:) :: ztra1_tmp
159  real, allocatable, dimension(:,:) :: xmass1_tmp
160
161! mp_redist_fract        Exchange particles between processes if relative difference
162!                        is larger. A good value is between 0.0 and 0.5
163! mp_maxpart_factor      Allocate more memory per process than strictly needed
164!                        (safety factor). Recommended value between 1.5 and 2.5
165! mp_min_redist          Do not redistribute particles if below this limit
166  real, parameter :: mp_redist_fract=0.2, mp_maxpart_factor=1.5
167  integer,parameter :: mp_min_redist=100000
168
169
170contains
171
172  subroutine mpif_init
173!***********************************************************************
174! DESCRIPTION
175!   Initialize MPI.
176!   
177!   Create the global communicator MPI_COMM_WORLD
178!   If using a separate MPI process for getfields/readwind, a subgroup
179!   is created for the other processes.
180!
181! VARIABLES
182!   mpi_mode    default 0, set to 2/3 if running MPI version
183!   mp_np       number of running processes, decided at run-time
184!***********************************************************************
185    use par_mod, only: maxpart, numwfmem, dep_prec
186    use com_mod, only: mpi_mode, verbosity
187
188    implicit none
189
190    integer :: i,j,s,addmaxpart=0
191
192! Each process gets an ID (mp_pid) in the range 0,..,mp_np-1
193    call MPI_INIT(mp_ierr)
194    if (mp_ierr /= 0) goto 100
195    call MPI_COMM_RANK(MPI_COMM_WORLD, mp_pid, mp_ierr)
196    if (mp_ierr /= 0) goto 100
197    call MPI_COMM_SIZE(MPI_COMM_WORLD, mp_np, mp_ierr)
198    if (mp_ierr /= 0) goto 100
199
200
201! Variable mpi_mode is used to handle subroutines common to parallel/serial version
202    if (lmp_sync) then
203      mpi_mode=2 ! hold 2 windfields in memory
204    else
205      mpi_mode=3 ! hold 3 windfields in memory
206    end if
207
208    if (mp_pid.ne.0) then
209      lroot = .false.
210    end if
211
212! Set MPI precision to use for transferring deposition fields
213!************************************************************
214    if (dep_prec==dp) then
215      mp_cp = MPI_REAL8
216! TODO: write info message for serial version as well
217      if (lroot.and.verbosity>0) write(*,*) 'Using double precision for deposition fields'
218    else if (dep_prec==sp) then
219      mp_cp = MPI_REAL4
220      if (lroot.and.verbosity>0) write(*,*) 'Using single precision for deposition fields'
221    else
222      write(*,*) 'ERROR: something went wrong setting MPI real precision'
223      stop
224    end if
225
226! Check for sensible combination of parameters
227!*********************************************
228
229    if (.not.lmp_sync.and.numwfmem.ne.3) then
230      if (lroot) then
231        write(*,FMT='(80("#"))')
232        write(*,*) '#### mpi_mod::mpif_init> ERROR: ', &
233             & 'numwfmem must be set to 3 for asyncronous reading ####'
234        write(*,FMT='(80("#"))')
235      end if
236      call MPI_FINALIZE(mp_ierr)
237      stop
238    else if (lmp_sync.and.numwfmem.ne.2.and.lroot) then
239      write(*,FMT='(80("#"))')
240      write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
241           & 'numwfmem should be set to 2 for syncronous'
242      write(*,*) ' reading. Results will still be valid, but unneccesary memory &
243           &is allocated.'
244      write(*,FMT='(80("#"))')
245! Force "syncronized" version if all processes will call getfields
246    else if ((.not.lmp_sync.and.mp_np.lt.read_grp_min).or.(mp_np.eq.1)) then
247      if (lroot) then
248        write(*,FMT='(80("#"))')
249        write(*,*) '#### mpi_mod::mpif_init> WARNING: ', &
250             & 'all procs call getfields. Setting lmp_sync=.true.'
251        write(*,FMT='(80("#"))')
252      end if
253      lmp_sync=.true.
254    end if
255
256! TODO: Add more warnings for unimplemented flexpart features
257
258! Set ID of process that calls getfield/readwind.
259! Using the last process in the group ensures statistical identical results
260! as running with one process less but not using separate read process
261!**********************************************************************
262
263!    id_read = min(mp_np-1, 1)
264    id_read = mp_np-1
265
266    if (mp_pid.eq.id_read) lmpreader=.true.
267
268    call MPI_Comm_group (MPI_COMM_WORLD, world_group_id, mp_ierr)
269
270! Create a MPI group/communicator that will calculate trajectories.
271! Skip this step if program is run with only a few processes
272!************************************************************************
273    allocate(mp_partgroup_rank(0:mp_np-2))
274
275! This allows checking for allocation of arrays when no subgroup is used
276    mp_partgroup_pid=0
277
278    if (read_grp_min.lt.2) then
279      write(*,*) '#### mpi_mod::mpif_init> ERROR ####', &
280           & 'read_grp_min should be at least 2. Exiting'
281      stop
282    end if
283
284    if (mp_np.ge.read_grp_min) then
285      lmp_use_reader = .true.
286
287! Map the subgroup IDs= 0,1,2,...,mp_np-2, skipping reader process
288      j=-1
289      do i=0, mp_np-2 !loop over all (subgroup) IDs
290        j=j+1
291        if (i.eq.id_read) then
292          j=j+1
293        end if
294        mp_partgroup_rank(i) = j
295      end do
296
297      call MPI_Group_incl (world_group_id, mp_np-1, mp_partgroup_rank, &
298           &mp_partgroup_pid, mp_ierr)
299      if (mp_ierr /= 0) goto 100
300      call MPI_Comm_create (MPI_COMM_WORLD, mp_partgroup_pid, mp_partgroup_comm, mp_ierr)
301      if (mp_ierr /= 0) goto 100
302
303      if (mp_pid.ne.id_read) then
304        call MPI_Comm_rank (mp_partgroup_comm, mp_partgroup_pid, mp_ierr)
305        if (mp_ierr /= 0) goto 100
306
307! Get size of new group (=mp_np-1)
308        call MPI_COMM_SIZE(mp_partgroup_comm, mp_partgroup_np, mp_ierr)
309        if (mp_ierr /= 0) goto 100
310        if (mp_partgroup_np.ne.mp_np-1) then
311          write(*,*)  '#### mpi_mod:: mpif_init> ERROR ####&
312               & mp_partgroup_np.ne.mp_np-1'
313          stop
314        endif
315
316      else
317        mp_partgroup_pid = -1
318      end if
319    end if
320
321    if (lmp_use_reader) then
322      mp_comm_used = mp_partgroup_comm
323      mp_partid = mp_partgroup_pid
324      mp_partgroup_np=mp_np-1
325    else
326      mp_comm_used = MPI_COMM_WORLD
327      mp_partgroup_np = mp_np
328      mp_partid = mp_pid
329    end if
330
331! Set maxpart per process
332! eso 08/2016: Increase maxpart per process, in case of unbalanced distribution
333    maxpart_mpi=int(mp_maxpart_factor*maxpart/mp_partgroup_np)
334
335! Do not allocate particle data arrays for readwind process
336    if (lmpreader.and.lmp_use_reader) then
337      maxpart_mpi=0
338    end if
339
340! Set random seed for each non-root process
341    if (mp_pid.gt.0) then
342      s = 244
343      mp_seed = -abs(mod((s*181)*((mp_pid-83)*359), 104729))
344    end if
345    if (mp_dbg_mode) then
346      mp_seed=0
347      if (lroot) write(*,*) 'MPI: setting seed=0'
348    end if
349
350! Allocate request array for asynchronous MPI
351    if (.not.lmp_sync) then
352      allocate(reqs(0:nvar_async*mp_np-1))
353      reqs(:)=MPI_REQUEST_NULL
354    else
355      allocate(reqs(0:1))
356      reqs(:)=MPI_REQUEST_NULL
357    end if
358
359    goto 101
360
361100 write(*,*) '#### mpi_mod::mpif_init> ERROR ####', mp_ierr
362    stop
363
364101 end subroutine mpif_init
365
366
367  subroutine mpif_mpi_barrier
368!***********************************************************************
369! Set MPI barrier (all processes are synchronized here), with the
370! possible exception of a separate 'readwind' process.
371! Optionally measure cpu/wall time.
372!
373!***********************************************************************
374    implicit none
375
376    if (mp_time_barrier) then
377      call cpu_time(mp_barrier_time_beg)
378      mp_barrier_wtime_beg = mpi_wtime()
379    endif
380
381    call MPI_BARRIER(mp_comm_used, mp_ierr)
382    if (mp_ierr /= 0) goto 100
383
384    if (mp_time_barrier) then
385      call cpu_time(mp_barrier_time_end)
386      mp_barrier_wtime_end = mpi_wtime()
387
388      mp_barrier_time_total=mp_barrier_time_total+&
389           &(mp_barrier_time_end-mp_barrier_time_beg)
390      mp_barrier_wtime_total=mp_barrier_wtime_total+&
391           &(mp_barrier_wtime_end-mp_barrier_wtime_beg)
392    end if
393
394    goto 101
395
396100 write(*,*) '#### mpi_mod::mpif_mpi_barrier> ERROR ####', mp_ierr
397    stop
398
399101 end subroutine mpif_mpi_barrier
400
401
402  subroutine mpif_com_mod_allocate
403!*******************************************************************************   
404! Dynamic allocation of arrays (in serial code these have size maxpart)
405!
406! Not in use anymore, moved to com_mod for interoperability with serial version
407!
408! TODO: remove
409!*******************************************************************************
410    use com_mod
411
412    implicit none
413
414    integer :: array_size
415
416    array_size = maxpart_mpi
417
418    allocate(itra1(array_size),npoint(array_size),&
419         & nclass(array_size),idt(array_size),&
420         & itramem(array_size),itrasplit(array_size),&
421         & xtra1(array_size),ytra1(array_size),&
422         & ztra1(array_size),xmass1(array_size, maxspec))
423
424    allocate(uap(array_size),ucp(array_size),&
425         & uzp(array_size),us(array_size),&
426         & vs(array_size),&
427         & ws(array_size),cbt(array_size))
428
429
430  end subroutine mpif_com_mod_allocate
431
432
433  subroutine mpif_tm_send_dims
434!***********************************************************************
435! Distribute array dimensions from pid0 to all processes.
436! numpart_mpi/numpart is sent to allow dynamic allocation
437!
438! Currently not in use.
439!
440! Variables of similar type (integer, real) are placed in an array
441! and sent collectively, to avoid the overhead associated with individual
442! MPI send/receives
443!
444!
445!***********************************************************************
446    use com_mod
447
448    implicit none
449
450    integer :: add_part=0
451
452    call MPI_Bcast(numpart, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
453
454! MPI subgroup does the particle-loop
455    if (lmp_use_reader) then
456      if (mod(numpart,mp_partgroup_np).ne.0) add_part=1
457      numpart_mpi=int(numpart/mp_partgroup_np)+add_part
458    else
459      if (mod(numpart,mp_np).ne.0) add_part=1
460      numpart_mpi=int(numpart/mp_np)+add_part
461    end if
462
463
464! redefine numpart as 'numpart per process' throughout the code
465!**************************************************************
466    numpart = numpart_mpi
467
468  end subroutine mpif_tm_send_dims
469
470
471  subroutine mpif_send_part_properties(num_part)
472!***********************************************************************
473! Distribute particle properties from root to all processes.
474
475! Used by init_domainfill_mpi
476!
477! Variables:
478!
479! num_part        input, number of particles per process (rounded up)
480!
481!***********************************************************************
482    use com_mod
483
484    implicit none
485
486    integer,intent(in) :: num_part
487    integer :: i,jj, addone
488
489! Exit if too many particles
490    if (num_part.gt.maxpart_mpi) then
491      write(*,*) '#####################################################'
492      write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED  ####'
493      write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE  ####'
494      write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
495      write(*,*) '#### OR INCREASE MAXPART.                        ####'
496      write(*,*) '#####################################################'
497!      call MPI_FINALIZE(mp_ierr)
498      stop
499    end if
500
501
502! Time for MPI communications
503!****************************
504    if (mp_measure_time) call mpif_mtime('commtime',0)
505
506! Distribute variables, send from pid 0 to other processes (including itself)
507!****************************************************************************
508
509    call MPI_SCATTER(nclass_tmp,num_part,MPI_INTEGER,nclass,&
510         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
511    if (mp_ierr /= 0) goto 600
512    call MPI_SCATTER(npoint_tmp,num_part,MPI_INTEGER,npoint,&
513         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
514    if (mp_ierr /= 0) goto 600
515    call MPI_SCATTER(itra1_tmp,num_part,MPI_INTEGER,itra1,&
516         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
517    if (mp_ierr /= 0) goto 600
518    call MPI_SCATTER(idt_tmp,num_part,MPI_INTEGER,idt,&
519         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
520    if (mp_ierr /= 0) goto 600
521    call MPI_SCATTER(itramem_tmp,num_part,MPI_INTEGER,itramem,&
522         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
523    if (mp_ierr /= 0) goto 600
524    call MPI_SCATTER(itrasplit_tmp,num_part,MPI_INTEGER,itrasplit,&
525         &num_part,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
526    if (mp_ierr /= 0) goto 600
527    call MPI_SCATTER(xtra1_tmp,num_part,mp_dp,xtra1,&
528         &num_part,mp_dp,id_root,mp_comm_used,mp_ierr)
529    if (mp_ierr /= 0) goto 600
530    call MPI_SCATTER(ytra1_tmp,num_part,mp_dp,ytra1,&
531         &num_part,mp_dp,id_root,mp_comm_used,mp_ierr)
532    if (mp_ierr /= 0) goto 600
533    call MPI_SCATTER(ztra1_tmp,num_part,mp_sp,ztra1,&
534         &num_part,mp_sp,id_root,mp_comm_used,mp_ierr)
535    if (mp_ierr /= 0) goto 600
536    do i=1,nspec
537      call MPI_SCATTER(xmass1_tmp(:,i),num_part,mp_sp,xmass1(:,i),&
538           &num_part,mp_sp,id_root,mp_comm_used,mp_ierr)
539      if (mp_ierr /= 0) goto 600
540    end do
541
542    if (mp_measure_time) call mpif_mtime('commtime',1)
543
544    goto 601
545
546600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
547    stop
548
549! After the transfer of particles it it possible that the value of
550! "numpart" is larger than the actual used, so we reduce it if there are
551! invalid particles at the end of the arrays
552601 do i=num_part, 1, -1
553      if (itra1(i).eq.-999999999) then
554        numpart=numpart-1
555      else
556        exit
557      end if
558    end do
559
560
561!601 end subroutine mpif_send_part_properties
562  end subroutine mpif_send_part_properties
563
564
565  subroutine mpif_calculate_part_redist(itime)
566!***********************************************************************
567! Calculate number of particles to redistribute between processes. This routine
568! can be called at regular time intervals to keep a level number of
569! particles on each process.
570!
571! First, all processes report their local "numpart" to each other, which is
572! stored in array "numpart_mpi(np)". This array is sorted from low to
573! high values, along with a corresponding process ID array "idx_arr(np)".
574! If the relative difference between the highest and lowest "numpart_mpi" value
575! is above a threshold, particles are transferred from process idx_arr(np-1)
576! to process idx_arr(0). Repeat for processes idx_arr(np-i) and idx_arr(i)
577! for all valid i.
578! Note: If np is an odd number, the process with the median
579! number of particles is left unchanged
580!
581! VARIABLES
582! itime        input, current time
583!***********************************************************************
584    use com_mod
585
586    implicit none
587   
588    integer, intent(in) :: itime
589    real :: pmin,z
590    integer :: i,jj,nn, num_part=1,m,imin, num_trans
591    logical :: first_iter
592    integer,allocatable,dimension(:) :: numparticles_mpi, idx_arr
593    real,allocatable,dimension(:) :: sorted ! TODO: we don't really need this
594
595! Exit if running with only 1 process
596!************************************************************************
597    if (mp_np.eq.1) return
598
599! All processes exchange information on number of particles
600!****************************************************************************
601    allocate(numparticles_mpi(0:mp_partgroup_np-1), &
602         &idx_arr(0:mp_partgroup_np-1), sorted(0:mp_partgroup_np-1))
603
604    call MPI_Allgather(numpart, 1, MPI_INTEGER, numparticles_mpi, &
605         & 1, MPI_INTEGER, mp_comm_used, mp_ierr)
606
607
608! Sort from lowest to highest
609! Initial guess: correct order
610    sorted(:) = numparticles_mpi(:)
611    do i=0, mp_partgroup_np-1
612      idx_arr(i) = i
613    end do
614
615! For each successive element in index array, see if a lower value exists
616    do i=0, mp_partgroup_np-2
617      pmin=sorted(i)
618      imin=idx_arr(i)
619      m=i+1
620      do jj=m, mp_partgroup_np-1
621        if (pmin.le.sorted(jj)) cycle
622        z=pmin
623        pmin=sorted(jj)
624        sorted(jj)=z
625
626        nn=imin
627        imin=idx_arr(jj)
628        idx_arr(jj)=nn
629      end do
630      sorted(i)=pmin
631      idx_arr(i)=imin
632    end do
633
634! For each pair of processes, transfer particles if the difference is above a
635! limit, and numpart at sending process large enough
636
637    m=mp_partgroup_np-1 ! index for last sorted process (most particles)
638    do i=0,mp_partgroup_np/2-1
639      num_trans = numparticles_mpi(idx_arr(m)) - numparticles_mpi(idx_arr(i))
640      if (mp_partid.eq.idx_arr(m).or.mp_partid.eq.idx_arr(i)) then
641        if ( numparticles_mpi(idx_arr(m)).gt.mp_min_redist.and.&
642             & real(num_trans)/real(numparticles_mpi(idx_arr(m))).gt.mp_redist_fract) then
643! DBG
644          ! write(*,*) 'mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi', &
645          !      &mp_partid, idx_arr(m), idx_arr(i), mp_min_redist, num_trans, numparticles_mpi
646! DBG
647          call mpif_redist_part(itime, idx_arr(m), idx_arr(i), num_trans/2)
648        end if
649      end if
650      m=m-1
651    end do
652
653    deallocate(numparticles_mpi, idx_arr, sorted)
654
655  end subroutine mpif_calculate_part_redist
656
657
658  subroutine mpif_redist_part(itime, src_proc, dest_proc, num_trans)
659!***********************************************************************
660! Transfer particle properties between two arbitrary processes.
661!
662! VARIABLES
663!
664! itime           input, current time
665! src_proc        input, ID of source process
666! dest_proc       input, ID of destination process
667! num_trans       input, number of particles to transfer
668!
669!************************************************************************
670    use com_mod !TODO:  ,only: nclass etc
671
672    implicit none
673
674    integer, intent(in) :: itime, src_proc, dest_proc, num_trans
675    integer :: ll, ul ! lower and upper indices in arrays
676    integer :: arr_sz ! size of temporary arrays
677    integer :: mtag   ! MPI message tag
678    integer :: i, j, minpart, ipart, maxnumpart
679 
680! Check for invalid input arguments
681!**********************************
682 if (src_proc.eq.dest_proc) then
683   write(*,*) '<mpi_mod::mpif_redist_part>: Error: &
684        &src_proc.eq.dest_proc'
685   stop
686 end if
687
688! Measure time for MPI communications
689!************************************
690    if (mp_measure_time) call mpif_mtime('commtime',0)
691
692! Send particles
693!***************
694    if (mp_partid.eq.src_proc) then
695      mtag = 2000
696
697! Array indices for the transferred particles
698      ll=numpart-num_trans+1
699      ul=numpart
700
701      if (mp_dev_mode) then
702        write(*,FMT='(72("#"))')
703        write(*,*) "Sending ", num_trans, "particles (from/to)", src_proc, dest_proc
704        write(*,*) "numpart before: ", numpart
705      end if
706
707      call MPI_SEND(nclass(ll:ul),num_trans,&
708           & MPI_INTEGER,dest_proc,mtag+1,mp_comm_used,mp_ierr)
709
710      call MPI_SEND(npoint(ll:ul),num_trans,&
711           & MPI_INTEGER,dest_proc,mtag+2,mp_comm_used,mp_ierr)
712
713      call MPI_SEND(itra1(ll:ul),num_trans, &
714           & MPI_INTEGER,dest_proc,mtag+3,mp_comm_used,mp_ierr)
715
716      call MPI_SEND(idt(ll:ul),num_trans, &
717           & MPI_INTEGER,dest_proc,mtag+4,mp_comm_used,mp_ierr)
718
719      call MPI_SEND(itramem(ll:ul),num_trans, &
720           & MPI_INTEGER,dest_proc,mtag+5,mp_comm_used,mp_ierr)
721
722      call MPI_SEND(itrasplit(ll:ul),num_trans,&
723           & MPI_INTEGER,dest_proc,mtag+6,mp_comm_used,mp_ierr)
724
725      call MPI_SEND(xtra1(ll:ul),num_trans, &
726           & mp_dp,dest_proc,mtag+7,mp_comm_used,mp_ierr)
727
728      call MPI_SEND(ytra1(ll:ul),num_trans,&
729           & mp_dp,dest_proc,mtag+8,mp_comm_used,mp_ierr)
730
731      call MPI_SEND(ztra1(ll:ul),num_trans,&
732           & mp_sp,dest_proc,mtag+9,mp_comm_used,mp_ierr)
733
734      do j=1,nspec
735        call MPI_SEND(xmass1(ll:ul,j),num_trans,&
736             & mp_sp,dest_proc,mtag+(9+j),mp_comm_used,mp_ierr)
737      end do
738
739! Terminate transferred particles, update numpart
740      itra1(ll:ul) = -999999999
741
742      numpart = numpart-num_trans
743
744      if (mp_dev_mode) then
745        write(*,*) "numpart after: ", numpart
746        write(*,FMT='(72("#"))')
747      end if
748
749    else if (mp_partid.eq.dest_proc) then
750
751! Receive particles
752!******************
753      mtag = 2000
754      ! Array indices for the transferred particles
755      ll=numpart+1
756      ul=numpart+num_trans
757
758      if (mp_dev_mode) then
759        write(*,FMT='(72("#"))')
760        write(*,*) "Receiving ", num_trans, "particles (from/to)", src_proc, dest_proc
761        write(*,*) "numpart before: ", numpart
762      end if
763
764! We could receive the data directly at nclass(ll:ul) etc., but this leaves
765! vacant spaces at lower indices. Using temporary arrays instead.
766      arr_sz = ul-ll+1
767      allocate(itra1_tmp(arr_sz),npoint_tmp(arr_sz),nclass_tmp(arr_sz),&
768           & idt_tmp(arr_sz),itramem_tmp(arr_sz),itrasplit_tmp(arr_sz),&
769           & xtra1_tmp(arr_sz),ytra1_tmp(arr_sz),ztra1_tmp(arr_sz),&
770           & xmass1_tmp(arr_sz, maxspec))
771     
772      call MPI_RECV(nclass_tmp,num_trans,&
773           & MPI_INTEGER,src_proc,mtag+1,mp_comm_used,mp_status,mp_ierr)
774
775      call MPI_RECV(npoint_tmp,num_trans,&
776           & MPI_INTEGER,src_proc,mtag+2,mp_comm_used,mp_status,mp_ierr)
777
778      call MPI_RECV(itra1_tmp,num_trans, &
779           & MPI_INTEGER,src_proc,mtag+3,mp_comm_used,mp_status,mp_ierr)
780
781      call MPI_RECV(idt_tmp,num_trans, &
782           & MPI_INTEGER,src_proc,mtag+4,mp_comm_used,mp_status,mp_ierr)
783
784      call MPI_RECV(itramem_tmp,num_trans, &
785           & MPI_INTEGER,src_proc,mtag+5,mp_comm_used,mp_status,mp_ierr)
786
787      call MPI_RECV(itrasplit_tmp,num_trans,&
788           & MPI_INTEGER,src_proc,mtag+6,mp_comm_used,mp_status,mp_ierr)
789
790      call MPI_RECV(xtra1_tmp,num_trans, &
791           & mp_dp,src_proc,mtag+7,mp_comm_used,mp_status,mp_ierr)
792
793      call MPI_RECV(ytra1_tmp,num_trans,&
794           & mp_dp,src_proc,mtag+8,mp_comm_used,mp_status,mp_ierr)
795
796      call MPI_RECV(ztra1_tmp,num_trans,&
797           & mp_sp,src_proc,mtag+9,mp_comm_used,mp_status,mp_ierr)
798
799      do j=1,nspec
800        call MPI_RECV(xmass1_tmp(:,j),num_trans,&
801             & mp_sp,src_proc,mtag+(9+j),mp_comm_used,mp_status,mp_ierr)
802      end do
803
804! This is the maximum value numpart can possibly have after the transfer
805      maxnumpart=numpart+num_trans
806
807! Search for vacant space and copy from temporary storage
808!********************************************************
809      minpart=1
810      do i=1, num_trans
811! Take into acount that we may have transferred invalid particles
812        if (itra1_tmp(i).ne.itime) cycle
813        do ipart=minpart,maxnumpart
814          if (itra1(ipart).ne.itime) then
815            itra1(ipart) = itra1_tmp(i)
816            npoint(ipart) = npoint_tmp(i)
817            nclass(ipart) = nclass_tmp(i)
818            idt(ipart) = idt_tmp(i)
819            itramem(ipart) = itramem_tmp(i)
820            itrasplit(ipart) = itrasplit_tmp(i)
821            xtra1(ipart) = xtra1_tmp(i)
822            ytra1(ipart) = ytra1_tmp(i)
823            ztra1(ipart) = ztra1_tmp(i)
824            xmass1(ipart,:) = xmass1_tmp(i,:)
825            goto 200 ! Storage space has been found, stop searching
826          end if
827! :TODO: add check for if too many particles requiried
828        end do
829200     minpart=ipart+1
830      end do
831! Increase numpart, if necessary
832      numpart=max(numpart,ipart)
833
834      deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,&
835           & xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp)
836
837      if (mp_dev_mode) then
838        write(*,*) "numpart after: ", numpart
839        write(*,FMT='(72("#"))')
840      end if
841
842    else
843! This routine should only be called by the two participating processes
844      write(*,*) "ERROR: wrong process has entered mpi_mod::mpif_redist_part"
845      stop
846      return
847    end if
848
849! Measure time for MPI communications
850!************************************
851    if (mp_measure_time) call mpif_mtime('commtime',1)
852
853  end subroutine mpif_redist_part
854
855
856  subroutine mpif_tm_send_vars
857!***********************************************************************
858! Distribute particle variables from pid0 to all processes.
859! Called from timemanager
860! *NOT IN USE* at the moment, but can be useful for debugging
861!
862!***********************************************************************
863    use com_mod
864
865    implicit none
866
867    integer :: i
868
869! Time for MPI communications
870!****************************
871    if (mp_measure_time) call mpif_mtime('commtime',0)
872
873! Distribute variables, send from pid 0 to other processes
874!**********************************************************************
875
876! integers
877    if (lroot) then
878      call MPI_SCATTER(npoint,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
879           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
880      if (mp_ierr /= 0) goto 600
881      call MPI_SCATTER(idt,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
882           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
883      if (mp_ierr /= 0) goto 600
884      call MPI_SCATTER(itra1,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
885           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
886      if (mp_ierr /= 0) goto 600
887      call MPI_SCATTER(nclass,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
888           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
889      if (mp_ierr /= 0) goto 600
890      call MPI_SCATTER(itramem,numpart_mpi,MPI_INTEGER,MPI_IN_PLACE,&
891           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
892      if (mp_ierr /= 0) goto 600
893
894! int2
895      call MPI_SCATTER(cbt,numpart_mpi,MPI_INTEGER2,MPI_IN_PLACE,&
896           & numpart_mpi,MPI_INTEGER2,id_root,mp_comm_used,mp_ierr)
897      if (mp_ierr /= 0) goto 600
898
899! real
900      call MPI_SCATTER(uap,numpart_mpi,mp_sp,MPI_IN_PLACE,&
901           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
902      if (mp_ierr /= 0) goto 600
903      call MPI_SCATTER(ucp,numpart_mpi,mp_sp,MPI_IN_PLACE,&
904           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
905      if (mp_ierr /= 0) goto 600
906      call MPI_SCATTER(uzp,numpart_mpi,mp_sp,MPI_IN_PLACE,&
907           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
908      if (mp_ierr /= 0) goto 600
909
910      call MPI_SCATTER(us,numpart_mpi,mp_sp,MPI_IN_PLACE,&
911           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
912      if (mp_ierr /= 0) goto 600
913      call MPI_SCATTER(vs,numpart_mpi,mp_sp,MPI_IN_PLACE,&
914           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
915      if (mp_ierr /= 0) goto 600
916      call MPI_SCATTER(ws,numpart_mpi,mp_sp,MPI_IN_PLACE,&
917           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
918      if (mp_ierr /= 0) goto 600
919
920      call MPI_SCATTER(xtra1,numpart_mpi,mp_dp,MPI_IN_PLACE,&
921           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
922      if (mp_ierr /= 0) goto 600
923      call MPI_SCATTER(ytra1,numpart_mpi,mp_dp,MPI_IN_PLACE,&
924           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
925      if (mp_ierr /= 0) goto 600
926      call MPI_SCATTER(ztra1,numpart_mpi,mp_sp,MPI_IN_PLACE,&
927           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
928      if (mp_ierr /= 0) goto 600
929
930      do i=1,nspec
931        call MPI_SCATTER(xmass1(:,i),numpart_mpi,mp_sp,MPI_IN_PLACE,&
932             & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
933        if (mp_ierr /= 0) goto 600
934      end do
935
936    else ! (mp_pid >= 1)
937! integers
938      call MPI_SCATTER(npoint,numpart_mpi,MPI_INTEGER,npoint,&
939           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
940      if (mp_ierr /= 0) goto 600
941      call MPI_SCATTER(idt,numpart_mpi,MPI_INTEGER,idt,&
942           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
943      if (mp_ierr /= 0) goto 600
944      call MPI_SCATTER(itra1,numpart_mpi,MPI_INTEGER,itra1,&
945           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
946      if (mp_ierr /= 0) goto 600
947      call MPI_SCATTER(nclass,numpart_mpi,MPI_INTEGER,nclass,&
948           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
949      if (mp_ierr /= 0) goto 600
950      call MPI_SCATTER(itramem,numpart_mpi,MPI_INTEGER,itramem,&
951           & numpart_mpi,MPI_INTEGER,id_root,mp_comm_used,mp_ierr)
952      if (mp_ierr /= 0) goto 600
953
954! int2
955      call MPI_SCATTER(cbt,numpart_mpi,MPI_INTEGER2,cbt,&
956           & numpart_mpi,MPI_INTEGER2,id_root,mp_comm_used,mp_ierr)
957      if (mp_ierr /= 0) goto 600
958
959! reals
960      call MPI_SCATTER(uap,numpart_mpi,mp_sp,uap,&
961           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
962      if (mp_ierr /= 0) goto 600
963      call MPI_SCATTER(ucp,numpart_mpi,mp_sp,ucp,&
964           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
965      if (mp_ierr /= 0) goto 600
966      call MPI_SCATTER(uzp,numpart_mpi,mp_sp,uzp,&
967           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
968      if (mp_ierr /= 0) goto 600
969
970      call MPI_SCATTER(us,numpart_mpi,mp_sp,us,&
971           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
972      if (mp_ierr /= 0) goto 600
973      call MPI_SCATTER(vs,numpart_mpi,mp_sp,vs,&
974           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
975      if (mp_ierr /= 0) goto 600
976      call MPI_SCATTER(ws,numpart_mpi,mp_sp,ws,&
977           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
978      if (mp_ierr /= 0) goto 600
979
980      call MPI_SCATTER(xtra1,numpart_mpi,mp_dp,xtra1,&
981           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
982      if (mp_ierr /= 0) goto 600
983      call MPI_SCATTER(ytra1,numpart_mpi,mp_dp,ytra1,&
984           & numpart_mpi,mp_dp,id_root,mp_comm_used,mp_ierr)
985      if (mp_ierr /= 0) goto 600
986      call MPI_SCATTER(ztra1,numpart_mpi,mp_sp,ztra1,&
987           & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
988      if (mp_ierr /= 0) goto 600
989
990      do i=1,nspec
991        call MPI_SCATTER(xmass1(:,i),numpart_mpi,mp_sp,xmass1(:,i),&
992             & numpart_mpi,mp_sp,id_root,mp_comm_used,mp_ierr)
993        if (mp_ierr /= 0) goto 600
994      end do
995
996    end if !(lroot)
997
998    if (mp_measure_time) call mpif_mtime('commtime',1)
999
1000    goto 601
1001
1002600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
1003    stop
1004601 end subroutine mpif_tm_send_vars
1005
1006
1007  subroutine mpif_tm_collect_vars
1008!*******************************************************************************
1009! Collect particle data to PID 0 from all processes.                           *
1010!                                                                              *
1011! This can be called from timemanager_mpi, before caclulating                  *
1012! concentrations/writing output grids.                                         *
1013!                                                                              *
1014! Currently *not in use*, as each process calculates concentrations            *
1015! separately, but can be useful for debugging                                  *
1016!                                                                              *
1017! To use this routine (together with mpif_tm_send_vars) to transfer data       *
1018! to the MPI root process (useful for debugging), insert code like this        *
1019! (in timemanager_mpi.f90)                                                     *
1020!                                                                              *
1021!      if (lroot) tot_numpart=numpart ! root stores total numpart              *
1022!      call mpif_tm_send_dims                                                  *
1023!      if (numpart>1) then                                                     *
1024!        call mpif_tm_send_vars                                                *
1025!      end if                                                                  *
1026!                                                                              *
1027!    To collect data from all processes to root process after a parallel       *
1028!    region:                                                                   *
1029!                                                                              *
1030!      if (numpart>0) then                                                     *
1031!          if (lroot) numpart = tot_numpart                                    *
1032!       call mpif_tm_collect_vars                                              *
1033!      end if                                                                  *
1034!                                                                              *
1035! Then a section that begins with "if (lroot) ..." will behave like            *
1036! serial flexpart                                                              *
1037!                                                                              *
1038!*******************************************************************************
1039    use com_mod !, only: numpart, nspec, itra1, npoint, nclass
1040
1041    implicit none
1042
1043    integer :: i
1044
1045    logical :: use_mp_vars = .false.! use mp_... type allocated vars
1046    logical :: use_in_place = .true.! use MPI_IN_PLACE to collect vars.
1047
1048
1049! Time for MPI communications
1050    if (mp_measure_time) call mpif_mtime('commtime',0)
1051
1052! Distribute variables, send from pid 0 to other processes
1053!**********************************************************************
1054
1055    if (.not. use_mp_vars.and..not.use_in_place) then
1056
1057! Integers:
1058      call MPI_GATHER(npoint, numpart_mpi, MPI_INTEGER, npoint, &
1059           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1060      if (mp_ierr /= 0) goto 600
1061      call MPI_GATHER(idt, numpart_mpi, MPI_INTEGER, idt, &
1062           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1063      if (mp_ierr /= 0) goto 600
1064      call MPI_GATHER(itra1, numpart_mpi, MPI_INTEGER, itra1, &
1065           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1066      if (mp_ierr /= 0) goto 600
1067      call MPI_GATHER(nclass, numpart_mpi, MPI_INTEGER, nclass, &
1068           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1069      if (mp_ierr /= 0) goto 600
1070      call MPI_GATHER(itramem, numpart_mpi, MPI_INTEGER, itramem, &
1071           &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1072      if (mp_ierr /= 0) goto 600
1073
1074! int2
1075      call MPI_GATHER(cbt, numpart_mpi, MPI_INTEGER2, cbt, &
1076           &numpart_mpi, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
1077      if (mp_ierr /= 0) goto 600
1078
1079! Reals:
1080      call MPI_GATHER(uap, numpart_mpi, mp_sp, uap, &
1081           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1082      if (mp_ierr /= 0) goto 600
1083      call MPI_GATHER(ucp, numpart_mpi, mp_sp, ucp, &
1084           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1085      if (mp_ierr /= 0) goto 600
1086      call MPI_GATHER(uzp, numpart_mpi, mp_sp, uzp, &
1087           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1088      if (mp_ierr /= 0) goto 600
1089
1090      call MPI_GATHER(us, numpart_mpi, mp_sp, us, &
1091           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1092      if (mp_ierr /= 0) goto 600
1093      call MPI_GATHER(vs, numpart_mpi, mp_sp, vs, &
1094           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1095      if (mp_ierr /= 0) goto 600
1096      call MPI_GATHER(ws, numpart_mpi, mp_sp, ws, &
1097           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1098      if (mp_ierr /= 0) goto 600
1099
1100
1101      call MPI_GATHER(xtra1, numpart_mpi, mp_dp, xtra1, &
1102           &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
1103      if (mp_ierr /= 0) goto 600
1104      call MPI_GATHER(ytra1, numpart_mpi, mp_dp, ytra1, &
1105           &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
1106      if (mp_ierr /= 0) goto 600
1107      call MPI_GATHER(ztra1, numpart_mpi, mp_sp, ztra1, &
1108           &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1109      if (mp_ierr /= 0) goto 600
1110
1111      do i=1, nspec
1112        call MPI_GATHER(xmass1(:,i), numpart_mpi, mp_sp, xmass1(:,i), &
1113             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1114        if (mp_ierr /= 0) goto 600
1115
1116      end do
1117
1118! Use variables npoint etc directly for communications
1119!***********************************************************************
1120    else if (use_in_place.and..not.use_mp_vars) then
1121      if (lroot) then
1122! Integers:
1123        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, npoint, &
1124             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1125        if (mp_ierr /= 0) goto 600
1126        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, idt, &
1127             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1128        if (mp_ierr /= 0) goto 600
1129        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, itra1, &
1130             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1131        if (mp_ierr /= 0) goto 600
1132        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, nclass, &
1133             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1134        if (mp_ierr /= 0) goto 600
1135        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER, itramem, &
1136             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1137        if (mp_ierr /= 0) goto 600
1138
1139! int2
1140        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, MPI_INTEGER2, cbt, &
1141             &numpart_mpi, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
1142        if (mp_ierr /= 0) goto 600
1143
1144! Reals:
1145        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, uap, &
1146             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1147        if (mp_ierr /= 0) goto 600
1148        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, ucp, &
1149             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1150        if (mp_ierr /= 0) goto 600
1151        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, uzp, &
1152             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1153        if (mp_ierr /= 0) goto 600
1154
1155        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, us, &
1156             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1157        if (mp_ierr /= 0) goto 600
1158        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, vs, &
1159             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1160        if (mp_ierr /= 0) goto 600
1161        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, ws, &
1162             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1163        if (mp_ierr /= 0) goto 600
1164
1165
1166        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_dp, xtra1, &
1167             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
1168        if (mp_ierr /= 0) goto 600
1169        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_dp, ytra1, &
1170             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
1171        if (mp_ierr /= 0) goto 600
1172        call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, ztra1, &
1173             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1174        if (mp_ierr /= 0) goto 600
1175
1176        do i=1, nspec
1177          call MPI_GATHER(MPI_IN_PLACE, numpart_mpi, mp_sp, xmass1(:,i), &
1178               &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1179          if (mp_ierr /= 0) goto 600
1180        end do
1181
1182      else ! (for mp_pid >= 1)
1183
1184! Integers:
1185        call MPI_GATHER(npoint, numpart_mpi, MPI_INTEGER, npoint, &
1186             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1187        if (mp_ierr /= 0) goto 600
1188        call MPI_GATHER(idt, numpart_mpi, MPI_INTEGER, idt, &
1189             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1190        if (mp_ierr /= 0) goto 600
1191        call MPI_GATHER(itra1, numpart_mpi, MPI_INTEGER, itra1, &
1192             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1193        if (mp_ierr /= 0) goto 600
1194        call MPI_GATHER(nclass, numpart_mpi, MPI_INTEGER, nclass, &
1195             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1196        if (mp_ierr /= 0) goto 600
1197        call MPI_GATHER(itramem, numpart_mpi, MPI_INTEGER, itramem, &
1198             &numpart_mpi, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
1199        if (mp_ierr /= 0) goto 600
1200
1201! int2
1202        call MPI_GATHER(cbt, numpart_mpi, MPI_INTEGER2, cbt, &
1203             &numpart_mpi, MPI_INTEGER2, id_root, mp_comm_used, mp_ierr)
1204        if (mp_ierr /= 0) goto 600
1205
1206! Reals:
1207        call MPI_GATHER(uap, numpart_mpi, mp_sp, uap, &
1208             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1209        if (mp_ierr /= 0) goto 600
1210        call MPI_GATHER(ucp, numpart_mpi, mp_sp, ucp, &
1211             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1212        if (mp_ierr /= 0) goto 600
1213        call MPI_GATHER(uzp, numpart_mpi, mp_sp, uzp, &
1214             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1215        if (mp_ierr /= 0) goto 600
1216
1217        call MPI_GATHER(us, numpart_mpi, mp_sp, us, &
1218             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1219        if (mp_ierr /= 0) goto 600
1220        call MPI_GATHER(vs, numpart_mpi, mp_sp, vs, &
1221             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1222        if (mp_ierr /= 0) goto 600
1223        call MPI_GATHER(ws, numpart_mpi, mp_sp, ws, &
1224             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1225        if (mp_ierr /= 0) goto 600
1226
1227
1228        call MPI_GATHER(xtra1, numpart_mpi, mp_dp, xtra1, &
1229             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
1230        if (mp_ierr /= 0) goto 600
1231        call MPI_GATHER(ytra1, numpart_mpi, mp_dp, ytra1, &
1232             &numpart_mpi, mp_dp, id_root, mp_comm_used, mp_ierr)
1233        if (mp_ierr /= 0) goto 600
1234        call MPI_GATHER(ztra1, numpart_mpi, mp_sp, ztra1, &
1235             &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1236        if (mp_ierr /= 0) goto 600
1237
1238        do i=1, nspec
1239          call MPI_GATHER(xmass1(:,i), numpart_mpi, mp_sp, xmass1(:,i), &
1240               &numpart_mpi, mp_sp, id_root, mp_comm_used, mp_ierr)
1241          if (mp_ierr /= 0) goto 600
1242        end do
1243      end if ! (lroot)
1244    end if !  (use_in_place)
1245
1246    if (mp_measure_time) call mpif_mtime('commtime',1)
1247
1248    goto 601
1249
1250600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
1251    stop
1252601 end subroutine mpif_tm_collect_vars
1253
1254
1255  subroutine mpif_gf_send_vars(memstat)
1256!*******************************************************************************
1257! DESCRIPTION
1258!   Distribute 'getfield' variables from reader process
1259!
1260!   Called from timemanager
1261!
1262! NOTE
1263!   This subroutine distributes windfields read from the reader process to
1264!   all other processes. Usually one set of fields are transfered, but at
1265!   first timestep when there are no fields in memory, both are transfered.
1266!   MPI_Bcast is used, so implicitly all processes are synchronized at this
1267!   step
1268!
1269!
1270!*******************************************************************************
1271    use com_mod
1272    use par_mod,only: numwfmem
1273
1274    implicit none
1275
1276    integer, intent(in) :: memstat
1277
1278! Common array sizes used for communications
1279    integer, parameter :: d3_size1 = nxmax*nymax*nzmax
1280    integer, parameter :: d3_size2 = nxmax*nymax*nuvzmax
1281    integer, parameter :: d2_size1 = nxmax*nymax
1282    integer, parameter :: d2_size2 = nxmax*nymax*maxspec
1283    integer, parameter :: d2_size3 = nxmax*nymax
1284    integer, parameter :: d1_size1 = maxwf
1285
1286    integer :: d3s1,d3s2,d2s1,d2s2
1287
1288! li,ui    lower,upper indices of fields to be transfered (1-3, ui>=li)
1289    integer :: li,ui
1290
1291! First time routine is called the unchangeable fields will be transfered   
1292    logical :: first_call=.true.
1293
1294!**********************************************************************
1295
1296! Sizes of arrays transfered
1297    d3s1=d3_size1
1298    d3s2=d3_size2
1299    d2s1=d2_size1
1300    d2s2=d2_size2
1301
1302! Decide which field(s) need to be transfered
1303    if (memstat.ge.32) then ! distribute 2 fields, to 1st/2nd indices
1304      li=1; ui=2
1305      d3s1=2*d3_size1
1306      d3s2=2*d3_size2
1307      d2s1=2*d2_size1
1308      d2s2=2*d2_size2
1309    else if (memstat.eq.17) then ! distribute 1 field, place on 1st index
1310      li=1; ui=1
1311    else if (memstat.eq.18) then ! distribute 1 field, place on 2nd index
1312      li=2; ui=2
1313    else if (memstat.eq.19) then ! distribute 1 field, place on 3nd index
1314      li=3; ui=3
1315    else
1316      write(*,*) '#### mpi_mod::mpif_gf_send_vars> ERROR: ', &
1317           & 'wrong value of memstat, exiting ####', memstat
1318      stop
1319    end if
1320
1321
1322! Time for MPI communications
1323    if (mp_measure_time) call mpif_mtime('commtime',0)
1324
1325! Send variables from getfield process (id_read) to other processes
1326!**********************************************************************
1327
1328! The non-reader processes need to know if cloud water was read.
1329    call MPI_Bcast(readclouds,1,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
1330    if (mp_ierr /= 0) goto 600
1331
1332! Static fields/variables sent only at startup
1333    if (first_call) then
1334      call MPI_Bcast(oro(:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1335      if (mp_ierr /= 0) goto 600
1336      call MPI_Bcast(excessoro(:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1337      if (mp_ierr /= 0) goto 600
1338      call MPI_Bcast(lsm(:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1339      if (mp_ierr /= 0) goto 600
1340      call MPI_Bcast(xlanduse(:,:,:),d2_size3*numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1341      if (mp_ierr /= 0) goto 600
1342      call MPI_Bcast(wftime,d1_size1,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
1343      if (mp_ierr /= 0) goto 600
1344      call MPI_Bcast(numbwf,1,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
1345      if (mp_ierr /= 0) goto 600
1346      call MPI_Bcast(nmixz,1,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
1347      if (mp_ierr /= 0) goto 600
1348      call MPI_Bcast(height,nzmax,MPI_INTEGER,id_read,MPI_COMM_WORLD,mp_ierr)
1349      if (mp_ierr /= 0) goto 600
1350
1351      first_call=.false.
1352    endif
1353
1354    call MPI_Bcast(uu(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1355    if (mp_ierr /= 0) goto 600
1356    call MPI_Bcast(vv(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1357    if (mp_ierr /= 0) goto 600
1358    call MPI_Bcast(uupol(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1359    if (mp_ierr /= 0) goto 600
1360    call MPI_Bcast(vvpol(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1361    if (mp_ierr /= 0) goto 600
1362    call MPI_Bcast(ww(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1363    if (mp_ierr /= 0) goto 600
1364    call MPI_Bcast(tt(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1365    if (mp_ierr /= 0) goto 600
1366    call MPI_Bcast(rho(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1367    if (mp_ierr /= 0) goto 600
1368    call MPI_Bcast(drhodz(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1369    if (mp_ierr /= 0) goto 600
1370    call MPI_Bcast(tth(:,:,:,li:ui),d3s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1371    if (mp_ierr /= 0) goto 600
1372    call MPI_Bcast(qvh(:,:,:,li:ui),d3s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1373    if (mp_ierr /= 0) goto 600
1374    call MPI_Bcast(qv(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1375    if (mp_ierr /= 0) goto 600
1376    call MPI_Bcast(pv(:,:,:,li:ui),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1377    if (mp_ierr /= 0) goto 600
1378    call MPI_Bcast(clouds(:,:,:,li:ui),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
1379    if (mp_ierr /= 0) goto 600
1380
1381! cloud water/ice:
1382    if (readclouds) then
1383      call MPI_Bcast(ctwc(:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1384      if (mp_ierr /= 0) goto 600
1385    end if
1386
1387! 2D fields
1388    call MPI_Bcast(cloudsh(:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1389    if (mp_ierr /= 0) goto 600
1390    call MPI_Bcast(vdep(:,:,:,li:ui),d2s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1391    if (mp_ierr /= 0) goto 600
1392    call MPI_Bcast(ps(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1393    if (mp_ierr /= 0) goto 600
1394    call MPI_Bcast(sd(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1395    if (mp_ierr /= 0) goto 600
1396    call MPI_Bcast(tcc(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1397    if (mp_ierr /= 0) goto 600
1398    call MPI_Bcast(tt2(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1399    if (mp_ierr /= 0) goto 600
1400    call MPI_Bcast(td2(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1401    if (mp_ierr /= 0) goto 600
1402    call MPI_Bcast(lsprec(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1403    if (mp_ierr /= 0) goto 600
1404    call MPI_Bcast(convprec(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1405    if (mp_ierr /= 0) goto 600
1406    call MPI_Bcast(ustar(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1407    if (mp_ierr /= 0) goto 600
1408    call MPI_Bcast(wstar(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1409    if (mp_ierr /= 0) goto 600
1410    call MPI_Bcast(hmix(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1411    if (mp_ierr /= 0) goto 600
1412    call MPI_Bcast(tropopause(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1413    if (mp_ierr /= 0) goto 600
1414    call MPI_Bcast(oli(:,:,:,li:ui),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1415    if (mp_ierr /= 0) goto 600
1416
1417    if (mp_measure_time) call mpif_mtime('commtime',1)
1418
1419    goto 601
1420
1421600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
1422    stop
1423
1424601 end subroutine mpif_gf_send_vars
1425
1426
1427  subroutine mpif_gf_send_vars_nest(memstat)
1428!***********************************************************************
1429! DESCRIPTION
1430!   Distribute 'getfield' variables from reader process to all processes.
1431!   For nested fields
1432!
1433!   Called from timemanager
1434!
1435! NOTE
1436!   This subroutine distributes nested windfields from the reader process to
1437!   all other processes. Usually one set of fields is transfered, but at
1438!   first timestep when there are no fields in memory, both are transfered.
1439!   MPI_Bcast is used, so implicitly all processes are synchronized at this
1440!   step
1441!
1442!
1443!***********************************************************************
1444    use com_mod
1445    use par_mod,only: numwfmem
1446
1447    implicit none
1448
1449    integer, intent(in) :: memstat
1450    integer :: i
1451! li,ui    lower,upper indices of fields to be transfered (1-3, ui>=li)
1452    integer :: li,ui
1453
1454! Common array sizes used for communications
1455    integer :: d3_size1 = nxmaxn*nymaxn*nzmax
1456    integer :: d3_size2 = nxmaxn*nymaxn*nuvzmax
1457    integer :: d2_size1 = nxmaxn*nymaxn
1458    integer :: d2_size2 = nxmaxn*nymaxn*maxspec
1459    integer :: d2_size3 = nxmaxn*nymaxn
1460
1461    integer :: d3s1,d3s2,d2s1,d2s2
1462
1463! First time routine is called the unchangeable fields will be transfered   
1464    logical :: first_call=.true.
1465
1466!**********************************************************************
1467
1468! Sizes of arrays transfered
1469    d3s1=d3_size1
1470    d3s2=d3_size2
1471    d2s1=d2_size1
1472    d2s2=d2_size2
1473
1474! Decide which field(s) need to be transfered
1475    if (memstat.ge.32) then ! distribute 2 fields
1476      li=1; ui=2
1477      d3s1=2*d3_size1
1478      d3s2=2*d3_size2
1479      d2s1=2*d2_size1
1480      d2s2=2*d2_size2
1481    else if (memstat.eq.17) then ! distribute 1 field, on 1st index
1482      li=1; ui=1
1483    else if (memstat.eq.18) then ! distribute 1 field, on 2nd index
1484      li=2; ui=2
1485    else if (memstat.eq.19) then ! distribute 1 field, on 3nd index
1486      li=3; ui=3
1487    else
1488      write(*,*) '#### mpi_mod::mpif_gf_send_vars_nest> ERROR: ', &
1489           & 'wrong value of memstat, exiting ####', memstat
1490      stop
1491    end if
1492
1493
1494! Time for MPI communications
1495    if (mp_measure_time) call mpif_mtime('commtime',0)
1496
1497! Distribute variables, send from getfield process (id_read) to other
1498! processes
1499!**********************************************************************
1500
1501! The non-reader processes need to know if cloud water was read.
1502    call MPI_Bcast(readclouds_nest,maxnests,MPI_LOGICAL,id_read,MPI_COMM_WORLD,mp_ierr)
1503    if (mp_ierr /= 0) goto 600
1504
1505! Static fields/variables sent only at startup
1506    if (first_call) then
1507      call MPI_Bcast(oron(:,:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1508      if (mp_ierr /= 0) goto 600
1509      call MPI_Bcast(excessoron(:,:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1510      if (mp_ierr /= 0) goto 600
1511      call MPI_Bcast(lsmn(:,:,:),d2_size3,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1512      if (mp_ierr /= 0) goto 600
1513      call MPI_Bcast(xlandusen(:,:,:,:),d2_size3*numclass,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1514      if (mp_ierr /= 0) goto 600
1515      first_call=.false.
1516    end if
1517
1518! MPI prefers contiguous arrays for sending (else a buffer is created),
1519! hence the loop over nests
1520!**********************************************************************
1521    do i=1, numbnests
1522! 3D fields
1523      call MPI_Bcast(uun(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1524      if (mp_ierr /= 0) goto 600
1525      call MPI_Bcast(vvn(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1526      if (mp_ierr /= 0) goto 600
1527      call MPI_Bcast(wwn(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1528      if (mp_ierr /= 0) goto 600
1529      call MPI_Bcast(ttn(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1530      if (mp_ierr /= 0) goto 600
1531      call MPI_Bcast(rhon(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1532      if (mp_ierr /= 0) goto 600
1533      call MPI_Bcast(drhodzn(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1534      if (mp_ierr /= 0) goto 600
1535      call MPI_Bcast(tthn(:,:,:,li:ui,i),d3s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1536      if (mp_ierr /= 0) goto 600
1537      call MPI_Bcast(qvhn(:,:,:,li:ui,i),d3s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1538      if (mp_ierr /= 0) goto 600
1539      call MPI_Bcast(qvn(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1540      if (mp_ierr /= 0) goto 600
1541      call MPI_Bcast(pvn(:,:,:,li:ui,i),d3s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1542      if (mp_ierr /= 0) goto 600
1543      call MPI_Bcast(cloudsn(:,:,:,li:ui,i),d3s1,MPI_INTEGER1,id_read,MPI_COMM_WORLD,mp_ierr)
1544      if (mp_ierr /= 0) goto 600
1545
1546! cloud water/ice:
1547      if (readclouds_nest(i)) then
1548! call MPI_Bcast(icloud_stats(:,:,:,li:ui),d2s1*5,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1549! if (mp_ierr /= 0) goto 600
1550        call MPI_Bcast(ctwcn(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1551        if (mp_ierr /= 0) goto 600
1552      end if
1553
1554! 2D fields
1555      call MPI_Bcast(cloudshn(:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1556      if (mp_ierr /= 0) goto 600
1557      call MPI_Bcast(vdepn(:,:,:,li:ui,i),d2s2,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1558      if (mp_ierr /= 0) goto 600
1559      call MPI_Bcast(psn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1560      if (mp_ierr /= 0) goto 600
1561      call MPI_Bcast(sdn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1562      if (mp_ierr /= 0) goto 600
1563      call MPI_Bcast(tccn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1564      if (mp_ierr /= 0) goto 600
1565      call MPI_Bcast(tt2n(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1566      if (mp_ierr /= 0) goto 600
1567      call MPI_Bcast(td2n(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1568      if (mp_ierr /= 0) goto 600
1569      call MPI_Bcast(lsprecn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1570      if (mp_ierr /= 0) goto 600
1571      call MPI_Bcast(convprecn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1572      if (mp_ierr /= 0) goto 600
1573      call MPI_Bcast(ustarn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1574      if (mp_ierr /= 0) goto 600
1575      call MPI_Bcast(wstarn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1576      if (mp_ierr /= 0) goto 600
1577      call MPI_Bcast(olin(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1578      if (mp_ierr /= 0) goto 600
1579      call MPI_Bcast(hmixn(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1580      if (mp_ierr /= 0) goto 600
1581      call MPI_Bcast(tropopausen(:,:,:,li:ui,i),d2s1,mp_sp,id_read,MPI_COMM_WORLD,mp_ierr)
1582      if (mp_ierr /= 0) goto 600
1583    end do
1584
1585    if (mp_measure_time) call mpif_mtime('commtime',1)
1586
1587    goto 601
1588
1589600 write(*,*) "mpi_mod> mp_ierr \= 0",mp_ierr
1590    stop
1591
1592601 end subroutine mpif_gf_send_vars_nest
1593
1594
1595  subroutine mpif_gf_send_vars_async(memstat)
1596!*******************************************************************************
1597! DESCRIPTION
1598!   Distribute 'getfield' variables from reader process to all processes.
1599!   Called from timemanager by reader process only
1600!
1601! NOTE
1602!   This version uses asynchronious sends. The newest fields are sent in the
1603!   background, so calculations can continue while
1604!   MPI communications are performed.
1605!
1606!   The way the MPI tags/requests are sequenced means that this routine must
1607!   carefully match routine 'mpif_gf_recv_vars_async'
1608!
1609! VARIABLES
1610!   memstat -- input, for resolving pointer to windfield index being read
1611!   mind    -- index where to place new fields
1612!
1613!
1614!*******************************************************************************
1615    use com_mod
1616
1617    implicit none
1618
1619    integer, intent(in) :: memstat
1620    integer :: mind
1621    integer :: dest,i
1622
1623! Common array sizes used for communications
1624    integer :: d3s1 = nxmax*nymax*nzmax
1625    integer :: d3s2 = nxmax*nymax*nuvzmax
1626    integer :: d2s1 = nxmax*nymax
1627    integer :: d2s2 = nxmax*nymax*maxspec
1628!    integer :: d1s1 = maxwf
1629
1630!*******************************************************************************
1631
1632! At the time the send is posted, the reader process is one step further
1633! in the permutation of memstat compared with the receiving processes
1634
1635    if (memstat.ge.32) then
1636! last read was synchronous, to indices 1 and 2, 3 is free
1637      write(*,*) "#### mpi_mod::mpif_gf_send_vars_async> ERROR: &
1638           & memstat>=32 should never happen here."
1639      stop
1640    else if (memstat.eq.17) then
1641! old fields on 1,2, send 3
1642      mind=3
1643    else if (memstat.eq.18) then
1644! old fields on 2,3, send 1
1645      mind=1
1646    else if (memstat.eq.19) then
1647! old fields on 3,1, send 2
1648      mind=2
1649    else
1650      write(*,*) "#### mpi_mod::mpif_gf_send_vars_async> ERROR: &
1651           & invalid memstat"
1652    end if
1653
1654    if (mp_dev_mode) then
1655      if (mp_pid.ne.id_read) then
1656        write(*,*) 'MPI_DEV_MODE: error calling mpif_gf_send_vars_async'
1657      end if
1658    end if
1659
1660    if (mp_dev_mode) write(*,*) '## in mpif_gf_send_vars_async, memstat:', memstat
1661
1662! Time for MPI communications
1663    if (mp_measure_time) call mpif_mtime('commtime',0)
1664
1665! Loop over receiving processes, initiate data sending
1666!*****************************************************
1667
1668    do dest=0,mp_np-1 ! mp_np-2 will also work if last proc reserved for reading
1669! TODO: use mp_partgroup_np here
1670      if (dest.eq.id_read) cycle
1671      i=dest*nvar_async
1672      call MPI_Isend(uu(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1673      if (mp_ierr /= 0) goto 600
1674      i=i+1
1675      call MPI_Isend(vv(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1676      if (mp_ierr /= 0) goto 600
1677      i=i+1
1678      call MPI_Isend(uupol(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1679      if (mp_ierr /= 0) goto 600
1680      i=i+1
1681      call MPI_Isend(vvpol(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1682      if (mp_ierr /= 0) goto 600
1683      i=i+1
1684      call MPI_Isend(ww(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1685      if (mp_ierr /= 0) goto 600
1686      i=i+1
1687      call MPI_Isend(tt(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1688      if (mp_ierr /= 0) goto 600
1689      i=i+1
1690      call MPI_Isend(rho(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1691      if (mp_ierr /= 0) goto 600
1692      i=i+1
1693      call MPI_Isend(drhodz(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1694      if (mp_ierr /= 0) goto 600
1695      i=i+1
1696      call MPI_Isend(tth(:,:,:,mind),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1697      if (mp_ierr /= 0) goto 600
1698      i=i+1
1699      call MPI_Isend(qvh(:,:,:,mind),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1700      if (mp_ierr /= 0) goto 600
1701      i=i+1
1702      call MPI_Isend(qv(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1703      if (mp_ierr /= 0) goto 600
1704      i=i+1
1705      call MPI_Isend(pv(:,:,:,mind),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1706      if (mp_ierr /= 0) goto 600
1707      i=i+1
1708      call MPI_Isend(clouds(:,:,:,mind),d3s1,MPI_INTEGER1,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1709      i=i+1
1710      if (mp_ierr /= 0) goto 600
1711      call MPI_Isend(cloudsh(:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1712      if (mp_ierr /= 0) goto 600
1713      i=i+1
1714      call MPI_Isend(vdep(:,:,:,mind),d2s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1715      if (mp_ierr /= 0) goto 600
1716      i=i+1
1717! 15
1718      call MPI_Isend(ps(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1719      if (mp_ierr /= 0) goto 600
1720      i=i+1
1721      call MPI_Isend(sd(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1722      if (mp_ierr /= 0) goto 600
1723      i=i+1
1724      call MPI_Isend(tcc(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1725      if (mp_ierr /= 0) goto 600
1726      i=i+1
1727      call MPI_Isend(tt2(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1728      if (mp_ierr /= 0) goto 600
1729      i=i+1
1730      call MPI_Isend(td2(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1731      if (mp_ierr /= 0) goto 600
1732      i=i+1
1733      call MPI_Isend(lsprec(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1734      if (mp_ierr /= 0) goto 600
1735      i=i+1
1736      call MPI_Isend(convprec(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1737      if (mp_ierr /= 0) goto 600
1738      i=i+1
1739      call MPI_Isend(ustar(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1740      if (mp_ierr /= 0) goto 600
1741      i=i+1
1742      call MPI_Isend(wstar(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1743      if (mp_ierr /= 0) goto 600
1744      i=i+1
1745      call MPI_Isend(hmix(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1746      if (mp_ierr /= 0) goto 600
1747      i=i+1
1748! 25
1749      call MPI_Isend(tropopause(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1750      if (mp_ierr /= 0) goto 600
1751      i=i+1
1752      call MPI_Isend(oli(:,:,:,mind),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
1753      if (mp_ierr /= 0) goto 600
1754
1755! Send cloud water if it exists. Increment counter always (as on receiving end)
1756      if (readclouds) then
1757        i=i+1
1758        call MPI_Isend(ctwc(:,:,mind),d2s1,mp_sp,dest,tm1,&
1759             &MPI_COMM_WORLD,reqs(i),mp_ierr)
1760        if (mp_ierr /= 0) goto 600
1761      end if
1762    end do
1763
1764    if (mp_measure_time) call mpif_mtime('commtime',1)
1765
1766    goto 601
1767
1768600 write(*,*) "#### mpi_mod::mpif_gf_send_vars_async> mp_ierr \= 0", mp_ierr
1769    stop
1770
1771601 end subroutine mpif_gf_send_vars_async
1772
1773
1774  subroutine mpif_gf_recv_vars_async(memstat)
1775!*******************************************************************************
1776! DESCRIPTION
1777!   Receive 'getfield' variables from reader process.
1778!   Called from timemanager by all processes except reader
1779!
1780! NOTE
1781!   This version uses asynchronious communications.
1782!
1783! VARIABLES
1784!   memstat -- input, used to resolve windfield index being received
1785!
1786!
1787!*******************************************************************************
1788    use com_mod
1789
1790    implicit none
1791
1792    integer, intent(in) :: memstat
1793    integer :: mind,j
1794
1795! Common array sizes used for communications
1796    integer :: d3s1 = nxmax*nymax*nzmax
1797    integer :: d3s2 = nxmax*nymax*nuvzmax
1798    integer :: d2s1 = nxmax*nymax
1799    integer :: d2s2 = nxmax*nymax*maxspec
1800!integer :: d1_size1 = maxwf
1801
1802!    integer :: d3s1,d3s2,d2s1,d2s2
1803!*******************************************************************************
1804
1805! At the time this immediate receive is posted, memstat is the state of
1806! windfield indices at the previous/current time. From this, the future
1807! state is deduced.
1808    if (memstat.eq.32) then
1809! last read was synchronous to indices 1 and 2, 3 is free
1810      mind=3
1811    else if (memstat.eq.17) then
1812! last read was asynchronous to index 3, 1 is free
1813      mind=1
1814    else if (memstat.eq.18) then
1815! last read was asynchronous to index 1, 2 is free
1816      mind=2
1817    else if (memstat.eq.19) then
1818! last read was asynchronous to index 2, 3 is free
1819      mind=3
1820    else
1821! illegal state
1822      write(*,*) 'mpi_mod> FLEXPART ERROR: Illegal memstat value. Exiting.'
1823      stop
1824    end if
1825
1826    if (mp_dev_mode) then
1827      if (mp_pid.eq.id_read) then
1828        write(*,*) 'MPI_DEV_MODE: error calling mpif_gf_recv_vars_async'
1829      end if
1830    end if
1831
1832! Time for MPI communications
1833    if (mp_measure_time) call mpif_mtime('commtime',0)
1834
1835    if (mp_dev_mode) write(*,*) '## in mpif_gf_send_vars_async, memstat:', memstat
1836
1837! Initiate receiving of data
1838!***************************
1839
1840! Get MPI tags/requests for communications
1841    j=mp_pid*nvar_async
1842    call MPI_Irecv(uu(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1843         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1844    if (mp_ierr /= 0) goto 600
1845    j=j+1
1846    call MPI_Irecv(vv(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1847         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1848    if (mp_ierr /= 0) goto 600
1849    j=j+1
1850    call MPI_Irecv(uupol(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1851         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1852    if (mp_ierr /= 0) goto 600
1853    j=j+1
1854    call MPI_Irecv(vvpol(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1855         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1856    if (mp_ierr /= 0) goto 600
1857    j=j+1
1858    call MPI_Irecv(ww(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1859         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1860    if (mp_ierr /= 0) goto 600
1861    j=j+1
1862    call MPI_Irecv(tt(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1863         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1864    if (mp_ierr /= 0) goto 600
1865    j=j+1
1866    call MPI_Irecv(rho(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1867         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1868    if (mp_ierr /= 0) goto 600
1869    j=j+1
1870    call MPI_Irecv(drhodz(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1871         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1872    if (mp_ierr /= 0) goto 600
1873    j=j+1
1874    call MPI_Irecv(tth(:,:,:,mind),d3s2,mp_sp,id_read,MPI_ANY_TAG,&
1875         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1876    if (mp_ierr /= 0) goto 600
1877    j=j+1
1878    call MPI_Irecv(qvh(:,:,:,mind),d3s2,mp_sp,id_read,MPI_ANY_TAG,&
1879         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1880    if (mp_ierr /= 0) goto 600
1881    j=j+1
1882    call MPI_Irecv(qv(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1883         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1884    if (mp_ierr /= 0) goto 600
1885    j=j+1
1886    call MPI_Irecv(pv(:,:,:,mind),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
1887         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1888    if (mp_ierr /= 0) goto 600
1889    j=j+1
1890    call MPI_Irecv(clouds(:,:,:,mind),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,&
1891         &MPI_COMM_WORLD,reqs(j),mp_ierr)   
1892    if (mp_ierr /= 0) goto 600
1893    j=j+1
1894    call MPI_Irecv(cloudsh(:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1895         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1896    if (mp_ierr /= 0) goto 600
1897    j=j+1
1898    call MPI_Irecv(vdep(:,:,:,mind),d2s2,mp_sp,id_read,MPI_ANY_TAG,&
1899         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1900    if (mp_ierr /= 0) goto 600
1901    j=j+1
1902    call MPI_Irecv(ps(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1903         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1904    if (mp_ierr /= 0) goto 600
1905    j=j+1
1906    call MPI_Irecv(sd(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1907         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1908    if (mp_ierr /= 0) goto 600
1909    j=j+1
1910    call MPI_Irecv(tcc(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1911         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1912    if (mp_ierr /= 0) goto 600
1913    j=j+1
1914    call MPI_Irecv(tt2(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1915         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1916    if (mp_ierr /= 0) goto 600
1917    j=j+1
1918    call MPI_Irecv(td2(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1919         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1920    if (mp_ierr /= 0) goto 600
1921    j=j+1
1922    call MPI_Irecv(lsprec(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1923         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1924    if (mp_ierr /= 0) goto 600
1925    j=j+1
1926    call MPI_Irecv(convprec(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1927         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1928    if (mp_ierr /= 0) goto 600
1929    call MPI_Irecv(ustar(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1930         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1931    if (mp_ierr /= 0) goto 600
1932    j=j+1
1933    call MPI_Irecv(wstar(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1934         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1935    if (mp_ierr /= 0) goto 600
1936    j=j+1
1937    call MPI_Irecv(hmix(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1938         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1939    if (mp_ierr /= 0) goto 600
1940    j=j+1
1941    call MPI_Irecv(tropopause(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1942         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1943    if (mp_ierr /= 0) goto 600
1944    j=j+1
1945    call MPI_Irecv(oli(:,:,:,mind),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
1946         &MPI_COMM_WORLD,reqs(j),mp_ierr)
1947    if (mp_ierr /= 0) goto 600
1948
1949! Post request for clwc. These data are possibly not sent, request must then be cancelled
1950! For now assume that data at all steps either have or do not have water
1951    if (readclouds) then
1952      j=j+1
1953      call MPI_Irecv(ctwc(:,:,mind),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
1954           &MPI_COMM_WORLD,reqs(j),mp_ierr)
1955      if (mp_ierr /= 0) goto 600
1956    end if
1957
1958
1959    if (mp_measure_time) call mpif_mtime('commtime',1)
1960
1961    goto 601
1962
1963600 write(*,*) "#### mpi_mod::mpif_gf_recv_vars_async> MPI ERROR ", mp_ierr
1964    stop
1965
1966601 end subroutine mpif_gf_recv_vars_async
1967
1968
1969  subroutine mpif_gf_send_vars_nest_async(memstat)
1970!*******************************************************************************
1971! DESCRIPTION
1972!   Distribute 'getfield' variables from reader process to all processes.
1973!   Called from timemanager by reader process only.
1974!   Version for nested wind fields
1975!
1976! NOTE
1977!   This version uses asynchronious sends. The newest fields are sent in the
1978!   background, so calculations can continue while
1979!   MPI communications are performed.
1980!
1981!   The way the MPI tags/requests are sequenced means that this routine must
1982!   carefully match routine 'mpif_gf_recv_vars_async'
1983!
1984! VARIABLES
1985!   memstat -- input, for resolving pointer to windfield index being read
1986!   mind    -- index where to place new fields
1987!
1988! TODO
1989!   Some unused arrays are currently sent (uupoln,..)
1990!*******************************************************************************
1991    use com_mod
1992
1993    implicit none
1994
1995    integer, intent(in) :: memstat
1996    integer :: mind
1997    integer :: dest,i,k
1998
1999! Common array sizes used for communications
2000    integer :: d3s1 = nxmaxn*nymaxn*nzmax
2001    integer :: d3s2 = nxmaxn*nymaxn*nuvzmax
2002    integer :: d2s1 = nxmaxn*nymaxn
2003    integer :: d2s2 = nxmaxn*nymaxn*maxspec
2004
2005!*******************************************************************************
2006
2007! At the time the send is posted, the reader process is one step further
2008! in the permutation of memstat compared with the receiving processes
2009
2010    if (memstat.ge.32) then
2011! last read was synchronous, to indices 1 and 2, 3 is free
2012      write(*,*) "#### mpi_mod::mpif_gf_send_vars_nest_async> ERROR: &
2013           & memstat>=32 should never happen here."
2014      stop
2015    else if (memstat.eq.17) then
2016! old fields on 1,2, send 3
2017      mind=3
2018    else if (memstat.eq.18) then
2019! old fields on 2,3, send 1
2020      mind=1
2021    else if (memstat.eq.19) then
2022! old fields on 3,1, send 2
2023      mind=2
2024    else
2025      write(*,*) "#### mpi_mod::mpif_gf_send_vars_nest_async> ERROR: &
2026           & invalid memstat"
2027    end if
2028
2029    if (mp_dev_mode) then
2030      if (mp_pid.ne.id_read) then
2031        write(*,*) 'MPI_DEV_MODE: error calling mpif_gf_send_vars_nest_async'
2032      end if
2033    end if
2034
2035    if (mp_dev_mode) write(*,*) '## in mpif_gf_send_vars_nest_async, memstat:', memstat
2036
2037! Time for MPI communications
2038    if (mp_measure_time) call mpif_mtime('commtime',0)
2039
2040! Loop over receiving processes, initiate data sending
2041!*****************************************************
2042
2043    do dest=0,mp_np-1 ! mp_np-2 will also work if last proc reserved for reading
2044! TODO: use mp_partgroup_np here
2045      if (dest.eq.id_read) cycle
2046      do k=1, numbnests
2047        i=dest*nvar_async
2048        call MPI_Isend(uun(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2049        if (mp_ierr /= 0) goto 600
2050        i=i+1
2051        call MPI_Isend(vvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2052        if (mp_ierr /= 0) goto 600
2053        i=i+1
2054        call MPI_Isend(wwn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2055        if (mp_ierr /= 0) goto 600
2056        i=i+1
2057        call MPI_Isend(ttn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2058        if (mp_ierr /= 0) goto 600
2059        i=i+1
2060        call MPI_Isend(rhon(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2061        if (mp_ierr /= 0) goto 600
2062        i=i+1
2063        call MPI_Isend(drhodzn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2064        if (mp_ierr /= 0) goto 600
2065        i=i+1
2066        call MPI_Isend(tthn(:,:,:,mind,k),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2067        if (mp_ierr /= 0) goto 600
2068        i=i+1
2069        call MPI_Isend(qvhn(:,:,:,mind,k),d3s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2070        if (mp_ierr /= 0) goto 600
2071        i=i+1
2072        call MPI_Isend(qvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2073        if (mp_ierr /= 0) goto 600
2074        i=i+1
2075        call MPI_Isend(pvn(:,:,:,mind,k),d3s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2076        if (mp_ierr /= 0) goto 600
2077        i=i+1
2078        call MPI_Isend(cloudsn(:,:,:,mind,k),d3s1,MPI_INTEGER1,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2079        i=i+1
2080        if (mp_ierr /= 0) goto 600
2081        call MPI_Isend(cloudshn(:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2082        if (mp_ierr /= 0) goto 600
2083        i=i+1
2084        call MPI_Isend(vdepn(:,:,:,mind,k),d2s2,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2085        if (mp_ierr /= 0) goto 600
2086        i=i+1
2087        call MPI_Isend(psn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2088        if (mp_ierr /= 0) goto 600
2089        i=i+1
2090        call MPI_Isend(sdn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2091        if (mp_ierr /= 0) goto 600
2092        i=i+1
2093! 15
2094        call MPI_Isend(tccn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2095        if (mp_ierr /= 0) goto 600
2096        i=i+1
2097        call MPI_Isend(tt2n(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2098        if (mp_ierr /= 0) goto 600
2099        i=i+1
2100        call MPI_Isend(td2n(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2101        if (mp_ierr /= 0) goto 600
2102        i=i+1
2103        call MPI_Isend(lsprecn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2104        if (mp_ierr /= 0) goto 600
2105        i=i+1
2106        call MPI_Isend(convprecn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2107        if (mp_ierr /= 0) goto 600
2108        i=i+1
2109        call MPI_Isend(ustarn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2110        if (mp_ierr /= 0) goto 600
2111        i=i+1
2112        call MPI_Isend(wstarn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2113        if (mp_ierr /= 0) goto 600
2114        i=i+1
2115        call MPI_Isend(hmixn(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2116        if (mp_ierr /= 0) goto 600
2117        i=i+1
2118        call MPI_Isend(tropopausen(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2119        if (mp_ierr /= 0) goto 600
2120        i=i+1
2121        call MPI_Isend(olin(:,:,:,mind,k),d2s1,mp_sp,dest,tm1,MPI_COMM_WORLD,reqs(i),mp_ierr)
2122        if (mp_ierr /= 0) goto 600
2123! 25
2124
2125! Send cloud water if it exists. Increment counter always (as on receiving end)
2126        if (readclouds) then
2127          i=i+1
2128          call MPI_Isend(ctwcn(:,:,mind,k),d2s1,mp_sp,dest,tm1,&
2129               &MPI_COMM_WORLD,reqs(i),mp_ierr)
2130          if (mp_ierr /= 0) goto 600
2131        end if
2132      end do
2133    end do
2134
2135    if (mp_measure_time) call mpif_mtime('commtime',1)
2136
2137    goto 601
2138
2139600 write(*,*) "#### mpi_mod::mpif_gf_send_vars_nest_async> mp_ierr \= 0", mp_ierr
2140    stop
2141
2142601 end subroutine mpif_gf_send_vars_nest_async
2143
2144
2145  subroutine mpif_gf_recv_vars_nest_async(memstat)
2146!*******************************************************************************
2147! DESCRIPTION
2148!   Receive 'getfield' variables from reader process.
2149!   Called from timemanager by all processes except reader
2150!   Version for nested wind fields
2151!
2152! NOTE
2153!   This version uses asynchronious communications.
2154!
2155! VARIABLES
2156!   memstat -- input, used to resolve windfield index being received
2157!
2158!
2159!*******************************************************************************
2160    use com_mod
2161
2162    implicit none
2163
2164    integer, intent(in) :: memstat
2165    integer :: mind,j,k
2166
2167! Common array sizes used for communications
2168    integer :: d3s1 = nxmaxn*nymaxn*nzmax
2169    integer :: d3s2 = nxmaxn*nymaxn*nuvzmax
2170    integer :: d2s1 = nxmaxn*nymaxn
2171    integer :: d2s2 = nxmaxn*nymaxn*maxspec
2172
2173!*******************************************************************************
2174
2175! At the time this immediate receive is posted, memstat is the state of
2176! windfield indices at the previous/current time. From this, the future
2177! state is deduced.
2178    if (memstat.eq.32) then
2179! last read was synchronous to indices 1 and 2, 3 is free
2180      mind=3
2181    else if (memstat.eq.17) then
2182! last read was asynchronous to index 3, 1 is free
2183      mind=1
2184    else if (memstat.eq.18) then
2185! last read was asynchronous to index 1, 2 is free
2186      mind=2
2187    else if (memstat.eq.19) then
2188! last read was asynchronous to index 2, 3 is free
2189      mind=3
2190    else
2191! illegal state
2192      write(*,*) 'mpi_mod> FLEXPART ERROR: Illegal memstat value. Exiting.'
2193      stop
2194    end if
2195
2196    if (mp_dev_mode) then
2197      if (mp_pid.eq.id_read) then
2198        write(*,*) 'MPI_DEV_MODE: error calling mpif_gf_recv_vars_async'
2199      end if
2200    end if
2201
2202! Time for MPI communications
2203    if (mp_measure_time) call mpif_mtime('commtime',0)
2204
2205    if (mp_dev_mode) write(*,*) '## in mpif_gf_send_vars_async, memstat:', memstat
2206
2207! Initiate receiving of data
2208!***************************
2209
2210    do k=1, numbnests
2211! Get MPI tags/requests for communications
2212      j=mp_pid*nvar_async
2213      call MPI_Irecv(uun(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2214           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2215      if (mp_ierr /= 0) goto 600
2216      j=j+1
2217      call MPI_Irecv(vvn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2218           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2219      if (mp_ierr /= 0) goto 600
2220      j=j+1
2221      call MPI_Irecv(wwn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2222           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2223      if (mp_ierr /= 0) goto 600
2224      j=j+1
2225      call MPI_Irecv(ttn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2226           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2227      if (mp_ierr /= 0) goto 600
2228      j=j+1
2229      call MPI_Irecv(rhon(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2230           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2231      if (mp_ierr /= 0) goto 600
2232      j=j+1
2233      call MPI_Irecv(drhodzn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2234           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2235      if (mp_ierr /= 0) goto 600
2236      j=j+1
2237      call MPI_Irecv(tthn(:,:,:,mind,k),d3s2,mp_sp,id_read,MPI_ANY_TAG,&
2238           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2239      if (mp_ierr /= 0) goto 600
2240      j=j+1
2241      call MPI_Irecv(qvhn(:,:,:,mind,k),d3s2,mp_sp,id_read,MPI_ANY_TAG,&
2242           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2243      if (mp_ierr /= 0) goto 600
2244      j=j+1
2245      call MPI_Irecv(qvn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2246           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2247      if (mp_ierr /= 0) goto 600
2248      j=j+1
2249      call MPI_Irecv(pvn(:,:,:,mind,k),d3s1,mp_sp,id_read,MPI_ANY_TAG,&
2250           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2251      if (mp_ierr /= 0) goto 600
2252      j=j+1
2253      call MPI_Irecv(cloudsn(:,:,:,mind,k),d3s1,MPI_INTEGER1,id_read,MPI_ANY_TAG,&
2254           &MPI_COMM_WORLD,reqs(j),mp_ierr)   
2255      if (mp_ierr /= 0) goto 600
2256      j=j+1
2257      call MPI_Irecv(cloudshn(:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2258           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2259      if (mp_ierr /= 0) goto 600
2260      j=j+1
2261      call MPI_Irecv(vdepn(:,:,:,mind,k),d2s2,mp_sp,id_read,MPI_ANY_TAG,&
2262           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2263      if (mp_ierr /= 0) goto 600
2264      j=j+1
2265      call MPI_Irecv(psn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2266           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2267      if (mp_ierr /= 0) goto 600
2268      j=j+1
2269      call MPI_Irecv(sdn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2270           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2271      if (mp_ierr /= 0) goto 600
2272      j=j+1
2273      call MPI_Irecv(tccn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2274           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2275      if (mp_ierr /= 0) goto 600
2276      j=j+1
2277      call MPI_Irecv(tt2n(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2278           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2279      if (mp_ierr /= 0) goto 600
2280      j=j+1
2281      call MPI_Irecv(td2n(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2282           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2283      if (mp_ierr /= 0) goto 600
2284      j=j+1
2285      call MPI_Irecv(lsprecn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2286           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2287      if (mp_ierr /= 0) goto 600
2288      j=j+1
2289      call MPI_Irecv(convprecn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2290           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2291      if (mp_ierr /= 0) goto 600
2292      call MPI_Irecv(ustarn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2293           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2294      if (mp_ierr /= 0) goto 600
2295      j=j+1
2296      call MPI_Irecv(wstarn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2297           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2298      if (mp_ierr /= 0) goto 600
2299      j=j+1
2300      call MPI_Irecv(hmixn(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2301           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2302      if (mp_ierr /= 0) goto 600
2303      j=j+1
2304      call MPI_Irecv(tropopausen(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2305           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2306      if (mp_ierr /= 0) goto 600
2307      j=j+1
2308      call MPI_Irecv(olin(:,:,:,mind,k),d2s1,mp_sp,id_read,MPI_ANY_TAG,&
2309           &MPI_COMM_WORLD,reqs(j),mp_ierr)
2310      if (mp_ierr /= 0) goto 600
2311
2312! Post request for clwc. These data are possibly not sent, request must then be cancelled
2313! For now assume that data at all steps either have or do not have water
2314      if (readclouds) then
2315        j=j+1
2316        call MPI_Irecv(ctwcn(:,:,mind,k),d2s1*5,mp_sp,id_read,MPI_ANY_TAG,&
2317             &MPI_COMM_WORLD,reqs(j),mp_ierr)
2318        if (mp_ierr /= 0) goto 600
2319      end if
2320    end do
2321
2322    if (mp_measure_time) call mpif_mtime('commtime',1)
2323
2324    goto 601
2325
2326600 write(*,*) "#### mpi_mod::mpif_gf_recv_vars_nest_async> MPI ERROR ", mp_ierr
2327    stop
2328
2329601 end subroutine mpif_gf_recv_vars_nest_async
2330
2331
2332  subroutine mpif_gf_request
2333!*******************************************************************************
2334! DESCRIPTION
2335!   Check for completion of MPI_Isend/MPI_Irecv data transfer.
2336!   
2337!   
2338! NOTE
2339!   implicit synchronisation between all processes takes place here
2340!
2341! TODO
2342!   take into account nested wind fields
2343!
2344!*******************************************************************************
2345!    use com_mod,only: readclouds
2346
2347    implicit none
2348
2349
2350    integer :: n_req !,i
2351
2352!***********************************************************************
2353
2354    n_req=nvar_async*mp_np
2355
2356    if (mp_measure_time) call mpif_mtime('commtime',0)
2357
2358!    call MPI_Wait(rm1,MPI_STATUS_IGNORE,mp_ierr)
2359
2360! TODO: cancel recv request if readclouds=.false.
2361!    if (readclouds) then
2362    call MPI_Waitall(n_req,reqs,MPI_STATUSES_IGNORE,mp_ierr)
2363!    endif
2364! else
2365!   do i = 0, nvar_async*mp_np-1
2366!     if (mod(i,27).eq.0 .or. mod(i,28).eq.0) then
2367!       call MPI_Cancel(reqs(i),mp_ierr)
2368!       cycle
2369!     end if
2370!     call MPI_Wait(reqs(i),MPI_STATUS_IGNORE,mp_ierr)
2371!   end do
2372! end if
2373
2374    if (mp_ierr /= 0) goto 600
2375
2376    if (mp_measure_time) call mpif_mtime('commtime',1)
2377
2378    goto 601
2379
2380600 write(*,*) "#### mpi_mod::mpif_gf_request> MPI ERROR ", mp_ierr
2381    stop
2382
2383601 end subroutine mpif_gf_request
2384
2385
2386  subroutine mpif_tm_reduce_grid
2387!***********************************************************************
2388! Collect grid variable to PID 0, adding from all processes.
2389!
2390! NOTE
2391!   - Older MPI implementations may lack the MPI_IN_PLACE operation.
2392!     If so use 1) below and change unc_mod
2393!
2394!***********************************************************************
2395    use com_mod
2396    use unc_mod
2397    use par_mod
2398
2399    implicit none
2400
2401    integer :: grid_size2d,grid_size3d
2402    integer, parameter :: rcpt_size=maxreceptor*maxspec
2403
2404!**********************************************************************
2405    grid_size2d=numxgrid*numygrid*maxspec* &
2406         & maxpointspec_act*nclassunc*maxageclass
2407    grid_size3d=numxgrid*numygrid*numzgrid*maxspec* &
2408         & maxpointspec_act*nclassunc*maxageclass
2409
2410
2411! Time for MPI communications
2412    if (mp_measure_time) call mpif_mtime('commtime',0)
2413
2414! 1) Using a separate grid (gridunc0) for received values
2415! call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2416!      & mp_comm_used, mp_ierr)
2417! if (mp_ierr /= 0) goto 600
2418
2419! 2) Using in-place reduction
2420    if (lroot) then
2421      call MPI_Reduce(MPI_IN_PLACE, gridunc, grid_size3d, mp_sp, MPI_SUM, id_root, &
2422           & mp_comm_used, mp_ierr)
2423      if (mp_ierr /= 0) goto 600
2424    else
2425      call MPI_Reduce(gridunc, 0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2426           & mp_comm_used, mp_ierr)
2427    end if
2428
2429    if ((WETDEP).and.(ldirect.gt.0)) then
2430      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2431           & mp_comm_used, mp_ierr)
2432      if (mp_ierr /= 0) goto 600
2433    end if
2434
2435    if ((DRYDEP).and.(ldirect.gt.0)) then
2436      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2437           & mp_comm_used, mp_ierr)
2438      if (mp_ierr /= 0) goto 600
2439    end if
2440
2441! Receptor concentrations   
2442    if (lroot) then
2443      call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
2444           & mp_comm_used,mp_ierr)
2445      if (mp_ierr /= 0) goto 600
2446    else
2447      call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &
2448           & mp_comm_used,mp_ierr)
2449    end if
2450
2451    if (mp_measure_time) call mpif_mtime('commtime',1)
2452
2453    goto 601
2454
2455600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
2456    stop
2457
2458601 end subroutine mpif_tm_reduce_grid
2459
2460
2461  subroutine mpif_tm_reduce_grid_nest
2462!***********************************************************************
2463! Collect nested grids to PID 0, adding from all processes.
2464!
2465! NOTE
2466!   - Compiling with 'fcheck=all' (gfortran) will cause a run-time error
2467!     as wetgriduncn0 is only allocated for root process. Possibly valid
2468!     MPI code but invalid gfortran.
2469!
2470!***********************************************************************
2471    use com_mod
2472    use unc_mod
2473    use par_mod
2474
2475    implicit none
2476
2477    integer :: grid_size2d,grid_size3d
2478
2479!**********************************************************************
2480
2481    grid_size3d=numxgridn*numygridn*numzgrid*maxspec* &
2482         & maxpointspec_act*nclassunc*maxageclass
2483    grid_size2d=numxgridn*numygridn*maxspec* &
2484         & maxpointspec_act*nclassunc*maxageclass
2485
2486
2487! Time for MPI communications
2488    if (mp_measure_time) call mpif_mtime('commtime',0)
2489
2490! Using a separate grid (gridunc0) for received values, for debugging
2491! call MPI_Reduce(griduncn, griduncn0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2492!      & mp_comm_used, mp_ierr)
2493! if (mp_ierr /= 0) goto 600
2494
2495! Using in-place reduction
2496    if (lroot) then
2497      call MPI_Reduce(MPI_IN_PLACE, griduncn, grid_size3d, mp_sp, MPI_SUM, id_root, &
2498           & mp_comm_used, mp_ierr)
2499      if (mp_ierr /= 0) goto 600
2500    else
2501      call MPI_Reduce(griduncn, 0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2502           & mp_comm_used, mp_ierr)
2503    end if
2504
2505    if ((WETDEP).and.(ldirect.gt.0)) then
2506      call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2507           & mp_comm_used, mp_ierr)
2508      if (mp_ierr /= 0) goto 600
2509    end if
2510
2511    if ((DRYDEP).and.(ldirect.gt.0)) then
2512      call MPI_Reduce(drygriduncn, drygriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2513           & mp_comm_used, mp_ierr)
2514      if (mp_ierr /= 0) goto 600
2515    end if
2516
2517    if (mp_measure_time) call mpif_mtime('commtime',1)
2518
2519    goto 601
2520
2521600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
2522    stop
2523
2524601 end subroutine mpif_tm_reduce_grid_nest
2525
2526
2527  subroutine mpif_mtime(ident,imode)
2528!***********************************************************************
2529! Measure CPU/Wall time in various parts of the code
2530!
2531! VARIABLES
2532!   ident        character, identifies routine to measure
2533!   imode        integer, 0:start clock(s)  1: stop clock(s)
2534!
2535!***********************************************************************
2536    implicit none
2537
2538    character(LEN=*), intent(in) :: ident
2539    integer, intent(in) :: imode
2540
2541!***********************************************************************
2542
2543    select case(ident)
2544
2545    case ('timemanager')
2546      if (imode.eq.0) then
2547        call cpu_time(tm_tot_beg)
2548      else
2549        call cpu_time(tm_tot_end)
2550        tm_tot_total = tm_tot_total + (tm_tot_end-tm_tot_beg)
2551      end if
2552
2553    case ('wetdepo')
2554      if (imode.eq.0) then
2555        mp_wetdepo_wtime_beg = mpi_wtime()
2556        call cpu_time(mp_wetdepo_time_beg)
2557      else
2558        mp_wetdepo_wtime_end = mpi_wtime()
2559        call cpu_time(mp_wetdepo_time_end)
2560
2561        mp_wetdepo_wtime_total = mp_wetdepo_wtime_total + &
2562             &(mp_wetdepo_wtime_end - mp_wetdepo_wtime_beg)
2563        mp_wetdepo_time_total = mp_wetdepo_time_total + &
2564             &(mp_wetdepo_time_end - mp_wetdepo_time_beg)
2565      end if
2566
2567    case ('advance')
2568      if (imode.eq.0) then
2569        mp_advance_wtime_beg = mpi_wtime()
2570      else
2571        mp_advance_wtime_end = mpi_wtime()
2572
2573        mp_advance_wtime_total = mp_advance_wtime_total + &
2574             &(mp_advance_wtime_end - mp_advance_wtime_beg)
2575      end if
2576
2577    case ('getfields')
2578      if (imode.eq.0) then
2579        mp_getfields_wtime_beg = mpi_wtime()
2580        call cpu_time(mp_getfields_time_beg)
2581      else
2582        mp_getfields_wtime_end = mpi_wtime()
2583        call cpu_time(mp_getfields_time_end)
2584
2585        mp_getfields_wtime_total = mp_getfields_wtime_total + &
2586             &(mp_getfields_wtime_end - mp_getfields_wtime_beg)
2587        mp_getfields_time_total = mp_getfields_time_total + &
2588             &(mp_getfields_time_end - mp_getfields_time_beg)
2589      end if
2590
2591    case ('partloop1')
2592      if (imode.eq.0) then
2593        call cpu_time(tm_nploop_beg)
2594      else
2595        call cpu_time(tm_nploop_end)
2596        tm_nploop_total = tm_nploop_total + (tm_nploop_end - tm_nploop_beg)
2597      end if
2598
2599    case ('conccalc')
2600      if (imode.eq.0) then
2601        call cpu_time(mp_conccalc_time_beg)
2602      else
2603        call cpu_time(mp_conccalc_time_end)
2604        mp_conccalc_time_total = mp_conccalc_time_total + mp_conccalc_time_end - &
2605             &mp_conccalc_time_beg
2606      end if
2607
2608    case ('rootonly')
2609      if (imode.eq.0) then
2610        call cpu_time(mp_root_time_beg)
2611        mp_root_wtime_beg = mpi_wtime()
2612      else
2613        call cpu_time(mp_root_time_end)
2614        mp_root_wtime_end = mpi_wtime()
2615        mp_root_time_total = mp_root_time_total + &
2616             &(mp_root_time_end - mp_root_time_beg)
2617        mp_root_wtime_total = mp_root_wtime_total + &
2618             &(mp_root_wtime_end - mp_root_wtime_beg)
2619      end if
2620
2621    case ('iotime')
2622      if (imode.eq.0) then
2623        mp_io_wtime_beg = mpi_wtime()
2624        call cpu_time(mp_io_time_beg)
2625      else
2626        mp_io_wtime_end = mpi_wtime()
2627        call cpu_time(mp_io_time_end)
2628
2629        mp_io_wtime_total = mp_io_wtime_total + (mp_io_wtime_end - &
2630             & mp_io_wtime_beg)
2631        mp_io_time_total = mp_io_time_total + (mp_io_time_end - &
2632             & mp_io_time_beg)
2633      end if
2634
2635    case ('verttransform')
2636      if (imode.eq.0) then
2637        mp_vt_wtime_beg = mpi_wtime()
2638        call cpu_time(mp_vt_time_beg)
2639      else
2640        mp_vt_wtime_end = mpi_wtime()
2641        call cpu_time(mp_vt_time_end)
2642
2643        mp_vt_wtime_total = mp_vt_wtime_total + (mp_vt_wtime_end - &
2644             & mp_vt_wtime_beg)
2645        mp_vt_time_total = mp_vt_time_total + (mp_vt_time_end - &
2646             & mp_vt_time_beg)
2647      end if
2648
2649    case ('readwind')
2650      if (imode.eq.0) then
2651        call cpu_time(mp_readwind_time_beg)
2652        mp_readwind_wtime_beg = mpi_wtime()
2653      else
2654        call cpu_time(mp_readwind_time_end)
2655        mp_readwind_wtime_end = mpi_wtime()
2656
2657        mp_readwind_time_total = mp_readwind_time_total + &
2658             &(mp_readwind_time_end - mp_readwind_time_beg)
2659        mp_readwind_wtime_total = mp_readwind_wtime_total + &
2660             &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
2661      end if
2662
2663    case ('commtime')
2664      if (imode.eq.0) then
2665        call cpu_time(mp_comm_time_beg)
2666        mp_comm_wtime_beg = mpi_wtime()
2667      else
2668        call cpu_time(mp_comm_time_end)
2669        mp_comm_wtime_end = mpi_wtime()
2670        mp_comm_time_total = mp_comm_time_total + &
2671             &(mp_comm_time_end - mp_comm_time_beg)
2672        mp_comm_wtime_total = mp_comm_wtime_total + &
2673             &(mp_comm_wtime_end - mp_comm_wtime_beg)
2674      end if
2675
2676    case ('flexpart')
2677      if (imode.eq.0) then
2678        mp_total_wtime_beg=mpi_wtime()
2679      else
2680        mp_total_wtime_end=mpi_wtime()
2681        mp_total_wtime_total = mp_total_wtime_end-mp_total_wtime_beg
2682      end if
2683
2684    case default
2685      write(*,*) 'mpi_mod::mpif_mtime> WARNING: unknown case identifier'
2686
2687    end select
2688
2689
2690  end subroutine mpif_mtime
2691
2692
2693  subroutine mpif_finalize
2694!***********************************************************************
2695! Finalize MPI
2696! Optionally print detailed time diagnostics
2697!***********************************************************************
2698    implicit none
2699
2700    integer :: ip !,j,r
2701
2702!***********************************************************************
2703
2704    IF (mp_measure_time) THEN
2705      do ip=0, mp_np-1
2706        call MPI_BARRIER(MPI_COMM_WORLD, mp_ierr)
2707
2708        if (ip == mp_pid) then
2709          write(*,FMT='(72("#"))')
2710          write(*,FMT='(A60,I3)') 'STATISTICS FOR MPI PROCESS:', mp_pid
2711          if (ip == 0) then
2712            write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR ROOT-ONLY (PID 0) &
2713                 &SECTIONS: ', mp_root_time_total
2714            write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR ROOT-ONLY (PID 0) &
2715                 &SECTIONS:', mp_root_wtime_total
2716          end if
2717          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FLEXPART: ' &
2718               &, mp_total_wtime_total
2719          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR MPI COMMUNICATIONS:',&
2720               & mp_comm_time_total
2721          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR MPI COMMUNICATIONS:',&
2722               & mp_comm_wtime_total
2723          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR MPI BARRIERS:',&
2724               & mp_barrier_time_total
2725          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR MPI BARRIERS:',&
2726               & mp_barrier_wtime_total
2727          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR TIMEMANAGER:',&
2728               & tm_tot_total
2729          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR TIMEMANAGER/NUMPART LOOP:',&
2730               & tm_nploop_total
2731          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR ADVANCE:',&
2732               & mp_advance_wtime_total
2733          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR GETFIELDS:',&
2734               & mp_getfields_wtime_total
2735          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR GETFIELDS:',&
2736               & mp_getfields_time_total
2737          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
2738               & mp_readwind_wtime_total
2739          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
2740               & mp_readwind_time_total
2741          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FILE IO:',&
2742               & mp_io_wtime_total
2743          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR FILE IO:',&
2744               & mp_io_time_total
2745          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR WETDEPO:',&
2746               & mp_wetdepo_wtime_total
2747          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR WETDEPO:',&
2748               & mp_wetdepo_time_total
2749          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR CONCCALC:',&
2750               & mp_conccalc_time_total
2751! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR VERTTRANSFORM:',&
2752!      & mp_vt_wtime_total
2753! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR VERTTRANSFORM:',&
2754!      & mp_vt_time_total
2755! NB: the 'flush' function is possibly a gfortran-specific extension,
2756! comment out if it gives problems
2757!          call flush()
2758        end if
2759      end do
2760    end if
2761
2762! This call to barrier is for correctly formatting output
2763    call MPI_BARRIER(MPI_COMM_WORLD, mp_ierr)
2764
2765    if (lroot.and.mp_measure_time) then
2766      write(*,FMT='(72("#"))')
2767      WRITE(*,*) "To turn off output of time measurements, set "
2768      WRITE(*,*) "    mp_measure_time=.false."
2769      WRITE(*,*) "in file mpi_mod.f90"
2770      write(*,FMT='(72("#"))')
2771    end if
2772
2773! j=mp_pid*nvar_async
2774! In the implementation with 3 fields, the processes may have posted
2775! MPI_Irecv requests that should be cancelled here
2776! if (.not.lmp_sync) then
2777!   r=mp_pid*nvar_async
2778!   do j=r,r+nvar_async-1
2779!     call MPI_Cancel(j,mp_ierr)
2780!     if (mp_ierr /= 0) write(*,*) '#### mpif_finalize::MPI_Cancel> ERROR ####'
2781!   end do
2782! end if
2783
2784    call MPI_FINALIZE(mp_ierr)
2785    if (mp_ierr /= 0) then
2786      write(*,*) '#### mpif_finalize::MPI_FINALIZE> MPI ERROR ', mp_ierr, ' ####'
2787      stop
2788    end if
2789
2790
2791  end subroutine mpif_finalize
2792
2793
2794  subroutine get_lun(my_lun)
2795!***********************************************************************
2796! get_lun:
2797!   Starting from 100, get next free logical unit number
2798!***********************************************************************
2799
2800    implicit none
2801
2802    integer, intent(inout) :: my_lun
2803    integer, save :: free_lun=100
2804    logical :: exists, iopen
2805
2806!***********************************************************************
2807
2808    loop1: do
2809      inquire(UNIT=free_lun, EXIST=exists, OPENED=iopen)
2810      if (exists .and. .not.iopen) exit loop1
2811      free_lun = free_lun+1
2812    end do loop1
2813    my_lun = free_lun
2814
2815  end subroutine get_lun
2816
2817
2818  subroutine write_data_dbg(array_in, array_name, tstep, ident)
2819!***********************************************************************
2820! Write one-dimensional arrays to file (for debugging purposes)
2821!***********************************************************************
2822    implicit none
2823
2824    real, intent(in), dimension(:) :: array_in
2825    integer, intent(in) :: tstep
2826    integer :: lios
2827    character(LEN=*), intent(in) :: ident, array_name
2828
2829    character(LEN=8) :: c_ts
2830    character(LEN=40) :: fn_1, fn_2
2831
2832!***********************************************************************
2833
2834    write(c_ts, FMT='(I8.8,BZ)') tstep
2835    fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)
2836    write(c_ts, FMT='(I2.2,BZ)') mp_np
2837    fn_2= trim(adjustl(array_name))//trim(adjustl(fn_1))//'-np'//trim(adjustl(c_ts))//'.dat'
2838
2839    call get_lun(dat_lun)
2840    open(UNIT=dat_lun, FILE=fn_2, IOSTAT=lios, ACTION='WRITE', &
2841         FORM='UNFORMATTED', STATUS='REPLACE')
2842    write(UNIT=dat_lun, IOSTAT=lios) array_in
2843    close(UNIT=dat_lun)
2844
2845  end subroutine write_data_dbg
2846
2847
2848  subroutine set_fields_synthetic()
2849!*******************************************************************************
2850! DESCRIPTION
2851!   Set all meteorological fields to synthetic (usually constant/homogeneous)
2852!   values.
2853!   Used for validation and error-checking
2854!
2855! NOTE
2856!   This version uses asynchronious communications.
2857!
2858! VARIABLES
2859!   
2860!
2861!
2862!*******************************************************************************
2863    use com_mod
2864
2865    implicit none
2866
2867    integer :: li=1, ui=2 ! wfmem indices (i.e, operate on entire field)
2868
2869    if (.not.lmp_sync) ui=3
2870
2871
2872! Variables transferred at initialization only
2873!*********************************************
2874!    readclouds=readclouds_
2875    oro(:,:)=0.0
2876    excessoro(:,:)=0.0
2877    lsm(:,:)=0.0
2878    xlanduse(:,:,:)=0.0
2879!    wftime
2880!    numbwf
2881!    nmixz
2882!    height
2883
2884! Time-varying fields:
2885    uu(:,:,:,li:ui) = 10.0
2886    vv(:,:,:,li:ui) = 0.0
2887    uupol(:,:,:,li:ui) = 10.0
2888    vvpol(:,:,:,li:ui)=0.0
2889    ww(:,:,:,li:ui)=0.
2890    tt(:,:,:,li:ui)=300.
2891    rho(:,:,:,li:ui)=1.3
2892    drhodz(:,:,:,li:ui)=.0
2893    tth(:,:,:,li:ui)=0.0
2894    qvh(:,:,:,li:ui)=1.0
2895    qv(:,:,:,li:ui)=1.0
2896
2897    pv(:,:,:,li:ui)=1.0
2898    clouds(:,:,:,li:ui)=0
2899
2900    clwc(:,:,:,li:ui)=0.0
2901    ciwc(:,:,:,li:ui)=0.0
2902
2903! 2D fields
2904
2905    cloudsh(:,:,li:ui)=0
2906    vdep(:,:,:,li:ui)=0.0
2907    ps(:,:,:,li:ui)=1.0e5
2908    sd(:,:,:,li:ui)=0.0
2909    tcc(:,:,:,li:ui)=0.0
2910    tt2(:,:,:,li:ui)=300.
2911    td2(:,:,:,li:ui)=250.
2912    lsprec(:,:,:,li:ui)=0.0
2913    convprec(:,:,:,li:ui)=0.0
2914    ustar(:,:,:,li:ui)=1.0
2915    wstar(:,:,:,li:ui)=1.0
2916    hmix(:,:,:,li:ui)=10000.
2917    tropopause(:,:,:,li:ui)=10000.
2918    oli(:,:,:,li:ui)=0.01
2919
2920  end subroutine set_fields_synthetic
2921
2922end module mpi_mod
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG