Home

Context Navigation

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

Last change on this file since 1319 was 1319, checked in by raasch, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 75.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |