Home

Context Navigation

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

Last change on this file since 4517 was 4514, checked in by suehring, 5 years ago
Bugfix in plant-canopy model for output of averaged transpiration rate after a restart; Revise check for output for plant heating rate and rename error message number; Surface-data output: enable output of mixing ratio and passive scalar concentration at the surface; Surface-data input: Add possibility to prescribe surface sensible and latent heat fluxes via static input file
Property svn:keywords set to `Id`
File size: 91.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |