Home

Context Navigation

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

Last change on this file since 2799 was 2798, checked in by suehring, 6 years ago
Bugfix initialization of %pt_surface array; Output of surface temperature also for default-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 2798 2018-02-09 17:16:39Z raasch $
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
916	ELSE
917	DO i = nxl, nxr
918	DO j = nys, nyn
919	local_pf(i,j,nzb+1) = r_a_av(j,i)
920	ENDDO
921	ENDDO
922	ENDIF
923	resorted = .TRUE.
924	two_d = .TRUE.
925	level_z(nzb+1) = zu(nzb+1)
926
927	CASE ( 'rho_ocean_xy', 'rho_ocean_xz', 'rho_ocean_yz' )
928	IF ( av == 0 ) THEN
929	to_be_resorted => rho_ocean
930	ELSE
931	to_be_resorted => rho_ocean_av
932	ENDIF
933
934	CASE ( 's_xy', 's_xz', 's_yz' )
935	IF ( av == 0 ) THEN
936	to_be_resorted => s
937	ELSE
938	to_be_resorted => s_av
939	ENDIF
940
941	CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
942	IF ( av == 0 ) THEN
943	to_be_resorted => sa
944	ELSE
945	to_be_resorted => sa_av
946	ENDIF
947
948	CASE ( 'shf*_xy' ) ! 2d-array
949	IF ( av == 0 ) THEN
950	!
951	!-- In case of default surfaces, clean-up flux by density.
952	!-- In case of land- and urban-surfaces, convert fluxes into
953	!-- dynamic units.
954	DO m = 1, surf_def_h(0)%ns
955	i = surf_def_h(0)%i(m)
956	j = surf_def_h(0)%j(m)
957	k = surf_def_h(0)%k(m)
958	local_pf(i,j,nzb+1) = surf_def_h(0)%shf(m) * &
959	heatflux_output_conversion(k)
960	ENDDO
961	DO m = 1, surf_lsm_h%ns
962	i = surf_lsm_h%i(m)
963	j = surf_lsm_h%j(m)
964	k = surf_lsm_h%k(m)
965	local_pf(i,j,nzb+1) = surf_lsm_h%shf(m) * cp
966	ENDDO
967	DO m = 1, surf_usm_h%ns
968	i = surf_usm_h%i(m)
969	j = surf_usm_h%j(m)
970	k = surf_usm_h%k(m)
971	local_pf(i,j,nzb+1) = surf_usm_h%shf(m) * cp
972	ENDDO
973	ELSE
974	DO i = nxl, nxr
975	DO j = nys, nyn
976	local_pf(i,j,nzb+1) = shf_av(j,i)
977	ENDDO
978	ENDDO
979	ENDIF
980	resorted = .TRUE.
981	two_d = .TRUE.
982	level_z(nzb+1) = zu(nzb+1)
983
984	CASE ( 'ssws*_xy' ) ! 2d-array
985	IF ( av == 0 ) THEN
986	DO m = 1, surf_def_h(0)%ns
987	i = surf_def_h(0)%i(m)
988	j = surf_def_h(0)%j(m)
989	local_pf(i,j,nzb+1) = surf_def_h(0)%ssws(m)
990	ENDDO
991	DO m = 1, surf_lsm_h%ns
992	i = surf_lsm_h%i(m)
993	j = surf_lsm_h%j(m)
994	local_pf(i,j,nzb+1) = surf_lsm_h%ssws(m)
995	ENDDO
996	DO m = 1, surf_usm_h%ns
997	i = surf_usm_h%i(m)
998	j = surf_usm_h%j(m)
999	local_pf(i,j,nzb+1) = surf_usm_h%ssws(m)
1000	ENDDO
1001	ELSE
1002	DO i = nxl, nxr
1003	DO j = nys, nyn
1004	local_pf(i,j,nzb+1) = ssws_av(j,i)
1005	ENDDO
1006	ENDDO
1007	ENDIF
1008	resorted = .TRUE.
1009	two_d = .TRUE.
1010	level_z(nzb+1) = zu(nzb+1)
1011
1012	CASE ( 't*_xy' ) ! 2d-array
1013	IF ( av == 0 ) THEN
1014	DO m = 1, surf_def_h(0)%ns
1015	i = surf_def_h(0)%i(m)
1016	j = surf_def_h(0)%j(m)
1017	local_pf(i,j,nzb+1) = surf_def_h(0)%ts(m)
1018	ENDDO
1019	DO m = 1, surf_lsm_h%ns
1020	i = surf_lsm_h%i(m)
1021	j = surf_lsm_h%j(m)
1022	local_pf(i,j,nzb+1) = surf_lsm_h%ts(m)
1023	ENDDO
1024	DO m = 1, surf_usm_h%ns
1025	i = surf_usm_h%i(m)
1026	j = surf_usm_h%j(m)
1027	local_pf(i,j,nzb+1) = surf_usm_h%ts(m)
1028	ENDDO
1029	ELSE
1030	DO i = nxl, nxr
1031	DO j = nys, nyn
1032	local_pf(i,j,nzb+1) = ts_av(j,i)
1033	ENDDO
1034	ENDDO
1035	ENDIF
1036	resorted = .TRUE.
1037	two_d = .TRUE.
1038	level_z(nzb+1) = zu(nzb+1)
1039
1040	CASE ( 'tsurf*_xy' ) ! 2d-array
1041	IF ( av == 0 ) THEN
1042	DO m = 1, surf_def_h(0)%ns
1043	i = surf_def_h(0)%i(m)
1044	j = surf_def_h(0)%j(m)
1045	local_pf(i,j,nzb+1) = surf_def_h(0)%pt_surface(m)
1046	ENDDO
1047
1048	DO m = 1, surf_lsm_h%ns
1049	i = surf_lsm_h%i(m)
1050	j = surf_lsm_h%j(m)
1051	local_pf(i,j,nzb+1) = surf_lsm_h%pt_surface(m)
1052	ENDDO
1053
1054	DO m = 1, surf_usm_h%ns
1055	i = surf_usm_h%i(m)
1056	j = surf_usm_h%j(m)
1057	local_pf(i,j,nzb+1) = surf_usm_h%pt_surface(m)
1058	ENDDO
1059
1060	ELSE
1061	DO i = nxl, nxr
1062	DO j = nys, nyn
1063	local_pf(i,j,nzb+1) = tsurf_av(j,i)
1064	ENDDO
1065	ENDDO
1066	ENDIF
1067	resorted = .TRUE.
1068	two_d = .TRUE.
1069	level_z(nzb+1) = zu(nzb+1)
1070
1071	CASE ( 'u_xy', 'u_xz', 'u_yz' )
1072	flag_nr = 1
1073	IF ( av == 0 ) THEN
1074	to_be_resorted => u
1075	ELSE
1076	to_be_resorted => u_av
1077	ENDIF
1078	IF ( mode == 'xy' ) level_z = zu
1079	!
1080	!-- Substitute the values generated by "mirror" boundary condition
1081	!-- at the bottom boundary by the real surface values.
1082	IF ( do2d(av,if) == 'u_xz' .OR. do2d(av,if) == 'u_yz' ) THEN
1083	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0_wp
1084	ENDIF
1085
1086	CASE ( 'u*_xy' ) ! 2d-array
1087	IF ( av == 0 ) THEN
1088	DO m = 1, surf_def_h(0)%ns
1089	i = surf_def_h(0)%i(m)
1090	j = surf_def_h(0)%j(m)
1091	local_pf(i,j,nzb+1) = surf_def_h(0)%us(m)
1092	ENDDO
1093	DO m = 1, surf_lsm_h%ns
1094	i = surf_lsm_h%i(m)
1095	j = surf_lsm_h%j(m)
1096	local_pf(i,j,nzb+1) = surf_lsm_h%us(m)
1097	ENDDO
1098	DO m = 1, surf_usm_h%ns
1099	i = surf_usm_h%i(m)
1100	j = surf_usm_h%j(m)
1101	local_pf(i,j,nzb+1) = surf_usm_h%us(m)
1102	ENDDO
1103	ELSE
1104	DO i = nxl, nxr
1105	DO j = nys, nyn
1106	local_pf(i,j,nzb+1) = us_av(j,i)
1107	ENDDO
1108	ENDDO
1109	ENDIF
1110	resorted = .TRUE.
1111	two_d = .TRUE.
1112	level_z(nzb+1) = zu(nzb+1)
1113
1114	CASE ( 'v_xy', 'v_xz', 'v_yz' )
1115	flag_nr = 2
1116	IF ( av == 0 ) THEN
1117	to_be_resorted => v
1118	ELSE
1119	to_be_resorted => v_av
1120	ENDIF
1121	IF ( mode == 'xy' ) level_z = zu
1122	!
1123	!-- Substitute the values generated by "mirror" boundary condition
1124	!-- at the bottom boundary by the real surface values.
1125	IF ( do2d(av,if) == 'v_xz' .OR. do2d(av,if) == 'v_yz' ) THEN
1126	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0_wp
1127	ENDIF
1128
1129	CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
1130	IF ( av == 0 ) THEN
1131	to_be_resorted => vpt
1132	ELSE
1133	to_be_resorted => vpt_av
1134	ENDIF
1135	IF ( mode == 'xy' ) level_z = zu
1136
1137	CASE ( 'w_xy', 'w_xz', 'w_yz' )
1138	flag_nr = 3
1139	IF ( av == 0 ) THEN
1140	to_be_resorted => w
1141	ELSE
1142	to_be_resorted => w_av
1143	ENDIF
1144	IF ( mode == 'xy' ) level_z = zw
1145
1146	CASE ( 'z0*_xy' ) ! 2d-array
1147	IF ( av == 0 ) THEN
1148	DO m = 1, surf_def_h(0)%ns
1149	i = surf_def_h(0)%i(m)
1150	j = surf_def_h(0)%j(m)
1151	local_pf(i,j,nzb+1) = surf_def_h(0)%z0(m)
1152	ENDDO
1153	DO m = 1, surf_lsm_h%ns
1154	i = surf_lsm_h%i(m)
1155	j = surf_lsm_h%j(m)
1156	local_pf(i,j,nzb+1) = surf_lsm_h%z0(m)
1157	ENDDO
1158	DO m = 1, surf_usm_h%ns
1159	i = surf_usm_h%i(m)
1160	j = surf_usm_h%j(m)
1161	local_pf(i,j,nzb+1) = surf_usm_h%z0(m)
1162	ENDDO
1163	ELSE
1164	DO i = nxl, nxr
1165	DO j = nys, nyn
1166	local_pf(i,j,nzb+1) = z0_av(j,i)
1167	ENDDO
1168	ENDDO
1169	ENDIF
1170	resorted = .TRUE.
1171	two_d = .TRUE.
1172	level_z(nzb+1) = zu(nzb+1)
1173
1174	CASE ( 'z0h*_xy' ) ! 2d-array
1175	IF ( av == 0 ) THEN
1176	DO m = 1, surf_def_h(0)%ns
1177	i = surf_def_h(0)%i(m)
1178	j = surf_def_h(0)%j(m)
1179	local_pf(i,j,nzb+1) = surf_def_h(0)%z0h(m)
1180	ENDDO
1181	DO m = 1, surf_lsm_h%ns
1182	i = surf_lsm_h%i(m)
1183	j = surf_lsm_h%j(m)
1184	local_pf(i,j,nzb+1) = surf_lsm_h%z0h(m)
1185	ENDDO
1186	DO m = 1, surf_usm_h%ns
1187	i = surf_usm_h%i(m)
1188	j = surf_usm_h%j(m)
1189	local_pf(i,j,nzb+1) = surf_usm_h%z0h(m)
1190	ENDDO
1191	ELSE
1192	DO i = nxl, nxr
1193	DO j = nys, nyn
1194	local_pf(i,j,nzb+1) = z0h_av(j,i)
1195	ENDDO
1196	ENDDO
1197	ENDIF
1198	resorted = .TRUE.
1199	two_d = .TRUE.
1200	level_z(nzb+1) = zu(nzb+1)
1201
1202	CASE ( 'z0q*_xy' ) ! 2d-array
1203	IF ( av == 0 ) THEN
1204	DO m = 1, surf_def_h(0)%ns
1205	i = surf_def_h(0)%i(m)
1206	j = surf_def_h(0)%j(m)
1207	local_pf(i,j,nzb+1) = surf_def_h(0)%z0q(m)
1208	ENDDO
1209	DO m = 1, surf_lsm_h%ns
1210	i = surf_lsm_h%i(m)
1211	j = surf_lsm_h%j(m)
1212	local_pf(i,j,nzb+1) = surf_lsm_h%z0q(m)
1213	ENDDO
1214	DO m = 1, surf_usm_h%ns
1215	i = surf_usm_h%i(m)
1216	j = surf_usm_h%j(m)
1217	local_pf(i,j,nzb+1) = surf_usm_h%z0q(m)
1218	ENDDO
1219	ELSE
1220	DO i = nxl, nxr
1221	DO j = nys, nyn
1222	local_pf(i,j,nzb+1) = z0q_av(j,i)
1223	ENDDO
1224	ENDDO
1225	ENDIF
1226	resorted = .TRUE.
1227	two_d = .TRUE.
1228	level_z(nzb+1) = zu(nzb+1)
1229
1230	CASE DEFAULT
1231
1232	!
1233	!-- Land surface model quantity
1234	IF ( land_surface ) THEN
1235	CALL lsm_data_output_2d( av, do2d(av,if), found, grid, mode,&
1236	local_pf, two_d, nzb_do, nzt_do )
1237	ENDIF
1238
1239	!
1240	!-- Turbulence closure variables
1241	IF ( .NOT. found ) THEN
1242	CALL tcm_data_output_2d( av, do2d(av,if), found, grid, mode,&
1243	local_pf, two_d, nzb_do, nzt_do )
1244	ENDIF
1245
1246	!
1247	!-- Radiation quantity
1248	IF ( .NOT. found .AND. radiation ) THEN
1249	CALL radiation_data_output_2d( av, do2d(av,if), found, grid,&
1250	mode, local_pf, two_d )
1251	ENDIF
1252
1253	!
1254	!-- UV exposure model quantity
1255	IF ( uv_exposure ) THEN
1256	CALL uvem_data_output_2d( av, do2d(av,if), found, grid, mode,&
1257	local_pf, two_d, nzb_do, nzt_do )
1258	ENDIF
1259
1260	!
1261	!-- User defined quantity
1262	IF ( .NOT. found ) THEN
1263	CALL user_data_output_2d( av, do2d(av,if), found, grid, &
1264	local_pf, two_d, nzb_do, nzt_do )
1265	ENDIF
1266
1267	resorted = .TRUE.
1268
1269	IF ( grid == 'zu' ) THEN
1270	IF ( mode == 'xy' ) level_z = zu
1271	ELSEIF ( grid == 'zw' ) THEN
1272	IF ( mode == 'xy' ) level_z = zw
1273	ELSEIF ( grid == 'zu1' ) THEN
1274	IF ( mode == 'xy' ) level_z(nzb+1) = zu(nzb+1)
1275	ELSEIF ( grid == 'zs' ) THEN
1276	IF ( mode == 'xy' ) level_z = zs
1277	ENDIF
1278
1279	IF ( .NOT. found ) THEN
1280	message_string = 'no output provided for: ' // &
1281	TRIM( do2d(av,if) )
1282	CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
1283	ENDIF
1284
1285	END SELECT
1286
1287	!
1288	!-- Resort the array to be output, if not done above. Flag topography
1289	!-- grid points with fill values, using the corresponding maksing flag.
1290	IF ( .NOT. resorted ) THEN
1291	DO i = nxl, nxr
1292	DO j = nys, nyn
1293	DO k = nzb_do, nzt_do
1294	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
1295	REAL( fill_value, KIND = wp ), &
1296	BTEST( wall_flags_0(k,j,i), &
1297	flag_nr ) )
1298	ENDDO
1299	ENDDO
1300	ENDDO
1301	ENDIF
1302
1303	!
1304	!-- Output of the individual cross-sections, depending on the cross-
1305	!-- section mode chosen.
1306	is = 1
1307	loop1: DO WHILE ( section(is,s_ind) /= -9999 .OR. two_d )
1308
1309	SELECT CASE ( mode )
1310
1311	CASE ( 'xy' )
1312	!
1313	!-- Determine the cross section index
1314	IF ( two_d ) THEN
1315	layer_xy = nzb+1
1316	ELSE
1317	layer_xy = section(is,s_ind)
1318	ENDIF
1319
1320	!
1321	!-- Exit the loop for layers beyond the data output domain
1322	!-- (used for soil model)
1323	IF ( layer_xy > nzt_do ) THEN
1324	EXIT loop1
1325	ENDIF
1326
1327	!
1328	!-- Update the netCDF xy cross section time axis.
1329	!-- In case of parallel output, this is only done by PE0
1330	!-- to increase the performance.
1331	IF ( simulated_time /= do2d_xy_last_time(av) ) THEN
1332	do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
1333	do2d_xy_last_time(av) = simulated_time
1334	IF ( myid == 0 ) THEN
1335	IF ( .NOT. data_output_2d_on_each_pe &
1336	.OR. netcdf_data_format > 4 ) &
1337	THEN
1338	#if defined( __netcdf )
1339	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1340	id_var_time_xy(av), &
1341	(/ time_since_reference_point /), &
1342	start = (/ do2d_xy_time_count(av) /), &
1343	count = (/ 1 /) )
1344	CALL netcdf_handle_error( 'data_output_2d', 53 )
1345	#endif
1346	ENDIF
1347	ENDIF
1348	ENDIF
1349	!
1350	!-- If required, carry out averaging along z
1351	IF ( section(is,s_ind) == -1 .AND. .NOT. two_d ) THEN
1352
1353	local_2d = 0.0_wp
1354	!
1355	!-- Carry out the averaging (all data are on the PE)
1356	DO k = nzb_do, nzt_do
1357	DO j = nys, nyn
1358	DO i = nxl, nxr
1359	local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
1360	ENDDO
1361	ENDDO
1362	ENDDO
1363
1364	local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
1365
1366	ELSE
1367	!
1368	!-- Just store the respective section on the local array
1369	local_2d = local_pf(:,:,layer_xy)
1370
1371	ENDIF
1372
1373	#if defined( __parallel )
1374	IF ( netcdf_data_format > 4 ) THEN
1375	!
1376	!-- Parallel output in netCDF4/HDF5 format.
1377	IF ( two_d ) THEN
1378	iis = 1
1379	ELSE
1380	iis = is
1381	ENDIF
1382
1383	#if defined( __netcdf )
1384	!
1385	!-- For parallel output, all cross sections are first stored
1386	!-- here on a local array and will be written to the output
1387	!-- file afterwards to increase the performance.
1388	DO i = nxl, nxr
1389	DO j = nys, nyn
1390	local_2d_sections(i,j,iis) = local_2d(i,j)
1391	ENDDO
1392	ENDDO
1393	#endif
1394	ELSE
1395
1396	IF ( data_output_2d_on_each_pe ) THEN
1397	!
1398	!-- Output of partial arrays on each PE
1399	#if defined( __netcdf )
1400	IF ( myid == 0 ) THEN
1401	WRITE ( 21 ) time_since_reference_point, &
1402	do2d_xy_time_count(av), av
1403	ENDIF
1404	#endif
1405	DO i = 0, io_blocks-1
1406	IF ( i == io_group ) THEN
1407	WRITE ( 21 ) nxl, nxr, nys, nyn, nys, nyn
1408	WRITE ( 21 ) local_2d
1409	ENDIF
1410	#if defined( __parallel )
1411	CALL MPI_BARRIER( comm2d, ierr )
1412	#endif
1413	ENDDO
1414
1415	ELSE
1416	!
1417	!-- PE0 receives partial arrays from all processors and
1418	!-- then outputs them. Here a barrier has to be set,
1419	!-- because otherwise "-MPI- FATAL: Remote protocol queue
1420	!-- full" may occur.
1421	CALL MPI_BARRIER( comm2d, ierr )
1422
1423	ngp = ( nxr-nxl+1 ) * ( nyn-nys+1 )
1424	IF ( myid == 0 ) THEN
1425	!
1426	!-- Local array can be relocated directly.
1427	total_2d(nxl:nxr,nys:nyn) = local_2d
1428	!
1429	!-- Receive data from all other PEs.
1430	DO n = 1, numprocs-1
1431	!
1432	!-- Receive index limits first, then array.
1433	!-- Index limits are received in arbitrary order from
1434	!-- the PEs.
1435	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1436	MPI_ANY_SOURCE, 0, comm2d, &
1437	status, ierr )
1438	sender = status(MPI_SOURCE)
1439	DEALLOCATE( local_2d )
1440	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1441	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1442	MPI_REAL, sender, 1, comm2d, &
1443	status, ierr )
1444	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1445	ENDDO
1446	!
1447	!-- Relocate the local array for the next loop increment
1448	DEALLOCATE( local_2d )
1449	ALLOCATE( local_2d(nxl:nxr,nys:nyn) )
1450
1451	#if defined( __netcdf )
1452	IF ( two_d ) THEN
1453	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1454	id_var_do2d(av,if), &
1455	total_2d(0:nx,0:ny), &
1456	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1457	count = (/ nx+1, ny+1, 1, 1 /) )
1458	ELSE
1459	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1460	id_var_do2d(av,if), &
1461	total_2d(0:nx,0:ny), &
1462	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1463	count = (/ nx+1, ny+1, 1, 1 /) )
1464	ENDIF
1465	CALL netcdf_handle_error( 'data_output_2d', 54 )
1466	#endif
1467
1468	ELSE
1469	!
1470	!-- First send the local index limits to PE0
1471	ind(1) = nxl; ind(2) = nxr
1472	ind(3) = nys; ind(4) = nyn
1473	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1474	comm2d, ierr )
1475	!
1476	!-- Send data to PE0
1477	CALL MPI_SEND( local_2d(nxl,nys), ngp, &
1478	MPI_REAL, 0, 1, comm2d, ierr )
1479	ENDIF
1480	!
1481	!-- A barrier has to be set, because otherwise some PEs may
1482	!-- proceed too fast so that PE0 may receive wrong data on
1483	!-- tag 0
1484	CALL MPI_BARRIER( comm2d, ierr )
1485	ENDIF
1486
1487	ENDIF
1488	#else
1489	#if defined( __netcdf )
1490	IF ( two_d ) THEN
1491	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1492	id_var_do2d(av,if), &
1493	local_2d(nxl:nxr,nys:nyn), &
1494	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1495	count = (/ nx+1, ny+1, 1, 1 /) )
1496	ELSE
1497	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1498	id_var_do2d(av,if), &
1499	local_2d(nxl:nxr,nys:nyn), &
1500	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1501	count = (/ nx+1, ny+1, 1, 1 /) )
1502	ENDIF
1503	CALL netcdf_handle_error( 'data_output_2d', 447 )
1504	#endif
1505	#endif
1506
1507	!
1508	!-- For 2D-arrays (e.g. u*) only one cross-section is available.
1509	!-- Hence exit loop of output levels.
1510	IF ( two_d ) THEN
1511	IF ( netcdf_data_format < 5 ) two_d = .FALSE.
1512	EXIT loop1
1513	ENDIF
1514
1515	CASE ( 'xz' )
1516	!
1517	!-- Update the netCDF xz cross section time axis.
1518	!-- In case of parallel output, this is only done by PE0
1519	!-- to increase the performance.
1520	IF ( simulated_time /= do2d_xz_last_time(av) ) THEN
1521	do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1522	do2d_xz_last_time(av) = simulated_time
1523	IF ( myid == 0 ) THEN
1524	IF ( .NOT. data_output_2d_on_each_pe &
1525	.OR. netcdf_data_format > 4 ) &
1526	THEN
1527	#if defined( __netcdf )
1528	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1529	id_var_time_xz(av), &
1530	(/ time_since_reference_point /), &
1531	start = (/ do2d_xz_time_count(av) /), &
1532	count = (/ 1 /) )
1533	CALL netcdf_handle_error( 'data_output_2d', 56 )
1534	#endif
1535	ENDIF
1536	ENDIF
1537	ENDIF
1538
1539	!
1540	!-- If required, carry out averaging along y
1541	IF ( section(is,s_ind) == -1 ) THEN
1542
1543	ALLOCATE( local_2d_l(nxl:nxr,nzb_do:nzt_do) )
1544	local_2d_l = 0.0_wp
1545	ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1546	!
1547	!-- First local averaging on the PE
1548	DO k = nzb_do, nzt_do
1549	DO j = nys, nyn
1550	DO i = nxl, nxr
1551	local_2d_l(i,k) = local_2d_l(i,k) + &
1552	local_pf(i,j,k)
1553	ENDDO
1554	ENDDO
1555	ENDDO
1556	#if defined( __parallel )
1557	!
1558	!-- Now do the averaging over all PEs along y
1559	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1560	CALL MPI_ALLREDUCE( local_2d_l(nxl,nzb_do), &
1561	local_2d(nxl,nzb_do), ngp, MPI_REAL, &
1562	MPI_SUM, comm1dy, ierr )
1563	#else
1564	local_2d = local_2d_l
1565	#endif
1566	local_2d = local_2d / ( ny + 1.0_wp )
1567
1568	DEALLOCATE( local_2d_l )
1569
1570	ELSE
1571	!
1572	!-- Just store the respective section on the local array
1573	!-- (but only if it is available on this PE!)
1574	IF ( section(is,s_ind) >= nys .AND. section(is,s_ind) <= nyn ) &
1575	THEN
1576	local_2d = local_pf(:,section(is,s_ind),nzb_do:nzt_do)
1577	ENDIF
1578
1579	ENDIF
1580
1581	#if defined( __parallel )
1582	IF ( netcdf_data_format > 4 ) THEN
1583	!
1584	!-- Output in netCDF4/HDF5 format.
1585	!-- Output only on those PEs where the respective cross
1586	!-- sections reside. Cross sections averaged along y are
1587	!-- output on the respective first PE along y (myidy=0).
1588	IF ( ( section(is,s_ind) >= nys .AND. &
1589	section(is,s_ind) <= nyn ) .OR. &
1590	( section(is,s_ind) == -1 .AND. myidy == 0 ) ) THEN
1591	#if defined( __netcdf )
1592	!
1593	!-- For parallel output, all cross sections are first
1594	!-- stored here on a local array and will be written to the
1595	!-- output file afterwards to increase the performance.
1596	DO i = nxl, nxr
1597	DO k = nzb_do, nzt_do
1598	local_2d_sections_l(i,is,k) = local_2d(i,k)
1599	ENDDO
1600	ENDDO
1601	#endif
1602	ENDIF
1603
1604	ELSE
1605
1606	IF ( data_output_2d_on_each_pe ) THEN
1607	!
1608	!-- Output of partial arrays on each PE. If the cross
1609	!-- section does not reside on the PE, output special
1610	!-- index values.
1611	#if defined( __netcdf )
1612	IF ( myid == 0 ) THEN
1613	WRITE ( 22 ) time_since_reference_point, &
1614	do2d_xz_time_count(av), av
1615	ENDIF
1616	#endif
1617	DO i = 0, io_blocks-1
1618	IF ( i == io_group ) THEN
1619	IF ( ( section(is,s_ind) >= nys .AND. &
1620	section(is,s_ind) <= nyn ) .OR. &
1621	( section(is,s_ind) == -1 .AND. &
1622	nys-1 == -1 ) ) &
1623	THEN
1624	WRITE (22) nxl, nxr, nzb_do, nzt_do, nzb, nzt+1
1625	WRITE (22) local_2d
1626	ELSE
1627	WRITE (22) -1, -1, -1, -1, -1, -1
1628	ENDIF
1629	ENDIF
1630	#if defined( __parallel )
1631	CALL MPI_BARRIER( comm2d, ierr )
1632	#endif
1633	ENDDO
1634
1635	ELSE
1636	!
1637	!-- PE0 receives partial arrays from all processors of the
1638	!-- respective cross section and outputs them. Here a
1639	!-- barrier has to be set, because otherwise
1640	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1641	CALL MPI_BARRIER( comm2d, ierr )
1642
1643	ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1644	IF ( myid == 0 ) THEN
1645	!
1646	!-- Local array can be relocated directly.
1647	IF ( ( section(is,s_ind) >= nys .AND. &
1648	section(is,s_ind) <= nyn ) .OR. &
1649	( section(is,s_ind) == -1 .AND. &
1650	nys-1 == -1 ) ) THEN
1651	total_2d(nxl:nxr,nzb_do:nzt_do) = local_2d
1652	ENDIF
1653	!
1654	!-- Receive data from all other PEs.
1655	DO n = 1, numprocs-1
1656	!
1657	!-- Receive index limits first, then array.
1658	!-- Index limits are received in arbitrary order from
1659	!-- the PEs.
1660	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1661	MPI_ANY_SOURCE, 0, comm2d, &
1662	status, ierr )
1663	!
1664	!-- Not all PEs have data for XZ-cross-section.
1665	IF ( ind(1) /= -9999 ) THEN
1666	sender = status(MPI_SOURCE)
1667	DEALLOCATE( local_2d )
1668	ALLOCATE( local_2d(ind(1):ind(2), &
1669	ind(3):ind(4)) )
1670	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1671	MPI_REAL, sender, 1, comm2d, &
1672	status, ierr )
1673	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1674	local_2d
1675	ENDIF
1676	ENDDO
1677	!
1678	!-- Relocate the local array for the next loop increment
1679	DEALLOCATE( local_2d )
1680	ALLOCATE( local_2d(nxl:nxr,nzb_do:nzt_do) )
1681
1682	#if defined( __netcdf )
1683	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1684	id_var_do2d(av,if), &
1685	total_2d(0:nx,nzb_do:nzt_do), &
1686	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1687	count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1688	CALL netcdf_handle_error( 'data_output_2d', 58 )
1689	#endif
1690
1691	ELSE
1692	!
1693	!-- If the cross section resides on the PE, send the
1694	!-- local index limits, otherwise send -9999 to PE0.
1695	IF ( ( section(is,s_ind) >= nys .AND. &
1696	section(is,s_ind) <= nyn ) .OR. &
1697	( section(is,s_ind) == -1 .AND. nys-1 == -1 ) ) &
1698	THEN
1699	ind(1) = nxl; ind(2) = nxr
1700	ind(3) = nzb_do; ind(4) = nzt_do
1701	ELSE
1702	ind(1) = -9999; ind(2) = -9999
1703	ind(3) = -9999; ind(4) = -9999
1704	ENDIF
1705	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1706	comm2d, ierr )
1707	!
1708	!-- If applicable, send data to PE0.
1709	IF ( ind(1) /= -9999 ) THEN
1710	CALL MPI_SEND( local_2d(nxl,nzb_do), ngp, &
1711	MPI_REAL, 0, 1, comm2d, ierr )
1712	ENDIF
1713	ENDIF
1714	!
1715	!-- A barrier has to be set, because otherwise some PEs may
1716	!-- proceed too fast so that PE0 may receive wrong data on
1717	!-- tag 0
1718	CALL MPI_BARRIER( comm2d, ierr )
1719	ENDIF
1720
1721	ENDIF
1722	#else
1723	#if defined( __netcdf )
1724	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1725	id_var_do2d(av,if), &
1726	local_2d(nxl:nxr,nzb_do:nzt_do), &
1727	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1728	count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1729	CALL netcdf_handle_error( 'data_output_2d', 451 )
1730	#endif
1731	#endif
1732
1733	CASE ( 'yz' )
1734	!
1735	!-- Update the netCDF yz cross section time axis.
1736	!-- In case of parallel output, this is only done by PE0
1737	!-- to increase the performance.
1738	IF ( simulated_time /= do2d_yz_last_time(av) ) THEN
1739	do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1740	do2d_yz_last_time(av) = simulated_time
1741	IF ( myid == 0 ) THEN
1742	IF ( .NOT. data_output_2d_on_each_pe &
1743	.OR. netcdf_data_format > 4 ) &
1744	THEN
1745	#if defined( __netcdf )
1746	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1747	id_var_time_yz(av), &
1748	(/ time_since_reference_point /), &
1749	start = (/ do2d_yz_time_count(av) /), &
1750	count = (/ 1 /) )
1751	CALL netcdf_handle_error( 'data_output_2d', 59 )
1752	#endif
1753	ENDIF
1754	ENDIF
1755	ENDIF
1756
1757	!
1758	!-- If required, carry out averaging along x
1759	IF ( section(is,s_ind) == -1 ) THEN
1760
1761	ALLOCATE( local_2d_l(nys:nyn,nzb_do:nzt_do) )
1762	local_2d_l = 0.0_wp
1763	ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1764	!
1765	!-- First local averaging on the PE
1766	DO k = nzb_do, nzt_do
1767	DO j = nys, nyn
1768	DO i = nxl, nxr
1769	local_2d_l(j,k) = local_2d_l(j,k) + &
1770	local_pf(i,j,k)
1771	ENDDO
1772	ENDDO
1773	ENDDO
1774	#if defined( __parallel )
1775	!
1776	!-- Now do the averaging over all PEs along x
1777	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1778	CALL MPI_ALLREDUCE( local_2d_l(nys,nzb_do), &
1779	local_2d(nys,nzb_do), ngp, MPI_REAL, &
1780	MPI_SUM, comm1dx, ierr )
1781	#else
1782	local_2d = local_2d_l
1783	#endif
1784	local_2d = local_2d / ( nx + 1.0_wp )
1785
1786	DEALLOCATE( local_2d_l )
1787
1788	ELSE
1789	!
1790	!-- Just store the respective section on the local array
1791	!-- (but only if it is available on this PE!)
1792	IF ( section(is,s_ind) >= nxl .AND. section(is,s_ind) <= nxr ) &
1793	THEN
1794	local_2d = local_pf(section(is,s_ind),:,nzb_do:nzt_do)
1795	ENDIF
1796
1797	ENDIF
1798
1799	#if defined( __parallel )
1800	IF ( netcdf_data_format > 4 ) THEN
1801	!
1802	!-- Output in netCDF4/HDF5 format.
1803	!-- Output only on those PEs where the respective cross
1804	!-- sections reside. Cross sections averaged along x are
1805	!-- output on the respective first PE along x (myidx=0).
1806	IF ( ( section(is,s_ind) >= nxl .AND. &
1807	section(is,s_ind) <= nxr ) .OR. &
1808	( section(is,s_ind) == -1 .AND. myidx == 0 ) ) THEN
1809	#if defined( __netcdf )
1810	!
1811	!-- For parallel output, all cross sections are first
1812	!-- stored here on a local array and will be written to the
1813	!-- output file afterwards to increase the performance.
1814	DO j = nys, nyn
1815	DO k = nzb_do, nzt_do
1816	local_2d_sections_l(is,j,k) = local_2d(j,k)
1817	ENDDO
1818	ENDDO
1819	#endif
1820	ENDIF
1821
1822	ELSE
1823
1824	IF ( data_output_2d_on_each_pe ) THEN
1825	!
1826	!-- Output of partial arrays on each PE. If the cross
1827	!-- section does not reside on the PE, output special
1828	!-- index values.
1829	#if defined( __netcdf )
1830	IF ( myid == 0 ) THEN
1831	WRITE ( 23 ) time_since_reference_point, &
1832	do2d_yz_time_count(av), av
1833	ENDIF
1834	#endif
1835	DO i = 0, io_blocks-1
1836	IF ( i == io_group ) THEN
1837	IF ( ( section(is,s_ind) >= nxl .AND. &
1838	section(is,s_ind) <= nxr ) .OR. &
1839	( section(is,s_ind) == -1 .AND. &
1840	nxl-1 == -1 ) ) &
1841	THEN
1842	WRITE (23) nys, nyn, nzb_do, nzt_do, nzb, nzt+1
1843	WRITE (23) local_2d
1844	ELSE
1845	WRITE (23) -1, -1, -1, -1, -1, -1
1846	ENDIF
1847	ENDIF
1848	#if defined( __parallel )
1849	CALL MPI_BARRIER( comm2d, ierr )
1850	#endif
1851	ENDDO
1852
1853	ELSE
1854	!
1855	!-- PE0 receives partial arrays from all processors of the
1856	!-- respective cross section and outputs them. Here a
1857	!-- barrier has to be set, because otherwise
1858	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1859	CALL MPI_BARRIER( comm2d, ierr )
1860
1861	ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1862	IF ( myid == 0 ) THEN
1863	!
1864	!-- Local array can be relocated directly.
1865	IF ( ( section(is,s_ind) >= nxl .AND. &
1866	section(is,s_ind) <= nxr ) .OR. &
1867	( section(is,s_ind) == -1 .AND. nxl-1 == -1 ) ) &
1868	THEN
1869	total_2d(nys:nyn,nzb_do:nzt_do) = local_2d
1870	ENDIF
1871	!
1872	!-- Receive data from all other PEs.
1873	DO n = 1, numprocs-1
1874	!
1875	!-- Receive index limits first, then array.
1876	!-- Index limits are received in arbitrary order from
1877	!-- the PEs.
1878	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1879	MPI_ANY_SOURCE, 0, comm2d, &
1880	status, ierr )
1881	!
1882	!-- Not all PEs have data for YZ-cross-section.
1883	IF ( ind(1) /= -9999 ) THEN
1884	sender = status(MPI_SOURCE)
1885	DEALLOCATE( local_2d )
1886	ALLOCATE( local_2d(ind(1):ind(2), &
1887	ind(3):ind(4)) )
1888	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1889	MPI_REAL, sender, 1, comm2d, &
1890	status, ierr )
1891	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1892	local_2d
1893	ENDIF
1894	ENDDO
1895	!
1896	!-- Relocate the local array for the next loop increment
1897	DEALLOCATE( local_2d )
1898	ALLOCATE( local_2d(nys:nyn,nzb_do:nzt_do) )
1899
1900	#if defined( __netcdf )
1901	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1902	id_var_do2d(av,if), &
1903	total_2d(0:ny,nzb_do:nzt_do), &
1904	start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1905	count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
1906	CALL netcdf_handle_error( 'data_output_2d', 61 )
1907	#endif
1908
1909	ELSE
1910	!
1911	!-- If the cross section resides on the PE, send the
1912	!-- local index limits, otherwise send -9999 to PE0.
1913	IF ( ( section(is,s_ind) >= nxl .AND. &
1914	section(is,s_ind) <= nxr ) .OR. &
1915	( section(is,s_ind) == -1 .AND. nxl-1 == -1 ) ) &
1916	THEN
1917	ind(1) = nys; ind(2) = nyn
1918	ind(3) = nzb_do; ind(4) = nzt_do
1919	ELSE
1920	ind(1) = -9999; ind(2) = -9999
1921	ind(3) = -9999; ind(4) = -9999
1922	ENDIF
1923	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1924	comm2d, ierr )
1925	!
1926	!-- If applicable, send data to PE0.
1927	IF ( ind(1) /= -9999 ) THEN
1928	CALL MPI_SEND( local_2d(nys,nzb_do), ngp, &
1929	MPI_REAL, 0, 1, comm2d, ierr )
1930	ENDIF
1931	ENDIF
1932	!
1933	!-- A barrier has to be set, because otherwise some PEs may
1934	!-- proceed too fast so that PE0 may receive wrong data on
1935	!-- tag 0
1936	CALL MPI_BARRIER( comm2d, ierr )
1937	ENDIF
1938
1939	ENDIF
1940	#else
1941	#if defined( __netcdf )
1942	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1943	id_var_do2d(av,if), &
1944	local_2d(nys:nyn,nzb_do:nzt_do), &
1945	start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1946	count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
1947	CALL netcdf_handle_error( 'data_output_2d', 452 )
1948	#endif
1949	#endif
1950
1951	END SELECT
1952
1953	is = is + 1
1954	ENDDO loop1
1955
1956	!
1957	!-- For parallel output, all data were collected before on a local array
1958	!-- and are written now to the netcdf file. This must be done to increase
1959	!-- the performance of the parallel output.
1960	#if defined( __netcdf )
1961	IF ( netcdf_data_format > 4 ) THEN
1962
1963	SELECT CASE ( mode )
1964
1965	CASE ( 'xy' )
1966	IF ( two_d ) THEN
1967	nis = 1
1968	two_d = .FALSE.
1969	ELSE
1970	nis = ns
1971	ENDIF
1972	!
1973	!-- Do not output redundant ghost point data except for the
1974	!-- boundaries of the total domain.
1975	! IF ( nxr == nx .AND. nyn /= ny ) THEN
1976	! nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1977	! id_var_do2d(av,if), &
1978	! local_2d_sections(nxl:nxr+1, &
1979	! nys:nyn,1:nis), &
1980	! start = (/ nxl+1, nys+1, 1, &
1981	! do2d_xy_time_count(av) /), &
1982	! count = (/ nxr-nxl+2, &
1983	! nyn-nys+1, nis, 1 &
1984	! /) )
1985	! ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
1986	! nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1987	! id_var_do2d(av,if), &
1988	! local_2d_sections(nxl:nxr, &
1989	! nys:nyn+1,1:nis), &
1990	! start = (/ nxl+1, nys+1, 1, &
1991	! do2d_xy_time_count(av) /), &
1992	! count = (/ nxr-nxl+1, &
1993	! nyn-nys+2, nis, 1 &
1994	! /) )
1995	! ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
1996	! nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1997	! id_var_do2d(av,if), &
1998	! local_2d_sections(nxl:nxr+1, &
1999	! nys:nyn+1,1:nis), &
2000	! start = (/ nxl+1, nys+1, 1, &
2001	! do2d_xy_time_count(av) /), &
2002	! count = (/ nxr-nxl+2, &
2003	! nyn-nys+2, nis, 1 &
2004	! /) )
2005	! ELSE
2006	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
2007	id_var_do2d(av,if), &
2008	local_2d_sections(nxl:nxr, &
2009	nys:nyn,1:nis), &
2010	start = (/ nxl+1, nys+1, 1, &
2011	do2d_xy_time_count(av) /), &
2012	count = (/ nxr-nxl+1, &
2013	nyn-nys+1, nis, 1 &
2014	/) )
2015	! ENDIF
2016
2017	CALL netcdf_handle_error( 'data_output_2d', 55 )
2018
2019	CASE ( 'xz' )
2020	!
2021	!-- First, all PEs get the information of all cross-sections.
2022	!-- Then the data are written to the output file by all PEs
2023	!-- while NF90_COLLECTIVE is set in subroutine
2024	!-- define_netcdf_header. Although redundant information are
2025	!-- written to the output file in that case, the performance
2026	!-- is significantly better compared to the case where only
2027	!-- the first row of PEs in x-direction (myidx = 0) is given
2028	!-- the output while NF90_INDEPENDENT is set.
2029	IF ( npey /= 1 ) THEN
2030
2031	#if defined( __parallel )
2032	!
2033	!-- Distribute data over all PEs along y
2034	ngp = ( nxr-nxl+1 ) * ( nzt_do-nzb_do+1 ) * ns
2035	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2036	CALL MPI_ALLREDUCE( local_2d_sections_l(nxl,1,nzb_do), &
2037	local_2d_sections(nxl,1,nzb_do), &
2038	ngp, MPI_REAL, MPI_SUM, comm1dy, &
2039	ierr )
2040	#else
2041	local_2d_sections = local_2d_sections_l
2042	#endif
2043	ENDIF
2044	!
2045	!-- Do not output redundant ghost point data except for the
2046	!-- boundaries of the total domain.
2047	! IF ( nxr == nx ) THEN
2048	! nc_stat = NF90_PUT_VAR( id_set_xz(av), &
2049	! id_var_do2d(av,if), &
2050	! local_2d_sections(nxl:nxr+1,1:ns, &
2051	! nzb_do:nzt_do), &
2052	! start = (/ nxl+1, 1, 1, &
2053	! do2d_xz_time_count(av) /), &
2054	! count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1, &
2055	! 1 /) )
2056	! ELSE
2057	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
2058	id_var_do2d(av,if), &
2059	local_2d_sections(nxl:nxr,1:ns, &
2060	nzb_do:nzt_do), &
2061	start = (/ nxl+1, 1, 1, &
2062	do2d_xz_time_count(av) /), &
2063	count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1, &
2064	1 /) )
2065	! ENDIF
2066
2067	CALL netcdf_handle_error( 'data_output_2d', 57 )
2068
2069	CASE ( 'yz' )
2070	!
2071	!-- First, all PEs get the information of all cross-sections.
2072	!-- Then the data are written to the output file by all PEs
2073	!-- while NF90_COLLECTIVE is set in subroutine
2074	!-- define_netcdf_header. Although redundant information are
2075	!-- written to the output file in that case, the performance
2076	!-- is significantly better compared to the case where only
2077	!-- the first row of PEs in y-direction (myidy = 0) is given
2078	!-- the output while NF90_INDEPENDENT is set.
2079	IF ( npex /= 1 ) THEN
2080
2081	#if defined( __parallel )
2082	!
2083	!-- Distribute data over all PEs along x
2084	ngp = ( nyn-nys+1 ) * ( nzt-nzb + 2 ) * ns
2085	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2086	CALL MPI_ALLREDUCE( local_2d_sections_l(1,nys,nzb_do), &
2087	local_2d_sections(1,nys,nzb_do), &
2088	ngp, MPI_REAL, MPI_SUM, comm1dx, &
2089	ierr )
2090	#else
2091	local_2d_sections = local_2d_sections_l
2092	#endif
2093	ENDIF
2094	!
2095	!-- Do not output redundant ghost point data except for the
2096	!-- boundaries of the total domain.
2097	! IF ( nyn == ny ) THEN
2098	! nc_stat = NF90_PUT_VAR( id_set_yz(av), &
2099	! id_var_do2d(av,if), &
2100	! local_2d_sections(1:ns, &
2101	! nys:nyn+1,nzb_do:nzt_do), &
2102	! start = (/ 1, nys+1, 1, &
2103	! do2d_yz_time_count(av) /), &
2104	! count = (/ ns, nyn-nys+2, &
2105	! nzt_do-nzb_do+1, 1 /) )
2106	! ELSE
2107	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
2108	id_var_do2d(av,if), &
2109	local_2d_sections(1:ns,nys:nyn, &
2110	nzb_do:nzt_do), &
2111	start = (/ 1, nys+1, 1, &
2112	do2d_yz_time_count(av) /), &
2113	count = (/ ns, nyn-nys+1, &
2114	nzt_do-nzb_do+1, 1 /) )
2115	! ENDIF
2116
2117	CALL netcdf_handle_error( 'data_output_2d', 60 )
2118
2119	CASE DEFAULT
2120	message_string = 'unknown cross-section: ' // TRIM( mode )
2121	CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
2122
2123	END SELECT
2124
2125	ENDIF
2126	#endif
2127	ENDIF
2128
2129	if = if + 1
2130	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
2131	do2d_mode = do2d(av,if)(l-1:l)
2132
2133	ENDDO
2134
2135	!
2136	!-- Deallocate temporary arrays.
2137	IF ( ALLOCATED( level_z ) ) DEALLOCATE( level_z )
2138	IF ( netcdf_data_format > 4 ) THEN
2139	DEALLOCATE( local_pf, local_2d, local_2d_sections )
2140	IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
2141	ENDIF
2142	#if defined( __parallel )
2143	IF ( .NOT. data_output_2d_on_each_pe .AND. myid == 0 ) THEN
2144	DEALLOCATE( total_2d )
2145	ENDIF
2146	#endif
2147
2148	!
2149	!-- Close plot output file.
2150	file_id = 20 + s_ind
2151
2152	IF ( data_output_2d_on_each_pe ) THEN
2153	DO i = 0, io_blocks-1
2154	IF ( i == io_group ) THEN
2155	CALL close_file( file_id )
2156	ENDIF
2157	#if defined( __parallel )
2158	CALL MPI_BARRIER( comm2d, ierr )
2159	#endif
2160	ENDDO
2161	ELSE
2162	IF ( myid == 0 ) CALL close_file( file_id )
2163	ENDIF
2164
2165	CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
2166
2167	END SUBROUTINE data_output_2d

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |