Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |