Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |