Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |