Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |