source: flexpart.git/src/mpi_mod.f90 @ 4c64400

devrelease-10univie
Last change on this file since 4c64400 was 4c64400, checked in by Espen Sollum ATMOS <eso@…>, 3 years ago

Bugfix for double precision dry deposition calculation. Added (hardcoded) option to not use output kernel. Version/date string updated.

  • Property mode set to 100644
File size: 106.9 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*real(maxpart)/real(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,&
835           &itrasplit_tmp,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!   NB: take into account nested wind fields by using a separate
2343!   subroutine mpif_gf_request_nest (and separate request arrays for the
2344!   nested variables)
2345!
2346!*******************************************************************************
2347!    use com_mod,only: readclouds
2348
2349    implicit none
2350
2351
2352    integer :: n_req !,i
2353
2354!***********************************************************************
2355
2356    n_req=nvar_async*mp_np
2357
2358    if (mp_measure_time) call mpif_mtime('commtime',0)
2359
2360!    call MPI_Wait(rm1,MPI_STATUS_IGNORE,mp_ierr)
2361
2362! TODO: cancel recv request if readclouds=.false.
2363!    if (readclouds) then
2364    call MPI_Waitall(n_req,reqs,MPI_STATUSES_IGNORE,mp_ierr)
2365!    endif
2366! else
2367!   do i = 0, nvar_async*mp_np-1
2368!     if (mod(i,27).eq.0 .or. mod(i,28).eq.0) then
2369!       call MPI_Cancel(reqs(i),mp_ierr)
2370!       cycle
2371!     end if
2372!     call MPI_Wait(reqs(i),MPI_STATUS_IGNORE,mp_ierr)
2373!   end do
2374! end if
2375
2376    if (mp_ierr /= 0) goto 600
2377
2378    if (mp_measure_time) call mpif_mtime('commtime',1)
2379
2380    goto 601
2381
2382600 write(*,*) "#### mpi_mod::mpif_gf_request> MPI ERROR ", mp_ierr
2383    stop
2384
2385601 end subroutine mpif_gf_request
2386
2387
2388  subroutine mpif_tm_reduce_grid
2389!***********************************************************************
2390! Collect grid variable to PID 0, adding from all processes.
2391!
2392! NOTE
2393!   - Older MPI implementations may lack the MPI_IN_PLACE operation.
2394!     If so use 1) below and change unc_mod
2395!
2396!***********************************************************************
2397    use com_mod
2398    use unc_mod
2399    use par_mod
2400
2401    implicit none
2402
2403    integer :: grid_size2d,grid_size3d
2404    integer, parameter :: rcpt_size=maxreceptor*maxspec
2405
2406!**********************************************************************
2407    grid_size2d=numxgrid*numygrid*maxspec* &
2408         & maxpointspec_act*nclassunc*maxageclass
2409    grid_size3d=numxgrid*numygrid*numzgrid*maxspec* &
2410         & maxpointspec_act*nclassunc*maxageclass
2411
2412
2413! Time for MPI communications
2414    if (mp_measure_time) call mpif_mtime('commtime',0)
2415
2416! 1) Using a separate grid (gridunc0) for received values
2417! call MPI_Reduce(gridunc, gridunc0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2418!      & mp_comm_used, mp_ierr)
2419! if (mp_ierr /= 0) goto 600
2420
2421! 2) Using in-place reduction
2422    if (lroot) then
2423      call MPI_Reduce(MPI_IN_PLACE, gridunc, grid_size3d, mp_sp, MPI_SUM, id_root, &
2424           & mp_comm_used, mp_ierr)
2425      if (mp_ierr /= 0) goto 600
2426    else
2427      call MPI_Reduce(gridunc, 0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2428           & mp_comm_used, mp_ierr)
2429    end if
2430
2431    if ((WETDEP).and.(ldirect.gt.0)) then
2432      call MPI_Reduce(wetgridunc, wetgridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2433           & mp_comm_used, mp_ierr)
2434      if (mp_ierr /= 0) goto 600
2435    end if
2436
2437    if ((DRYDEP).and.(ldirect.gt.0)) then
2438      call MPI_Reduce(drygridunc, drygridunc0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2439           & mp_comm_used, mp_ierr)
2440      if (mp_ierr /= 0) goto 600
2441    end if
2442
2443! Receptor concentrations   
2444    if (lroot) then
2445      call MPI_Reduce(MPI_IN_PLACE,creceptor,rcpt_size,mp_sp,MPI_SUM,id_root, &
2446           & mp_comm_used,mp_ierr)
2447      if (mp_ierr /= 0) goto 600
2448    else
2449      call MPI_Reduce(creceptor,0,rcpt_size,mp_sp,MPI_SUM,id_root, &
2450           & mp_comm_used,mp_ierr)
2451    end if
2452
2453    if (mp_measure_time) call mpif_mtime('commtime',1)
2454
2455    goto 601
2456
2457600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
2458    stop
2459
2460601 end subroutine mpif_tm_reduce_grid
2461
2462
2463  subroutine mpif_tm_reduce_grid_nest
2464!***********************************************************************
2465! Collect nested grids to PID 0, adding from all processes.
2466!
2467! NOTE
2468!   - Compiling with 'fcheck=all' (gfortran) will cause a run-time error
2469!     as wetgriduncn0 is only allocated for root process. Possibly valid
2470!     MPI code but invalid gfortran.
2471!
2472!***********************************************************************
2473    use com_mod
2474    use unc_mod
2475    use par_mod
2476
2477    implicit none
2478
2479    integer :: grid_size2d,grid_size3d
2480
2481!**********************************************************************
2482
2483    grid_size3d=numxgridn*numygridn*numzgrid*maxspec* &
2484         & maxpointspec_act*nclassunc*maxageclass
2485    grid_size2d=numxgridn*numygridn*maxspec* &
2486         & maxpointspec_act*nclassunc*maxageclass
2487
2488
2489! Time for MPI communications
2490    if (mp_measure_time) call mpif_mtime('commtime',0)
2491
2492! Using a separate grid (gridunc0) for received values, for debugging
2493! call MPI_Reduce(griduncn, griduncn0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2494!      & mp_comm_used, mp_ierr)
2495! if (mp_ierr /= 0) goto 600
2496
2497! Using in-place reduction
2498    if (lroot) then
2499      call MPI_Reduce(MPI_IN_PLACE, griduncn, grid_size3d, mp_sp, MPI_SUM, id_root, &
2500           & mp_comm_used, mp_ierr)
2501      if (mp_ierr /= 0) goto 600
2502    else
2503      call MPI_Reduce(griduncn, 0, grid_size3d, mp_sp, MPI_SUM, id_root, &
2504           & mp_comm_used, mp_ierr)
2505    end if
2506
2507    if ((WETDEP).and.(ldirect.gt.0)) then
2508      call MPI_Reduce(wetgriduncn, wetgriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2509           & mp_comm_used, mp_ierr)
2510      if (mp_ierr /= 0) goto 600
2511    end if
2512
2513    if ((DRYDEP).and.(ldirect.gt.0)) then
2514      call MPI_Reduce(drygriduncn, drygriduncn0, grid_size2d, mp_cp, MPI_SUM, id_root, &
2515           & mp_comm_used, mp_ierr)
2516      if (mp_ierr /= 0) goto 600
2517    end if
2518
2519    if (mp_measure_time) call mpif_mtime('commtime',1)
2520
2521    goto 601
2522
2523600 write(*,*) "mpi_mod> mp_ierr \= 0", mp_ierr
2524    stop
2525
2526601 end subroutine mpif_tm_reduce_grid_nest
2527
2528
2529  subroutine mpif_mtime(ident,imode)
2530!***********************************************************************
2531! Measure CPU/Wall time in various parts of the code
2532!
2533! VARIABLES
2534!   ident        character, identifies routine to measure
2535!   imode        integer, 0:start clock(s)  1: stop clock(s)
2536!
2537!***********************************************************************
2538    implicit none
2539
2540    character(LEN=*), intent(in) :: ident
2541    integer, intent(in) :: imode
2542
2543!***********************************************************************
2544
2545    select case(ident)
2546
2547    case ('timemanager')
2548      if (imode.eq.0) then
2549        call cpu_time(tm_tot_beg)
2550      else
2551        call cpu_time(tm_tot_end)
2552        tm_tot_total = tm_tot_total + (tm_tot_end-tm_tot_beg)
2553      end if
2554
2555    case ('wetdepo')
2556      if (imode.eq.0) then
2557        mp_wetdepo_wtime_beg = mpi_wtime()
2558        call cpu_time(mp_wetdepo_time_beg)
2559      else
2560        mp_wetdepo_wtime_end = mpi_wtime()
2561        call cpu_time(mp_wetdepo_time_end)
2562
2563        mp_wetdepo_wtime_total = mp_wetdepo_wtime_total + &
2564             &(mp_wetdepo_wtime_end - mp_wetdepo_wtime_beg)
2565        mp_wetdepo_time_total = mp_wetdepo_time_total + &
2566             &(mp_wetdepo_time_end - mp_wetdepo_time_beg)
2567      end if
2568
2569    case ('advance')
2570      if (imode.eq.0) then
2571        mp_advance_wtime_beg = mpi_wtime()
2572      else
2573        mp_advance_wtime_end = mpi_wtime()
2574
2575        mp_advance_wtime_total = mp_advance_wtime_total + &
2576             &(mp_advance_wtime_end - mp_advance_wtime_beg)
2577      end if
2578
2579    case ('getfields')
2580      if (imode.eq.0) then
2581        mp_getfields_wtime_beg = mpi_wtime()
2582        call cpu_time(mp_getfields_time_beg)
2583      else
2584        mp_getfields_wtime_end = mpi_wtime()
2585        call cpu_time(mp_getfields_time_end)
2586
2587        mp_getfields_wtime_total = mp_getfields_wtime_total + &
2588             &(mp_getfields_wtime_end - mp_getfields_wtime_beg)
2589        mp_getfields_time_total = mp_getfields_time_total + &
2590             &(mp_getfields_time_end - mp_getfields_time_beg)
2591      end if
2592
2593    case ('partloop1')
2594      if (imode.eq.0) then
2595        call cpu_time(tm_nploop_beg)
2596      else
2597        call cpu_time(tm_nploop_end)
2598        tm_nploop_total = tm_nploop_total + (tm_nploop_end - tm_nploop_beg)
2599      end if
2600
2601    case ('conccalc')
2602      if (imode.eq.0) then
2603        call cpu_time(mp_conccalc_time_beg)
2604      else
2605        call cpu_time(mp_conccalc_time_end)
2606        mp_conccalc_time_total = mp_conccalc_time_total + mp_conccalc_time_end - &
2607             &mp_conccalc_time_beg
2608      end if
2609
2610    case ('rootonly')
2611      if (imode.eq.0) then
2612        call cpu_time(mp_root_time_beg)
2613        mp_root_wtime_beg = mpi_wtime()
2614      else
2615        call cpu_time(mp_root_time_end)
2616        mp_root_wtime_end = mpi_wtime()
2617        mp_root_time_total = mp_root_time_total + &
2618             &(mp_root_time_end - mp_root_time_beg)
2619        mp_root_wtime_total = mp_root_wtime_total + &
2620             &(mp_root_wtime_end - mp_root_wtime_beg)
2621      end if
2622
2623    case ('iotime')
2624      if (imode.eq.0) then
2625        mp_io_wtime_beg = mpi_wtime()
2626        call cpu_time(mp_io_time_beg)
2627      else
2628        mp_io_wtime_end = mpi_wtime()
2629        call cpu_time(mp_io_time_end)
2630
2631        mp_io_wtime_total = mp_io_wtime_total + (mp_io_wtime_end - &
2632             & mp_io_wtime_beg)
2633        mp_io_time_total = mp_io_time_total + (mp_io_time_end - &
2634             & mp_io_time_beg)
2635      end if
2636
2637    case ('verttransform')
2638      if (imode.eq.0) then
2639        mp_vt_wtime_beg = mpi_wtime()
2640        call cpu_time(mp_vt_time_beg)
2641      else
2642        mp_vt_wtime_end = mpi_wtime()
2643        call cpu_time(mp_vt_time_end)
2644
2645        mp_vt_wtime_total = mp_vt_wtime_total + (mp_vt_wtime_end - &
2646             & mp_vt_wtime_beg)
2647        mp_vt_time_total = mp_vt_time_total + (mp_vt_time_end - &
2648             & mp_vt_time_beg)
2649      end if
2650
2651    case ('readwind')
2652      if (imode.eq.0) then
2653        call cpu_time(mp_readwind_time_beg)
2654        mp_readwind_wtime_beg = mpi_wtime()
2655      else
2656        call cpu_time(mp_readwind_time_end)
2657        mp_readwind_wtime_end = mpi_wtime()
2658
2659        mp_readwind_time_total = mp_readwind_time_total + &
2660             &(mp_readwind_time_end - mp_readwind_time_beg)
2661        mp_readwind_wtime_total = mp_readwind_wtime_total + &
2662             &(mp_readwind_wtime_end - mp_readwind_wtime_beg)
2663      end if
2664
2665    case ('commtime')
2666      if (imode.eq.0) then
2667        call cpu_time(mp_comm_time_beg)
2668        mp_comm_wtime_beg = mpi_wtime()
2669      else
2670        call cpu_time(mp_comm_time_end)
2671        mp_comm_wtime_end = mpi_wtime()
2672        mp_comm_time_total = mp_comm_time_total + &
2673             &(mp_comm_time_end - mp_comm_time_beg)
2674        mp_comm_wtime_total = mp_comm_wtime_total + &
2675             &(mp_comm_wtime_end - mp_comm_wtime_beg)
2676      end if
2677
2678    case ('flexpart')
2679      if (imode.eq.0) then
2680        mp_total_wtime_beg=mpi_wtime()
2681      else
2682        mp_total_wtime_end=mpi_wtime()
2683        mp_total_wtime_total = mp_total_wtime_end-mp_total_wtime_beg
2684      end if
2685
2686    case default
2687      write(*,*) 'mpi_mod::mpif_mtime> WARNING: unknown case identifier'
2688
2689    end select
2690
2691
2692  end subroutine mpif_mtime
2693
2694
2695  subroutine mpif_finalize
2696!***********************************************************************
2697! Finalize MPI
2698! Optionally print detailed time diagnostics
2699!***********************************************************************
2700    implicit none
2701
2702    integer :: ip !,j,r
2703
2704!***********************************************************************
2705
2706    IF (mp_measure_time) THEN
2707      do ip=0, mp_np-1
2708        call MPI_BARRIER(MPI_COMM_WORLD, mp_ierr)
2709
2710        if (ip == mp_pid) then
2711          write(*,FMT='(72("#"))')
2712          write(*,FMT='(A60,I3)') 'STATISTICS FOR MPI PROCESS:', mp_pid
2713          if (ip == 0) then
2714            write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR ROOT-ONLY (PID 0) &
2715                 &SECTIONS: ', mp_root_time_total
2716            write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR ROOT-ONLY (PID 0) &
2717                 &SECTIONS:', mp_root_wtime_total
2718          end if
2719          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FLEXPART: ' &
2720               &, mp_total_wtime_total
2721          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR MPI COMMUNICATIONS:',&
2722               & mp_comm_time_total
2723          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR MPI COMMUNICATIONS:',&
2724               & mp_comm_wtime_total
2725          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR MPI BARRIERS:',&
2726               & mp_barrier_time_total
2727          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR MPI BARRIERS:',&
2728               & mp_barrier_wtime_total
2729          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR TIMEMANAGER:',&
2730               & tm_tot_total
2731          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR TIMEMANAGER/NUMPART LOOP:',&
2732               & tm_nploop_total
2733          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR ADVANCE:',&
2734               & mp_advance_wtime_total
2735          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR GETFIELDS:',&
2736               & mp_getfields_wtime_total
2737          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR GETFIELDS:',&
2738               & mp_getfields_time_total
2739          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR READWIND:',&
2740               & mp_readwind_wtime_total
2741          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR READWIND:',&
2742               & mp_readwind_time_total
2743          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR FILE IO:',&
2744               & mp_io_wtime_total
2745          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR FILE IO:',&
2746               & mp_io_time_total
2747          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR WETDEPO:',&
2748               & mp_wetdepo_wtime_total
2749          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR WETDEPO:',&
2750               & mp_wetdepo_time_total
2751          write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR CONCCALC:',&
2752               & mp_conccalc_time_total
2753! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL WALL TIME FOR VERTTRANSFORM:',&
2754!      & mp_vt_wtime_total
2755! write(*,FMT='(A60,TR1,F9.2)') 'TOTAL CPU TIME FOR VERTTRANSFORM:',&
2756!      & mp_vt_time_total
2757! NB: the 'flush' function is possibly a gfortran-specific extension,
2758! comment out if it gives problems
2759!          call flush()
2760        end if
2761      end do
2762    end if
2763
2764! This call to barrier is for correctly formatting output
2765    call MPI_BARRIER(MPI_COMM_WORLD, mp_ierr)
2766
2767    if (lroot.and.mp_measure_time) then
2768      write(*,FMT='(72("#"))')
2769      WRITE(*,*) "To turn off output of time measurements, set "
2770      WRITE(*,*) "    mp_measure_time=.false."
2771      WRITE(*,*) "in file mpi_mod.f90"
2772      write(*,FMT='(72("#"))')
2773    end if
2774
2775! j=mp_pid*nvar_async
2776! In the implementation with 3 fields, the processes may have posted
2777! MPI_Irecv requests that should be cancelled here
2778! if (.not.lmp_sync) then
2779!   r=mp_pid*nvar_async
2780!   do j=r,r+nvar_async-1
2781!     call MPI_Cancel(j,mp_ierr)
2782!     if (mp_ierr /= 0) write(*,*) '#### mpif_finalize::MPI_Cancel> ERROR ####'
2783!   end do
2784! end if
2785
2786    call MPI_FINALIZE(mp_ierr)
2787    if (mp_ierr /= 0) then
2788      write(*,*) '#### mpif_finalize::MPI_FINALIZE> MPI ERROR ', mp_ierr, ' ####'
2789      stop
2790    end if
2791
2792
2793  end subroutine mpif_finalize
2794
2795
2796  subroutine get_lun(my_lun)
2797!***********************************************************************
2798! get_lun:
2799!   Starting from 100, get next free logical unit number
2800!***********************************************************************
2801
2802    implicit none
2803
2804    integer, intent(inout) :: my_lun
2805    integer, save :: free_lun=100
2806    logical :: exists, iopen
2807
2808!***********************************************************************
2809
2810    loop1: do
2811      inquire(UNIT=free_lun, EXIST=exists, OPENED=iopen)
2812      if (exists .and. .not.iopen) exit loop1
2813      free_lun = free_lun+1
2814    end do loop1
2815    my_lun = free_lun
2816
2817  end subroutine get_lun
2818
2819
2820  subroutine write_data_dbg(array_in, array_name, tstep, ident)
2821!***********************************************************************
2822! Write one-dimensional arrays to file (for debugging purposes)
2823!***********************************************************************
2824    implicit none
2825
2826    real, intent(in), dimension(:) :: array_in
2827    integer, intent(in) :: tstep
2828    integer :: lios
2829    character(LEN=*), intent(in) :: ident, array_name
2830
2831    character(LEN=8) :: c_ts
2832    character(LEN=40) :: fn_1, fn_2
2833
2834!***********************************************************************
2835
2836    write(c_ts, FMT='(I8.8,BZ)') tstep
2837    fn_1='-'//trim(adjustl(c_ts))//'-'//trim(ident)
2838    write(c_ts, FMT='(I2.2,BZ)') mp_np
2839    fn_2= trim(adjustl(array_name))//trim(adjustl(fn_1))//'-np'//trim(adjustl(c_ts))//'.dat'
2840
2841    call get_lun(dat_lun)
2842    open(UNIT=dat_lun, FILE=fn_2, IOSTAT=lios, ACTION='WRITE', &
2843         FORM='UNFORMATTED', STATUS='REPLACE')
2844    write(UNIT=dat_lun, IOSTAT=lios) array_in
2845    close(UNIT=dat_lun)
2846
2847  end subroutine write_data_dbg
2848
2849
2850  subroutine set_fields_synthetic()
2851!*******************************************************************************
2852! DESCRIPTION
2853!   Set all meteorological fields to synthetic (usually constant/homogeneous)
2854!   values.
2855!   Used for validation and error-checking
2856!
2857! NOTE
2858!   This version uses asynchronious communications.
2859!
2860! VARIABLES
2861!   
2862!
2863!
2864!*******************************************************************************
2865    use com_mod
2866
2867    implicit none
2868
2869    integer :: li=1, ui=2 ! wfmem indices (i.e, operate on entire field)
2870
2871    if (.not.lmp_sync) ui=3
2872
2873
2874! Variables transferred at initialization only
2875!*********************************************
2876!    readclouds=readclouds_
2877    oro(:,:)=0.0
2878    excessoro(:,:)=0.0
2879    lsm(:,:)=0.0
2880    xlanduse(:,:,:)=0.0
2881!    wftime
2882!    numbwf
2883!    nmixz
2884!    height
2885
2886! Time-varying fields:
2887    uu(:,:,:,li:ui) = 10.0
2888    vv(:,:,:,li:ui) = 0.0
2889    uupol(:,:,:,li:ui) = 10.0
2890    vvpol(:,:,:,li:ui)=0.0
2891    ww(:,:,:,li:ui)=0.
2892    tt(:,:,:,li:ui)=300.
2893    rho(:,:,:,li:ui)=1.3
2894    drhodz(:,:,:,li:ui)=.0
2895    tth(:,:,:,li:ui)=0.0
2896    qvh(:,:,:,li:ui)=1.0
2897    qv(:,:,:,li:ui)=1.0
2898
2899    pv(:,:,:,li:ui)=1.0
2900    clouds(:,:,:,li:ui)=0
2901
2902    clwc(:,:,:,li:ui)=0.0
2903    ciwc(:,:,:,li:ui)=0.0
2904
2905! 2D fields
2906
2907    cloudsh(:,:,li:ui)=0
2908    vdep(:,:,:,li:ui)=0.0
2909    ps(:,:,:,li:ui)=1.0e5
2910    sd(:,:,:,li:ui)=0.0
2911    tcc(:,:,:,li:ui)=0.0
2912    tt2(:,:,:,li:ui)=300.
2913    td2(:,:,:,li:ui)=250.
2914    lsprec(:,:,:,li:ui)=0.0
2915    convprec(:,:,:,li:ui)=0.0
2916    ustar(:,:,:,li:ui)=1.0
2917    wstar(:,:,:,li:ui)=1.0
2918    hmix(:,:,:,li:ui)=10000.
2919    tropopause(:,:,:,li:ui)=10000.
2920    oli(:,:,:,li:ui)=0.01
2921
2922  end subroutine set_fields_synthetic
2923
2924end module mpi_mod
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG