Home

Context Navigation

source: palm/trunk/SOURCE/data_output_2d.f90 @ 1429

Last change on this file since 1429 was 1360, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 72.1 KB

Line
1	SUBROUTINE data_output_2d( mode, av )
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: data_output_2d.f90 1360 2014-04-11 17:20:32Z knoop $
27	!
28	! 1359 2014-04-11 17:15:14Z hoffmann
29	! New particle structure integrated.
30	!
31	! 1353 2014-04-08 15:21:23Z heinze
32	! REAL constants provided with KIND-attribute
33	!
34	! 1327 2014-03-21 11:00:16Z raasch
35	! parts concerning iso2d output removed,
36	! -netcdf output queries
37	!
38	! 1320 2014-03-20 08:40:49Z raasch
39	! ONLY-attribute added to USE-statements,
40	! kind-parameters added to all INTEGER and REAL declaration statements,
41	! kinds are defined in new module kinds,
42	! revision history before 2012 removed,
43	! comment fields (!:) to be used for variable explanations added to
44	! all variable declaration statements
45	!
46	! 1318 2014-03-17 13:35:16Z raasch
47	! barrier argument removed from cpu_log.
48	! module interfaces removed
49	!
50	! 1311 2014-03-14 12:13:39Z heinze
51	! bugfix: close #if defined( __netcdf )
52	!
53	! 1308 2014-03-13 14:58:42Z fricke
54	! +local_2d_sections, local_2d_sections_l, ns
55	! Check, if the limit of the time dimension is exceeded for parallel output
56	! To increase the performance for parallel output, the following is done:
57	! - Update of time axis is only done by PE0
58	! - Cross sections are first stored on a local array and are written
59	! collectively to the output file by all PEs.
60	!
61	! 1115 2013-03-26 18:16:16Z hoffmann
62	! ql is calculated by calc_liquid_water_content
63	!
64	! 1076 2012-12-05 08:30:18Z hoffmann
65	! Bugfix in output of ql
66	!
67	! 1065 2012-11-22 17:42:36Z hoffmann
68	! Bugfix: Output of cross sections of ql
69	!
70	! 1053 2012-11-13 17:11:03Z hoffmann
71	! +qr, nr, qc and cross sections
72	!
73	! 1036 2012-10-22 13:43:42Z raasch
74	! code put under GPL (PALM 3.9)
75	!
76	! 1031 2012-10-19 14:35:30Z raasch
77	! netCDF4 without parallel file support implemented
78	!
79	! 1007 2012-09-19 14:30:36Z franke
80	! Bugfix: missing calculation of ql_vp added
81	!
82	! 978 2012-08-09 08:28:32Z fricke
83	! +z0h
84	!
85	! Revision 1.1 1997/08/11 06:24:09 raasch
86	! Initial revision
87	!
88	!
89	! Description:
90	! ------------
91	! Data output of horizontal cross-sections in netCDF format or binary format
92	! compatible to old graphic software iso2d.
93	! Attention: The position of the sectional planes is still not always computed
94	! --------- correctly. (zu is used always)!
95	!------------------------------------------------------------------------------!
96
97	USE arrays_3d, &
98	ONLY: dzw, e, nr, p, pt, q, qc, ql, ql_c, ql_v, ql_vp, qr, qsws, &
99	rho, sa, shf, tend, ts, u, us, v, vpt, w, z0, z0h, zu, zw
100
101	USE averaging
102
103	USE cloud_parameters, &
104	ONLY: hyrho, l_d_cp, precipitation_amount, precipitation_rate, prr, &
105	pt_d_t
106
107	USE control_parameters, &
108	ONLY: cloud_physics, data_output_2d_on_each_pe, data_output_xy, &
109	data_output_xz, data_output_yz, do2d, &
110	do2d_xy_last_time, do2d_xy_n, do2d_xy_time_count, &
111	do2d_xz_last_time, do2d_xz_n, do2d_xz_time_count, &
112	do2d_yz_last_time, do2d_yz_n, do2d_yz_time_count, &
113	ibc_uv_b, icloud_scheme, io_blocks, io_group, &
114	message_string, netcdf_data_format, &
115	ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz, psolver, section, &
116	simulated_time, simulated_time_chr, time_since_reference_point
117
118	USE cpulog, &
119	ONLY: cpu_log, log_point
120
121	USE grid_variables, &
122	ONLY: dx, dy
123
124	USE indices, &
125	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
126	nz, nzb, nzt
127
128	USE kinds
129
130	USE netcdf_control
131
132	USE particle_attributes, &
133	ONLY: grid_particles, number_of_particles, particle_advection_start, &
134	particles, prt_count
135
136	USE pegrid
137
138	IMPLICIT NONE
139
140	CHARACTER (LEN=2) :: do2d_mode !:
141	CHARACTER (LEN=2) :: mode !:
142	CHARACTER (LEN=4) :: grid !:
143	CHARACTER (LEN=25) :: section_chr !:
144	CHARACTER (LEN=50) :: rtext !:
145
146	INTEGER(iwp) :: av !:
147	INTEGER(iwp) :: ngp !:
148	INTEGER(iwp) :: file_id !:
149	INTEGER(iwp) :: i !:
150	INTEGER(iwp) :: if !:
151	INTEGER(iwp) :: is !:
152	INTEGER(iwp) :: iis !:
153	INTEGER(iwp) :: j !:
154	INTEGER(iwp) :: k !:
155	INTEGER(iwp) :: l !:
156	INTEGER(iwp) :: layer_xy !:
157	INTEGER(iwp) :: n !:
158	INTEGER(iwp) :: ns !:
159	INTEGER(iwp) :: psi !:
160	INTEGER(iwp) :: s !:
161	INTEGER(iwp) :: sender !:
162	INTEGER(iwp) :: ind(4) !:
163
164	LOGICAL :: found !:
165	LOGICAL :: resorted !:
166	LOGICAL :: two_d !:
167
168	REAL(wp) :: mean_r !:
169	REAL(wp) :: s_r2 !:
170	REAL(wp) :: s_r3 !:
171
172	REAL(wp), DIMENSION(:), ALLOCATABLE :: level_z !:
173	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: local_2d !:
174	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: local_2d_l !:
175	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !:
176	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_2d_sections !:
177	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_2d_sections_l !:
178
179	#if defined( __parallel )
180	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: total_2d !:
181	#endif
182	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !:
183
184	NAMELIST /LOCAL/ rtext
185
186	!
187	!-- Immediate return, if no output is requested (no respective sections
188	!-- found in parameter data_output)
189	IF ( mode == 'xy' .AND. .NOT. data_output_xy(av) ) RETURN
190	IF ( mode == 'xz' .AND. .NOT. data_output_xz(av) ) RETURN
191	IF ( mode == 'yz' .AND. .NOT. data_output_yz(av) ) RETURN
192	!
193	!-- For parallel netcdf output the time axis must be limited. Return, if this
194	!-- limit is exceeded. This could be the case, if the simulated time exceeds
195	!-- the given end time by the length of the given output interval.
196	IF ( netcdf_data_format > 4 ) THEN
197	IF ( mode == 'xy' .AND. do2d_xy_time_count(av) + 1 > &
198	ntdim_2d_xy(av) ) THEN
199	WRITE ( message_string, * ) 'Output of xy cross-sections is not ', &
200	'given at t=', simulated_time, '&because the', &
201	' maximum number of output time levels is exceeded.'
202	CALL message( 'data_output_2d', 'PA0384', 0, 1, 0, 6, 0 )
203	RETURN
204	ENDIF
205	IF ( mode == 'xz' .AND. do2d_xz_time_count(av) + 1 > &
206	ntdim_2d_xz(av) ) THEN
207	WRITE ( message_string, * ) 'Output of xz cross-sections is not ', &
208	'given at t=', simulated_time, '&because the', &
209	' maximum number of output time levels is exceeded.'
210	CALL message( 'data_output_2d', 'PA0385', 0, 1, 0, 6, 0 )
211	RETURN
212	ENDIF
213	IF ( mode == 'yz' .AND. do2d_yz_time_count(av) + 1 > &
214	ntdim_2d_yz(av) ) THEN
215	WRITE ( message_string, * ) 'Output of yz cross-sections is not ', &
216	'given at t=', simulated_time, '&because the', &
217	' maximum number of output time levels is exceeded.'
218	CALL message( 'data_output_2d', 'PA0386', 0, 1, 0, 6, 0 )
219	RETURN
220	ENDIF
221	ENDIF
222
223	CALL cpu_log (log_point(3),'data_output_2d','start')
224
225	two_d = .FALSE. ! local variable to distinguish between output of pure 2D
226	! arrays and cross-sections of 3D arrays.
227
228	!
229	!-- Depending on the orientation of the cross-section, the respective output
230	!-- files have to be opened.
231	SELECT CASE ( mode )
232
233	CASE ( 'xy' )
234	s = 1
235	ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
236
237	IF ( netcdf_data_format > 4 ) THEN
238	ns = 1
239	DO WHILE ( section(ns,s) /= -9999 .AND. ns <= 100 )
240	ns = ns + 1
241	ENDDO
242	ns = ns - 1
243	ALLOCATE( local_2d_sections(nxlg:nxrg,nysg:nyng,1:ns) )
244	local_2d_sections = 0.0_wp
245	ENDIF
246
247	!
248	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
249	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
250	CALL check_open( 101+av*10 )
251	ENDIF
252
253	IF ( data_output_2d_on_each_pe ) THEN
254	CALL check_open( 21 )
255	ELSE
256	IF ( myid == 0 ) THEN
257	#if defined( __parallel )
258	ALLOCATE( total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp) )
259	#endif
260	ENDIF
261	ENDIF
262
263	CASE ( 'xz' )
264	s = 2
265	ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
266
267	IF ( netcdf_data_format > 4 ) THEN
268	ns = 1
269	DO WHILE ( section(ns,s) /= -9999 .AND. ns <= 100 )
270	ns = ns + 1
271	ENDDO
272	ns = ns - 1
273	ALLOCATE( local_2d_sections(nxlg:nxrg,1:ns,nzb:nzt+1) )
274	ALLOCATE( local_2d_sections_l(nxlg:nxrg,1:ns,nzb:nzt+1) )
275	local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
276	ENDIF
277
278	!
279	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
280	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
281	CALL check_open( 102+av*10 )
282	ENDIF
283
284	IF ( data_output_2d_on_each_pe ) THEN
285	CALL check_open( 22 )
286	ELSE
287	IF ( myid == 0 ) THEN
288	#if defined( __parallel )
289	ALLOCATE( total_2d(-nbgp:nx+nbgp,nzb:nzt+1) )
290	#endif
291	ENDIF
292	ENDIF
293
294	CASE ( 'yz' )
295	s = 3
296	ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
297
298	IF ( netcdf_data_format > 4 ) THEN
299	ns = 1
300	DO WHILE ( section(ns,s) /= -9999 .AND. ns <= 100 )
301	ns = ns + 1
302	ENDDO
303	ns = ns - 1
304	ALLOCATE( local_2d_sections(1:ns,nysg:nyng,nzb:nzt+1) )
305	ALLOCATE( local_2d_sections_l(1:ns,nysg:nyng,nzb:nzt+1) )
306	local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
307	ENDIF
308
309	!
310	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
311	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
312	CALL check_open( 103+av*10 )
313	ENDIF
314
315	IF ( data_output_2d_on_each_pe ) THEN
316	CALL check_open( 23 )
317	ELSE
318	IF ( myid == 0 ) THEN
319	#if defined( __parallel )
320	ALLOCATE( total_2d(-nbgp:ny+nbgp,nzb:nzt+1) )
321	#endif
322	ENDIF
323	ENDIF
324
325	CASE DEFAULT
326	message_string = 'unknown cross-section: ' // TRIM( mode )
327	CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
328
329	END SELECT
330
331	!
332	!-- Allocate a temporary array for resorting (kji -> ijk).
333	ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nzt+1) )
334
335	!
336	!-- Loop of all variables to be written.
337	!-- Output dimensions chosen
338	if = 1
339	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
340	do2d_mode = do2d(av,if)(l-1:l)
341
342	DO WHILE ( do2d(av,if)(1:1) /= ' ' )
343
344	IF ( do2d_mode == mode ) THEN
345	!
346	!-- Store the array chosen on the temporary array.
347	resorted = .FALSE.
348	SELECT CASE ( TRIM( do2d(av,if) ) )
349
350	CASE ( 'e_xy', 'e_xz', 'e_yz' )
351	IF ( av == 0 ) THEN
352	to_be_resorted => e
353	ELSE
354	to_be_resorted => e_av
355	ENDIF
356	IF ( mode == 'xy' ) level_z = zu
357
358	CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
359	IF ( av == 0 ) THEN
360	to_be_resorted => pt
361	ELSE
362	to_be_resorted => lpt_av
363	ENDIF
364	IF ( mode == 'xy' ) level_z = zu
365
366	CASE ( 'lwp*_xy' ) ! 2d-array
367	IF ( av == 0 ) THEN
368	DO i = nxlg, nxrg
369	DO j = nysg, nyng
370	local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
371	dzw(1:nzt+1) )
372	ENDDO
373	ENDDO
374	ELSE
375	DO i = nxlg, nxrg
376	DO j = nysg, nyng
377	local_pf(i,j,nzb+1) = lwp_av(j,i)
378	ENDDO
379	ENDDO
380	ENDIF
381	resorted = .TRUE.
382	two_d = .TRUE.
383	level_z(nzb+1) = zu(nzb+1)
384
385	CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
386	IF ( av == 0 ) THEN
387	to_be_resorted => nr
388	ELSE
389	to_be_resorted => nr_av
390	ENDIF
391	IF ( mode == 'xy' ) level_z = zu
392
393	CASE ( 'p_xy', 'p_xz', 'p_yz' )
394	IF ( av == 0 ) THEN
395	IF ( psolver /= 'sor' ) CALL exchange_horiz( p, nbgp )
396	to_be_resorted => p
397	ELSE
398	IF ( psolver /= 'sor' ) CALL exchange_horiz( p_av, nbgp )
399	to_be_resorted => p_av
400	ENDIF
401	IF ( mode == 'xy' ) level_z = zu
402
403	CASE ( 'pc_xy', 'pc_xz', 'pc_yz' ) ! particle concentration
404	IF ( av == 0 ) THEN
405	IF ( simulated_time >= particle_advection_start ) THEN
406	tend = prt_count
407	CALL exchange_horiz( tend, nbgp )
408	ELSE
409	tend = 0.0_wp
410	ENDIF
411	DO i = nxlg, nxrg
412	DO j = nysg, nyng
413	DO k = nzb, nzt+1
414	local_pf(i,j,k) = tend(k,j,i)
415	ENDDO
416	ENDDO
417	ENDDO
418	resorted = .TRUE.
419	ELSE
420	CALL exchange_horiz( pc_av, nbgp )
421	to_be_resorted => pc_av
422	ENDIF
423
424	CASE ( 'pr_xy', 'pr_xz', 'pr_yz' ) ! mean particle radius (effective radius)
425	IF ( av == 0 ) THEN
426	IF ( simulated_time >= particle_advection_start ) THEN
427	DO i = nxl, nxr
428	DO j = nys, nyn
429	DO k = nzb, nzt+1
430	number_of_particles = prt_count(k,j,i)
431	IF (number_of_particles <= 0) CYCLE
432	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
433	s_r2 = 0.0_wp
434	s_r3 = 0.0_wp
435	DO n = 1, number_of_particles
436	IF ( particles(n)%particle_mask ) THEN
437	s_r2 = s_r2 + particles(n)%radius*2 &
438	particles(n)%weight_factor
439	s_r3 = s_r3 + particles(n)%radius*3 &
440	particles(n)%weight_factor
441	ENDIF
442	ENDDO
443	IF ( s_r2 > 0.0_wp ) THEN
444	mean_r = s_r3 / s_r2
445	ELSE
446	mean_r = 0.0_wp
447	ENDIF
448	tend(k,j,i) = mean_r
449	ENDDO
450	ENDDO
451	ENDDO
452	CALL exchange_horiz( tend, nbgp )
453	ELSE
454	tend = 0.0_wp
455	ENDIF
456	DO i = nxlg, nxrg
457	DO j = nysg, nyng
458	DO k = nzb, nzt+1
459	local_pf(i,j,k) = tend(k,j,i)
460	ENDDO
461	ENDDO
462	ENDDO
463	resorted = .TRUE.
464	ELSE
465	CALL exchange_horiz( pr_av, nbgp )
466	to_be_resorted => pr_av
467	ENDIF
468
469	CASE ( 'pra*_xy' ) ! 2d-array / integral quantity => no av
470	CALL exchange_horiz_2d( precipitation_amount )
471	DO i = nxlg, nxrg
472	DO j = nysg, nyng
473	local_pf(i,j,nzb+1) = precipitation_amount(j,i)
474	ENDDO
475	ENDDO
476	precipitation_amount = 0.0_wp ! reset for next integ. interval
477	resorted = .TRUE.
478	two_d = .TRUE.
479	level_z(nzb+1) = zu(nzb+1)
480
481	CASE ( 'prr*_xy' ) ! 2d-array
482	IF ( icloud_scheme == 1 ) THEN
483	IF ( av == 0 ) THEN
484	CALL exchange_horiz_2d( precipitation_rate )
485	DO i = nxlg, nxrg
486	DO j = nysg, nyng
487	local_pf(i,j,nzb+1) = precipitation_rate(j,i)
488	ENDDO
489	ENDDO
490	ELSE
491	CALL exchange_horiz_2d( precipitation_rate_av )
492	DO i = nxlg, nxrg
493	DO j = nysg, nyng
494	local_pf(i,j,nzb+1) = precipitation_rate_av(j,i)
495	ENDDO
496	ENDDO
497	ENDIF
498	ELSE
499	IF ( av == 0 ) THEN
500	CALL exchange_horiz_2d( prr(nzb+1,:,:) )
501	DO i = nxlg, nxrg
502	DO j = nysg, nyng
503	local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
504	ENDDO
505	ENDDO
506	ELSE
507	CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
508	DO i = nxlg, nxrg
509	DO j = nysg, nyng
510	local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) * &
511	hyrho(nzb+1)
512	ENDDO
513	ENDDO
514	ENDIF
515	ENDIF
516	resorted = .TRUE.
517	two_d = .TRUE.
518	level_z(nzb+1) = zu(nzb+1)
519
520	CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
521	IF ( av == 0 ) THEN
522	CALL exchange_horiz( prr, nbgp )
523	DO i = nxlg, nxrg
524	DO j = nysg, nyng
525	DO k = nzb, nzt+1
526	local_pf(i,j,k) = prr(k,j,i)
527	ENDDO
528	ENDDO
529	ENDDO
530	ELSE
531	CALL exchange_horiz( prr_av, nbgp )
532	DO i = nxlg, nxrg
533	DO j = nysg, nyng
534	DO k = nzb, nzt+1
535	local_pf(i,j,k) = prr_av(k,j,i)
536	ENDDO
537	ENDDO
538	ENDDO
539	ENDIF
540	resorted = .TRUE.
541	IF ( mode == 'xy' ) level_z = zu
542
543	CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
544	IF ( av == 0 ) THEN
545	IF ( .NOT. cloud_physics ) THEN
546	to_be_resorted => pt
547	ELSE
548	DO i = nxlg, nxrg
549	DO j = nysg, nyng
550	DO k = nzb, nzt+1
551	local_pf(i,j,k) = pt(k,j,i) + l_d_cp * &
552	pt_d_t(k) * &
553	ql(k,j,i)
554	ENDDO
555	ENDDO
556	ENDDO
557	resorted = .TRUE.
558	ENDIF
559	ELSE
560	to_be_resorted => pt_av
561	ENDIF
562	IF ( mode == 'xy' ) level_z = zu
563
564	CASE ( 'q_xy', 'q_xz', 'q_yz' )
565	IF ( av == 0 ) THEN
566	to_be_resorted => q
567	ELSE
568	to_be_resorted => q_av
569	ENDIF
570	IF ( mode == 'xy' ) level_z = zu
571
572	CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
573	IF ( av == 0 ) THEN
574	to_be_resorted => qc
575	ELSE
576	to_be_resorted => qc_av
577	ENDIF
578	IF ( mode == 'xy' ) level_z = zu
579
580	CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
581	IF ( av == 0 ) THEN
582	to_be_resorted => ql
583	ELSE
584	to_be_resorted => ql_av
585	ENDIF
586	IF ( mode == 'xy' ) level_z = zu
587
588	CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
589	IF ( av == 0 ) THEN
590	to_be_resorted => ql_c
591	ELSE
592	to_be_resorted => ql_c_av
593	ENDIF
594	IF ( mode == 'xy' ) level_z = zu
595
596	CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
597	IF ( av == 0 ) THEN
598	to_be_resorted => ql_v
599	ELSE
600	to_be_resorted => ql_v_av
601	ENDIF
602	IF ( mode == 'xy' ) level_z = zu
603
604	CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
605	IF ( av == 0 ) THEN
606	IF ( simulated_time >= particle_advection_start ) THEN
607	DO i = nxl, nxr
608	DO j = nys, nyn
609	DO k = nzb, nzt+1
610	number_of_particles = prt_count(k,j,i)
611	IF (number_of_particles <= 0) CYCLE
612	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
613	DO n = 1, number_of_particles
614	IF ( particles(n)%particle_mask ) THEN
615	tend(k,j,i) = tend(k,j,i) + &
616	particles(n)%weight_factor / &
617	prt_count(k,j,i)
618	ENDIF
619	ENDDO
620	ENDDO
621	ENDDO
622	ENDDO
623	CALL exchange_horiz( tend, nbgp )
624	ELSE
625	tend = 0.0_wp
626	ENDIF
627	DO i = nxlg, nxrg
628	DO j = nysg, nyng
629	DO k = nzb, nzt+1
630	local_pf(i,j,k) = tend(k,j,i)
631	ENDDO
632	ENDDO
633	ENDDO
634	resorted = .TRUE.
635	ELSE
636	CALL exchange_horiz( ql_vp_av, nbgp )
637	to_be_resorted => ql_vp
638	ENDIF
639	IF ( mode == 'xy' ) level_z = zu
640
641	CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
642	IF ( av == 0 ) THEN
643	to_be_resorted => qr
644	ELSE
645	to_be_resorted => qr_av
646	ENDIF
647	IF ( mode == 'xy' ) level_z = zu
648
649	CASE ( 'qsws*_xy' ) ! 2d-array
650	IF ( av == 0 ) THEN
651	DO i = nxlg, nxrg
652	DO j = nysg, nyng
653	local_pf(i,j,nzb+1) = qsws(j,i)
654	ENDDO
655	ENDDO
656	ELSE
657	DO i = nxlg, nxrg
658	DO j = nysg, nyng
659	local_pf(i,j,nzb+1) = qsws_av(j,i)
660	ENDDO
661	ENDDO
662	ENDIF
663	resorted = .TRUE.
664	two_d = .TRUE.
665	level_z(nzb+1) = zu(nzb+1)
666
667	CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
668	IF ( av == 0 ) THEN
669	DO i = nxlg, nxrg
670	DO j = nysg, nyng
671	DO k = nzb, nzt+1
672	local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
673	ENDDO
674	ENDDO
675	ENDDO
676	resorted = .TRUE.
677	ELSE
678	to_be_resorted => qv_av
679	ENDIF
680	IF ( mode == 'xy' ) level_z = zu
681
682	CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
683	IF ( av == 0 ) THEN
684	to_be_resorted => rho
685	ELSE
686	to_be_resorted => rho_av
687	ENDIF
688
689	CASE ( 's_xy', 's_xz', 's_yz' )
690	IF ( av == 0 ) THEN
691	to_be_resorted => q
692	ELSE
693	to_be_resorted => s_av
694	ENDIF
695
696	CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
697	IF ( av == 0 ) THEN
698	to_be_resorted => sa
699	ELSE
700	to_be_resorted => sa_av
701	ENDIF
702
703	CASE ( 'shf*_xy' ) ! 2d-array
704	IF ( av == 0 ) THEN
705	DO i = nxlg, nxrg
706	DO j = nysg, nyng
707	local_pf(i,j,nzb+1) = shf(j,i)
708	ENDDO
709	ENDDO
710	ELSE
711	DO i = nxlg, nxrg
712	DO j = nysg, nyng
713	local_pf(i,j,nzb+1) = shf_av(j,i)
714	ENDDO
715	ENDDO
716	ENDIF
717	resorted = .TRUE.
718	two_d = .TRUE.
719	level_z(nzb+1) = zu(nzb+1)
720
721	CASE ( 't*_xy' ) ! 2d-array
722	IF ( av == 0 ) THEN
723	DO i = nxlg, nxrg
724	DO j = nysg, nyng
725	local_pf(i,j,nzb+1) = ts(j,i)
726	ENDDO
727	ENDDO
728	ELSE
729	DO i = nxlg, nxrg
730	DO j = nysg, nyng
731	local_pf(i,j,nzb+1) = ts_av(j,i)
732	ENDDO
733	ENDDO
734	ENDIF
735	resorted = .TRUE.
736	two_d = .TRUE.
737	level_z(nzb+1) = zu(nzb+1)
738
739	CASE ( 'u_xy', 'u_xz', 'u_yz' )
740	IF ( av == 0 ) THEN
741	to_be_resorted => u
742	ELSE
743	to_be_resorted => u_av
744	ENDIF
745	IF ( mode == 'xy' ) level_z = zu
746	!
747	!-- Substitute the values generated by "mirror" boundary condition
748	!-- at the bottom boundary by the real surface values.
749	IF ( do2d(av,if) == 'u_xz' .OR. do2d(av,if) == 'u_yz' ) THEN
750	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0_wp
751	ENDIF
752
753	CASE ( 'u*_xy' ) ! 2d-array
754	IF ( av == 0 ) THEN
755	DO i = nxlg, nxrg
756	DO j = nysg, nyng
757	local_pf(i,j,nzb+1) = us(j,i)
758	ENDDO
759	ENDDO
760	ELSE
761	DO i = nxlg, nxrg
762	DO j = nysg, nyng
763	local_pf(i,j,nzb+1) = us_av(j,i)
764	ENDDO
765	ENDDO
766	ENDIF
767	resorted = .TRUE.
768	two_d = .TRUE.
769	level_z(nzb+1) = zu(nzb+1)
770
771	CASE ( 'v_xy', 'v_xz', 'v_yz' )
772	IF ( av == 0 ) THEN
773	to_be_resorted => v
774	ELSE
775	to_be_resorted => v_av
776	ENDIF
777	IF ( mode == 'xy' ) level_z = zu
778	!
779	!-- Substitute the values generated by "mirror" boundary condition
780	!-- at the bottom boundary by the real surface values.
781	IF ( do2d(av,if) == 'v_xz' .OR. do2d(av,if) == 'v_yz' ) THEN
782	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0_wp
783	ENDIF
784
785	CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
786	IF ( av == 0 ) THEN
787	to_be_resorted => vpt
788	ELSE
789	to_be_resorted => vpt_av
790	ENDIF
791	IF ( mode == 'xy' ) level_z = zu
792
793	CASE ( 'w_xy', 'w_xz', 'w_yz' )
794	IF ( av == 0 ) THEN
795	to_be_resorted => w
796	ELSE
797	to_be_resorted => w_av
798	ENDIF
799	IF ( mode == 'xy' ) level_z = zw
800
801	CASE ( 'z0*_xy' ) ! 2d-array
802	IF ( av == 0 ) THEN
803	DO i = nxlg, nxrg
804	DO j = nysg, nyng
805	local_pf(i,j,nzb+1) = z0(j,i)
806	ENDDO
807	ENDDO
808	ELSE
809	DO i = nxlg, nxrg
810	DO j = nysg, nyng
811	local_pf(i,j,nzb+1) = z0_av(j,i)
812	ENDDO
813	ENDDO
814	ENDIF
815	resorted = .TRUE.
816	two_d = .TRUE.
817	level_z(nzb+1) = zu(nzb+1)
818
819	CASE ( 'z0h*_xy' ) ! 2d-array
820	IF ( av == 0 ) THEN
821	DO i = nxlg, nxrg
822	DO j = nysg, nyng
823	local_pf(i,j,nzb+1) = z0h(j,i)
824	ENDDO
825	ENDDO
826	ELSE
827	DO i = nxlg, nxrg
828	DO j = nysg, nyng
829	local_pf(i,j,nzb+1) = z0h_av(j,i)
830	ENDDO
831	ENDDO
832	ENDIF
833	resorted = .TRUE.
834	two_d = .TRUE.
835	level_z(nzb+1) = zu(nzb+1)
836
837	CASE DEFAULT
838	!
839	!-- User defined quantity
840	CALL user_data_output_2d( av, do2d(av,if), found, grid, &
841	local_pf, two_d )
842	resorted = .TRUE.
843
844	IF ( grid == 'zu' ) THEN
845	IF ( mode == 'xy' ) level_z = zu
846	ELSEIF ( grid == 'zw' ) THEN
847	IF ( mode == 'xy' ) level_z = zw
848	ELSEIF ( grid == 'zu1' ) THEN
849	IF ( mode == 'xy' ) level_z(nzb+1) = zu(nzb+1)
850	ENDIF
851
852	IF ( .NOT. found ) THEN
853	message_string = 'no output provided for: ' // &
854	TRIM( do2d(av,if) )
855	CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
856	ENDIF
857
858	END SELECT
859
860	!
861	!-- Resort the array to be output, if not done above
862	IF ( .NOT. resorted ) THEN
863	DO i = nxlg, nxrg
864	DO j = nysg, nyng
865	DO k = nzb, nzt+1
866	local_pf(i,j,k) = to_be_resorted(k,j,i)
867	ENDDO
868	ENDDO
869	ENDDO
870	ENDIF
871
872	!
873	!-- Output of the individual cross-sections, depending on the cross-
874	!-- section mode chosen.
875	is = 1
876	loop1: DO WHILE ( section(is,s) /= -9999 .OR. two_d )
877
878	SELECT CASE ( mode )
879
880	CASE ( 'xy' )
881	!
882	!-- Determine the cross section index
883	IF ( two_d ) THEN
884	layer_xy = nzb+1
885	ELSE
886	layer_xy = section(is,s)
887	ENDIF
888
889	!
890	!-- Update the netCDF xy cross section time axis.
891	!-- In case of parallel output, this is only done by PE0
892	!-- to increase the performance.
893	IF ( simulated_time /= do2d_xy_last_time(av) ) THEN
894	do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
895	do2d_xy_last_time(av) = simulated_time
896	IF ( myid == 0 ) THEN
897	IF ( .NOT. data_output_2d_on_each_pe &
898	.OR. netcdf_data_format > 4 ) &
899	THEN
900	#if defined( __netcdf )
901	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
902	id_var_time_xy(av), &
903	(/ time_since_reference_point /), &
904	start = (/ do2d_xy_time_count(av) /), &
905	count = (/ 1 /) )
906	CALL handle_netcdf_error( 'data_output_2d', 53 )
907	#endif
908	ENDIF
909	ENDIF
910	ENDIF
911	!
912	!-- If required, carry out averaging along z
913	IF ( section(is,s) == -1 .AND. .NOT. two_d ) THEN
914
915	local_2d = 0.0_wp
916	!
917	!-- Carry out the averaging (all data are on the PE)
918	DO k = nzb, nzt+1
919	DO j = nysg, nyng
920	DO i = nxlg, nxrg
921	local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
922	ENDDO
923	ENDDO
924	ENDDO
925
926	local_2d = local_2d / ( nzt -nzb + 2.0_wp)
927
928	ELSE
929	!
930	!-- Just store the respective section on the local array
931	local_2d = local_pf(:,:,layer_xy)
932
933	ENDIF
934
935	#if defined( __parallel )
936	IF ( netcdf_data_format > 4 ) THEN
937	!
938	!-- Parallel output in netCDF4/HDF5 format.
939	IF ( two_d ) THEN
940	iis = 1
941	ELSE
942	iis = is
943	ENDIF
944
945	#if defined( __netcdf )
946	!
947	!-- For parallel output, all cross sections are first stored
948	!-- here on a local array and will be written to the output
949	!-- file afterwards to increase the performance.
950	DO i = nxlg, nxrg
951	DO j = nysg, nyng
952	local_2d_sections(i,j,iis) = local_2d(i,j)
953	ENDDO
954	ENDDO
955	#endif
956	ELSE
957
958	IF ( data_output_2d_on_each_pe ) THEN
959	!
960	!-- Output of partial arrays on each PE
961	#if defined( __netcdf )
962	IF ( myid == 0 ) THEN
963	WRITE ( 21 ) time_since_reference_point, &
964	do2d_xy_time_count(av), av
965	ENDIF
966	#endif
967	DO i = 0, io_blocks-1
968	IF ( i == io_group ) THEN
969	WRITE ( 21 ) nxlg, nxrg, nysg, nyng
970	WRITE ( 21 ) local_2d
971	ENDIF
972	#if defined( __parallel )
973	CALL MPI_BARRIER( comm2d, ierr )
974	#endif
975	ENDDO
976
977	ELSE
978	!
979	!-- PE0 receives partial arrays from all processors and
980	!-- then outputs them. Here a barrier has to be set,
981	!-- because otherwise "-MPI- FATAL: Remote protocol queue
982	!-- full" may occur.
983	CALL MPI_BARRIER( comm2d, ierr )
984
985	ngp = ( nxrg-nxlg+1 ) * ( nyng-nysg+1 )
986	IF ( myid == 0 ) THEN
987	!
988	!-- Local array can be relocated directly.
989	total_2d(nxlg:nxrg,nysg:nyng) = local_2d
990	!
991	!-- Receive data from all other PEs.
992	DO n = 1, numprocs-1
993	!
994	!-- Receive index limits first, then array.
995	!-- Index limits are received in arbitrary order from
996	!-- the PEs.
997	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
998	MPI_ANY_SOURCE, 0, comm2d, &
999	status, ierr )
1000	sender = status(MPI_SOURCE)
1001	DEALLOCATE( local_2d )
1002	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1003	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1004	MPI_REAL, sender, 1, comm2d, &
1005	status, ierr )
1006	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1007	ENDDO
1008	!
1009	!-- Relocate the local array for the next loop increment
1010	DEALLOCATE( local_2d )
1011	ALLOCATE( local_2d(nxlg:nxrg,nysg:nyng) )
1012
1013	#if defined( __netcdf )
1014	IF ( two_d ) THEN
1015	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1016	id_var_do2d(av,if), &
1017	total_2d(0:nx+1,0:ny+1), &
1018	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1019	count = (/ nx+2, ny+2, 1, 1 /) )
1020	ELSE
1021	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1022	id_var_do2d(av,if), &
1023	total_2d(0:nx+1,0:ny+1), &
1024	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1025	count = (/ nx+2, ny+2, 1, 1 /) )
1026	ENDIF
1027	CALL handle_netcdf_error( 'data_output_2d', 54 )
1028	#endif
1029
1030	ELSE
1031	!
1032	!-- First send the local index limits to PE0
1033	ind(1) = nxlg; ind(2) = nxrg
1034	ind(3) = nysg; ind(4) = nyng
1035	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1036	comm2d, ierr )
1037	!
1038	!-- Send data to PE0
1039	CALL MPI_SEND( local_2d(nxlg,nysg), ngp, &
1040	MPI_REAL, 0, 1, comm2d, ierr )
1041	ENDIF
1042	!
1043	!-- A barrier has to be set, because otherwise some PEs may
1044	!-- proceed too fast so that PE0 may receive wrong data on
1045	!-- tag 0
1046	CALL MPI_BARRIER( comm2d, ierr )
1047	ENDIF
1048
1049	ENDIF
1050	#else
1051	#if defined( __netcdf )
1052	IF ( two_d ) THEN
1053	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1054	id_var_do2d(av,if), &
1055	local_2d(nxl:nxr+1,nys:nyn+1), &
1056	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1057	count = (/ nx+2, ny+2, 1, 1 /) )
1058	ELSE
1059	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1060	id_var_do2d(av,if), &
1061	local_2d(nxl:nxr+1,nys:nyn+1), &
1062	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1063	count = (/ nx+2, ny+2, 1, 1 /) )
1064	ENDIF
1065	CALL handle_netcdf_error( 'data_output_2d', 447 )
1066	#endif
1067	#endif
1068	do2d_xy_n = do2d_xy_n + 1
1069	!
1070	!-- For 2D-arrays (e.g. u*) only one cross-section is available.
1071	!-- Hence exit loop of output levels.
1072	IF ( two_d ) THEN
1073	two_d = .FALSE.
1074	EXIT loop1
1075	ENDIF
1076
1077	CASE ( 'xz' )
1078	!
1079	!-- Update the netCDF xz cross section time axis.
1080	!-- In case of parallel output, this is only done by PE0
1081	!-- to increase the performance.
1082	IF ( simulated_time /= do2d_xz_last_time(av) ) THEN
1083	do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1084	do2d_xz_last_time(av) = simulated_time
1085	IF ( myid == 0 ) THEN
1086	IF ( .NOT. data_output_2d_on_each_pe &
1087	.OR. netcdf_data_format > 4 ) &
1088	THEN
1089	#if defined( __netcdf )
1090	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1091	id_var_time_xz(av), &
1092	(/ time_since_reference_point /), &
1093	start = (/ do2d_xz_time_count(av) /), &
1094	count = (/ 1 /) )
1095	CALL handle_netcdf_error( 'data_output_2d', 56 )
1096	#endif
1097	ENDIF
1098	ENDIF
1099	ENDIF
1100
1101	!
1102	!-- If required, carry out averaging along y
1103	IF ( section(is,s) == -1 ) THEN
1104
1105	ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
1106	local_2d_l = 0.0_wp
1107	ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
1108	!
1109	!-- First local averaging on the PE
1110	DO k = nzb, nzt+1
1111	DO j = nys, nyn
1112	DO i = nxlg, nxrg
1113	local_2d_l(i,k) = local_2d_l(i,k) + &
1114	local_pf(i,j,k)
1115	ENDDO
1116	ENDDO
1117	ENDDO
1118	#if defined( __parallel )
1119	!
1120	!-- Now do the averaging over all PEs along y
1121	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1122	CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb), &
1123	local_2d(nxlg,nzb), ngp, MPI_REAL, &
1124	MPI_SUM, comm1dy, ierr )
1125	#else
1126	local_2d = local_2d_l
1127	#endif
1128	local_2d = local_2d / ( ny + 1.0_wp )
1129
1130	DEALLOCATE( local_2d_l )
1131
1132	ELSE
1133	!
1134	!-- Just store the respective section on the local array
1135	!-- (but only if it is available on this PE!)
1136	IF ( section(is,s) >= nys .AND. section(is,s) <= nyn ) &
1137	THEN
1138	local_2d = local_pf(:,section(is,s),nzb:nzt+1)
1139	ENDIF
1140
1141	ENDIF
1142
1143	#if defined( __parallel )
1144	IF ( netcdf_data_format > 4 ) THEN
1145	!
1146	!-- Output in netCDF4/HDF5 format.
1147	!-- Output only on those PEs where the respective cross
1148	!-- sections reside. Cross sections averaged along y are
1149	!-- output on the respective first PE along y (myidy=0).
1150	IF ( ( section(is,s) >= nys .AND. &
1151	section(is,s) <= nyn ) .OR. &
1152	( section(is,s) == -1 .AND. myidy == 0 ) ) THEN
1153	#if defined( __netcdf )
1154	!
1155	!-- For parallel output, all cross sections are first
1156	!-- stored here on a local array and will be written to the
1157	!-- output file afterwards to increase the performance.
1158	DO i = nxlg, nxrg
1159	DO k = nzb, nzt+1
1160	local_2d_sections_l(i,is,k) = local_2d(i,k)
1161	ENDDO
1162	ENDDO
1163	#endif
1164	ENDIF
1165
1166	ELSE
1167
1168	IF ( data_output_2d_on_each_pe ) THEN
1169	!
1170	!-- Output of partial arrays on each PE. If the cross
1171	!-- section does not reside on the PE, output special
1172	!-- index values.
1173	#if defined( __netcdf )
1174	IF ( myid == 0 ) THEN
1175	WRITE ( 22 ) time_since_reference_point, &
1176	do2d_xz_time_count(av), av
1177	ENDIF
1178	#endif
1179	DO i = 0, io_blocks-1
1180	IF ( i == io_group ) THEN
1181	IF ( ( section(is,s) >= nys .AND. &
1182	section(is,s) <= nyn ) .OR. &
1183	( section(is,s) == -1 .AND. &
1184	nys-1 == -1 ) ) &
1185	THEN
1186	WRITE (22) nxlg, nxrg, nzb, nzt+1
1187	WRITE (22) local_2d
1188	ELSE
1189	WRITE (22) -1, -1, -1, -1
1190	ENDIF
1191	ENDIF
1192	#if defined( __parallel )
1193	CALL MPI_BARRIER( comm2d, ierr )
1194	#endif
1195	ENDDO
1196
1197	ELSE
1198	!
1199	!-- PE0 receives partial arrays from all processors of the
1200	!-- respective cross section and outputs them. Here a
1201	!-- barrier has to be set, because otherwise
1202	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1203	CALL MPI_BARRIER( comm2d, ierr )
1204
1205	ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
1206	IF ( myid == 0 ) THEN
1207	!
1208	!-- Local array can be relocated directly.
1209	IF ( ( section(is,s) >= nys .AND. &
1210	section(is,s) <= nyn ) .OR. &
1211	( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
1212	THEN
1213	total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
1214	ENDIF
1215	!
1216	!-- Receive data from all other PEs.
1217	DO n = 1, numprocs-1
1218	!
1219	!-- Receive index limits first, then array.
1220	!-- Index limits are received in arbitrary order from
1221	!-- the PEs.
1222	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1223	MPI_ANY_SOURCE, 0, comm2d, &
1224	status, ierr )
1225	!
1226	!-- Not all PEs have data for XZ-cross-section.
1227	IF ( ind(1) /= -9999 ) THEN
1228	sender = status(MPI_SOURCE)
1229	DEALLOCATE( local_2d )
1230	ALLOCATE( local_2d(ind(1):ind(2), &
1231	ind(3):ind(4)) )
1232	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1233	MPI_REAL, sender, 1, comm2d, &
1234	status, ierr )
1235	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1236	local_2d
1237	ENDIF
1238	ENDDO
1239	!
1240	!-- Relocate the local array for the next loop increment
1241	DEALLOCATE( local_2d )
1242	ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
1243
1244	#if defined( __netcdf )
1245	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1246	id_var_do2d(av,if), &
1247	total_2d(0:nx+1,nzb:nzt+1),&
1248	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1249	count = (/ nx+2, 1, nz+2, 1 /) )
1250	CALL handle_netcdf_error( 'data_output_2d', 58 )
1251	#endif
1252
1253	ELSE
1254	!
1255	!-- If the cross section resides on the PE, send the
1256	!-- local index limits, otherwise send -9999 to PE0.
1257	IF ( ( section(is,s) >= nys .AND. &
1258	section(is,s) <= nyn ) .OR. &
1259	( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
1260	THEN
1261	ind(1) = nxlg; ind(2) = nxrg
1262	ind(3) = nzb; ind(4) = nzt+1
1263	ELSE
1264	ind(1) = -9999; ind(2) = -9999
1265	ind(3) = -9999; ind(4) = -9999
1266	ENDIF
1267	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1268	comm2d, ierr )
1269	!
1270	!-- If applicable, send data to PE0.
1271	IF ( ind(1) /= -9999 ) THEN
1272	CALL MPI_SEND( local_2d(nxlg,nzb), ngp, &
1273	MPI_REAL, 0, 1, comm2d, ierr )
1274	ENDIF
1275	ENDIF
1276	!
1277	!-- A barrier has to be set, because otherwise some PEs may
1278	!-- proceed too fast so that PE0 may receive wrong data on
1279	!-- tag 0
1280	CALL MPI_BARRIER( comm2d, ierr )
1281	ENDIF
1282
1283	ENDIF
1284	#else
1285	#if defined( __netcdf )
1286	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1287	id_var_do2d(av,if), &
1288	local_2d(nxl:nxr+1,nzb:nzt+1), &
1289	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1290	count = (/ nx+2, 1, nz+2, 1 /) )
1291	CALL handle_netcdf_error( 'data_output_2d', 451 )
1292	#endif
1293	#endif
1294	do2d_xz_n = do2d_xz_n + 1
1295
1296	CASE ( 'yz' )
1297	!
1298	!-- Update the netCDF yz cross section time axis.
1299	!-- In case of parallel output, this is only done by PE0
1300	!-- to increase the performance.
1301	IF ( simulated_time /= do2d_yz_last_time(av) ) THEN
1302	do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1303	do2d_yz_last_time(av) = simulated_time
1304	IF ( myid == 0 ) THEN
1305	IF ( .NOT. data_output_2d_on_each_pe &
1306	.OR. netcdf_data_format > 4 ) &
1307	THEN
1308	#if defined( __netcdf )
1309	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1310	id_var_time_yz(av), &
1311	(/ time_since_reference_point /), &
1312	start = (/ do2d_yz_time_count(av) /), &
1313	count = (/ 1 /) )
1314	CALL handle_netcdf_error( 'data_output_2d', 59 )
1315	#endif
1316	ENDIF
1317	ENDIF
1318	ENDIF
1319
1320	!
1321	!-- If required, carry out averaging along x
1322	IF ( section(is,s) == -1 ) THEN
1323
1324	ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
1325	local_2d_l = 0.0_wp
1326	ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
1327	!
1328	!-- First local averaging on the PE
1329	DO k = nzb, nzt+1
1330	DO j = nysg, nyng
1331	DO i = nxl, nxr
1332	local_2d_l(j,k) = local_2d_l(j,k) + &
1333	local_pf(i,j,k)
1334	ENDDO
1335	ENDDO
1336	ENDDO
1337	#if defined( __parallel )
1338	!
1339	!-- Now do the averaging over all PEs along x
1340	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1341	CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb), &
1342	local_2d(nysg,nzb), ngp, MPI_REAL, &
1343	MPI_SUM, comm1dx, ierr )
1344	#else
1345	local_2d = local_2d_l
1346	#endif
1347	local_2d = local_2d / ( nx + 1.0_wp )
1348
1349	DEALLOCATE( local_2d_l )
1350
1351	ELSE
1352	!
1353	!-- Just store the respective section on the local array
1354	!-- (but only if it is available on this PE!)
1355	IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr ) &
1356	THEN
1357	local_2d = local_pf(section(is,s),:,nzb:nzt+1)
1358	ENDIF
1359
1360	ENDIF
1361
1362	#if defined( __parallel )
1363	IF ( netcdf_data_format > 4 ) THEN
1364	!
1365	!-- Output in netCDF4/HDF5 format.
1366	!-- Output only on those PEs where the respective cross
1367	!-- sections reside. Cross sections averaged along x are
1368	!-- output on the respective first PE along x (myidx=0).
1369	IF ( ( section(is,s) >= nxl .AND. &
1370	section(is,s) <= nxr ) .OR. &
1371	( section(is,s) == -1 .AND. myidx == 0 ) ) THEN
1372	#if defined( __netcdf )
1373	!
1374	!-- For parallel output, all cross sections are first
1375	!-- stored here on a local array and will be written to the
1376	!-- output file afterwards to increase the performance.
1377	DO j = nysg, nyng
1378	DO k = nzb, nzt+1
1379	local_2d_sections_l(is,j,k) = local_2d(j,k)
1380	ENDDO
1381	ENDDO
1382	#endif
1383	ENDIF
1384
1385	ELSE
1386
1387	IF ( data_output_2d_on_each_pe ) THEN
1388	!
1389	!-- Output of partial arrays on each PE. If the cross
1390	!-- section does not reside on the PE, output special
1391	!-- index values.
1392	#if defined( __netcdf )
1393	IF ( myid == 0 ) THEN
1394	WRITE ( 23 ) time_since_reference_point, &
1395	do2d_yz_time_count(av), av
1396	ENDIF
1397	#endif
1398	DO i = 0, io_blocks-1
1399	IF ( i == io_group ) THEN
1400	IF ( ( section(is,s) >= nxl .AND. &
1401	section(is,s) <= nxr ) .OR. &
1402	( section(is,s) == -1 .AND. &
1403	nxl-1 == -1 ) ) &
1404	THEN
1405	WRITE (23) nysg, nyng, nzb, nzt+1
1406	WRITE (23) local_2d
1407	ELSE
1408	WRITE (23) -1, -1, -1, -1
1409	ENDIF
1410	ENDIF
1411	#if defined( __parallel )
1412	CALL MPI_BARRIER( comm2d, ierr )
1413	#endif
1414	ENDDO
1415
1416	ELSE
1417	!
1418	!-- PE0 receives partial arrays from all processors of the
1419	!-- respective cross section and outputs them. Here a
1420	!-- barrier has to be set, because otherwise
1421	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1422	CALL MPI_BARRIER( comm2d, ierr )
1423
1424	ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
1425	IF ( myid == 0 ) THEN
1426	!
1427	!-- Local array can be relocated directly.
1428	IF ( ( section(is,s) >= nxl .AND. &
1429	section(is,s) <= nxr ) .OR. &
1430	( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1431	THEN
1432	total_2d(nysg:nyng,nzb:nzt+1) = local_2d
1433	ENDIF
1434	!
1435	!-- Receive data from all other PEs.
1436	DO n = 1, numprocs-1
1437	!
1438	!-- Receive index limits first, then array.
1439	!-- Index limits are received in arbitrary order from
1440	!-- the PEs.
1441	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1442	MPI_ANY_SOURCE, 0, comm2d, &
1443	status, ierr )
1444	!
1445	!-- Not all PEs have data for YZ-cross-section.
1446	IF ( ind(1) /= -9999 ) THEN
1447	sender = status(MPI_SOURCE)
1448	DEALLOCATE( local_2d )
1449	ALLOCATE( local_2d(ind(1):ind(2), &
1450	ind(3):ind(4)) )
1451	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1452	MPI_REAL, sender, 1, comm2d, &
1453	status, ierr )
1454	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1455	local_2d
1456	ENDIF
1457	ENDDO
1458	!
1459	!-- Relocate the local array for the next loop increment
1460	DEALLOCATE( local_2d )
1461	ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
1462
1463	#if defined( __netcdf )
1464	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1465	id_var_do2d(av,if), &
1466	total_2d(0:ny+1,nzb:nzt+1),&
1467	start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1468	count = (/ 1, ny+2, nz+2, 1 /) )
1469	CALL handle_netcdf_error( 'data_output_2d', 61 )
1470	#endif
1471
1472	ELSE
1473	!
1474	!-- If the cross section resides on the PE, send the
1475	!-- local index limits, otherwise send -9999 to PE0.
1476	IF ( ( section(is,s) >= nxl .AND. &
1477	section(is,s) <= nxr ) .OR. &
1478	( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1479	THEN
1480	ind(1) = nysg; ind(2) = nyng
1481	ind(3) = nzb; ind(4) = nzt+1
1482	ELSE
1483	ind(1) = -9999; ind(2) = -9999
1484	ind(3) = -9999; ind(4) = -9999
1485	ENDIF
1486	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1487	comm2d, ierr )
1488	!
1489	!-- If applicable, send data to PE0.
1490	IF ( ind(1) /= -9999 ) THEN
1491	CALL MPI_SEND( local_2d(nysg,nzb), ngp, &
1492	MPI_REAL, 0, 1, comm2d, ierr )
1493	ENDIF
1494	ENDIF
1495	!
1496	!-- A barrier has to be set, because otherwise some PEs may
1497	!-- proceed too fast so that PE0 may receive wrong data on
1498	!-- tag 0
1499	CALL MPI_BARRIER( comm2d, ierr )
1500	ENDIF
1501
1502	ENDIF
1503	#else
1504	#if defined( __netcdf )
1505	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1506	id_var_do2d(av,if), &
1507	local_2d(nys:nyn+1,nzb:nzt+1), &
1508	start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1509	count = (/ 1, ny+2, nz+2, 1 /) )
1510	CALL handle_netcdf_error( 'data_output_2d', 452 )
1511	#endif
1512	#endif
1513	do2d_yz_n = do2d_yz_n + 1
1514
1515	END SELECT
1516
1517	is = is + 1
1518	ENDDO loop1
1519
1520	!
1521	!-- For parallel output, all data were collected before on a local array
1522	!-- and are written now to the netcdf file. This must be done to increase
1523	!-- the performance of the parallel output.
1524	#if defined( __netcdf )
1525	IF ( netcdf_data_format > 4 ) THEN
1526
1527	SELECT CASE ( mode )
1528
1529	CASE ( 'xy' )
1530	IF ( two_d ) THEN
1531	iis = 1
1532	ELSE
1533	iis = is-1
1534	ENDIF
1535	!
1536	!-- Do not output redundant ghost point data except for the
1537	!-- boundaries of the total domain.
1538	IF ( nxr == nx .AND. nyn /= ny ) THEN
1539	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1540	id_var_do2d(av,if), &
1541	local_2d_sections(nxl:nxr+1, &
1542	nys:nyn,1:ns), &
1543	start = (/ nxl+1, nys+1, 1, &
1544	do2d_xy_time_count(av) /), &
1545	count = (/ nxr-nxl+2, &
1546	nyn-nys+1, ns, 1 &
1547	/) )
1548	ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
1549	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1550	id_var_do2d(av,if), &
1551	local_2d_sections(nxl:nxr, &
1552	nys:nyn+1,1:ns), &
1553	start = (/ nxl+1, nys+1, 1, &
1554	do2d_xy_time_count(av) /), &
1555	count = (/ nxr-nxl+1, &
1556	nyn-nys+2, ns, 1 &
1557	/) )
1558	ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
1559	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1560	id_var_do2d(av,if), &
1561	local_2d_sections(nxl:nxr+1, &
1562	nys:nyn+1,1:ns), &
1563	start = (/ nxl+1, nys+1, 1, &
1564	do2d_xy_time_count(av) /), &
1565	count = (/ nxr-nxl+2, &
1566	nyn-nys+2, ns, 1 &
1567	/) )
1568	ELSE
1569	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1570	id_var_do2d(av,if), &
1571	local_2d_sections(nxl:nxr, &
1572	nys:nyn,1:ns), &
1573	start = (/ nxl+1, nys+1, 1, &
1574	do2d_xy_time_count(av) /), &
1575	count = (/ nxr-nxl+1, &
1576	nyn-nys+1, ns, 1 &
1577	/) )
1578	ENDIF
1579
1580	CALL handle_netcdf_error( 'data_output_2d', 55 )
1581
1582	CASE ( 'xz' )
1583	!
1584	!-- First, all PEs get the information of all cross-sections.
1585	!-- Then the data are written to the output file by all PEs
1586	!-- while NF90_COLLECTIVE is set in subroutine
1587	!-- define_netcdf_header. Although redundant information are
1588	!-- written to the output file in that case, the performance
1589	!-- is significantly better compared to the case where only
1590	!-- the first row of PEs in x-direction (myidx = 0) is given
1591	!-- the output while NF90_INDEPENDENT is set.
1592	IF ( npey /= 1 ) THEN
1593
1594	#if defined( __parallel )
1595	!
1596	!-- Distribute data over all PEs along y
1597	ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) * ns
1598	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1599	CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb), &
1600	local_2d_sections(nxlg,1,nzb), &
1601	ngp, MPI_REAL, MPI_SUM, comm1dy, &
1602	ierr )
1603	#else
1604	local_2d_sections = local_2d_sections_l
1605	#endif
1606	ENDIF
1607	!
1608	!-- Do not output redundant ghost point data except for the
1609	!-- boundaries of the total domain.
1610	IF ( nxr == nx ) THEN
1611	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1612	id_var_do2d(av,if), &
1613	local_2d_sections(nxl:nxr+1,1:ns, &
1614	nzb:nzt+1), &
1615	start = (/ nxl+1, 1, 1, &
1616	do2d_xz_time_count(av) /), &
1617	count = (/ nxr-nxl+2, ns, nzt+2, &
1618	1 /) )
1619	ELSE
1620	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1621	id_var_do2d(av,if), &
1622	local_2d_sections(nxl:nxr,1:ns, &
1623	nzb:nzt+1), &
1624	start = (/ nxl+1, 1, 1, &
1625	do2d_xz_time_count(av) /), &
1626	count = (/ nxr-nxl+1, ns, nzt+2, &
1627	1 /) )
1628	ENDIF
1629
1630	CALL handle_netcdf_error( 'data_output_2d', 57 )
1631
1632	CASE ( 'yz' )
1633	!
1634	!-- First, all PEs get the information of all cross-sections.
1635	!-- Then the data are written to the output file by all PEs
1636	!-- while NF90_COLLECTIVE is set in subroutine
1637	!-- define_netcdf_header. Although redundant information are
1638	!-- written to the output file in that case, the performance
1639	!-- is significantly better compared to the case where only
1640	!-- the first row of PEs in y-direction (myidy = 0) is given
1641	!-- the output while NF90_INDEPENDENT is set.
1642	IF ( npex /= 1 ) THEN
1643
1644	#if defined( __parallel )
1645	!
1646	!-- Distribute data over all PEs along x
1647	ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) * ns
1648	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1649	CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb), &
1650	local_2d_sections(1,nysg,nzb), &
1651	ngp, MPI_REAL, MPI_SUM, comm1dx, &
1652	ierr )
1653	#else
1654	local_2d_sections = local_2d_sections_l
1655	#endif
1656	ENDIF
1657	!
1658	!-- Do not output redundant ghost point data except for the
1659	!-- boundaries of the total domain.
1660	IF ( nyn == ny ) THEN
1661	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1662	id_var_do2d(av,if), &
1663	local_2d_sections(1:ns, &
1664	nys:nyn+1,nzb:nzt+1), &
1665	start = (/ 1, nys+1, 1, &
1666	do2d_yz_time_count(av) /), &
1667	count = (/ ns, nyn-nys+2, &
1668	nzt+2, 1 /) )
1669	ELSE
1670	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1671	id_var_do2d(av,if), &
1672	local_2d_sections(1:ns,nys:nyn, &
1673	nzb:nzt+1), &
1674	start = (/ 1, nys+1, 1, &
1675	do2d_yz_time_count(av) /), &
1676	count = (/ ns, nyn-nys+1, &
1677	nzt+2, 1 /) )
1678	ENDIF
1679
1680	CALL handle_netcdf_error( 'data_output_2d', 60 )
1681
1682	CASE DEFAULT
1683	message_string = 'unknown cross-section: ' // TRIM( mode )
1684	CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
1685
1686	END SELECT
1687
1688	ENDIF
1689	#endif
1690	ENDIF
1691
1692	if = if + 1
1693	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
1694	do2d_mode = do2d(av,if)(l-1:l)
1695
1696	ENDDO
1697
1698	!
1699	!-- Deallocate temporary arrays.
1700	IF ( ALLOCATED( level_z ) ) DEALLOCATE( level_z )
1701	IF ( netcdf_data_format > 4 ) THEN
1702	DEALLOCATE( local_pf, local_2d, local_2d_sections )
1703	IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
1704	ENDIF
1705	#if defined( __parallel )
1706	IF ( .NOT. data_output_2d_on_each_pe .AND. myid == 0 ) THEN
1707	DEALLOCATE( total_2d )
1708	ENDIF
1709	#endif
1710
1711	!
1712	!-- Close plot output file.
1713	file_id = 20 + s
1714
1715	IF ( data_output_2d_on_each_pe ) THEN
1716	DO i = 0, io_blocks-1
1717	IF ( i == io_group ) THEN
1718	CALL close_file( file_id )
1719	ENDIF
1720	#if defined( __parallel )
1721	CALL MPI_BARRIER( comm2d, ierr )
1722	#endif
1723	ENDDO
1724	ELSE
1725	IF ( myid == 0 ) CALL close_file( file_id )
1726	ENDIF
1727
1728	CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
1729
1730	END SUBROUTINE data_output_2d

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |