Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |