Home

Context Navigation

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

Last change on this file since 2813 was 2805, checked in by suehring, 6 years ago
Bugfix in re-mapping surface-element data in case of restarts; bugfix in initialization of water-surface roughness; bugfix - uninitialized resistance at urban-type surfaces
Property svn:keywords set to `Id`
File size: 89.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |