Context Navigation

mpi_mod.f90 @ 4c64400

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

Last change on this file since 4c64400 was 4c64400, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago
Bugfix for double precision dry deposition calculation. Added (hardcoded) option to not use output kernel. Version/date string updated.
Property mode set to `100644`
File size: 106.9 KB

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

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

Download in other formats:

Original Format