Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |