Home

Context Navigation

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

Last change on this file since 4467 was 4457, checked in by raasch, 5 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 89.0 KB

Line
1	!> @file data_output_2d.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! 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-2020 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: data_output_2d.f90 4457 2020-03-11 14:20:43Z maronga $
27	! use statement for exchange horiz added
28	!
29	! 4444 2020-03-05 15:59:50Z raasch
30	! bugfix: cpp-directives for serial mode added
31	!
32	! 4442 2020-03-04 19:21:13Z suehring
33	! Change order of dimension in surface array %frac to allow for better
34	! vectorization.
35	!
36	! 4441 2020-03-04 19:20:35Z suehring
37	! Introduction of wall_flags_total_0, which currently sets bits based on static
38	! topography information used in wall_flags_static_0
39	!
40	! 4331 2019-12-10 18:25:02Z suehring
41	! Move 2-m potential temperature output to diagnostic_output_quantities
42	!
43	! 4329 2019-12-10 15:46:36Z motisi
44	! Renamed wall_flags_0 to wall_flags_static_0
45	!
46	! 4182 2019-08-22 15:20:23Z scharf
47	! Corrected "Former revisions" section
48	!
49	! 4048 2019-06-21 21:00:21Z knoop
50	! Removed turbulence_closure_mod dependency
51	!
52	! 4039 2019-06-18 10:32:41Z suehring
53	! modularize diagnostic output
54	!
55	! 3994 2019-05-22 18:08:09Z suehring
56	! output of turbulence intensity added
57	!
58	! 3987 2019-05-22 09:52:13Z kanani
59	! Introduce alternative switch for debug output during timestepping
60	!
61	! 3943 2019-05-02 09:50:41Z maronga
62	! Added output of qsws for green roofs.
63	!
64	! 3885 2019-04-11 11:29:34Z kanani
65	! Changes related to global restructuring of location messages and introduction
66	! of additional debug messages
67	!
68	! 3766 2019-02-26 16:23:41Z raasch
69	! unused variables removed
70	!
71	! 3655 2019-01-07 16:51:22Z knoop
72	! Bugfix: use time_since_reference_point instead of simulated_time (relevant
73	! when using wall/soil spinup)
74	!
75	! Revision 1.1 1997/08/11 06:24:09 raasch
76	! Initial revision
77	!
78	!
79	! Description:
80	! ------------
81	!> Data output of cross-sections in netCDF format or binary format
82	!> to be later converted to NetCDF by helper routine combine_plot_fields.
83	!> Attention: The position of the sectional planes is still not always computed
84	!> --------- correctly. (zu is used always)!
85	!------------------------------------------------------------------------------!
86	SUBROUTINE data_output_2d( mode, av )
87
88
89	USE arrays_3d, &
90	ONLY: dzw, d_exner, e, heatflux_output_conversion, p, pt, q, ql, ql_c, ql_v, s, tend, u, &
91	v, vpt, w, waterflux_output_conversion, zu, zw
92
93	USE averaging
94
95	USE basic_constants_and_equations_mod, &
96	ONLY: c_p, lv_d_cp, l_v
97
98	USE bulk_cloud_model_mod, &
99	ONLY: bulk_cloud_model
100
101	USE control_parameters, &
102	ONLY: data_output_2d_on_each_pe, &
103	data_output_xy, data_output_xz, data_output_yz, &
104	debug_output_timestep, &
105	do2d, &
106	do2d_xy_last_time, do2d_xy_time_count, &
107	do2d_xz_last_time, do2d_xz_time_count, &
108	do2d_yz_last_time, do2d_yz_time_count, &
109	ibc_uv_b, io_blocks, io_group, message_string, &
110	ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz, &
111	psolver, section, &
112	time_since_reference_point
113
114	USE cpulog, &
115	ONLY: cpu_log, log_point
116
117	USE exchange_horiz_mod, &
118	ONLY: exchange_horiz
119
120	USE indices, &
121	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
122	nzb, nzt, wall_flags_total_0
123
124	USE kinds
125
126	USE land_surface_model_mod, &
127	ONLY: zs
128
129	USE module_interface, &
130	ONLY: module_interface_data_output_2d
131
132	#if defined( __netcdf )
133	USE NETCDF
134	#endif
135
136	USE netcdf_interface, &
137	ONLY: fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d, &
138	id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat, &
139	netcdf_data_format, netcdf_handle_error
140
141	USE particle_attributes, &
142	ONLY: grid_particles, number_of_particles, particle_advection_start, &
143	particles, prt_count
144
145	USE pegrid
146
147	USE surface_mod, &
148	ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h, &
149	surf_lsm_h, surf_usm_h
150
151
152	IMPLICIT NONE
153
154	CHARACTER (LEN=2) :: do2d_mode !< output mode of variable ('xy', 'xz', 'yz')
155	CHARACTER (LEN=2) :: mode !< mode with which the routine is called ('xy', 'xz', 'yz')
156	CHARACTER (LEN=4) :: grid !< string defining the vertical grid
157
158	INTEGER(iwp) :: av !< flag for (non-)average output
159	INTEGER(iwp) :: ngp !< number of grid points of an output slice
160	INTEGER(iwp) :: file_id !< id of output files
161	INTEGER(iwp) :: flag_nr !< number of masking flag
162	INTEGER(iwp) :: i !< loop index
163	INTEGER(iwp) :: is !< slice index
164	INTEGER(iwp) :: ivar !< variable index
165	INTEGER(iwp) :: j !< loop index
166	INTEGER(iwp) :: k !< loop index
167	INTEGER(iwp) :: l !< loop index
168	INTEGER(iwp) :: layer_xy !< vertical index of a xy slice in array 'local_pf'
169	INTEGER(iwp) :: m !< loop index
170	INTEGER(iwp) :: n !< loop index
171	INTEGER(iwp) :: nis !< number of vertical slices to be written via parallel NetCDF output
172	INTEGER(iwp) :: ns !< number of output slices
173	INTEGER(iwp) :: nzb_do !< lower limit of the data field (usually nzb)
174	INTEGER(iwp) :: nzt_do !< upper limit of the data field (usually nzt+1)
175	INTEGER(iwp) :: s_ind !< index of slice types (xy=1, xz=2, yz=3)
176	#if defined( __parallel )
177	INTEGER(iwp) :: iis !< vertical index of a xy slice in array 'local_2d_sections'
178	INTEGER(iwp) :: sender !< PE id of sending PE
179	INTEGER(iwp) :: ind(4) !< index limits (lower/upper bounds) of array 'local_2d'
180	#endif
181
182	LOGICAL :: found !< true if output variable was found
183	LOGICAL :: resorted !< true if variable is resorted
184	LOGICAL :: two_d !< true if variable is only two dimensional
185
186	REAL(wp) :: mean_r !< mean particle radius
187	REAL(wp) :: s_r2 !< sum( particle-radius**2 )
188	REAL(wp) :: s_r3 !< sum( particle-radius**3 )
189
190	REAL(wp), DIMENSION(:), ALLOCATABLE :: level_z !< z levels for output array
191	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: local_2d !< local 2-dimensional array containing output values
192	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: local_2d_l !< local 2-dimensional array containing output values
193
194	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !< output array
195	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_2d_sections !< local array containing values at all slices
196	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_2d_sections_l !< local array containing values at all slices
197
198	#if defined( __parallel )
199	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: total_2d !< same as local_2d
200	#endif
201	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which shall be output
202
203
204	IF ( debug_output_timestep ) CALL debug_message( 'data_output_2d', 'start' )
205	!
206	!-- Immediate return, if no output is requested (no respective sections
207	!-- found in parameter data_output)
208	IF ( mode == 'xy' .AND. .NOT. data_output_xy(av) ) RETURN
209	IF ( mode == 'xz' .AND. .NOT. data_output_xz(av) ) RETURN
210	IF ( mode == 'yz' .AND. .NOT. data_output_yz(av) ) RETURN
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_ind = 1
224	ALLOCATE( level_z(nzb:nzt+1), local_2d(nxl:nxr,nys:nyn) )
225
226	IF ( netcdf_data_format > 4 ) THEN
227	ns = 1
228	DO WHILE ( section(ns,s_ind) /= -9999 .AND. ns <= 100 )
229	ns = ns + 1
230	ENDDO
231	ns = ns - 1
232	ALLOCATE( local_2d_sections(nxl:nxr,nys:nyn,1:ns) )
233	local_2d_sections = 0.0_wp
234	ENDIF
235
236	!
237	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
238	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
239	CALL check_open( 101+av*10 )
240	ENDIF
241	IF ( data_output_2d_on_each_pe .AND. netcdf_data_format < 5 ) THEN
242	CALL check_open( 21 )
243	ELSE
244	IF ( myid == 0 ) THEN
245	#if defined( __parallel )
246	ALLOCATE( total_2d(0:nx,0:ny) )
247	#endif
248	ENDIF
249	ENDIF
250
251	CASE ( 'xz' )
252	s_ind = 2
253	ALLOCATE( local_2d(nxl:nxr,nzb:nzt+1) )
254
255	IF ( netcdf_data_format > 4 ) THEN
256	ns = 1
257	DO WHILE ( section(ns,s_ind) /= -9999 .AND. ns <= 100 )
258	ns = ns + 1
259	ENDDO
260	ns = ns - 1
261	ALLOCATE( local_2d_sections(nxl:nxr,1:ns,nzb:nzt+1) )
262	ALLOCATE( local_2d_sections_l(nxl:nxr,1:ns,nzb:nzt+1) )
263	local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
264	ENDIF
265
266	!
267	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
268	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
269	CALL check_open( 102+av*10 )
270	ENDIF
271
272	IF ( data_output_2d_on_each_pe .AND. netcdf_data_format < 5 ) THEN
273	CALL check_open( 22 )
274	ELSE
275	IF ( myid == 0 ) THEN
276	#if defined( __parallel )
277	ALLOCATE( total_2d(0:nx,nzb:nzt+1) )
278	#endif
279	ENDIF
280	ENDIF
281
282	CASE ( 'yz' )
283	s_ind = 3
284	ALLOCATE( local_2d(nys:nyn,nzb:nzt+1) )
285
286	IF ( netcdf_data_format > 4 ) THEN
287	ns = 1
288	DO WHILE ( section(ns,s_ind) /= -9999 .AND. ns <= 100 )
289	ns = ns + 1
290	ENDDO
291	ns = ns - 1
292	ALLOCATE( local_2d_sections(1:ns,nys:nyn,nzb:nzt+1) )
293	ALLOCATE( local_2d_sections_l(1:ns,nys:nyn,nzb:nzt+1) )
294	local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
295	ENDIF
296
297	!
298	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
299	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
300	CALL check_open( 103+av*10 )
301	ENDIF
302
303	IF ( data_output_2d_on_each_pe .AND. netcdf_data_format < 5 ) THEN
304	CALL check_open( 23 )
305	ELSE
306	IF ( myid == 0 ) THEN
307	#if defined( __parallel )
308	ALLOCATE( total_2d(0:ny,nzb:nzt+1) )
309	#endif
310	ENDIF
311	ENDIF
312
313	CASE DEFAULT
314	message_string = 'unknown cross-section: ' // TRIM( mode )
315	CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
316
317	END SELECT
318
319	!
320	!-- For parallel netcdf output the time axis must be limited. Return, if this
321	!-- limit is exceeded. This could be the case, if the simulated time exceeds
322	!-- the given end time by the length of the given output interval.
323	IF ( netcdf_data_format > 4 ) THEN
324	IF ( mode == 'xy' .AND. do2d_xy_time_count(av) + 1 > &
325	ntdim_2d_xy(av) ) THEN
326	WRITE ( message_string, * ) 'Output of xy cross-sections is not ', &
327	'given at t=', time_since_reference_point, 's because the', &
328	' maximum number of output time levels is exceeded.'
329	CALL message( 'data_output_2d', 'PA0384', 0, 1, 0, 6, 0 )
330	CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
331	RETURN
332	ENDIF
333	IF ( mode == 'xz' .AND. do2d_xz_time_count(av) + 1 > &
334	ntdim_2d_xz(av) ) THEN
335	WRITE ( message_string, * ) 'Output of xz cross-sections is not ', &
336	'given at t=', time_since_reference_point, 's because the', &
337	' maximum number of output time levels is exceeded.'
338	CALL message( 'data_output_2d', 'PA0385', 0, 1, 0, 6, 0 )
339	CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
340	RETURN
341	ENDIF
342	IF ( mode == 'yz' .AND. do2d_yz_time_count(av) + 1 > &
343	ntdim_2d_yz(av) ) THEN
344	WRITE ( message_string, * ) 'Output of yz cross-sections is not ', &
345	'given at t=', time_since_reference_point, 's because the', &
346	' maximum number of output time levels is exceeded.'
347	CALL message( 'data_output_2d', 'PA0386', 0, 1, 0, 6, 0 )
348	CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
349	RETURN
350	ENDIF
351	ENDIF
352
353	!
354	!-- Allocate a temporary array for resorting (kji -> ijk).
355	ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb:nzt+1) )
356	local_pf = 0.0
357
358	!
359	!-- Loop of all variables to be written.
360	!-- Output dimensions chosen
361	ivar = 1
362	l = MAX( 2, LEN_TRIM( do2d(av,ivar) ) )
363	do2d_mode = do2d(av,ivar)(l-1:l)
364
365	DO WHILE ( do2d(av,ivar)(1:1) /= ' ' )
366
367	IF ( do2d_mode == mode ) THEN
368	!
369	!-- Set flag to steer output of radiation, land-surface, or user-defined
370	!-- quantities
371	found = .FALSE.
372
373	nzb_do = nzb
374	nzt_do = nzt+1
375	!
376	!-- Before each output, set array local_pf to fill value
377	local_pf = fill_value
378	!
379	!-- Set masking flag for topography for not resorted arrays
380	flag_nr = 0
381
382	!
383	!-- Store the array chosen on the temporary array.
384	resorted = .FALSE.
385	SELECT CASE ( TRIM( do2d(av,ivar) ) )
386	CASE ( 'e_xy', 'e_xz', 'e_yz' )
387	IF ( av == 0 ) THEN
388	to_be_resorted => e
389	ELSE
390	IF ( .NOT. ALLOCATED( e_av ) ) THEN
391	ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
392	e_av = REAL( fill_value, KIND = wp )
393	ENDIF
394	to_be_resorted => e_av
395	ENDIF
396	IF ( mode == 'xy' ) level_z = zu
397
398	CASE ( 'thetal_xy', 'thetal_xz', 'thetal_yz' )
399	IF ( av == 0 ) THEN
400	to_be_resorted => pt
401	ELSE
402	IF ( .NOT. ALLOCATED( lpt_av ) ) THEN
403	ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
404	lpt_av = REAL( fill_value, KIND = wp )
405	ENDIF
406	to_be_resorted => lpt_av
407	ENDIF
408	IF ( mode == 'xy' ) level_z = zu
409
410	CASE ( 'lwp*_xy' ) ! 2d-array
411	IF ( av == 0 ) THEN
412	DO i = nxl, nxr
413	DO j = nys, nyn
414	local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
415	dzw(1:nzt+1) )
416	ENDDO
417	ENDDO
418	ELSE
419	IF ( .NOT. ALLOCATED( lwp_av ) ) THEN
420	ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
421	lwp_av = REAL( fill_value, KIND = wp )
422	ENDIF
423	DO i = nxl, nxr
424	DO j = nys, nyn
425	local_pf(i,j,nzb+1) = lwp_av(j,i)
426	ENDDO
427	ENDDO
428	ENDIF
429	resorted = .TRUE.
430	two_d = .TRUE.
431	level_z(nzb+1) = zu(nzb+1)
432
433	CASE ( 'ghf*_xy' ) ! 2d-array
434	IF ( av == 0 ) THEN
435	DO m = 1, surf_lsm_h%ns
436	i = surf_lsm_h%i(m)
437	j = surf_lsm_h%j(m)
438	local_pf(i,j,nzb+1) = surf_lsm_h%ghf(m)
439	ENDDO
440	DO m = 1, surf_usm_h%ns
441	i = surf_usm_h%i(m)
442	j = surf_usm_h%j(m)
443	local_pf(i,j,nzb+1) = surf_usm_h%frac(m,ind_veg_wall) * &
444	surf_usm_h%wghf_eb(m) + &
445	surf_usm_h%frac(m,ind_pav_green) * &
446	surf_usm_h%wghf_eb_green(m) + &
447	surf_usm_h%frac(m,ind_wat_win) * &
448	surf_usm_h%wghf_eb_window(m)
449	ENDDO
450	ELSE
451	IF ( .NOT. ALLOCATED( ghf_av ) ) THEN
452	ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
453	ghf_av = REAL( fill_value, KIND = wp )
454	ENDIF
455	DO i = nxl, nxr
456	DO j = nys, nyn
457	local_pf(i,j,nzb+1) = ghf_av(j,i)
458	ENDDO
459	ENDDO
460	ENDIF
461
462	resorted = .TRUE.
463	two_d = .TRUE.
464	level_z(nzb+1) = zu(nzb+1)
465
466	CASE ( 'ol*_xy' ) ! 2d-array
467	IF ( av == 0 ) THEN
468	DO m = 1, surf_def_h(0)%ns
469	i = surf_def_h(0)%i(m)
470	j = surf_def_h(0)%j(m)
471	local_pf(i,j,nzb+1) = surf_def_h(0)%ol(m)
472	ENDDO
473	DO m = 1, surf_lsm_h%ns
474	i = surf_lsm_h%i(m)
475	j = surf_lsm_h%j(m)
476	local_pf(i,j,nzb+1) = surf_lsm_h%ol(m)
477	ENDDO
478	DO m = 1, surf_usm_h%ns
479	i = surf_usm_h%i(m)
480	j = surf_usm_h%j(m)
481	local_pf(i,j,nzb+1) = surf_usm_h%ol(m)
482	ENDDO
483	ELSE
484	IF ( .NOT. ALLOCATED( ol_av ) ) THEN
485	ALLOCATE( ol_av(nysg:nyng,nxlg:nxrg) )
486	ol_av = REAL( fill_value, KIND = wp )
487	ENDIF
488	DO i = nxl, nxr
489	DO j = nys, nyn
490	local_pf(i,j,nzb+1) = ol_av(j,i)
491	ENDDO
492	ENDDO
493	ENDIF
494	resorted = .TRUE.
495	two_d = .TRUE.
496	level_z(nzb+1) = zu(nzb+1)
497
498	CASE ( 'p_xy', 'p_xz', 'p_yz' )
499	IF ( av == 0 ) THEN
500	IF ( psolver /= 'sor' ) CALL exchange_horiz( p, nbgp )
501	to_be_resorted => p
502	ELSE
503	IF ( .NOT. ALLOCATED( p_av ) ) THEN
504	ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
505	p_av = REAL( fill_value, KIND = wp )
506	ENDIF
507	IF ( psolver /= 'sor' ) CALL exchange_horiz( p_av, nbgp )
508	to_be_resorted => p_av
509	ENDIF
510	IF ( mode == 'xy' ) level_z = zu
511
512	CASE ( 'pc_xy', 'pc_xz', 'pc_yz' ) ! particle concentration
513	IF ( av == 0 ) THEN
514	IF ( time_since_reference_point >= particle_advection_start ) THEN
515	tend = prt_count
516	! CALL exchange_horiz( tend, nbgp )
517	ELSE
518	tend = 0.0_wp
519	ENDIF
520	DO i = nxl, nxr
521	DO j = nys, nyn
522	DO k = nzb, nzt+1
523	local_pf(i,j,k) = tend(k,j,i)
524	ENDDO
525	ENDDO
526	ENDDO
527	resorted = .TRUE.
528	ELSE
529	IF ( .NOT. ALLOCATED( pc_av ) ) THEN
530	ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
531	pc_av = REAL( fill_value, KIND = wp )
532	ENDIF
533	! CALL exchange_horiz( pc_av, nbgp )
534	to_be_resorted => pc_av
535	ENDIF
536
537	CASE ( 'pr_xy', 'pr_xz', 'pr_yz' ) ! mean particle radius (effective radius)
538	IF ( av == 0 ) THEN
539	IF ( time_since_reference_point >= particle_advection_start ) THEN
540	DO i = nxl, nxr
541	DO j = nys, nyn
542	DO k = nzb, nzt+1
543	number_of_particles = prt_count(k,j,i)
544	IF (number_of_particles <= 0) CYCLE
545	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
546	s_r2 = 0.0_wp
547	s_r3 = 0.0_wp
548	DO n = 1, number_of_particles
549	IF ( particles(n)%particle_mask ) THEN
550	s_r2 = s_r2 + particles(n)%radius*2 &
551	particles(n)%weight_factor
552	s_r3 = s_r3 + particles(n)%radius*3 &
553	particles(n)%weight_factor
554	ENDIF
555	ENDDO
556	IF ( s_r2 > 0.0_wp ) THEN
557	mean_r = s_r3 / s_r2
558	ELSE
559	mean_r = 0.0_wp
560	ENDIF
561	tend(k,j,i) = mean_r
562	ENDDO
563	ENDDO
564	ENDDO
565	! CALL exchange_horiz( tend, nbgp )
566	ELSE
567	tend = 0.0_wp
568	ENDIF
569	DO i = nxl, nxr
570	DO j = nys, nyn
571	DO k = nzb, nzt+1
572	local_pf(i,j,k) = tend(k,j,i)
573	ENDDO
574	ENDDO
575	ENDDO
576	resorted = .TRUE.
577	ELSE
578	IF ( .NOT. ALLOCATED( pr_av ) ) THEN
579	ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
580	pr_av = REAL( fill_value, KIND = wp )
581	ENDIF
582	! CALL exchange_horiz( pr_av, nbgp )
583	to_be_resorted => pr_av
584	ENDIF
585
586	CASE ( 'theta_xy', 'theta_xz', 'theta_yz' )
587	IF ( av == 0 ) THEN
588	IF ( .NOT. bulk_cloud_model ) THEN
589	to_be_resorted => pt
590	ELSE
591	DO i = nxl, nxr
592	DO j = nys, nyn
593	DO k = nzb, nzt+1
594	local_pf(i,j,k) = pt(k,j,i) + lv_d_cp * &
595	d_exner(k) * &
596	ql(k,j,i)
597	ENDDO
598	ENDDO
599	ENDDO
600	resorted = .TRUE.
601	ENDIF
602	ELSE
603	IF ( .NOT. ALLOCATED( pt_av ) ) THEN
604	ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
605	pt_av = REAL( fill_value, KIND = wp )
606	ENDIF
607	to_be_resorted => pt_av
608	ENDIF
609	IF ( mode == 'xy' ) level_z = zu
610
611	CASE ( 'q_xy', 'q_xz', 'q_yz' )
612	IF ( av == 0 ) THEN
613	to_be_resorted => q
614	ELSE
615	IF ( .NOT. ALLOCATED( q_av ) ) THEN
616	ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
617	q_av = REAL( fill_value, KIND = wp )
618	ENDIF
619	to_be_resorted => q_av
620	ENDIF
621	IF ( mode == 'xy' ) level_z = zu
622
623	CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
624	IF ( av == 0 ) THEN
625	to_be_resorted => ql
626	ELSE
627	IF ( .NOT. ALLOCATED( ql_av ) ) THEN
628	ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
629	ql_av = REAL( fill_value, KIND = wp )
630	ENDIF
631	to_be_resorted => ql_av
632	ENDIF
633	IF ( mode == 'xy' ) level_z = zu
634
635	CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
636	IF ( av == 0 ) THEN
637	to_be_resorted => ql_c
638	ELSE
639	IF ( .NOT. ALLOCATED( ql_c_av ) ) THEN
640	ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
641	ql_c_av = REAL( fill_value, KIND = wp )
642	ENDIF
643	to_be_resorted => ql_c_av
644	ENDIF
645	IF ( mode == 'xy' ) level_z = zu
646
647	CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
648	IF ( av == 0 ) THEN
649	to_be_resorted => ql_v
650	ELSE
651	IF ( .NOT. ALLOCATED( ql_v_av ) ) THEN
652	ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
653	ql_v_av = REAL( fill_value, KIND = wp )
654	ENDIF
655	to_be_resorted => ql_v_av
656	ENDIF
657	IF ( mode == 'xy' ) level_z = zu
658
659	CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
660	IF ( av == 0 ) THEN
661	IF ( time_since_reference_point >= particle_advection_start ) THEN
662	DO i = nxl, nxr
663	DO j = nys, nyn
664	DO k = nzb, nzt+1
665	number_of_particles = prt_count(k,j,i)
666	IF (number_of_particles <= 0) CYCLE
667	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
668	DO n = 1, number_of_particles
669	IF ( particles(n)%particle_mask ) THEN
670	tend(k,j,i) = tend(k,j,i) + &
671	particles(n)%weight_factor / &
672	prt_count(k,j,i)
673	ENDIF
674	ENDDO
675	ENDDO
676	ENDDO
677	ENDDO
678	! CALL exchange_horiz( tend, nbgp )
679	ELSE
680	tend = 0.0_wp
681	ENDIF
682	DO i = nxl, nxr
683	DO j = nys, nyn
684	DO k = nzb, nzt+1
685	local_pf(i,j,k) = tend(k,j,i)
686	ENDDO
687	ENDDO
688	ENDDO
689	resorted = .TRUE.
690	ELSE
691	IF ( .NOT. ALLOCATED( ql_vp_av ) ) THEN
692	ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
693	ql_vp_av = REAL( fill_value, KIND = wp )
694	ENDIF
695	! CALL exchange_horiz( ql_vp_av, nbgp )
696	to_be_resorted => ql_vp_av
697	ENDIF
698	IF ( mode == 'xy' ) level_z = zu
699
700	CASE ( 'qsws*_xy' ) ! 2d-array
701	IF ( av == 0 ) THEN
702	local_pf(:,:,nzb+1) = REAL( fill_value, KIND = wp )
703	!
704	!-- In case of default surfaces, clean-up flux by density.
705	!-- In case of land-surfaces, convert fluxes into
706	!-- dynamic units
707	DO m = 1, surf_def_h(0)%ns
708	i = surf_def_h(0)%i(m)
709	j = surf_def_h(0)%j(m)
710	k = surf_def_h(0)%k(m)
711	local_pf(i,j,nzb+1) = surf_def_h(0)%qsws(m) * &
712	waterflux_output_conversion(k)
713	ENDDO
714	DO m = 1, surf_lsm_h%ns
715	i = surf_lsm_h%i(m)
716	j = surf_lsm_h%j(m)
717	k = surf_lsm_h%k(m)
718	local_pf(i,j,nzb+1) = surf_lsm_h%qsws(m) * l_v
719	ENDDO
720	DO m = 1, surf_usm_h%ns
721	i = surf_usm_h%i(m)
722	j = surf_usm_h%j(m)
723	k = surf_usm_h%k(m)
724	local_pf(i,j,nzb+1) = surf_usm_h%qsws(m) * l_v
725	ENDDO
726	ELSE
727	IF ( .NOT. ALLOCATED( qsws_av ) ) THEN
728	ALLOCATE( qsws_av(nysg:nyng,nxlg:nxrg) )
729	qsws_av = REAL( fill_value, KIND = wp )
730	ENDIF
731	DO i = nxl, nxr
732	DO j = nys, nyn
733	local_pf(i,j,nzb+1) = qsws_av(j,i)
734	ENDDO
735	ENDDO
736	ENDIF
737	resorted = .TRUE.
738	two_d = .TRUE.
739	level_z(nzb+1) = zu(nzb+1)
740
741	CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
742	IF ( av == 0 ) THEN
743	DO i = nxl, nxr
744	DO j = nys, nyn
745	DO k = nzb, nzt+1
746	local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
747	ENDDO
748	ENDDO
749	ENDDO
750	resorted = .TRUE.
751	ELSE
752	IF ( .NOT. ALLOCATED( qv_av ) ) THEN
753	ALLOCATE( qv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
754	qv_av = REAL( fill_value, KIND = wp )
755	ENDIF
756	to_be_resorted => qv_av
757	ENDIF
758	IF ( mode == 'xy' ) level_z = zu
759
760	CASE ( 'r_a*_xy' ) ! 2d-array
761	IF ( av == 0 ) THEN
762	DO m = 1, surf_lsm_h%ns
763	i = surf_lsm_h%i(m)
764	j = surf_lsm_h%j(m)
765	local_pf(i,j,nzb+1) = surf_lsm_h%r_a(m)
766	ENDDO
767
768	DO m = 1, surf_usm_h%ns
769	i = surf_usm_h%i(m)
770	j = surf_usm_h%j(m)
771	local_pf(i,j,nzb+1) = &
772	( surf_usm_h%frac(m,ind_veg_wall) * &
773	surf_usm_h%r_a(m) + &
774	surf_usm_h%frac(m,ind_pav_green) * &
775	surf_usm_h%r_a_green(m) + &
776	surf_usm_h%frac(m,ind_wat_win) * &
777	surf_usm_h%r_a_window(m) )
778	ENDDO
779	ELSE
780	IF ( .NOT. ALLOCATED( r_a_av ) ) THEN
781	ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
782	r_a_av = REAL( fill_value, KIND = wp )
783	ENDIF
784	DO i = nxl, nxr
785	DO j = nys, nyn
786	local_pf(i,j,nzb+1) = r_a_av(j,i)
787	ENDDO
788	ENDDO
789	ENDIF
790	resorted = .TRUE.
791	two_d = .TRUE.
792	level_z(nzb+1) = zu(nzb+1)
793
794	CASE ( 's_xy', 's_xz', 's_yz' )
795	IF ( av == 0 ) THEN
796	to_be_resorted => s
797	ELSE
798	IF ( .NOT. ALLOCATED( s_av ) ) THEN
799	ALLOCATE( s_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
800	s_av = REAL( fill_value, KIND = wp )
801	ENDIF
802	to_be_resorted => s_av
803	ENDIF
804
805	CASE ( 'shf*_xy' ) ! 2d-array
806	IF ( av == 0 ) THEN
807	!
808	!-- In case of default surfaces, clean-up flux by density.
809	!-- In case of land- and urban-surfaces, convert fluxes into
810	!-- dynamic units.
811	DO m = 1, surf_def_h(0)%ns
812	i = surf_def_h(0)%i(m)
813	j = surf_def_h(0)%j(m)
814	k = surf_def_h(0)%k(m)
815	local_pf(i,j,nzb+1) = surf_def_h(0)%shf(m) * &
816	heatflux_output_conversion(k)
817	ENDDO
818	DO m = 1, surf_lsm_h%ns
819	i = surf_lsm_h%i(m)
820	j = surf_lsm_h%j(m)
821	k = surf_lsm_h%k(m)
822	local_pf(i,j,nzb+1) = surf_lsm_h%shf(m) * c_p
823	ENDDO
824	DO m = 1, surf_usm_h%ns
825	i = surf_usm_h%i(m)
826	j = surf_usm_h%j(m)
827	k = surf_usm_h%k(m)
828	local_pf(i,j,nzb+1) = surf_usm_h%shf(m) * c_p
829	ENDDO
830	ELSE
831	IF ( .NOT. ALLOCATED( shf_av ) ) THEN
832	ALLOCATE( shf_av(nysg:nyng,nxlg:nxrg) )
833	shf_av = REAL( fill_value, KIND = wp )
834	ENDIF
835	DO i = nxl, nxr
836	DO j = nys, nyn
837	local_pf(i,j,nzb+1) = shf_av(j,i)
838	ENDDO
839	ENDDO
840	ENDIF
841	resorted = .TRUE.
842	two_d = .TRUE.
843	level_z(nzb+1) = zu(nzb+1)
844
845	CASE ( 'ssws*_xy' ) ! 2d-array
846	IF ( av == 0 ) THEN
847	DO m = 1, surf_def_h(0)%ns
848	i = surf_def_h(0)%i(m)
849	j = surf_def_h(0)%j(m)
850	local_pf(i,j,nzb+1) = surf_def_h(0)%ssws(m)
851	ENDDO
852	DO m = 1, surf_lsm_h%ns
853	i = surf_lsm_h%i(m)
854	j = surf_lsm_h%j(m)
855	local_pf(i,j,nzb+1) = surf_lsm_h%ssws(m)
856	ENDDO
857	DO m = 1, surf_usm_h%ns
858	i = surf_usm_h%i(m)
859	j = surf_usm_h%j(m)
860	local_pf(i,j,nzb+1) = surf_usm_h%ssws(m)
861	ENDDO
862	ELSE
863	IF ( .NOT. ALLOCATED( ssws_av ) ) THEN
864	ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
865	ssws_av = REAL( fill_value, KIND = wp )
866	ENDIF
867	DO i = nxl, nxr
868	DO j = nys, nyn
869	local_pf(i,j,nzb+1) = ssws_av(j,i)
870	ENDDO
871	ENDDO
872	ENDIF
873	resorted = .TRUE.
874	two_d = .TRUE.
875	level_z(nzb+1) = zu(nzb+1)
876
877	CASE ( 't*_xy' ) ! 2d-array
878	IF ( av == 0 ) THEN
879	DO m = 1, surf_def_h(0)%ns
880	i = surf_def_h(0)%i(m)
881	j = surf_def_h(0)%j(m)
882	local_pf(i,j,nzb+1) = surf_def_h(0)%ts(m)
883	ENDDO
884	DO m = 1, surf_lsm_h%ns
885	i = surf_lsm_h%i(m)
886	j = surf_lsm_h%j(m)
887	local_pf(i,j,nzb+1) = surf_lsm_h%ts(m)
888	ENDDO
889	DO m = 1, surf_usm_h%ns
890	i = surf_usm_h%i(m)
891	j = surf_usm_h%j(m)
892	local_pf(i,j,nzb+1) = surf_usm_h%ts(m)
893	ENDDO
894	ELSE
895	IF ( .NOT. ALLOCATED( ts_av ) ) THEN
896	ALLOCATE( ts_av(nysg:nyng,nxlg:nxrg) )
897	ts_av = REAL( fill_value, KIND = wp )
898	ENDIF
899	DO i = nxl, nxr
900	DO j = nys, nyn
901	local_pf(i,j,nzb+1) = ts_av(j,i)
902	ENDDO
903	ENDDO
904	ENDIF
905	resorted = .TRUE.
906	two_d = .TRUE.
907	level_z(nzb+1) = zu(nzb+1)
908
909	CASE ( 'tsurf*_xy' ) ! 2d-array
910	IF ( av == 0 ) THEN
911	DO m = 1, surf_def_h(0)%ns
912	i = surf_def_h(0)%i(m)
913	j = surf_def_h(0)%j(m)
914	local_pf(i,j,nzb+1) = surf_def_h(0)%pt_surface(m)
915	ENDDO
916
917	DO m = 1, surf_lsm_h%ns
918	i = surf_lsm_h%i(m)
919	j = surf_lsm_h%j(m)
920	local_pf(i,j,nzb+1) = surf_lsm_h%pt_surface(m)
921	ENDDO
922
923	DO m = 1, surf_usm_h%ns
924	i = surf_usm_h%i(m)
925	j = surf_usm_h%j(m)
926	local_pf(i,j,nzb+1) = surf_usm_h%pt_surface(m)
927	ENDDO
928
929	ELSE
930	IF ( .NOT. ALLOCATED( tsurf_av ) ) THEN
931	ALLOCATE( tsurf_av(nysg:nyng,nxlg:nxrg) )
932	tsurf_av = REAL( fill_value, KIND = wp )
933	ENDIF
934	DO i = nxl, nxr
935	DO j = nys, nyn
936	local_pf(i,j,nzb+1) = tsurf_av(j,i)
937	ENDDO
938	ENDDO
939	ENDIF
940	resorted = .TRUE.
941	two_d = .TRUE.
942	level_z(nzb+1) = zu(nzb+1)
943
944	CASE ( 'u_xy', 'u_xz', 'u_yz' )
945	flag_nr = 1
946	IF ( av == 0 ) THEN
947	to_be_resorted => u
948	ELSE
949	IF ( .NOT. ALLOCATED( u_av ) ) THEN
950	ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
951	u_av = REAL( fill_value, KIND = wp )
952	ENDIF
953	to_be_resorted => u_av
954	ENDIF
955	IF ( mode == 'xy' ) level_z = zu
956	!
957	!-- Substitute the values generated by "mirror" boundary condition
958	!-- at the bottom boundary by the real surface values.
959	IF ( do2d(av,ivar) == 'u_xz' .OR. do2d(av,ivar) == 'u_yz' ) THEN
960	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0_wp
961	ENDIF
962
963	CASE ( 'us*_xy' ) ! 2d-array
964	IF ( av == 0 ) THEN
965	DO m = 1, surf_def_h(0)%ns
966	i = surf_def_h(0)%i(m)
967	j = surf_def_h(0)%j(m)
968	local_pf(i,j,nzb+1) = surf_def_h(0)%us(m)
969	ENDDO
970	DO m = 1, surf_lsm_h%ns
971	i = surf_lsm_h%i(m)
972	j = surf_lsm_h%j(m)
973	local_pf(i,j,nzb+1) = surf_lsm_h%us(m)
974	ENDDO
975	DO m = 1, surf_usm_h%ns
976	i = surf_usm_h%i(m)
977	j = surf_usm_h%j(m)
978	local_pf(i,j,nzb+1) = surf_usm_h%us(m)
979	ENDDO
980	ELSE
981	IF ( .NOT. ALLOCATED( us_av ) ) THEN
982	ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
983	us_av = REAL( fill_value, KIND = wp )
984	ENDIF
985	DO i = nxl, nxr
986	DO j = nys, nyn
987	local_pf(i,j,nzb+1) = us_av(j,i)
988	ENDDO
989	ENDDO
990	ENDIF
991	resorted = .TRUE.
992	two_d = .TRUE.
993	level_z(nzb+1) = zu(nzb+1)
994
995	CASE ( 'v_xy', 'v_xz', 'v_yz' )
996	flag_nr = 2
997	IF ( av == 0 ) THEN
998	to_be_resorted => v
999	ELSE
1000	IF ( .NOT. ALLOCATED( v_av ) ) THEN
1001	ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1002	v_av = REAL( fill_value, KIND = wp )
1003	ENDIF
1004	to_be_resorted => v_av
1005	ENDIF
1006	IF ( mode == 'xy' ) level_z = zu
1007	!
1008	!-- Substitute the values generated by "mirror" boundary condition
1009	!-- at the bottom boundary by the real surface values.
1010	IF ( do2d(av,ivar) == 'v_xz' .OR. do2d(av,ivar) == 'v_yz' ) THEN
1011	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0_wp
1012	ENDIF
1013
1014	CASE ( 'thetav_xy', 'thetav_xz', 'thetav_yz' )
1015	IF ( av == 0 ) THEN
1016	to_be_resorted => vpt
1017	ELSE
1018	IF ( .NOT. ALLOCATED( vpt_av ) ) THEN
1019	ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1020	vpt_av = REAL( fill_value, KIND = wp )
1021	ENDIF
1022	to_be_resorted => vpt_av
1023	ENDIF
1024	IF ( mode == 'xy' ) level_z = zu
1025
1026	CASE ( 'w_xy', 'w_xz', 'w_yz' )
1027	flag_nr = 3
1028	IF ( av == 0 ) THEN
1029	to_be_resorted => w
1030	ELSE
1031	IF ( .NOT. ALLOCATED( w_av ) ) THEN
1032	ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1033	w_av = REAL( fill_value, KIND = wp )
1034	ENDIF
1035	to_be_resorted => w_av
1036	ENDIF
1037	IF ( mode == 'xy' ) level_z = zw
1038
1039	CASE ( 'z0*_xy' ) ! 2d-array
1040	IF ( av == 0 ) THEN
1041	DO m = 1, surf_def_h(0)%ns
1042	i = surf_def_h(0)%i(m)
1043	j = surf_def_h(0)%j(m)
1044	local_pf(i,j,nzb+1) = surf_def_h(0)%z0(m)
1045	ENDDO
1046	DO m = 1, surf_lsm_h%ns
1047	i = surf_lsm_h%i(m)
1048	j = surf_lsm_h%j(m)
1049	local_pf(i,j,nzb+1) = surf_lsm_h%z0(m)
1050	ENDDO
1051	DO m = 1, surf_usm_h%ns
1052	i = surf_usm_h%i(m)
1053	j = surf_usm_h%j(m)
1054	local_pf(i,j,nzb+1) = surf_usm_h%z0(m)
1055	ENDDO
1056	ELSE
1057	IF ( .NOT. ALLOCATED( z0_av ) ) THEN
1058	ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
1059	z0_av = REAL( fill_value, KIND = wp )
1060	ENDIF
1061	DO i = nxl, nxr
1062	DO j = nys, nyn
1063	local_pf(i,j,nzb+1) = z0_av(j,i)
1064	ENDDO
1065	ENDDO
1066	ENDIF
1067	resorted = .TRUE.
1068	two_d = .TRUE.
1069	level_z(nzb+1) = zu(nzb+1)
1070
1071	CASE ( 'z0h*_xy' ) ! 2d-array
1072	IF ( av == 0 ) THEN
1073	DO m = 1, surf_def_h(0)%ns
1074	i = surf_def_h(0)%i(m)
1075	j = surf_def_h(0)%j(m)
1076	local_pf(i,j,nzb+1) = surf_def_h(0)%z0h(m)
1077	ENDDO
1078	DO m = 1, surf_lsm_h%ns
1079	i = surf_lsm_h%i(m)
1080	j = surf_lsm_h%j(m)
1081	local_pf(i,j,nzb+1) = surf_lsm_h%z0h(m)
1082	ENDDO
1083	DO m = 1, surf_usm_h%ns
1084	i = surf_usm_h%i(m)
1085	j = surf_usm_h%j(m)
1086	local_pf(i,j,nzb+1) = surf_usm_h%z0h(m)
1087	ENDDO
1088	ELSE
1089	IF ( .NOT. ALLOCATED( z0h_av ) ) THEN
1090	ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
1091	z0h_av = REAL( fill_value, KIND = wp )
1092	ENDIF
1093	DO i = nxl, nxr
1094	DO j = nys, nyn
1095	local_pf(i,j,nzb+1) = z0h_av(j,i)
1096	ENDDO
1097	ENDDO
1098	ENDIF
1099	resorted = .TRUE.
1100	two_d = .TRUE.
1101	level_z(nzb+1) = zu(nzb+1)
1102
1103	CASE ( 'z0q*_xy' ) ! 2d-array
1104	IF ( av == 0 ) THEN
1105	DO m = 1, surf_def_h(0)%ns
1106	i = surf_def_h(0)%i(m)
1107	j = surf_def_h(0)%j(m)
1108	local_pf(i,j,nzb+1) = surf_def_h(0)%z0q(m)
1109	ENDDO
1110	DO m = 1, surf_lsm_h%ns
1111	i = surf_lsm_h%i(m)
1112	j = surf_lsm_h%j(m)
1113	local_pf(i,j,nzb+1) = surf_lsm_h%z0q(m)
1114	ENDDO
1115	DO m = 1, surf_usm_h%ns
1116	i = surf_usm_h%i(m)
1117	j = surf_usm_h%j(m)
1118	local_pf(i,j,nzb+1) = surf_usm_h%z0q(m)
1119	ENDDO
1120	ELSE
1121	IF ( .NOT. ALLOCATED( z0q_av ) ) THEN
1122	ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
1123	z0q_av = REAL( fill_value, KIND = wp )
1124	ENDIF
1125	DO i = nxl, nxr
1126	DO j = nys, nyn
1127	local_pf(i,j,nzb+1) = z0q_av(j,i)
1128	ENDDO
1129	ENDDO
1130	ENDIF
1131	resorted = .TRUE.
1132	two_d = .TRUE.
1133	level_z(nzb+1) = zu(nzb+1)
1134
1135	CASE DEFAULT
1136
1137	!
1138	!-- Quantities of other modules
1139	IF ( .NOT. found ) THEN
1140	CALL module_interface_data_output_2d( &
1141	av, do2d(av,ivar), found, grid, mode, &
1142	local_pf, two_d, nzb_do, nzt_do, &
1143	fill_value &
1144	)
1145	ENDIF
1146
1147	resorted = .TRUE.
1148
1149	IF ( grid == 'zu' ) THEN
1150	IF ( mode == 'xy' ) level_z = zu
1151	ELSEIF ( grid == 'zw' ) THEN
1152	IF ( mode == 'xy' ) level_z = zw
1153	ELSEIF ( grid == 'zu1' ) THEN
1154	IF ( mode == 'xy' ) level_z(nzb+1) = zu(nzb+1)
1155	ELSEIF ( grid == 'zs' ) THEN
1156	IF ( mode == 'xy' ) level_z = zs
1157	ENDIF
1158
1159	IF ( .NOT. found ) THEN
1160	message_string = 'no output provided for: ' // &
1161	TRIM( do2d(av,ivar) )
1162	CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
1163	ENDIF
1164
1165	END SELECT
1166
1167	!
1168	!-- Resort the array to be output, if not done above. Flag topography
1169	!-- grid points with fill values, using the corresponding maksing flag.
1170	IF ( .NOT. resorted ) THEN
1171	DO i = nxl, nxr
1172	DO j = nys, nyn
1173	DO k = nzb_do, nzt_do
1174	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
1175	REAL( fill_value, KIND = wp ), &
1176	BTEST( wall_flags_total_0(k,j,i), &
1177	flag_nr ) )
1178	ENDDO
1179	ENDDO
1180	ENDDO
1181	ENDIF
1182
1183	!
1184	!-- Output of the individual cross-sections, depending on the cross-
1185	!-- section mode chosen.
1186	is = 1
1187	loop1: DO WHILE ( section(is,s_ind) /= -9999 .OR. two_d )
1188
1189	SELECT CASE ( mode )
1190
1191	CASE ( 'xy' )
1192	!
1193	!-- Determine the cross section index
1194	IF ( two_d ) THEN
1195	layer_xy = nzb+1
1196	ELSE
1197	layer_xy = section(is,s_ind)
1198	ENDIF
1199
1200	!
1201	!-- Exit the loop for layers beyond the data output domain
1202	!-- (used for soil model)
1203	IF ( layer_xy > nzt_do ) THEN
1204	EXIT loop1
1205	ENDIF
1206
1207	!
1208	!-- Update the netCDF xy cross section time axis.
1209	!-- In case of parallel output, this is only done by PE0
1210	!-- to increase the performance.
1211	IF ( time_since_reference_point /= do2d_xy_last_time(av) ) THEN
1212	do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
1213	do2d_xy_last_time(av) = time_since_reference_point
1214	IF ( myid == 0 ) THEN
1215	IF ( .NOT. data_output_2d_on_each_pe &
1216	.OR. netcdf_data_format > 4 ) &
1217	THEN
1218	#if defined( __netcdf )
1219	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1220	id_var_time_xy(av), &
1221	(/ time_since_reference_point /), &
1222	start = (/ do2d_xy_time_count(av) /), &
1223	count = (/ 1 /) )
1224	CALL netcdf_handle_error( 'data_output_2d', 53 )
1225	#endif
1226	ENDIF
1227	ENDIF
1228	ENDIF
1229	!
1230	!-- If required, carry out averaging along z
1231	IF ( section(is,s_ind) == -1 .AND. .NOT. two_d ) THEN
1232
1233	local_2d = 0.0_wp
1234	!
1235	!-- Carry out the averaging (all data are on the PE)
1236	DO k = nzb_do, nzt_do
1237	DO j = nys, nyn
1238	DO i = nxl, nxr
1239	local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
1240	ENDDO
1241	ENDDO
1242	ENDDO
1243
1244	local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
1245
1246	ELSE
1247	!
1248	!-- Just store the respective section on the local array
1249	local_2d = local_pf(:,:,layer_xy)
1250
1251	ENDIF
1252
1253	#if defined( __parallel )
1254	IF ( netcdf_data_format > 4 ) THEN
1255	!
1256	!-- Parallel output in netCDF4/HDF5 format.
1257	IF ( two_d ) THEN
1258	iis = 1
1259	ELSE
1260	iis = is
1261	ENDIF
1262
1263	#if defined( __netcdf )
1264	!
1265	!-- For parallel output, all cross sections are first stored
1266	!-- here on a local array and will be written to the output
1267	!-- file afterwards to increase the performance.
1268	DO i = nxl, nxr
1269	DO j = nys, nyn
1270	local_2d_sections(i,j,iis) = local_2d(i,j)
1271	ENDDO
1272	ENDDO
1273	#endif
1274	ELSE
1275
1276	IF ( data_output_2d_on_each_pe ) THEN
1277	!
1278	!-- Output of partial arrays on each PE
1279	#if defined( __netcdf )
1280	IF ( myid == 0 ) THEN
1281	WRITE ( 21 ) time_since_reference_point, &
1282	do2d_xy_time_count(av), av
1283	ENDIF
1284	#endif
1285	DO i = 0, io_blocks-1
1286	IF ( i == io_group ) THEN
1287	WRITE ( 21 ) nxl, nxr, nys, nyn, nys, nyn
1288	WRITE ( 21 ) local_2d
1289	ENDIF
1290	#if defined( __parallel )
1291	CALL MPI_BARRIER( comm2d, ierr )
1292	#endif
1293	ENDDO
1294
1295	ELSE
1296	!
1297	!-- PE0 receives partial arrays from all processors and
1298	!-- then outputs them. Here a barrier has to be set,
1299	!-- because otherwise "-MPI- FATAL: Remote protocol queue
1300	!-- full" may occur.
1301	CALL MPI_BARRIER( comm2d, ierr )
1302
1303	ngp = ( nxr-nxl+1 ) * ( nyn-nys+1 )
1304	IF ( myid == 0 ) THEN
1305	!
1306	!-- Local array can be relocated directly.
1307	total_2d(nxl:nxr,nys:nyn) = local_2d
1308	!
1309	!-- Receive data from all other PEs.
1310	DO n = 1, numprocs-1
1311	!
1312	!-- Receive index limits first, then array.
1313	!-- Index limits are received in arbitrary order from
1314	!-- the PEs.
1315	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1316	MPI_ANY_SOURCE, 0, comm2d, &
1317	status, ierr )
1318	sender = status(MPI_SOURCE)
1319	DEALLOCATE( local_2d )
1320	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1321	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1322	MPI_REAL, sender, 1, comm2d, &
1323	status, ierr )
1324	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1325	ENDDO
1326	!
1327	!-- Relocate the local array for the next loop increment
1328	DEALLOCATE( local_2d )
1329	ALLOCATE( local_2d(nxl:nxr,nys:nyn) )
1330
1331	#if defined( __netcdf )
1332	IF ( two_d ) THEN
1333	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1334	id_var_do2d(av,ivar), &
1335	total_2d(0:nx,0:ny), &
1336	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1337	count = (/ nx+1, ny+1, 1, 1 /) )
1338	ELSE
1339	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1340	id_var_do2d(av,ivar), &
1341	total_2d(0:nx,0:ny), &
1342	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1343	count = (/ nx+1, ny+1, 1, 1 /) )
1344	ENDIF
1345	CALL netcdf_handle_error( 'data_output_2d', 54 )
1346	#endif
1347
1348	ELSE
1349	!
1350	!-- First send the local index limits to PE0
1351	ind(1) = nxl; ind(2) = nxr
1352	ind(3) = nys; ind(4) = nyn
1353	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1354	comm2d, ierr )
1355	!
1356	!-- Send data to PE0
1357	CALL MPI_SEND( local_2d(nxl,nys), ngp, &
1358	MPI_REAL, 0, 1, comm2d, ierr )
1359	ENDIF
1360	!
1361	!-- A barrier has to be set, because otherwise some PEs may
1362	!-- proceed too fast so that PE0 may receive wrong data on
1363	!-- tag 0
1364	CALL MPI_BARRIER( comm2d, ierr )
1365	ENDIF
1366
1367	ENDIF
1368	#else
1369	#if defined( __netcdf )
1370	IF ( two_d ) THEN
1371	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1372	id_var_do2d(av,ivar), &
1373	local_2d(nxl:nxr,nys:nyn), &
1374	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1375	count = (/ nx+1, ny+1, 1, 1 /) )
1376	ELSE
1377	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1378	id_var_do2d(av,ivar), &
1379	local_2d(nxl:nxr,nys:nyn), &
1380	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1381	count = (/ nx+1, ny+1, 1, 1 /) )
1382	ENDIF
1383	CALL netcdf_handle_error( 'data_output_2d', 447 )
1384	#endif
1385	#endif
1386
1387	!
1388	!-- For 2D-arrays (e.g. u*) only one cross-section is available.
1389	!-- Hence exit loop of output levels.
1390	IF ( two_d ) THEN
1391	IF ( netcdf_data_format < 5 ) two_d = .FALSE.
1392	EXIT loop1
1393	ENDIF
1394
1395	CASE ( 'xz' )
1396	!
1397	!-- Update the netCDF xz cross section time axis.
1398	!-- In case of parallel output, this is only done by PE0
1399	!-- to increase the performance.
1400	IF ( time_since_reference_point /= do2d_xz_last_time(av) ) THEN
1401	do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1402	do2d_xz_last_time(av) = time_since_reference_point
1403	IF ( myid == 0 ) THEN
1404	IF ( .NOT. data_output_2d_on_each_pe &
1405	.OR. netcdf_data_format > 4 ) &
1406	THEN
1407	#if defined( __netcdf )
1408	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1409	id_var_time_xz(av), &
1410	(/ time_since_reference_point /), &
1411	start = (/ do2d_xz_time_count(av) /), &
1412	count = (/ 1 /) )
1413	CALL netcdf_handle_error( 'data_output_2d', 56 )
1414	#endif
1415	ENDIF
1416	ENDIF
1417	ENDIF
1418
1419	!
1420	!-- If required, carry out averaging along y
1421	IF ( section(is,s_ind) == -1 ) THEN
1422
1423	ALLOCATE( local_2d_l(nxl:nxr,nzb_do:nzt_do) )
1424	local_2d_l = 0.0_wp
1425	ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1426	!
1427	!-- First local averaging on the PE
1428	DO k = nzb_do, nzt_do
1429	DO j = nys, nyn
1430	DO i = nxl, nxr
1431	local_2d_l(i,k) = local_2d_l(i,k) + &
1432	local_pf(i,j,k)
1433	ENDDO
1434	ENDDO
1435	ENDDO
1436	#if defined( __parallel )
1437	!
1438	!-- Now do the averaging over all PEs along y
1439	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1440	CALL MPI_ALLREDUCE( local_2d_l(nxl,nzb_do), &
1441	local_2d(nxl,nzb_do), ngp, MPI_REAL, &
1442	MPI_SUM, comm1dy, ierr )
1443	#else
1444	local_2d = local_2d_l
1445	#endif
1446	local_2d = local_2d / ( ny + 1.0_wp )
1447
1448	DEALLOCATE( local_2d_l )
1449
1450	ELSE
1451	!
1452	!-- Just store the respective section on the local array
1453	!-- (but only if it is available on this PE!)
1454	IF ( section(is,s_ind) >= nys .AND. section(is,s_ind) <= nyn ) &
1455	THEN
1456	local_2d = local_pf(:,section(is,s_ind),nzb_do:nzt_do)
1457	ENDIF
1458
1459	ENDIF
1460
1461	#if defined( __parallel )
1462	IF ( netcdf_data_format > 4 ) THEN
1463	!
1464	!-- Output in netCDF4/HDF5 format.
1465	!-- Output only on those PEs where the respective cross
1466	!-- sections reside. Cross sections averaged along y are
1467	!-- output on the respective first PE along y (myidy=0).
1468	IF ( ( section(is,s_ind) >= nys .AND. &
1469	section(is,s_ind) <= nyn ) .OR. &
1470	( section(is,s_ind) == -1 .AND. myidy == 0 ) ) THEN
1471	#if defined( __netcdf )
1472	!
1473	!-- For parallel output, all cross sections are first
1474	!-- stored here on a local array and will be written to the
1475	!-- output file afterwards to increase the performance.
1476	DO i = nxl, nxr
1477	DO k = nzb_do, nzt_do
1478	local_2d_sections_l(i,is,k) = local_2d(i,k)
1479	ENDDO
1480	ENDDO
1481	#endif
1482	ENDIF
1483
1484	ELSE
1485
1486	IF ( data_output_2d_on_each_pe ) THEN
1487	!
1488	!-- Output of partial arrays on each PE. If the cross
1489	!-- section does not reside on the PE, output special
1490	!-- index values.
1491	#if defined( __netcdf )
1492	IF ( myid == 0 ) THEN
1493	WRITE ( 22 ) time_since_reference_point, &
1494	do2d_xz_time_count(av), av
1495	ENDIF
1496	#endif
1497	DO i = 0, io_blocks-1
1498	IF ( i == io_group ) THEN
1499	IF ( ( section(is,s_ind) >= nys .AND. &
1500	section(is,s_ind) <= nyn ) .OR. &
1501	( section(is,s_ind) == -1 .AND. &
1502	nys-1 == -1 ) ) &
1503	THEN
1504	WRITE (22) nxl, nxr, nzb_do, nzt_do, nzb, nzt+1
1505	WRITE (22) local_2d
1506	ELSE
1507	WRITE (22) -1, -1, -1, -1, -1, -1
1508	ENDIF
1509	ENDIF
1510	#if defined( __parallel )
1511	CALL MPI_BARRIER( comm2d, ierr )
1512	#endif
1513	ENDDO
1514
1515	ELSE
1516	!
1517	!-- PE0 receives partial arrays from all processors of the
1518	!-- respective cross section and outputs them. Here a
1519	!-- barrier has to be set, because otherwise
1520	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1521	CALL MPI_BARRIER( comm2d, ierr )
1522
1523	ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1524	IF ( myid == 0 ) THEN
1525	!
1526	!-- Local array can be relocated directly.
1527	IF ( ( section(is,s_ind) >= nys .AND. &
1528	section(is,s_ind) <= nyn ) .OR. &
1529	( section(is,s_ind) == -1 .AND. &
1530	nys-1 == -1 ) ) THEN
1531	total_2d(nxl:nxr,nzb_do:nzt_do) = local_2d
1532	ENDIF
1533	!
1534	!-- Receive data from all other PEs.
1535	DO n = 1, numprocs-1
1536	!
1537	!-- Receive index limits first, then array.
1538	!-- Index limits are received in arbitrary order from
1539	!-- the PEs.
1540	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1541	MPI_ANY_SOURCE, 0, comm2d, &
1542	status, ierr )
1543	!
1544	!-- Not all PEs have data for XZ-cross-section.
1545	IF ( ind(1) /= -9999 ) THEN
1546	sender = status(MPI_SOURCE)
1547	DEALLOCATE( local_2d )
1548	ALLOCATE( local_2d(ind(1):ind(2), &
1549	ind(3):ind(4)) )
1550	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1551	MPI_REAL, sender, 1, comm2d, &
1552	status, ierr )
1553	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1554	local_2d
1555	ENDIF
1556	ENDDO
1557	!
1558	!-- Relocate the local array for the next loop increment
1559	DEALLOCATE( local_2d )
1560	ALLOCATE( local_2d(nxl:nxr,nzb_do:nzt_do) )
1561
1562	#if defined( __netcdf )
1563	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1564	id_var_do2d(av,ivar), &
1565	total_2d(0:nx,nzb_do:nzt_do), &
1566	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1567	count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1568	CALL netcdf_handle_error( 'data_output_2d', 58 )
1569	#endif
1570
1571	ELSE
1572	!
1573	!-- If the cross section resides on the PE, send the
1574	!-- local index limits, otherwise send -9999 to PE0.
1575	IF ( ( section(is,s_ind) >= nys .AND. &
1576	section(is,s_ind) <= nyn ) .OR. &
1577	( section(is,s_ind) == -1 .AND. nys-1 == -1 ) ) &
1578	THEN
1579	ind(1) = nxl; ind(2) = nxr
1580	ind(3) = nzb_do; ind(4) = nzt_do
1581	ELSE
1582	ind(1) = -9999; ind(2) = -9999
1583	ind(3) = -9999; ind(4) = -9999
1584	ENDIF
1585	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1586	comm2d, ierr )
1587	!
1588	!-- If applicable, send data to PE0.
1589	IF ( ind(1) /= -9999 ) THEN
1590	CALL MPI_SEND( local_2d(nxl,nzb_do), ngp, &
1591	MPI_REAL, 0, 1, comm2d, ierr )
1592	ENDIF
1593	ENDIF
1594	!
1595	!-- A barrier has to be set, because otherwise some PEs may
1596	!-- proceed too fast so that PE0 may receive wrong data on
1597	!-- tag 0
1598	CALL MPI_BARRIER( comm2d, ierr )
1599	ENDIF
1600
1601	ENDIF
1602	#else
1603	#if defined( __netcdf )
1604	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1605	id_var_do2d(av,ivar), &
1606	local_2d(nxl:nxr,nzb_do:nzt_do), &
1607	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1608	count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1609	CALL netcdf_handle_error( 'data_output_2d', 451 )
1610	#endif
1611	#endif
1612
1613	CASE ( 'yz' )
1614	!
1615	!-- Update the netCDF yz cross section time axis.
1616	!-- In case of parallel output, this is only done by PE0
1617	!-- to increase the performance.
1618	IF ( time_since_reference_point /= do2d_yz_last_time(av) ) THEN
1619	do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1620	do2d_yz_last_time(av) = time_since_reference_point
1621	IF ( myid == 0 ) THEN
1622	IF ( .NOT. data_output_2d_on_each_pe &
1623	.OR. netcdf_data_format > 4 ) &
1624	THEN
1625	#if defined( __netcdf )
1626	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1627	id_var_time_yz(av), &
1628	(/ time_since_reference_point /), &
1629	start = (/ do2d_yz_time_count(av) /), &
1630	count = (/ 1 /) )
1631	CALL netcdf_handle_error( 'data_output_2d', 59 )
1632	#endif
1633	ENDIF
1634	ENDIF
1635	ENDIF
1636
1637	!
1638	!-- If required, carry out averaging along x
1639	IF ( section(is,s_ind) == -1 ) THEN
1640
1641	ALLOCATE( local_2d_l(nys:nyn,nzb_do:nzt_do) )
1642	local_2d_l = 0.0_wp
1643	ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1644	!
1645	!-- First local averaging on the PE
1646	DO k = nzb_do, nzt_do
1647	DO j = nys, nyn
1648	DO i = nxl, nxr
1649	local_2d_l(j,k) = local_2d_l(j,k) + &
1650	local_pf(i,j,k)
1651	ENDDO
1652	ENDDO
1653	ENDDO
1654	#if defined( __parallel )
1655	!
1656	!-- Now do the averaging over all PEs along x
1657	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1658	CALL MPI_ALLREDUCE( local_2d_l(nys,nzb_do), &
1659	local_2d(nys,nzb_do), ngp, MPI_REAL, &
1660	MPI_SUM, comm1dx, ierr )
1661	#else
1662	local_2d = local_2d_l
1663	#endif
1664	local_2d = local_2d / ( nx + 1.0_wp )
1665
1666	DEALLOCATE( local_2d_l )
1667
1668	ELSE
1669	!
1670	!-- Just store the respective section on the local array
1671	!-- (but only if it is available on this PE!)
1672	IF ( section(is,s_ind) >= nxl .AND. section(is,s_ind) <= nxr ) &
1673	THEN
1674	local_2d = local_pf(section(is,s_ind),:,nzb_do:nzt_do)
1675	ENDIF
1676
1677	ENDIF
1678
1679	#if defined( __parallel )
1680	IF ( netcdf_data_format > 4 ) THEN
1681	!
1682	!-- Output in netCDF4/HDF5 format.
1683	!-- Output only on those PEs where the respective cross
1684	!-- sections reside. Cross sections averaged along x are
1685	!-- output on the respective first PE along x (myidx=0).
1686	IF ( ( section(is,s_ind) >= nxl .AND. &
1687	section(is,s_ind) <= nxr ) .OR. &
1688	( section(is,s_ind) == -1 .AND. myidx == 0 ) ) THEN
1689	#if defined( __netcdf )
1690	!
1691	!-- For parallel output, all cross sections are first
1692	!-- stored here on a local array and will be written to the
1693	!-- output file afterwards to increase the performance.
1694	DO j = nys, nyn
1695	DO k = nzb_do, nzt_do
1696	local_2d_sections_l(is,j,k) = local_2d(j,k)
1697	ENDDO
1698	ENDDO
1699	#endif
1700	ENDIF
1701
1702	ELSE
1703
1704	IF ( data_output_2d_on_each_pe ) THEN
1705	!
1706	!-- Output of partial arrays on each PE. If the cross
1707	!-- section does not reside on the PE, output special
1708	!-- index values.
1709	#if defined( __netcdf )
1710	IF ( myid == 0 ) THEN
1711	WRITE ( 23 ) time_since_reference_point, &
1712	do2d_yz_time_count(av), av
1713	ENDIF
1714	#endif
1715	DO i = 0, io_blocks-1
1716	IF ( i == io_group ) THEN
1717	IF ( ( section(is,s_ind) >= nxl .AND. &
1718	section(is,s_ind) <= nxr ) .OR. &
1719	( section(is,s_ind) == -1 .AND. &
1720	nxl-1 == -1 ) ) &
1721	THEN
1722	WRITE (23) nys, nyn, nzb_do, nzt_do, nzb, nzt+1
1723	WRITE (23) local_2d
1724	ELSE
1725	WRITE (23) -1, -1, -1, -1, -1, -1
1726	ENDIF
1727	ENDIF
1728	#if defined( __parallel )
1729	CALL MPI_BARRIER( comm2d, ierr )
1730	#endif
1731	ENDDO
1732
1733	ELSE
1734	!
1735	!-- PE0 receives partial arrays from all processors of the
1736	!-- respective cross section and outputs them. Here a
1737	!-- barrier has to be set, because otherwise
1738	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1739	CALL MPI_BARRIER( comm2d, ierr )
1740
1741	ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1742	IF ( myid == 0 ) THEN
1743	!
1744	!-- Local array can be relocated directly.
1745	IF ( ( section(is,s_ind) >= nxl .AND. &
1746	section(is,s_ind) <= nxr ) .OR. &
1747	( section(is,s_ind) == -1 .AND. nxl-1 == -1 ) ) &
1748	THEN
1749	total_2d(nys:nyn,nzb_do:nzt_do) = local_2d
1750	ENDIF
1751	!
1752	!-- Receive data from all other PEs.
1753	DO n = 1, numprocs-1
1754	!
1755	!-- Receive index limits first, then array.
1756	!-- Index limits are received in arbitrary order from
1757	!-- the PEs.
1758	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1759	MPI_ANY_SOURCE, 0, comm2d, &
1760	status, ierr )
1761	!
1762	!-- Not all PEs have data for YZ-cross-section.
1763	IF ( ind(1) /= -9999 ) THEN
1764	sender = status(MPI_SOURCE)
1765	DEALLOCATE( local_2d )
1766	ALLOCATE( local_2d(ind(1):ind(2), &
1767	ind(3):ind(4)) )
1768	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1769	MPI_REAL, sender, 1, comm2d, &
1770	status, ierr )
1771	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1772	local_2d
1773	ENDIF
1774	ENDDO
1775	!
1776	!-- Relocate the local array for the next loop increment
1777	DEALLOCATE( local_2d )
1778	ALLOCATE( local_2d(nys:nyn,nzb_do:nzt_do) )
1779
1780	#if defined( __netcdf )
1781	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1782	id_var_do2d(av,ivar), &
1783	total_2d(0:ny,nzb_do:nzt_do), &
1784	start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1785	count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
1786	CALL netcdf_handle_error( 'data_output_2d', 61 )
1787	#endif
1788
1789	ELSE
1790	!
1791	!-- If the cross section resides on the PE, send the
1792	!-- local index limits, otherwise send -9999 to PE0.
1793	IF ( ( section(is,s_ind) >= nxl .AND. &
1794	section(is,s_ind) <= nxr ) .OR. &
1795	( section(is,s_ind) == -1 .AND. nxl-1 == -1 ) ) &
1796	THEN
1797	ind(1) = nys; ind(2) = nyn
1798	ind(3) = nzb_do; ind(4) = nzt_do
1799	ELSE
1800	ind(1) = -9999; ind(2) = -9999
1801	ind(3) = -9999; ind(4) = -9999
1802	ENDIF
1803	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1804	comm2d, ierr )
1805	!
1806	!-- If applicable, send data to PE0.
1807	IF ( ind(1) /= -9999 ) THEN
1808	CALL MPI_SEND( local_2d(nys,nzb_do), ngp, &
1809	MPI_REAL, 0, 1, comm2d, ierr )
1810	ENDIF
1811	ENDIF
1812	!
1813	!-- A barrier has to be set, because otherwise some PEs may
1814	!-- proceed too fast so that PE0 may receive wrong data on
1815	!-- tag 0
1816	CALL MPI_BARRIER( comm2d, ierr )
1817	ENDIF
1818
1819	ENDIF
1820	#else
1821	#if defined( __netcdf )
1822	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1823	id_var_do2d(av,ivar), &
1824	local_2d(nys:nyn,nzb_do:nzt_do), &
1825	start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1826	count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
1827	CALL netcdf_handle_error( 'data_output_2d', 452 )
1828	#endif
1829	#endif
1830
1831	END SELECT
1832
1833	is = is + 1
1834	ENDDO loop1
1835
1836	!
1837	!-- For parallel output, all data were collected before on a local array
1838	!-- and are written now to the netcdf file. This must be done to increase
1839	!-- the performance of the parallel output.
1840	#if defined( __netcdf )
1841	IF ( netcdf_data_format > 4 ) THEN
1842
1843	SELECT CASE ( mode )
1844
1845	CASE ( 'xy' )
1846	IF ( two_d ) THEN
1847	nis = 1
1848	two_d = .FALSE.
1849	ELSE
1850	nis = ns
1851	ENDIF
1852	!
1853	!-- Do not output redundant ghost point data except for the
1854	!-- boundaries of the total domain.
1855	! IF ( nxr == nx .AND. nyn /= ny ) THEN
1856	! nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1857	! id_var_do2d(av,ivar), &
1858	! local_2d_sections(nxl:nxr+1, &
1859	! nys:nyn,1:nis), &
1860	! start = (/ nxl+1, nys+1, 1, &
1861	! do2d_xy_time_count(av) /), &
1862	! count = (/ nxr-nxl+2, &
1863	! nyn-nys+1, nis, 1 &
1864	! /) )
1865	! ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
1866	! nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1867	! id_var_do2d(av,ivar), &
1868	! local_2d_sections(nxl:nxr, &
1869	! nys:nyn+1,1:nis), &
1870	! start = (/ nxl+1, nys+1, 1, &
1871	! do2d_xy_time_count(av) /), &
1872	! count = (/ nxr-nxl+1, &
1873	! nyn-nys+2, nis, 1 &
1874	! /) )
1875	! ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
1876	! nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1877	! id_var_do2d(av,ivar), &
1878	! local_2d_sections(nxl:nxr+1, &
1879	! nys:nyn+1,1:nis), &
1880	! start = (/ nxl+1, nys+1, 1, &
1881	! do2d_xy_time_count(av) /), &
1882	! count = (/ nxr-nxl+2, &
1883	! nyn-nys+2, nis, 1 &
1884	! /) )
1885	! ELSE
1886	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1887	id_var_do2d(av,ivar), &
1888	local_2d_sections(nxl:nxr, &
1889	nys:nyn,1:nis), &
1890	start = (/ nxl+1, nys+1, 1, &
1891	do2d_xy_time_count(av) /), &
1892	count = (/ nxr-nxl+1, &
1893	nyn-nys+1, nis, 1 &
1894	/) )
1895	! ENDIF
1896
1897	CALL netcdf_handle_error( 'data_output_2d', 55 )
1898
1899	CASE ( 'xz' )
1900	!
1901	!-- First, all PEs get the information of all cross-sections.
1902	!-- Then the data are written to the output file by all PEs
1903	!-- while NF90_COLLECTIVE is set in subroutine
1904	!-- define_netcdf_header. Although redundant information are
1905	!-- written to the output file in that case, the performance
1906	!-- is significantly better compared to the case where only
1907	!-- the first row of PEs in x-direction (myidx = 0) is given
1908	!-- the output while NF90_INDEPENDENT is set.
1909	IF ( npey /= 1 ) THEN
1910
1911	#if defined( __parallel )
1912	!
1913	!-- Distribute data over all PEs along y
1914	ngp = ( nxr-nxl+1 ) * ( nzt_do-nzb_do+1 ) * ns
1915	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1916	CALL MPI_ALLREDUCE( local_2d_sections_l(nxl,1,nzb_do), &
1917	local_2d_sections(nxl,1,nzb_do), &
1918	ngp, MPI_REAL, MPI_SUM, comm1dy, &
1919	ierr )
1920	#else
1921	local_2d_sections = local_2d_sections_l
1922	#endif
1923	ENDIF
1924	!
1925	!-- Do not output redundant ghost point data except for the
1926	!-- boundaries of the total domain.
1927	! IF ( nxr == nx ) THEN
1928	! nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1929	! id_var_do2d(av,ivar), &
1930	! local_2d_sections(nxl:nxr+1,1:ns, &
1931	! nzb_do:nzt_do), &
1932	! start = (/ nxl+1, 1, 1, &
1933	! do2d_xz_time_count(av) /), &
1934	! count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1, &
1935	! 1 /) )
1936	! ELSE
1937	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1938	id_var_do2d(av,ivar), &
1939	local_2d_sections(nxl:nxr,1:ns, &
1940	nzb_do:nzt_do), &
1941	start = (/ nxl+1, 1, 1, &
1942	do2d_xz_time_count(av) /), &
1943	count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1, &
1944	1 /) )
1945	! ENDIF
1946
1947	CALL netcdf_handle_error( 'data_output_2d', 57 )
1948
1949	CASE ( 'yz' )
1950	!
1951	!-- First, all PEs get the information of all cross-sections.
1952	!-- Then the data are written to the output file by all PEs
1953	!-- while NF90_COLLECTIVE is set in subroutine
1954	!-- define_netcdf_header. Although redundant information are
1955	!-- written to the output file in that case, the performance
1956	!-- is significantly better compared to the case where only
1957	!-- the first row of PEs in y-direction (myidy = 0) is given
1958	!-- the output while NF90_INDEPENDENT is set.
1959	IF ( npex /= 1 ) THEN
1960
1961	#if defined( __parallel )
1962	!
1963	!-- Distribute data over all PEs along x
1964	ngp = ( nyn-nys+1 ) * ( nzt-nzb + 2 ) * ns
1965	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1966	CALL MPI_ALLREDUCE( local_2d_sections_l(1,nys,nzb_do), &
1967	local_2d_sections(1,nys,nzb_do), &
1968	ngp, MPI_REAL, MPI_SUM, comm1dx, &
1969	ierr )
1970	#else
1971	local_2d_sections = local_2d_sections_l
1972	#endif
1973	ENDIF
1974	!
1975	!-- Do not output redundant ghost point data except for the
1976	!-- boundaries of the total domain.
1977	! IF ( nyn == ny ) THEN
1978	! nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1979	! id_var_do2d(av,ivar), &
1980	! local_2d_sections(1:ns, &
1981	! nys:nyn+1,nzb_do:nzt_do), &
1982	! start = (/ 1, nys+1, 1, &
1983	! do2d_yz_time_count(av) /), &
1984	! count = (/ ns, nyn-nys+2, &
1985	! nzt_do-nzb_do+1, 1 /) )
1986	! ELSE
1987	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1988	id_var_do2d(av,ivar), &
1989	local_2d_sections(1:ns,nys:nyn, &
1990	nzb_do:nzt_do), &
1991	start = (/ 1, nys+1, 1, &
1992	do2d_yz_time_count(av) /), &
1993	count = (/ ns, nyn-nys+1, &
1994	nzt_do-nzb_do+1, 1 /) )
1995	! ENDIF
1996
1997	CALL netcdf_handle_error( 'data_output_2d', 60 )
1998
1999	CASE DEFAULT
2000	message_string = 'unknown cross-section: ' // TRIM( mode )
2001	CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
2002
2003	END SELECT
2004
2005	ENDIF
2006	#endif
2007	ENDIF
2008
2009	ivar = ivar + 1
2010	l = MAX( 2, LEN_TRIM( do2d(av,ivar) ) )
2011	do2d_mode = do2d(av,ivar)(l-1:l)
2012
2013	ENDDO
2014
2015	!
2016	!-- Deallocate temporary arrays.
2017	IF ( ALLOCATED( level_z ) ) DEALLOCATE( level_z )
2018	IF ( netcdf_data_format > 4 ) THEN
2019	DEALLOCATE( local_pf, local_2d, local_2d_sections )
2020	IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
2021	ENDIF
2022	#if defined( __parallel )
2023	IF ( .NOT. data_output_2d_on_each_pe .AND. myid == 0 ) THEN
2024	DEALLOCATE( total_2d )
2025	ENDIF
2026	#endif
2027
2028	!
2029	!-- Close plot output file.
2030	file_id = 20 + s_ind
2031
2032	IF ( data_output_2d_on_each_pe ) THEN
2033	DO i = 0, io_blocks-1
2034	IF ( i == io_group ) THEN
2035	CALL close_file( file_id )
2036	ENDIF
2037	#if defined( __parallel )
2038	CALL MPI_BARRIER( comm2d, ierr )
2039	#endif
2040	ENDDO
2041	ELSE
2042	IF ( myid == 0 ) CALL close_file( file_id )
2043	ENDIF
2044
2045	CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
2046
2047	IF ( debug_output_timestep ) CALL debug_message( 'data_output_2d', 'end' )
2048
2049
2050	END SUBROUTINE data_output_2d

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |