Context Navigation

mpi_mod.f90 @ aa8c34a

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

Last change on this file since aa8c34a was aa8c34a, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago
Minor edits
Property mode set to `100644`
File size: 106.9 KB

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

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

Download in other formats:

Original Format