Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |