Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |