Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |