Context Navigation

mpi_mod.f90 @ 861805a

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format