Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |