Home

Context Navigation

source: palm/trunk/SOURCE/surface_data_output_mod.f90 @ 4682

Last change on this file since 4682 was 4671, checked in by pavelkrc, 5 years ago
Radiative transfer model RTM version 4.1
Property svn:keywords set to `Id`
File size: 244.9 KB

Line
1	!> @file surface_data_output_mod.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 terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2020 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: surface_data_output_mod.f90 4671 2020-09-09 20:27:58Z pavelkrc $
27	! Implementation of downward facing USM and LSM surfaces
28	!
29	! 4601 2020-07-14 12:06:09Z suehring
30	! Minor simplification in name creation for IO variables in restart files.
31	!
32	! 4600 2020-07-13 18:50:12Z suehring
33	! - Change: adjustmens for mpi-io - surface data is transformed to a 2D-based surface array
34	! before writing.
35	! - Bugfix in counting of surface elements
36	! - Bugfix in data-output of averaged surface data in case of restarts
37	!
38	! 4577 2020-06-25 09:53:58Z raasch
39	! File re-formatted to follow the PALM coding standard
40	!
41	!
42	! 4547 2020-05-27 09:05:24Z moh.hefny
43	! Added surface albedo and emissivity, which are defined using the tile approach
44	!
45	! 4535 2020-05-15 12:07:23Z raasch
46	! Bugfix for restart data format query
47	!
48	! 4535 2020-05-15 12:07:23Z raasch
49	! Bugfix for restart data format query
50	!
51	! 4517 2020-05-03 14:29:30Z raasch
52	! Added restart with MPI-IO for reading local arrays
53	!
54	! 4502 2020-04-17 16:14:16Z schwenkel
55	! Implementation of ice microphysics
56	!
57	! 4500 2020-04-17 10:12:45Z suehring
58	! - Correct output of ground/wall heat flux at USM surfaces
59	! - Add conversion factor to heat and momentum-flux output
60	!
61	! 4495 2020-04-13 20:11:20Z raasch
62	! Restart data handling with MPI-IO added
63	!
64	! 4444 2020-03-05 15:59:50Z raasch
65	! Bugfix: cpp-directives for serial mode added
66	!
67	! 4360 2020-01-07 11:25:50Z suehring
68	! Fix wrongly declared nc_stat variable in surface_data_output_mod
69	!
70	! 4205 2019-08-30 13:25:00Z suehring
71	! - Correct x,y-coordinates of vertical surfaces in netcdf output
72	! - Change definition of azimuth angle, reference is north 0 degree
73	! - Zenith angle is always defined, also for vertical surfaces where it is 90 degree, while azimuth
74	! angle is only defined for vertical surfaces, not for horizontal ones
75	!
76	! 4182 2019-08-22 15:20:23Z scharf
77	! Corrected "Former revisions" section
78	!
79	! 4129 2019-07-31 12:56:07Z gronemeier
80	! - Bugfix: corrected loop over horizontal default surfaces
81	! - Change default setting of to_vtk and to_netcdf
82	!
83	! 4029 2019-06-14 14:04:35Z raasch
84	! Netcdf variable NF90_NOFILL is used as argument instead of "1" in call to NF90_DEF_VAR_FILL
85	!
86	! 3881 2019-04-10 09:31:22Z suehring
87	! Check for zero output timestep (not allowed in parallel NetCDF output mode)
88	!
89	! 3817 2019-03-26 13:53:57Z suehring
90	! Correct output coordinates of vertical surface elements
91	!
92	! 3766 2019-02-26 16:23:41Z raasch
93	! Bugfix in surface_data_output_rrd_local (variable k removed)
94	!
95	! 3762 2019-02-25 16:54:16Z suehring
96	! Remove unused variables and add preprocessor directives for variables that are used only when
97	! netcdf4 is defined
98	!
99	! 3745 2019-02-15 18:57:56Z suehring
100	! Output of waste_heat and innermost wall flux from indoor model
101	!
102	! 3744 2019-02-15 18:38:58Z suehring
103	! Add azimuth and zenith to output file; set long-name attributes; clean-up coding layout
104	!
105	! 3735 2019-02-12 09:52:40Z suehring
106	! - Split initialization into initialization of arrays and further initialization in order to enable
107	! reading of restart data.
108	! - Consider restarts in surface data averaging.
109	! - Correct error message numbers
110	!
111	! 3731 2019-02-11 13:06:27Z suehring
112	! Bugfix: add cpp options
113	!
114	! 3727 2019-02-08 14:52:10Z gronemeier
115	! Enable NetCDF output for surface data (suehring, gronemeier)
116	!
117	! 3691 2019-01-23 09:57:04Z suehring
118	! Add output of surface-parallel flow speed
119	!
120	! 3648 2019-01-02 16:35:46Z suehring
121	! Rename module and subroutines
122	! 3420 2018-10-24 17:30:08Z gronemeier
123	! Initial implementation from Klaus Ketelsen and Matthias Suehring
124	!
125	!
126	! Authors:
127	! --------
128	! @author Klaus Ketelsen, Matthias Suehring, Tobias Gronemeier
129	!
130	!--------------------------------------------------------------------------------------------------!
131	! Description:
132	! ------------
133	!> Generate output for surface data.
134	!>
135	!> @todo Create namelist file for post-processing tool.
136	!--------------------------------------------------------------------------------------------------!
137
138	MODULE surface_data_output_mod
139
140	USE kinds
141
142	USE arrays_3d, &
143	ONLY: heatflux_output_conversion, &
144	momentumflux_output_conversion, &
145	waterflux_output_conversion, &
146	zu, &
147	zw
148
149	USE control_parameters, &
150	ONLY: coupling_char, &
151	data_output_during_spinup, &
152	end_time, &
153	message_string, &
154	restart_data_format_output, &
155	run_description_header, &
156	simulated_time_at_begin, &
157	spinup_time, &
158	surface_output
159
160	USE grid_variables, &
161	ONLY: dx, &
162	dy
163
164	USE indices, &
165	ONLY: nxl, &
166	nxr, &
167	nys, &
168	nyn, &
169	nzb, &
170	nzt
171
172	#if defined( __netcdf )
173	USE NETCDF
174	#endif
175
176	USE netcdf_data_input_mod, &
177	ONLY: init_model
178
179	USE netcdf_interface, &
180	ONLY: nc_stat, &
181	netcdf_create_att, &
182	netcdf_create_dim, &
183	netcdf_create_file, &
184	netcdf_create_global_atts, &
185	netcdf_create_var, &
186	netcdf_data_format, &
187	netcdf_handle_error
188
189	USE pegrid
190
191	USE restart_data_mpi_io_mod, &
192	ONLY: rrd_mpi_io, &
193	rd_mpi_io_check_array, &
194	rrd_mpi_io_surface, &
195	rd_mpi_io_surface_filetypes, &
196	wrd_mpi_io, &
197	wrd_mpi_io_surface
198
199	USE surface_mod, &
200	ONLY: ind_pav_green, &
201	ind_veg_wall, &
202	ind_wat_win, &
203	surf_def_h, &
204	surf_def_v, &
205	surf_lsm_h, &
206	surf_lsm_v, &
207	surf_usm_h, &
208	surf_usm_v
209
210	IMPLICIT NONE
211
212	TYPE surf_out !< data structure which contains all surfaces elements of all types on subdomain
213
214	INTEGER(iwp) :: ns !< number of surface elements on subdomain
215	INTEGER(iwp) :: ns_total !< total number of surface elements
216	INTEGER(iwp) :: npoints !< number of points / vertices which define a surface element (on subdomain)
217	INTEGER(iwp) :: npoints_total !< total number of points / vertices which define a surface element
218
219	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: s !< coordinate for NetCDF output, number of the surface element
220
221	REAL(wp) :: fillvalue = -9999.0_wp !< fillvalue for surface elements which are not defined
222
223	REAL(wp), DIMENSION(:), ALLOCATABLE :: azimuth !< azimuth orientation coordinate for NetCDF output
224	REAL(wp), DIMENSION(:), ALLOCATABLE :: es_utm !< E-UTM coordinate for NetCDF output
225	REAL(wp), DIMENSION(:), ALLOCATABLE :: ns_utm !< E-UTM coordinate for NetCDF output
226	REAL(wp), DIMENSION(:), ALLOCATABLE :: xs !< x-coordinate for NetCDF output
227	REAL(wp), DIMENSION(:), ALLOCATABLE :: ys !< y-coordinate for NetCDF output
228	REAL(wp), DIMENSION(:), ALLOCATABLE :: zs !< z-coordinate for NetCDF output
229	REAL(wp), DIMENSION(:), ALLOCATABLE :: zenith !< zenith orientation coordinate for NetCDF output
230	REAL(wp), DIMENSION(:), ALLOCATABLE :: var_out !< output variable
231	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: var_av !< variable used for averaging
232	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: points !< points / vertices of a surface element
233	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: polygons !< polygon data of a surface element
234	END TYPE surf_out
235
236	CHARACTER(LEN=100), DIMENSION(300) :: data_output_surf = ' ' !< namelist variable which describes the output variables
237	CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf = ' ' !< internal variable which describes the output variables
238	!< and separates averaged from non-averaged output
239	CHARACTER(LEN=100), DIMENSION(0:1,300) :: dosurf_unit = ' ' !< internal variable which holds the unit of the given output
240	!< variable
241
242	INTEGER(iwp) :: average_count_surf = 0 !< number of ensemble members used for averaging
243	INTEGER(iwp) :: dosurf_no(0:1) = 0 !< number of surface output quantities
244	#if defined( __netcdf4_parallel )
245	INTEGER(iwp) :: oldmode !< save old set-fill-mode of netcdf file (not needed, but required for routine call)
246
247	INTEGER(iwp), DIMENSION(0:1) :: dosurf_time_count = 0 !< count of output time steps
248	INTEGER(iwp), DIMENSION(0:1) :: id_dim_s_surf !< netcdf ID for dimension s
249	INTEGER(iwp), DIMENSION(0:1) :: id_dim_time_surf !< netcdf ID for dimension time
250	INTEGER(iwp), DIMENSION(0:1) :: id_set_surf !< netcdf ID for file
251	INTEGER(iwp), DIMENSION(0:1) :: id_var_azimuth_surf !< netcdf ID for variable azimuth
252	INTEGER(iwp), DIMENSION(0:1) :: id_var_etum_surf !< netcdf ID for variable Es_UTM
253	INTEGER(iwp), DIMENSION(0:1) :: id_var_nutm_surf !< netcdf ID for variable Ns_UTM
254	INTEGER(iwp), DIMENSION(0:1) :: id_var_time_surf !< netcdf ID for variable time
255	INTEGER(iwp), DIMENSION(0:1) :: id_var_s_surf !< netcdf ID for variable s
256	INTEGER(iwp), DIMENSION(0:1) :: id_var_xs_surf !< netcdf ID for variable xs
257	INTEGER(iwp), DIMENSION(0:1) :: id_var_ys_surf !< netcdf ID for variable ys
258	INTEGER(iwp), DIMENSION(0:1) :: id_var_zenith_surf !< netcdf ID for variable zenith
259	INTEGER(iwp), DIMENSION(0:1) :: id_var_zs_surf !< netcdf ID for variable zs
260	INTEGER(iwp), DIMENSION(0:1) :: ntdim_surf !< number of output time steps
261
262	INTEGER(iwp), DIMENSION(0:1,300) :: id_var_dosurf !< netcdf ID for output variables
263	#endif
264
265	LOGICAL :: first_output(0:1) = .FALSE. !< true if first output was already called
266	LOGICAL :: to_netcdf = .FALSE. !< flag indicating parallel NetCDF output
267	LOGICAL :: to_vtk = .FALSE. !< flag indicating binary surface-data output that can be further
268	!< processed to VTK format
269
270	REAL(wp) :: averaging_interval_surf = 9999999.9_wp !< averaging interval
271	REAL(wp) :: dt_dosurf = 9999999.9_wp !< time interval for instantaneous data output
272	REAL(wp) :: dt_dosurf_av = 9999999.9_wp !< time interval for averaged data output
273	REAL(wp) :: skip_time_dosurf = 0.0_wp !< skip time for instantaneous data output
274	REAL(wp) :: skip_time_dosurf_av = 0.0_wp !< skip time for averaged data output
275	REAL(wp) :: time_dosurf = 0.0_wp !< internal counter variable to check for instantaneous data output
276	REAL(wp) :: time_dosurf_av = 0.0_wp !< internal counter variable to check for averaged data output
277
278	TYPE(surf_out) :: surfaces !< variable which contains all required output information
279
280	SAVE
281
282	PRIVATE
283
284	INTERFACE surface_data_output
285	MODULE PROCEDURE surface_data_output
286	END INTERFACE surface_data_output
287
288	INTERFACE surface_data_output_averaging
289	MODULE PROCEDURE surface_data_output_averaging
290	END INTERFACE surface_data_output_averaging
291
292	INTERFACE surface_data_output_check_parameters
293	MODULE PROCEDURE surface_data_output_check_parameters
294	END INTERFACE surface_data_output_check_parameters
295
296	INTERFACE surface_data_output_init
297	MODULE PROCEDURE surface_data_output_init
298	END INTERFACE surface_data_output_init
299
300	INTERFACE surface_data_output_init_arrays
301	MODULE PROCEDURE surface_data_output_init_arrays
302	END INTERFACE surface_data_output_init_arrays
303
304	INTERFACE surface_data_output_last_action
305	MODULE PROCEDURE surface_data_output_last_action
306	END INTERFACE surface_data_output_last_action
307
308	INTERFACE surface_data_output_parin
309	MODULE PROCEDURE surface_data_output_parin
310	END INTERFACE surface_data_output_parin
311
312	INTERFACE surface_data_output_rrd_global
313	MODULE PROCEDURE surface_data_output_rrd_global_ftn
314	MODULE PROCEDURE surface_data_output_rrd_global_mpi
315	END INTERFACE surface_data_output_rrd_global
316
317	INTERFACE surface_data_output_rrd_local
318	MODULE PROCEDURE surface_data_output_rrd_local_ftn
319	MODULE PROCEDURE surface_data_output_rrd_local_mpi
320	END INTERFACE surface_data_output_rrd_local
321
322	INTERFACE surface_data_output_wrd_global
323	MODULE PROCEDURE surface_data_output_wrd_global
324	END INTERFACE surface_data_output_wrd_global
325
326	INTERFACE surface_data_output_wrd_local
327	MODULE PROCEDURE surface_data_output_wrd_local
328	END INTERFACE surface_data_output_wrd_local
329
330	INTERFACE surface_data_output_sum_up
331	MODULE PROCEDURE surface_data_output_sum_up_1d
332	MODULE PROCEDURE surface_data_output_sum_up_2d
333	END INTERFACE surface_data_output_sum_up
334
335	INTERFACE surface_data_output_collect
336	MODULE PROCEDURE surface_data_output_collect_1d
337	MODULE PROCEDURE surface_data_output_collect_2d
338	END INTERFACE surface_data_output_collect
339
340	!
341	!--Public subroutines
342	PUBLIC surface_data_output, &
343	surface_data_output_averaging, &
344	surface_data_output_check_parameters, &
345	surface_data_output_init, &
346	surface_data_output_init_arrays, &
347	surface_data_output_last_action, &
348	surface_data_output_parin, &
349	surface_data_output_rrd_global, &
350	surface_data_output_rrd_local, &
351	surface_data_output_wrd_local, &
352	surface_data_output_wrd_global
353	!
354	!--Public variables
355	PUBLIC average_count_surf, &
356	averaging_interval_surf, &
357	dt_dosurf, &
358	dt_dosurf_av, &
359	skip_time_dosurf, &
360	skip_time_dosurf_av, &
361	time_dosurf, &
362	time_dosurf_av
363
364	CONTAINS
365
366	!--------------------------------------------------------------------------------------------------!
367	! Description:
368	! ------------
369	!> This routine counts the number of surfaces on each core and allocates arrays.
370	!--------------------------------------------------------------------------------------------------!
371	SUBROUTINE surface_data_output_init_arrays
372
373	IMPLICIT NONE
374
375	!
376	!-- Determine the number of surface elements on subdomain
377	surfaces%ns = surf_def_h(0)%ns + surf_lsm_h(0)%ns + surf_usm_h(0)%ns & !horizontal upward-facing
378	+ surf_def_h(1)%ns + surf_lsm_h(1)%ns + surf_usm_h(1)%ns & !horizontal downard-facing
379	+ surf_def_v(0)%ns + surf_lsm_v(0)%ns + surf_usm_v(0)%ns & !northward-facing
380	+ surf_def_v(1)%ns + surf_lsm_v(1)%ns + surf_usm_v(1)%ns & !southward-facing
381	+ surf_def_v(2)%ns + surf_lsm_v(2)%ns + surf_usm_v(2)%ns & !westward-facing
382	+ surf_def_v(3)%ns + surf_lsm_v(3)%ns + surf_usm_v(3)%ns !eastward-facing
383	!
384	!-- Determine the total number of surfaces in the model domain
385	#if defined( __parallel )
386	CALL MPI_ALLREDUCE( surfaces%ns, surfaces%ns_total, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
387	#else
388	surfaces%ns_total = surfaces%ns
389	#endif
390	!
391	!-- Allocate output variable and set to _FillValue attribute
392	ALLOCATE ( surfaces%var_out(1:surfaces%ns) )
393	surfaces%var_out = surfaces%fillvalue
394	!
395	!-- If there is an output of time average output variables, allocate the required array.
396	IF ( dosurf_no(1) > 0 ) THEN
397	ALLOCATE ( surfaces%var_av(1:surfaces%ns,1:dosurf_no(1)) )
398	surfaces%var_av = 0.0_wp
399	ENDIF
400
401	END SUBROUTINE surface_data_output_init_arrays
402
403
404	!--------------------------------------------------------------------------------------------------!
405	! Description:
406	! ------------
407	!> Initialization surface-data output data structure: calculation of vertices and polygon data for
408	!> the surface elements, allocation of required arrays.
409	!--------------------------------------------------------------------------------------------------!
410	SUBROUTINE surface_data_output_init
411
412	IMPLICIT NONE
413
414	#if defined( __netcdf4_parallel )
415	CHARACTER (LEN=100) :: filename !< name of output file
416	CHARACTER (LEN=80) :: time_average_text !< string written to file attribute time_avg
417	CHARACTER (LEN=4000) :: var_list !< list of variables written to NetCDF file
418
419	INTEGER(iwp) :: av !< flag for averaged (=1) and non-averaged (=0) data
420	#endif
421	INTEGER(iwp) :: i !< grid index in x-direction, also running variable for counting non-average data output
422	INTEGER(iwp) :: j !< grid index in y-direction, also running variable for counting average data output
423	INTEGER(iwp) :: k !< grid index in z-direction
424	INTEGER(iwp) :: l !< running index for surface-element orientation
425	INTEGER(iwp) :: m !< running index for surface elements
426	INTEGER(iwp) :: mm !< local counting variable for surface elements
427	INTEGER(iwp) :: npg !< counter variable for all surface elements ( or polygons )
428	INTEGER(iwp) :: point_index_count !< local counter variable for point index
429	INTEGER(iwp) :: start_count !< local start counter for the surface index
430
431	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_points_on_pe !< array which contains the number of points on all mpi ranks
432	INTEGER(iwp), DIMENSION(0:numprocs-1) :: num_surfaces_on_pe !< array which contains the number of surfaces on all mpi ranks
433	INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:,:) :: point_index !< dummy array used to check where the reference points for
434	!< surface polygons are located
435
436	REAL(wp) :: az !< azimuth angle, indicated the vertical orientation of a surface element
437	REAL(wp) :: off_x !< grid offset in x-direction between the stored grid index and the actual wall
438	REAL(wp) :: off_y !< grid offset in y-direction between the stored grid index and the actual wall
439	#if defined( __netcdf4_parallel )
440	REAL(wp), DIMENSION(:), ALLOCATABLE :: netcdf_data_1d !< dummy array to output 1D data into netcdf file
441	#endif
442
443	!
444	!-- If output to VTK format is enabled, initialize point and polygon data.
445	!-- In a first step, count the number of points which are defining the surfaces and the polygons.
446	IF ( to_vtk ) THEN
447	ALLOCATE( point_index(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
448	point_index = -1
449	!
450	!-- Horizontal default surfaces
451	surfaces%npoints = 0
452	DO l = 0, 1
453	DO m = 1, surf_def_h(l)%ns
454	!
455	!-- Determine the indices of the respective grid cell inside the topography
456	i = surf_def_h(l)%i(m) + surf_def_h(l)%ioff
457	j = surf_def_h(l)%j(m) + surf_def_h(l)%joff
458	k = surf_def_h(l)%k(m) + surf_def_h(l)%koff
459	!
460	!-- Check if the vertices that define the surface element are already defined, if not,
461	!-- increment the counter.
462	IF ( point_index(k,j,i) < 0 ) THEN
463	surfaces%npoints = surfaces%npoints + 1
464	point_index(k,j,i) = surfaces%npoints - 1
465	ENDIF
466	IF ( point_index(k,j,i+1) < 0 ) THEN
467	surfaces%npoints = surfaces%npoints + 1
468	point_index(k,j,i+1) = surfaces%npoints - 1
469	ENDIF
470	IF ( point_index(k,j+1,i+1) < 0 ) THEN
471	surfaces%npoints = surfaces%npoints + 1
472	point_index(k,j+1,i+1) = surfaces%npoints - 1
473	ENDIF
474	IF ( point_index(k,j+1,i) < 0 ) THEN
475	surfaces%npoints = surfaces%npoints + 1
476	point_index(k,j+1,i) = surfaces%npoints - 1
477	ENDIF
478	ENDDO
479	DO m = 1, surf_lsm_h(l)%ns
480	i = surf_lsm_h(l)%i(m) + surf_lsm_h(l)%ioff
481	j = surf_lsm_h(l)%j(m) + surf_lsm_h(l)%joff
482	k = surf_lsm_h(l)%k(m) + surf_lsm_h(l)%koff
483
484	IF ( point_index(k,j,i) < 0 ) THEN
485	surfaces%npoints = surfaces%npoints + 1
486	point_index(k,j,i) = surfaces%npoints - 1
487	ENDIF
488	IF ( point_index(k,j,i+1) < 0 ) THEN
489	surfaces%npoints = surfaces%npoints + 1
490	point_index(k,j,i+1) = surfaces%npoints - 1
491	ENDIF
492	IF ( point_index(k,j+1,i+1) < 0 ) THEN
493	surfaces%npoints = surfaces%npoints + 1
494	point_index(k,j+1,i+1) = surfaces%npoints - 1
495	ENDIF
496	IF ( point_index(k,j+1,i) < 0 ) THEN
497	surfaces%npoints = surfaces%npoints + 1
498	point_index(k,j+1,i) = surfaces%npoints - 1
499	ENDIF
500	ENDDO
501	DO m = 1, surf_usm_h(l)%ns
502	i = surf_usm_h(l)%i(m) + surf_usm_h(l)%ioff
503	j = surf_usm_h(l)%j(m) + surf_usm_h(l)%joff
504	k = surf_usm_h(l)%k(m) + surf_usm_h(l)%koff
505
506	IF ( point_index(k,j,i) < 0 ) THEN
507	surfaces%npoints = surfaces%npoints + 1
508	point_index(k,j,i) = surfaces%npoints - 1
509	ENDIF
510	IF ( point_index(k,j,i+1) < 0 ) THEN
511	surfaces%npoints = surfaces%npoints + 1
512	point_index(k,j,i+1) = surfaces%npoints - 1
513	ENDIF
514	IF ( point_index(k,j+1,i+1) < 0 ) THEN
515	surfaces%npoints = surfaces%npoints + 1
516	point_index(k,j+1,i+1) = surfaces%npoints - 1
517	ENDIF
518	IF ( point_index(k,j+1,i) < 0 ) THEN
519	surfaces%npoints = surfaces%npoints + 1
520	point_index(k,j+1,i) = surfaces%npoints - 1
521	ENDIF
522	ENDDO
523	ENDDO
524	!
525	!-- Vertical surfaces
526	DO l = 0, 3
527	DO m = 1, surf_def_v(l)%ns
528	!
529	!-- Determine the indices of the respective grid cell inside the topography. Please note,
530	!-- j-index for north-facing surfaces ( l==0 ) is identical to the reference j-index
531	!-- outside the grid box. Equivalent for east-facing surfaces and i-index.
532	i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
533	j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
534	k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
535	!
536	!-- Lower left /front vertex
537	IF ( point_index(k,j,i) < 0 ) THEN
538	surfaces%npoints = surfaces%npoints + 1
539	point_index(k,j,i) = surfaces%npoints - 1
540	ENDIF
541	!
542	!-- Upper / lower right index for north- and south-facing surfaces
543	IF ( l == 0 .OR. l == 1 ) THEN
544	IF ( point_index(k,j,i+1) < 0 ) THEN
545	surfaces%npoints = surfaces%npoints + 1
546	point_index(k,j,i+1) = surfaces%npoints - 1
547	ENDIF
548	IF ( point_index(k+1,j,i+1) < 0 ) THEN
549	surfaces%npoints = surfaces%npoints + 1
550	point_index(k+1,j,i+1) = surfaces%npoints - 1
551	ENDIF
552	!
553	!-- Upper / lower front index for east- and west-facing surfaces
554	ELSEIF ( l == 2 .OR. l == 3 ) THEN
555	IF ( point_index(k,j+1,i) < 0 ) THEN
556	surfaces%npoints = surfaces%npoints + 1
557	point_index(k,j+1,i) = surfaces%npoints - 1
558	ENDIF
559	IF ( point_index(k+1,j+1,i) < 0 ) THEN
560	surfaces%npoints = surfaces%npoints + 1
561	point_index(k+1,j+1,i) = surfaces%npoints - 1
562	ENDIF
563	ENDIF
564	!
565	!-- Upper left / front vertex
566	IF ( point_index(k+1,j,i) < 0 ) THEN
567	surfaces%npoints = surfaces%npoints + 1
568	point_index(k+1,j,i) = surfaces%npoints - 1
569	ENDIF
570	ENDDO
571	DO m = 1, surf_lsm_v(l)%ns
572	i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
573	j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
574	k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
575	!
576	!-- Lower left /front vertex
577	IF ( point_index(k,j,i) < 0 ) THEN
578	surfaces%npoints = surfaces%npoints + 1
579	point_index(k,j,i) = surfaces%npoints - 1
580	ENDIF
581	!
582	!-- Upper / lower right index for north- and south-facing surfaces
583	IF ( l == 0 .OR. l == 1 ) THEN
584	IF ( point_index(k,j,i+1) < 0 ) THEN
585	surfaces%npoints = surfaces%npoints + 1
586	point_index(k,j,i+1) = surfaces%npoints - 1
587	ENDIF
588	IF ( point_index(k+1,j,i+1) < 0 ) THEN
589	surfaces%npoints = surfaces%npoints + 1
590	point_index(k+1,j,i+1) = surfaces%npoints - 1
591	ENDIF
592	!
593	!-- Upper / lower front index for east- and west-facing surfaces
594	ELSEIF ( l == 2 .OR. l == 3 ) THEN
595	IF ( point_index(k,j+1,i) < 0 ) THEN
596	surfaces%npoints = surfaces%npoints + 1
597	point_index(k,j+1,i) = surfaces%npoints - 1
598	ENDIF
599	IF ( point_index(k+1,j+1,i) < 0 ) THEN
600	surfaces%npoints = surfaces%npoints + 1
601	point_index(k+1,j+1,i) = surfaces%npoints - 1
602	ENDIF
603	ENDIF
604	!
605	!-- Upper left / front vertex
606	IF ( point_index(k+1,j,i) < 0 ) THEN
607	surfaces%npoints = surfaces%npoints + 1
608	point_index(k+1,j,i) = surfaces%npoints - 1
609	ENDIF
610	ENDDO
611
612	DO m = 1, surf_usm_v(l)%ns
613	i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
614	j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
615	k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
616	!
617	!-- Lower left /front vertex
618	IF ( point_index(k,j,i) < 0 ) THEN
619	surfaces%npoints = surfaces%npoints + 1
620	point_index(k,j,i) = surfaces%npoints - 1
621	ENDIF
622	!
623	!-- Upper / lower right index for north- and south-facing surfaces
624	IF ( l == 0 .OR. l == 1 ) THEN
625	IF ( point_index(k,j,i+1) < 0 ) THEN
626	surfaces%npoints = surfaces%npoints + 1
627	point_index(k,j,i+1) = surfaces%npoints - 1
628	ENDIF
629	IF ( point_index(k+1,j,i+1) < 0 ) THEN
630	surfaces%npoints = surfaces%npoints + 1
631	point_index(k+1,j,i+1) = surfaces%npoints - 1
632	ENDIF
633	!
634	!-- Upper / lower front index for east- and west-facing surfaces
635	ELSEIF ( l == 2 .OR. l == 3 ) THEN
636	IF ( point_index(k,j+1,i) < 0 ) THEN
637	surfaces%npoints = surfaces%npoints + 1
638	point_index(k,j+1,i) = surfaces%npoints - 1
639	ENDIF
640	IF ( point_index(k+1,j+1,i) < 0 ) THEN
641	surfaces%npoints = surfaces%npoints + 1
642	point_index(k+1,j+1,i) = surfaces%npoints - 1
643	ENDIF
644	ENDIF
645	!
646	!-- Upper left / front vertex
647	IF ( point_index(k+1,j,i) < 0 ) THEN
648	surfaces%npoints = surfaces%npoints + 1
649	point_index(k+1,j,i) = surfaces%npoints - 1
650	ENDIF
651	ENDDO
652
653	ENDDO
654	!
655	!-- Allocate the number of points and polygons. Note, the number of polygons is identical to the
656	!-- number of surfaces elements, whereas the number of points (vertices), which define the
657	!-- polygons, can be larger.
658	ALLOCATE( surfaces%points(3,1:surfaces%npoints) )
659	ALLOCATE( surfaces%polygons(5,1:surfaces%ns) )
660	!
661	!-- Note, PARAVIEW expects consecutively ordered points, in order to unambiguously identify
662	!-- surfaces. Hence, all PEs should know where they start counting, depending on the number of
663	!-- points on the other PE's with lower MPI rank.
664	#if defined( __parallel )
665	CALL MPI_ALLGATHER( surfaces%npoints, 1, MPI_INTEGER, num_points_on_pe, 1, MPI_INTEGER, &
666	comm2d, ierr )
667	#else
668	num_points_on_pe = surfaces%npoints
669	#endif
670
671	!
672	!-- After the number of vertices is counted, repeat the loops and define the vertices. Start with
673	!-- the horizontal default surfaces. First, however, determine the offset where couting of points
674	!-- should be started, which is the sum of points of all PE's with lower MPI rank.
675	i = 0
676	point_index_count = 0
677	DO WHILE ( i < myid .AND. i <= SIZE( num_points_on_pe ) )
678	point_index_count = point_index_count + num_points_on_pe(i)
679	i = i + 1
680	ENDDO
681
682	surfaces%npoints = 0
683	point_index = -1
684	npg = 0
685
686	DO l = 0, 1
687	DO m = 1, surf_def_h(l)%ns
688	!
689	!-- Determine the indices of the respective grid cell inside the topography.
690	i = surf_def_h(l)%i(m) + surf_def_h(l)%ioff
691	j = surf_def_h(l)%j(m) + surf_def_h(l)%joff
692	k = surf_def_h(l)%k(m) + surf_def_h(l)%koff
693	!
694	!-- Check if the vertices that define the surface element are already defined, if not,
695	!-- increment the counter.
696	IF ( point_index(k,j,i) < 0 ) THEN
697	surfaces%npoints = surfaces%npoints + 1
698	point_index(k,j,i) = point_index_count
699	point_index_count = point_index_count + 1
700	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
701	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
702	surfaces%points(3,surfaces%npoints) = zw(k)
703	ENDIF
704	IF ( point_index(k,j,i+1) < 0 ) THEN
705	surfaces%npoints = surfaces%npoints + 1
706	point_index(k,j,i+1) = point_index_count
707	point_index_count = point_index_count + 1
708	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
709	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
710	surfaces%points(3,surfaces%npoints) = zw(k)
711	ENDIF
712	IF ( point_index(k,j+1,i+1) < 0 ) THEN
713	surfaces%npoints = surfaces%npoints + 1
714	point_index(k,j+1,i+1) = point_index_count
715	point_index_count = point_index_count + 1
716	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
717	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
718	surfaces%points(3,surfaces%npoints) = zw(k)
719	ENDIF
720	IF ( point_index(k,j+1,i) < 0 ) THEN
721	surfaces%npoints = surfaces%npoints + 1
722	point_index(k,j+1,i) = point_index_count
723	point_index_count = point_index_count + 1
724	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
725	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
726	surfaces%points(3,surfaces%npoints) = zw(k)
727	ENDIF
728
729	npg = npg + 1
730	surfaces%polygons(1,npg) = 4
731	surfaces%polygons(2,npg) = point_index(k,j,i)
732	surfaces%polygons(3,npg) = point_index(k,j,i+1)
733	surfaces%polygons(4,npg) = point_index(k,j+1,i+1)
734	surfaces%polygons(5,npg) = point_index(k,j+1,i)
735	ENDDO
736	DO m = 1, surf_lsm_h(l)%ns
737	i = surf_lsm_h(l)%i(m) + surf_lsm_h(l)%ioff
738	j = surf_lsm_h(l)%j(m) + surf_lsm_h(l)%joff
739	k = surf_lsm_h(l)%k(m) + surf_lsm_h(l)%koff
740	IF ( point_index(k,j,i) < 0 ) THEN
741	surfaces%npoints = surfaces%npoints + 1
742	point_index(k,j,i) = point_index_count
743	point_index_count = point_index_count + 1
744	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
745	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
746	surfaces%points(3,surfaces%npoints) = zw(k)
747	ENDIF
748	IF ( point_index(k,j,i+1) < 0 ) THEN
749	surfaces%npoints = surfaces%npoints + 1
750	point_index(k,j,i+1) = point_index_count
751	point_index_count = point_index_count + 1
752	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
753	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
754	surfaces%points(3,surfaces%npoints) = zw(k)
755	ENDIF
756	IF ( point_index(k,j+1,i+1) < 0 ) THEN
757	surfaces%npoints = surfaces%npoints + 1
758	point_index(k,j+1,i+1) = point_index_count
759	point_index_count = point_index_count + 1
760	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
761	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
762	surfaces%points(3,surfaces%npoints) = zw(k)
763	ENDIF
764	IF ( point_index(k,j+1,i) < 0 ) THEN
765	surfaces%npoints = surfaces%npoints + 1
766	point_index(k,j+1,i) = point_index_count
767	point_index_count = point_index_count + 1
768	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
769	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
770	surfaces%points(3,surfaces%npoints) = zw(k)
771	ENDIF
772
773	npg = npg + 1
774	surfaces%polygons(1,npg) = 4
775	surfaces%polygons(2,npg) = point_index(k,j,i)
776	surfaces%polygons(3,npg) = point_index(k,j,i+1)
777	surfaces%polygons(4,npg) = point_index(k,j+1,i+1)
778	surfaces%polygons(5,npg) = point_index(k,j+1,i)
779	ENDDO
780
781	DO m = 1, surf_usm_h(l)%ns
782	i = surf_usm_h(l)%i(m) + surf_usm_h(l)%ioff
783	j = surf_usm_h(l)%j(m) + surf_usm_h(l)%joff
784	k = surf_usm_h(l)%k(m) + surf_usm_h(l)%koff
785
786	IF ( point_index(k,j,i) < 0 ) THEN
787	surfaces%npoints = surfaces%npoints + 1
788	point_index(k,j,i) = point_index_count
789	point_index_count = point_index_count + 1
790	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
791	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
792	surfaces%points(3,surfaces%npoints) = zw(k)
793	ENDIF
794	IF ( point_index(k,j,i+1) < 0 ) THEN
795	surfaces%npoints = surfaces%npoints + 1
796	point_index(k,j,i+1) = point_index_count
797	point_index_count = point_index_count + 1
798	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
799	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
800	surfaces%points(3,surfaces%npoints) = zw(k)
801	ENDIF
802	IF ( point_index(k,j+1,i+1) < 0 ) THEN
803	surfaces%npoints = surfaces%npoints + 1
804	point_index(k,j+1,i+1) = point_index_count
805	point_index_count = point_index_count + 1
806	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
807	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
808	surfaces%points(3,surfaces%npoints) = zw(k)
809	ENDIF
810	IF ( point_index(k,j+1,i) < 0 ) THEN
811	surfaces%npoints = surfaces%npoints + 1
812	point_index(k,j+1,i) = point_index_count
813	point_index_count = point_index_count + 1
814	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
815	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
816	surfaces%points(3,surfaces%npoints) = zw(k)
817	ENDIF
818
819	npg = npg + 1
820	surfaces%polygons(1,npg) = 4
821	surfaces%polygons(2,npg) = point_index(k,j,i)
822	surfaces%polygons(3,npg) = point_index(k,j,i+1)
823	surfaces%polygons(4,npg) = point_index(k,j+1,i+1)
824	surfaces%polygons(5,npg) = point_index(k,j+1,i)
825	ENDDO
826	ENDDO
827
828	DO l = 0, 3
829	DO m = 1, surf_def_v(l)%ns
830	!
831	!-- Determine the indices of the respective grid cell inside the topography.
832	!-- NOTE, j-index for north-facing surfaces ( l==0 ) is identical to the reference j-index
833	!-- outside the grid box. Equivalent for east-facing surfaces and i-index.
834	i = surf_def_v(l)%i(m) + MERGE( surf_def_v(l)%ioff, 0, l == 3 )
835	j = surf_def_v(l)%j(m) + MERGE( surf_def_v(l)%joff, 0, l == 1 )
836	k = surf_def_v(l)%k(m) + surf_def_v(l)%koff
837	!
838	!-- Lower left /front vertex
839	IF ( point_index(k,j,i) < 0 ) THEN
840	surfaces%npoints = surfaces%npoints + 1
841	point_index(k,j,i) = point_index_count
842	point_index_count = point_index_count + 1
843	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
844	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
845	surfaces%points(3,surfaces%npoints) = zw(k-1)
846	ENDIF
847	!
848	!-- Upper / lower right index for north- and south-facing surfaces
849	IF ( l == 0 .OR. l == 1 ) THEN
850	IF ( point_index(k,j,i+1) < 0 ) THEN
851	surfaces%npoints = surfaces%npoints + 1
852	point_index(k,j,i+1) = point_index_count
853	point_index_count = point_index_count + 1
854	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
855	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
856	surfaces%points(3,surfaces%npoints) = zw(k-1)
857	ENDIF
858	IF ( point_index(k+1,j,i+1) < 0 ) THEN
859	surfaces%npoints = surfaces%npoints + 1
860	point_index(k+1,j,i+1) = point_index_count
861	point_index_count = point_index_count + 1
862	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
863	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
864	surfaces%points(3,surfaces%npoints) = zw(k)
865	ENDIF
866	!
867	!-- Upper / lower front index for east- and west-facing surfaces
868	ELSEIF ( l == 2 .OR. l == 3 ) THEN
869	IF ( point_index(k,j+1,i) < 0 ) THEN
870	surfaces%npoints = surfaces%npoints + 1
871	point_index(k,j+1,i) = point_index_count
872	point_index_count = point_index_count + 1
873	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
874	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
875	surfaces%points(3,surfaces%npoints) = zw(k-1)
876	ENDIF
877	IF ( point_index(k+1,j+1,i) < 0 ) THEN
878	surfaces%npoints = surfaces%npoints + 1
879	point_index(k+1,j+1,i) = point_index_count
880	point_index_count = point_index_count + 1
881	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
882	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
883	surfaces%points(3,surfaces%npoints) = zw(k)
884	ENDIF
885	ENDIF
886	!
887	!-- Upper left / front vertex
888	IF ( point_index(k+1,j,i) < 0 ) THEN
889	surfaces%npoints = surfaces%npoints + 1
890	point_index(k+1,j,i) = point_index_count
891	point_index_count = point_index_count + 1
892	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
893	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
894	surfaces%points(3,surfaces%npoints) = zw(k)
895	ENDIF
896
897	npg = npg + 1
898	IF ( l == 0 .OR. l == 1 ) THEN
899	surfaces%polygons(1,npg) = 4
900	surfaces%polygons(2,npg) = point_index(k,j,i)
901	surfaces%polygons(3,npg) = point_index(k,j,i+1)
902	surfaces%polygons(4,npg) = point_index(k+1,j,i+1)
903	surfaces%polygons(5,npg) = point_index(k+1,j,i)
904	ELSE
905	surfaces%polygons(1,npg) = 4
906	surfaces%polygons(2,npg) = point_index(k,j,i)
907	surfaces%polygons(3,npg) = point_index(k,j+1,i)
908	surfaces%polygons(4,npg) = point_index(k+1,j+1,i)
909	surfaces%polygons(5,npg) = point_index(k+1,j,i)
910	ENDIF
911
912	ENDDO
913
914	DO m = 1, surf_lsm_v(l)%ns
915	i = surf_lsm_v(l)%i(m) + MERGE( surf_lsm_v(l)%ioff, 0, l == 3 )
916	j = surf_lsm_v(l)%j(m) + MERGE( surf_lsm_v(l)%joff, 0, l == 1 )
917	k = surf_lsm_v(l)%k(m) + surf_lsm_v(l)%koff
918	!
919	!-- Lower left /front vertex
920	IF ( point_index(k,j,i) < 0 ) THEN
921	surfaces%npoints = surfaces%npoints + 1
922	point_index(k,j,i) = point_index_count
923	point_index_count = point_index_count + 1
924	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
925	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
926	surfaces%points(3,surfaces%npoints) = zw(k-1)
927	ENDIF
928	!
929	!-- Upper / lower right index for north- and south-facing surfaces
930	IF ( l == 0 .OR. l == 1 ) THEN
931	IF ( point_index(k,j,i+1) < 0 ) THEN
932	surfaces%npoints = surfaces%npoints + 1
933	point_index(k,j,i+1) = point_index_count
934	point_index_count = point_index_count + 1
935	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
936	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
937	surfaces%points(3,surfaces%npoints) = zw(k-1)
938	ENDIF
939	IF ( point_index(k+1,j,i+1) < 0 ) THEN
940	surfaces%npoints = surfaces%npoints + 1
941	point_index(k+1,j,i+1) = point_index_count
942	point_index_count = point_index_count + 1
943	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
944	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
945	surfaces%points(3,surfaces%npoints) = zw(k)
946	ENDIF
947	!
948	!-- Upper / lower front index for east- and west-facing surfaces
949	ELSEIF ( l == 2 .OR. l == 3 ) THEN
950	IF ( point_index(k,j+1,i) < 0 ) THEN
951	surfaces%npoints = surfaces%npoints + 1
952	point_index(k,j+1,i) = point_index_count
953	point_index_count = point_index_count + 1
954	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
955	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
956	surfaces%points(3,surfaces%npoints) = zw(k-1)
957	ENDIF
958	IF ( point_index(k+1,j+1,i) < 0 ) THEN
959	surfaces%npoints = surfaces%npoints + 1
960	point_index(k+1,j+1,i) = point_index_count
961	point_index_count = point_index_count + 1
962	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
963	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
964	surfaces%points(3,surfaces%npoints) = zw(k)
965	ENDIF
966	ENDIF
967	!
968	!-- Upper left / front vertex
969	IF ( point_index(k+1,j,i) < 0 ) THEN
970	surfaces%npoints = surfaces%npoints + 1
971	point_index(k+1,j,i) = point_index_count
972	point_index_count = point_index_count + 1
973	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
974	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
975	surfaces%points(3,surfaces%npoints) = zw(k)
976	ENDIF
977
978	npg = npg + 1
979	IF ( l == 0 .OR. l == 1 ) THEN
980	surfaces%polygons(1,npg) = 4
981	surfaces%polygons(2,npg) = point_index(k,j,i)
982	surfaces%polygons(3,npg) = point_index(k,j,i+1)
983	surfaces%polygons(4,npg) = point_index(k+1,j,i+1)
984	surfaces%polygons(5,npg) = point_index(k+1,j,i)
985	ELSE
986	surfaces%polygons(1,npg) = 4
987	surfaces%polygons(2,npg) = point_index(k,j,i)
988	surfaces%polygons(3,npg) = point_index(k,j+1,i)
989	surfaces%polygons(4,npg) = point_index(k+1,j+1,i)
990	surfaces%polygons(5,npg) = point_index(k+1,j,i)
991	ENDIF
992	ENDDO
993	DO m = 1, surf_usm_v(l)%ns
994	i = surf_usm_v(l)%i(m) + MERGE( surf_usm_v(l)%ioff, 0, l == 3 )
995	j = surf_usm_v(l)%j(m) + MERGE( surf_usm_v(l)%joff, 0, l == 1 )
996	k = surf_usm_v(l)%k(m) + surf_usm_v(l)%koff
997	!
998	!-- Lower left /front vertex
999	IF ( point_index(k,j,i) < 0 ) THEN
1000	surfaces%npoints = surfaces%npoints + 1
1001	point_index(k,j,i) = point_index_count
1002	point_index_count = point_index_count + 1
1003	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
1004	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
1005	surfaces%points(3,surfaces%npoints) = zw(k-1)
1006	ENDIF
1007	!
1008	!-- Upper / lower right index for north- and south-facing surfaces
1009	IF ( l == 0 .OR. l == 1 ) THEN
1010	IF ( point_index(k,j,i+1) < 0 ) THEN
1011	surfaces%npoints = surfaces%npoints + 1
1012	point_index(k,j,i+1) = point_index_count
1013	point_index_count = point_index_count + 1
1014	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
1015	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
1016	surfaces%points(3,surfaces%npoints) = zw(k-1)
1017	ENDIF
1018	IF ( point_index(k+1,j,i+1) < 0 ) THEN
1019	surfaces%npoints = surfaces%npoints + 1
1020	point_index(k+1,j,i+1) = point_index_count
1021	point_index_count = point_index_count + 1
1022	surfaces%points(1,surfaces%npoints) = ( i + 1 - 0.5_wp ) * dx
1023	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
1024	surfaces%points(3,surfaces%npoints) = zw(k)
1025	ENDIF
1026	!
1027	!-- Upper / lower front index for east- and west-facing surfaces
1028	ELSEIF ( l == 2 .OR. l == 3 ) THEN
1029	IF ( point_index(k,j+1,i) < 0 ) THEN
1030	surfaces%npoints = surfaces%npoints + 1
1031	point_index(k,j+1,i) = point_index_count
1032	point_index_count = point_index_count + 1
1033	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
1034	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
1035	surfaces%points(3,surfaces%npoints) = zw(k-1)
1036	ENDIF
1037	IF ( point_index(k+1,j+1,i) < 0 ) THEN
1038	surfaces%npoints = surfaces%npoints + 1
1039	point_index(k+1,j+1,i) = point_index_count
1040	point_index_count = point_index_count + 1
1041	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
1042	surfaces%points(2,surfaces%npoints) = ( j + 1 - 0.5_wp ) * dy
1043	surfaces%points(3,surfaces%npoints) = zw(k)
1044	ENDIF
1045	ENDIF
1046	!
1047	!-- Upper left / front vertex
1048	IF ( point_index(k+1,j,i) < 0 ) THEN
1049	surfaces%npoints = surfaces%npoints + 1
1050	point_index(k+1,j,i) = point_index_count
1051	point_index_count = point_index_count + 1
1052	surfaces%points(1,surfaces%npoints) = ( i - 0.5_wp ) * dx
1053	surfaces%points(2,surfaces%npoints) = ( j - 0.5_wp ) * dy
1054	surfaces%points(3,surfaces%npoints) = zw(k)
1055	ENDIF
1056
1057	npg = npg + 1
1058	IF ( l == 0 .OR. l == 1 ) THEN
1059	surfaces%polygons(1,npg) = 4
1060	surfaces%polygons(2,npg) = point_index(k,j,i)
1061	surfaces%polygons(3,npg) = point_index(k,j,i+1)
1062	surfaces%polygons(4,npg) = point_index(k+1,j,i+1)
1063	surfaces%polygons(5,npg) = point_index(k+1,j,i)
1064	ELSE
1065	surfaces%polygons(1,npg) = 4
1066	surfaces%polygons(2,npg) = point_index(k,j,i)
1067	surfaces%polygons(3,npg) = point_index(k,j+1,i)
1068	surfaces%polygons(4,npg) = point_index(k+1,j+1,i)
1069	surfaces%polygons(5,npg) = point_index(k+1,j,i)
1070	ENDIF
1071	ENDDO
1072
1073	ENDDO
1074	!
1075	!-- Deallocate temporary dummy variable
1076	DEALLOCATE ( point_index )
1077	!
1078	!-- Sum-up total number of vertices on domain. This will be needed for post-processing.
1079	surfaces%npoints_total = 0
1080	#if defined( __parallel )
1081	CALL MPI_ALLREDUCE( surfaces%npoints, surfaces%npoints_total, 1, MPI_INTEGER, MPI_SUM, &
1082	comm2d, ierr )
1083	#else
1084	surfaces%npoints_total = surfaces%npoints
1085	#endif
1086	ENDIF
1087	!
1088	!-- If output to netcdf is enabled, set-up the coordinate arrays that unambiguously describe the
1089	!-- position and orientation of each surface element.
1090	IF ( to_netcdf ) THEN
1091	!
1092	!-- Allocate local coordinate arrays
1093	ALLOCATE( surfaces%s(1:surfaces%ns) )
1094	ALLOCATE( surfaces%xs(1:surfaces%ns) )
1095	ALLOCATE( surfaces%ys(1:surfaces%ns) )
1096	ALLOCATE( surfaces%zs(1:surfaces%ns) )
1097	ALLOCATE( surfaces%azimuth(1:surfaces%ns) )
1098	ALLOCATE( surfaces%zenith(1:surfaces%ns) )
1099	ALLOCATE( surfaces%es_utm(1:surfaces%ns) )
1100	ALLOCATE( surfaces%ns_utm(1:surfaces%ns) )
1101	!
1102	!-- Gather the number of surface on each processor, in order to number the surface elements in
1103	!-- ascending order with respect to the total number of surfaces in the domain.
1104	#if defined( __parallel )
1105	CALL MPI_ALLGATHER( surfaces%ns, 1, MPI_INTEGER, num_surfaces_on_pe, 1, MPI_INTEGER, &
1106	comm2d, ierr )
1107	#else
1108	num_surfaces_on_pe = surfaces%ns
1109	#endif
1110	!
1111	!-- First, however, determine the offset where couting of the surfaces should start (the sum of
1112	!-- surfaces on all PE's with lower MPI rank).
1113	i = 0
1114	start_count = 1
1115	DO WHILE ( i < myid .AND. i <= SIZE( num_surfaces_on_pe ) )
1116	start_count = start_count + num_surfaces_on_pe(i)
1117	i = i + 1
1118	ENDDO
1119	!
1120	!-- Set coordinate arrays. For horizontal surfaces, azimuth angles are not defined (fill value).
1121	!-- Zenith angle is 0 (180) for upward (downward)-facing surfaces.
1122	i = start_count
1123	mm = 1
1124	DO l = 0, 1
1125	DO m = 1, surf_def_h(l)%ns
1126	surfaces%s(mm) = i
1127	surfaces%xs(mm) = ( surf_def_h(l)%i(m) + 0.5_wp ) * dx
1128	surfaces%ys(mm) = ( surf_def_h(l)%j(m) + 0.5_wp ) * dy
1129	surfaces%zs(mm) = zw(surf_def_h(l)%k(m)+surf_def_h(l)%koff)
1130	surfaces%azimuth(mm) = surfaces%fillvalue
1131	surfaces%zenith(mm) = 180.0_wp * l
1132	i = i + 1
1133	mm = mm + 1
1134	ENDDO
1135	DO m = 1, surf_lsm_h(l)%ns
1136	surfaces%s(mm) = i
1137	surfaces%xs(mm) = ( surf_lsm_h(l)%i(m) + 0.5_wp ) * dx
1138	surfaces%ys(mm) = ( surf_lsm_h(l)%j(m) + 0.5_wp ) * dy
1139	surfaces%zs(mm) = zw(surf_lsm_h(l)%k(m)+surf_lsm_h(l)%koff)
1140	surfaces%azimuth(mm) = surfaces%fillvalue
1141	surfaces%zenith(mm) = 180.0_wp * l
1142	i = i + 1
1143	mm = mm + 1
1144	ENDDO
1145	DO m = 1, surf_usm_h(l)%ns
1146	surfaces%s(mm) = i
1147	surfaces%xs(mm) = ( surf_usm_h(l)%i(m) + 0.5_wp ) * dx
1148	surfaces%ys(mm) = ( surf_usm_h(l)%j(m) + 0.5_wp ) * dy
1149	surfaces%zs(mm) = zw(surf_usm_h(l)%k(m)+surf_usm_h(l)%koff)
1150	surfaces%azimuth(mm) = surfaces%fillvalue
1151	surfaces%zenith(mm) = 180.0_wp * l
1152	i = i + 1
1153	mm = mm + 1
1154	ENDDO
1155	ENDDO
1156	!
1157	!-- For vertical surfaces, zenith angles are not defined (fill value).
1158	!-- Azimuth angle: northward (0), eastward (90), southward (180), westward (270).
1159	!-- Note, for vertical surfaces, zenith angles are 90.0_wp.
1160	DO l = 0, 3
1161	IF ( l == 0 ) THEN
1162	az = 0.0_wp
1163	off_x = 0.5_wp
1164	off_y = 0.0_wp
1165	ELSEIF ( l == 1 ) THEN
1166	az = 180.0_wp
1167	off_x = 0.5_wp
1168	off_y = 1.0_wp
1169	ELSEIF ( l == 2 ) THEN
1170	az = 90.0_wp
1171	off_x = 0.0_wp
1172	off_y = 0.5_wp
1173	ELSEIF ( l == 3 ) THEN
1174	az = 270.0_wp
1175	off_x = 1.0_wp
1176	off_y = 0.5_wp
1177	ENDIF
1178
1179	DO m = 1, surf_def_v(l)%ns
1180	surfaces%s(mm) = i
1181	surfaces%xs(mm) = ( surf_def_v(l)%i(m) + off_x ) * dx
1182	surfaces%ys(mm) = ( surf_def_v(l)%j(m) + off_y ) * dy
1183	surfaces%zs(mm) = zu(surf_def_v(l)%k(m))
1184	surfaces%azimuth(mm) = az
1185	surfaces%zenith(mm) = 90.0_wp
1186	i = i + 1
1187	mm = mm + 1
1188	ENDDO
1189	DO m = 1, surf_lsm_v(l)%ns
1190	surfaces%s(mm) = i
1191	surfaces%xs(mm) = ( surf_lsm_v(l)%i(m) + off_x ) * dx
1192	surfaces%ys(mm) = ( surf_lsm_v(l)%j(m) + off_y ) * dy
1193	surfaces%zs(mm) = zu(surf_lsm_v(l)%k(m))
1194	surfaces%azimuth(mm) = az
1195	surfaces%zenith(mm) = 90.0_wp
1196	i = i + 1
1197	mm = mm + 1
1198	ENDDO
1199	DO m = 1, surf_usm_v(l)%ns
1200	surfaces%s(mm) = i
1201	surfaces%xs(mm) = ( surf_usm_v(l)%i(m) + off_x ) * dx
1202	surfaces%ys(mm) = ( surf_usm_v(l)%j(m) + off_y ) * dy
1203	surfaces%zs(mm) = zu(surf_usm_v(l)%k(m))
1204	surfaces%azimuth(mm) = az
1205	surfaces%zenith(mm) = 90.0_wp
1206	i = i + 1
1207	mm = mm + 1
1208	ENDDO
1209	ENDDO
1210	!
1211	!-- Finally, define UTM coordinates, which are the x/y-coordinates plus the origin (lower-left
1212	!-- coordinate of the model domain).
1213	surfaces%es_utm = surfaces%xs + init_model%origin_x
1214	surfaces%ns_utm = surfaces%ys + init_model%origin_y
1215	!
1216	!-- Initialize NetCDF data output. Please note, local start position for the surface elements in
1217	!-- the NetCDF file is surfaces%s(1), while the number of surfaces on the subdomain is given by
1218	!-- surfaces%ns.
1219	#if defined( __netcdf4_parallel )
1220
1221	!
1222	!-- Calculate number of time steps to be output
1223	ntdim_surf(0) = dosurf_time_count(0) + CEILING( ( end_time - MAX( &
1224	MERGE( skip_time_dosurf, skip_time_dosurf + spinup_time, &
1225	data_output_during_spinup ), simulated_time_at_begin ) &
1226	) / dt_dosurf )
1227
1228	ntdim_surf(1) = dosurf_time_count(1) + CEILING( ( end_time - MAX( &
1229	MERGE( skip_time_dosurf_av, skip_time_dosurf_av + spinup_time, &
1230	data_output_during_spinup ), simulated_time_at_begin ) &
1231	) / dt_dosurf_av )
1232
1233	!
1234	!-- Create NetCDF4 files for parallel writing
1235	DO av = 0, 1
1236	!
1237	!-- If there is no instantaneous data (av=0) or averaged data (av=1) requested, do not create
1238	!-- the corresponding NetCDF file
1239	IF ( dosurf_no(av) == 0 ) CYCLE
1240
1241	IF ( av == 0 ) THEN
1242	filename = 'SURFACE_DATA_NETCDF' // TRIM( coupling_char )
1243	ELSE
1244	filename = 'SURFACE_DATA_AV_NETCDF' // TRIM( coupling_char )
1245	ENDIF
1246	!
1247	!-- Open file using netCDF4/HDF5 format, parallel
1248	nc_stat = NF90_CREATE( TRIM(filename), &
1249	IOR( NF90_NOCLOBBER, IOR( NF90_NETCDF4, NF90_MPIIO ) ), &
1250	id_set_surf(av), COMM = comm2d, INFO = MPI_INFO_NULL )
1251	CALL netcdf_handle_error( 'surface_data_output_mod', 5550 )
1252
1253	!- Write some global attributes
1254	IF ( av == 0 ) THEN
1255	CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data', &
1256	TRIM( run_description_header ), 5551 )
1257	time_average_text = ' '
1258	ELSE
1259	CALL netcdf_create_global_atts( id_set_surf(av), 'surface-data_av', &
1260	TRIM( run_description_header ), 5552 )
1261	WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_surf
1262	nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'time_avg', &
1263	TRIM( time_average_text ) )
1264	CALL netcdf_handle_error( 'surface_data_output_mod', 5553 )
1265	ENDIF
1266
1267
1268	!
1269	!-- Define time coordinate for surface data.
1270	!-- For parallel output the time dimension has to be limited (ntdim_surf), otherwise the
1271	!-- performance drops significantly.
1272	CALL netcdf_create_dim( id_set_surf(av), 'time', ntdim_surf(av), id_dim_time_surf(av), &
1273	5554 )
1274
1275	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_time_surf(av) /), &
1276	'time', NF90_DOUBLE, id_var_time_surf(av), &
1277	'seconds since '// TRIM(init_model%origin_time), &
1278	'time', 5555, 5555, 5555 )
1279
1280	CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'standard_name', 'time', &
1281	5556)
1282
1283	CALL netcdf_create_att( id_set_surf(av), id_var_time_surf(av), 'axis', 'T', 5557)
1284	!
1285	!-- Define spatial dimensions and coordinates:
1286	!-- Define index of surface element
1287	CALL netcdf_create_dim( id_set_surf(av), 's', surfaces%ns_total, id_dim_s_surf(av), &
1288	5558 )
1289	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 's', NF90_DOUBLE, &
1290	id_var_s_surf(av), '1', 'number of surface element', 5559, &
1291	5559, 5559 )
1292	!
1293	!-- Define x coordinate
1294	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'xs', NF90_DOUBLE, &
1295	id_var_xs_surf(av), 'meters', &
1296	'distance to origin in x-direction', 5561, 5561, 5561 )
1297	!
1298	!-- Define y coordinate
1299	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'ys', NF90_DOUBLE, &
1300	id_var_ys_surf(av), 'meters', &
1301	'distance to origin in y-direction', 5562, 5562, 5562 )
1302	!
1303	!-- Define z coordinate
1304	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'zs', NF90_DOUBLE, &
1305	id_var_zs_surf(av), 'meters', 'height', 5560, 5560, 5560 )
1306	CALL netcdf_create_att( id_set_surf(av), id_var_zs_surf(av), 'standard_name', 'height', &
1307	5583 )
1308
1309	!
1310	!-- Define UTM coordinates
1311	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'Es_UTM', &
1312	NF90_DOUBLE, id_var_etum_surf(av), 'meters', '', 5563, 5563, &
1313	5563 )
1314	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'Ns_UTM', &
1315	NF90_DOUBLE, id_var_nutm_surf(av), 'meters', '', 5564, 5564, &
1316	5564 )
1317
1318	!
1319	!-- Define angles
1320	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'azimuth', &
1321	NF90_DOUBLE, id_var_azimuth_surf(av), 'degree', &
1322	'azimuth angle', 5577, 5578, 5579, fill = .TRUE. )
1323	CALL netcdf_create_att( id_set_surf(av), id_var_azimuth_surf(av), 'standard_name', &
1324	'surface_azimuth_angle', 5584 )
1325
1326	CALL netcdf_create_var( id_set_surf(av), (/ id_dim_s_surf(av) /), 'zenith', &
1327	NF90_DOUBLE, id_var_zenith_surf(av), 'degree', '', 5580, 5581, &
1328	5582, fill = .TRUE. )
1329	!
1330	!-- Define the variables
1331	var_list = ';'
1332	i = 1
1333
1334	DO WHILE ( dosurf(av,i)(1:1) /= ' ' )
1335
1336	CALL netcdf_create_var( id_set_surf(av), &
1337	(/ id_dim_s_surf(av), id_dim_time_surf(av) /), dosurf(av,i), &
1338	NF90_REAL4, id_var_dosurf(av,i), dosurf_unit(av,i), &
1339	dosurf(av,i), 5565, 5565, 5565, .TRUE. )
1340	!
1341	!-- Set no fill for every variable to increase performance.
1342	nc_stat = NF90_DEF_VAR_FILL( id_set_surf(av), id_var_dosurf(av,i), NF90_NOFILL, 0 )
1343	CALL netcdf_handle_error( 'surface_data_output_init', 5566 )
1344	!
1345	!-- Set collective io operations for parallel io
1346	nc_stat = NF90_VAR_PAR_ACCESS( id_set_surf(av), id_var_dosurf(av,i), NF90_COLLECTIVE )
1347	CALL netcdf_handle_error( 'surface_data_output_init', 5567 )
1348	var_list = TRIM( var_list ) // TRIM( dosurf(av,i) ) // ';'
1349
1350	i = i + 1
1351
1352	ENDDO
1353	!
1354	!-- Write the list of variables as global attribute (this is used by restart runs and by
1355	!-- combine_plot_fields)
1356	nc_stat = NF90_PUT_ATT( id_set_surf(av), NF90_GLOBAL, 'VAR_LIST', var_list )
1357	CALL netcdf_handle_error( 'surface_data_output_init', 5568 )
1358
1359	!
1360	!-- Set general no fill, otherwise the performance drops significantly for parallel output.
1361	nc_stat = NF90_SET_FILL( id_set_surf(av), NF90_NOFILL, oldmode )
1362	CALL netcdf_handle_error( 'surface_data_output_init', 5569 )
1363
1364	!
1365	!-- Leave netCDF define mode
1366	nc_stat = NF90_ENDDEF( id_set_surf(av) )
1367	CALL netcdf_handle_error( 'surface_data_output_init', 5570 )
1368
1369	!
1370	!-- These data are only written by PE0 for parallel output to increase the performance.
1371	IF ( myid == 0 ) THEN
1372	!
1373	!-- Write data for surface indices
1374	ALLOCATE( netcdf_data_1d(1:surfaces%ns_total) )
1375
1376	DO i = 1, surfaces%ns_total
1377	netcdf_data_1d(i) = i
1378	ENDDO
1379
1380	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_s_surf(av), netcdf_data_1d, &
1381	start = (/ 1 /), count = (/ surfaces%ns_total /) )
1382	CALL netcdf_handle_error( 'surface_data_output_init', 5571 )
1383
1384	DEALLOCATE( netcdf_data_1d )
1385
1386	ENDIF
1387
1388	!
1389	!-- Write surface positions to file
1390	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_xs_surf(av), &
1391	surfaces%xs, start = (/ surfaces%s(1) /), &
1392	count = (/ surfaces%ns /) )
1393	CALL netcdf_handle_error( 'surface_data_output_init', 5572 )
1394
1395	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_ys_surf(av), surfaces%ys, &
1396	start = (/ surfaces%s(1) /), count = (/ surfaces%ns /) )
1397	CALL netcdf_handle_error( 'surface_data_output_init', 5573 )
1398
1399	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zs_surf(av), surfaces%zs, &
1400	start = (/ surfaces%s(1) /), count = (/ surfaces%ns /) )
1401	CALL netcdf_handle_error( 'surface_data_output_init', 5574 )
1402
1403	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_etum_surf(av), surfaces%es_utm, &
1404	start = (/ surfaces%s(1) /), count = (/ surfaces%ns /) )
1405	CALL netcdf_handle_error( 'surface_data_output_init', 5575 )
1406
1407	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_nutm_surf(av), surfaces%ns_utm, &
1408	start = (/ surfaces%s(1) /), count = (/ surfaces%ns /) )
1409	CALL netcdf_handle_error( 'surface_data_output_init', 5576 )
1410
1411	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_azimuth_surf(av), surfaces%azimuth, &
1412	start = (/ surfaces%s(1) /), count = (/ surfaces%ns /) )
1413	CALL netcdf_handle_error( 'surface_data_output_init', 5585 )
1414
1415	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_zenith_surf(av), surfaces%zenith, &
1416	start = (/ surfaces%s(1) /), count = (/ surfaces%ns /) )
1417	CALL netcdf_handle_error( 'surface_data_output_init', 5586 )
1418
1419	ENDDO
1420	#endif
1421
1422	ENDIF
1423
1424	END SUBROUTINE surface_data_output_init
1425
1426	!--------------------------------------------------------------------------------------------------!
1427	! Description:
1428	! ------------
1429	!> Routine for controlling the data output. Surface data is collected from different types of
1430	!> surfaces (default, natural, urban) and different orientation and written to one 1D-output array.
1431	!> Further, NetCDF routines are called to write the surface data in the respective NetCDF files.
1432	!--------------------------------------------------------------------------------------------------!
1433	SUBROUTINE surface_data_output( av )
1434
1435	USE control_parameters, &
1436	ONLY: io_blocks, &
1437	io_group, &
1438	time_since_reference_point
1439
1440	#if defined( __parallel )
1441	USE pegrid, &
1442	ONLY: comm2d, &
1443	ierr
1444	#endif
1445
1446
1447	IMPLICIT NONE
1448
1449	CHARACTER(LEN=100) :: trimvar = ' ' !< dummy for single output variable
1450
1451	INTEGER(iwp) :: av !< id indicating average or non-average data output
1452	INTEGER(iwp) :: i !< loop index
1453	INTEGER(iwp) :: l !< running index for surface orientation
1454	INTEGER(iwp) :: m !< running index for surface elements
1455	INTEGER(iwp) :: n_out !< counter variables for surface output
1456
1457	!
1458	!-- Return, if nothing to output
1459	IF ( dosurf_no(av) == 0 ) RETURN
1460	!
1461	!-- In case of VTK output, check if binary files are open and write coordinates.
1462	IF ( to_vtk ) THEN
1463
1464	CALL check_open( 25 + av )
1465
1466	IF ( .NOT. first_output(av) ) THEN
1467	DO i = 0, io_blocks - 1
1468	IF ( i == io_group ) THEN
1469	WRITE ( 25 + av ) surfaces%npoints
1470	WRITE ( 25 + av ) surfaces%npoints_total
1471	WRITE ( 25 + av ) surfaces%ns
1472	WRITE ( 25 + av ) surfaces%ns_total
1473	WRITE ( 25 + av ) surfaces%points
1474	WRITE ( 25 + av ) surfaces%polygons
1475	ENDIF
1476	#if defined( __parallel )
1477	CALL MPI_BARRIER( comm2d, ierr )
1478	#endif
1479	first_output(av) = .TRUE.
1480	ENDDO
1481	ENDIF
1482	ENDIF
1483	!
1484	!-- In case of NetCDF output, check if enough time steps are available in file and update time axis.
1485	IF ( to_netcdf ) THEN
1486	#if defined( __netcdf4_parallel )
1487	IF ( dosurf_time_count(av) + 1 > ntdim_surf(av) ) THEN
1488	WRITE ( message_string, * ) 'Output of surface data is not given at t=', &
1489	time_since_reference_point, 's because the maximum ', &
1490	'number of output time levels is exceeded.'
1491	CALL message( 'surface_data_output', 'PA0539', 0, 1, 0, 6, 0 )
1492
1493	RETURN
1494
1495	ENDIF
1496	!
1497	!-- Update the netCDF time axis
1498	!-- In case of parallel output, this is only done by PE0 to increase the performance.
1499	dosurf_time_count(av) = dosurf_time_count(av) + 1
1500	IF ( myid == 0 ) THEN
1501	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_time_surf(av), &
1502	(/ time_since_reference_point /), &
1503	start = (/ dosurf_time_count(av) /), count = (/ 1 /) )
1504	CALL netcdf_handle_error( 'surface_data_output', 6666 )
1505	ENDIF
1506	#endif
1507	ENDIF
1508
1509	!
1510	!-- Cycle through output quantities and write them to file.
1511	n_out = 0
1512	DO WHILE ( dosurf(av,n_out+1)(1:1) /= ' ' )
1513
1514	n_out = n_out + 1
1515	trimvar = TRIM( dosurf(av,n_out) )
1516	!
1517	!-- Set the output array to the _FillValue in case it is not defined for each type of surface.
1518	surfaces%var_out = surfaces%fillvalue
1519	SELECT CASE ( trimvar )
1520
1521	CASE ( 'us' )
1522	!
1523	!-- Output of instantaneous data
1524	IF ( av == 0 ) THEN
1525	CALL surface_data_output_collect( &
1526	surf_def_h(0)%us, surf_lsm_h(0)%us, surf_usm_h(0)%us, &
1527	surf_def_h(1)%us, surf_lsm_h(1)%us, surf_usm_h(1)%us, &
1528	surf_def_v(0)%us, surf_lsm_v(0)%us, surf_usm_v(0)%us, &
1529	surf_def_v(1)%us, surf_lsm_v(1)%us, surf_usm_v(1)%us, &
1530	surf_def_v(2)%us, surf_lsm_v(2)%us, surf_usm_v(2)%us, &
1531	surf_def_v(3)%us, surf_lsm_v(3)%us, surf_usm_v(3)%us )
1532	ELSE
1533	!
1534	!-- Output of averaged data
1535	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1536	surfaces%var_av(:,n_out) = 0.0_wp
1537
1538	ENDIF
1539
1540	CASE ( 'ts' )
1541	!
1542	!-- Output of instantaneous data
1543	IF ( av == 0 ) THEN
1544	CALL surface_data_output_collect( &
1545	surf_def_h(0)%ts, surf_lsm_h(0)%ts, surf_usm_h(0)%ts, &
1546	surf_def_h(1)%ts, surf_lsm_h(1)%ts, surf_usm_h(1)%ts, &
1547	surf_def_v(0)%ts, surf_lsm_v(0)%ts, surf_usm_v(0)%ts, &
1548	surf_def_v(1)%ts, surf_lsm_v(1)%ts, surf_usm_v(1)%ts, &
1549	surf_def_v(2)%ts, surf_lsm_v(2)%ts, surf_usm_v(2)%ts, &
1550	surf_def_v(3)%ts, surf_lsm_v(3)%ts, surf_usm_v(3)%ts )
1551	ELSE
1552	!
1553	!-- Output of averaged data
1554	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1555	surfaces%var_av(:,n_out) = 0.0_wp
1556
1557	ENDIF
1558
1559	CASE ( 'qs' )
1560	!
1561	!-- Output of instantaneous data
1562	IF ( av == 0 ) THEN
1563	CALL surface_data_output_collect( &
1564	surf_def_h(0)%qs, surf_lsm_h(0)%qs, surf_usm_h(0)%qs, &
1565	surf_def_h(1)%qs, surf_lsm_h(1)%qs, surf_usm_h(1)%qs, &
1566	surf_def_v(0)%qs, surf_lsm_v(0)%qs, surf_usm_v(0)%qs, &
1567	surf_def_v(1)%qs, surf_lsm_v(1)%qs, surf_usm_v(1)%qs, &
1568	surf_def_v(2)%qs, surf_lsm_v(2)%qs, surf_usm_v(2)%qs, &
1569	surf_def_v(3)%qs, surf_lsm_v(3)%qs, surf_usm_v(3)%qs )
1570	ELSE
1571	!
1572	!-- Output of averaged data
1573	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1574	surfaces%var_av(:,n_out) = 0.0_wp
1575
1576	ENDIF
1577
1578	CASE ( 'ss' )
1579	!
1580	!-- Output of instantaneous data
1581	IF ( av == 0 ) THEN
1582	CALL surface_data_output_collect( &
1583	surf_def_h(0)%ss, surf_lsm_h(0)%ss, surf_usm_h(0)%ss, &
1584	surf_def_h(1)%ss, surf_lsm_h(1)%ss, surf_usm_h(1)%ss, &
1585	surf_def_v(0)%ss, surf_lsm_v(0)%ss, surf_usm_v(0)%ss, &
1586	surf_def_v(1)%ss, surf_lsm_v(1)%ss, surf_usm_v(1)%ss, &
1587	surf_def_v(2)%ss, surf_lsm_v(2)%ss, surf_usm_v(2)%ss, &
1588	surf_def_v(3)%ss, surf_lsm_v(3)%ss, surf_usm_v(3)%ss )
1589	ELSE
1590	!
1591	!-- Output of averaged data
1592	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1593	surfaces%var_av(:,n_out) = 0.0_wp
1594
1595	ENDIF
1596
1597	CASE ( 'qcs' )
1598	!
1599	!-- Output of instantaneous data
1600	IF ( av == 0 ) THEN
1601	CALL surface_data_output_collect( &
1602	surf_def_h(0)%qcs, surf_lsm_h(0)%qcs, surf_usm_h(0)%qcs, &
1603	surf_def_h(1)%qcs, surf_lsm_h(1)%qcs, surf_usm_h(1)%qcs, &
1604	surf_def_v(0)%qcs, surf_lsm_v(0)%qcs, surf_usm_v(0)%qcs, &
1605	surf_def_v(1)%qcs, surf_lsm_v(1)%qcs, surf_usm_v(1)%qcs, &
1606	surf_def_v(2)%qcs, surf_lsm_v(2)%qcs, surf_usm_v(2)%qcs, &
1607	surf_def_v(3)%qcs, surf_lsm_v(3)%qcs, surf_usm_v(3)%qcs )
1608	ELSE
1609	!
1610	!-- Output of averaged data
1611	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1612	surfaces%var_av(:,n_out) = 0.0_wp
1613
1614	ENDIF
1615
1616	CASE ( 'ncs' )
1617	!
1618	!-- Output of instantaneous data
1619	IF ( av == 0 ) THEN
1620	CALL surface_data_output_collect( &
1621	surf_def_h(0)%ncs, surf_lsm_h(0)%ncs, surf_usm_h(0)%ncs, &
1622	surf_def_h(1)%ncs, surf_lsm_h(1)%ncs, surf_usm_h(1)%ncs, &
1623	surf_def_v(0)%ncs, surf_lsm_v(0)%ncs, surf_usm_v(0)%ncs, &
1624	surf_def_v(1)%ncs, surf_lsm_v(1)%ncs, surf_usm_v(1)%ncs, &
1625	surf_def_v(2)%ncs, surf_lsm_v(2)%ncs, surf_usm_v(2)%ncs, &
1626	surf_def_v(3)%ncs, surf_lsm_v(3)%ncs, surf_usm_v(3)%ncs )
1627	ELSE
1628	!
1629	!-- Output of averaged data
1630	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1631	surfaces%var_av(:,n_out) = 0.0_wp
1632
1633	ENDIF
1634
1635	CASE ( 'qis' )
1636	!
1637	!-- Output of instantaneous data
1638	IF ( av == 0 ) THEN
1639	CALL surface_data_output_collect( &
1640	surf_def_h(0)%qis, surf_lsm_h(0)%qis, surf_usm_h(0)%qis, &
1641	surf_def_h(1)%qis, surf_lsm_h(1)%qis, surf_usm_h(1)%qis, &
1642	surf_def_v(0)%qis, surf_lsm_v(0)%qis, surf_usm_v(0)%qis, &
1643	surf_def_v(1)%qis, surf_lsm_v(1)%qis, surf_usm_v(1)%qis, &
1644	surf_def_v(2)%qis, surf_lsm_v(2)%qis, surf_usm_v(2)%qis, &
1645	surf_def_v(3)%qis, surf_lsm_v(3)%qis, surf_usm_v(3)%qis )
1646	ELSE
1647	!
1648	!-- Output of averaged data
1649	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1650	surfaces%var_av(:,n_out) = 0.0_wp
1651
1652	ENDIF
1653
1654	CASE ( 'nis' )
1655	!
1656	!-- Output of instantaneous data
1657	IF ( av == 0 ) THEN
1658	CALL surface_data_output_collect( &
1659	surf_def_h(0)%nis, surf_lsm_h(0)%nis, surf_usm_h(0)%nis, &
1660	surf_def_h(1)%nis, surf_lsm_h(1)%nis, surf_usm_h(1)%nis, &
1661	surf_def_v(0)%nis, surf_lsm_v(0)%nis, surf_usm_v(0)%nis, &
1662	surf_def_v(1)%nis, surf_lsm_v(1)%nis, surf_usm_v(1)%nis, &
1663	surf_def_v(2)%nis, surf_lsm_v(2)%nis, surf_usm_v(2)%nis, &
1664	surf_def_v(3)%nis, surf_lsm_v(3)%nis, surf_usm_v(3)%nis )
1665	ELSE
1666	!
1667	!-- Output of averaged data
1668	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1669	surfaces%var_av(:,n_out) = 0.0_wp
1670
1671	ENDIF
1672
1673	CASE ( 'qrs' )
1674	!
1675	!-- Output of instantaneous data
1676	IF ( av == 0 ) THEN
1677	CALL surface_data_output_collect( &
1678	surf_def_h(0)%qrs, surf_lsm_h(0)%qrs, surf_usm_h(0)%qrs, &
1679	surf_def_h(1)%qrs, surf_lsm_h(1)%qrs, surf_usm_h(1)%qrs, &
1680	surf_def_v(0)%qrs, surf_lsm_v(0)%qrs, surf_usm_v(0)%qrs, &
1681	surf_def_v(1)%qrs, surf_lsm_v(1)%qrs, surf_usm_v(1)%qrs, &
1682	surf_def_v(2)%qrs, surf_lsm_v(2)%qrs, surf_usm_v(2)%qrs, &
1683	surf_def_v(3)%qrs, surf_lsm_v(3)%qrs, surf_usm_v(3)%qrs )
1684	ELSE
1685	!
1686	!-- Output of averaged data
1687	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1688	surfaces%var_av(:,n_out) = 0.0_wp
1689
1690	ENDIF
1691
1692	CASE ( 'nrs' )
1693	!
1694	!-- Output of instantaneous data
1695	IF ( av == 0 ) THEN
1696	CALL surface_data_output_collect( &
1697	surf_def_h(0)%nrs, surf_lsm_h(0)%nrs, surf_usm_h(0)%nrs, &
1698	surf_def_h(1)%nrs, surf_lsm_h(1)%nrs, surf_usm_h(1)%nrs, &
1699	surf_def_v(0)%nrs, surf_lsm_v(0)%nrs, surf_usm_v(0)%nrs, &
1700	surf_def_v(1)%nrs, surf_lsm_v(1)%nrs, surf_usm_v(1)%nrs, &
1701	surf_def_v(2)%nrs, surf_lsm_v(2)%nrs, surf_usm_v(2)%nrs, &
1702	surf_def_v(3)%nrs, surf_lsm_v(3)%nrs, surf_usm_v(3)%nrs )
1703	ELSE
1704	!
1705	!-- Output of averaged data
1706	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1707	surfaces%var_av(:,n_out) = 0.0_wp
1708
1709	ENDIF
1710
1711	CASE ( 'ol' )
1712	!
1713	!-- Output of instantaneous data
1714	IF ( av == 0 ) THEN
1715	CALL surface_data_output_collect( &
1716	surf_def_h(0)%ol, surf_lsm_h(0)%ol, surf_usm_h(0)%ol, &
1717	surf_def_h(1)%ol, surf_lsm_h(1)%ol, surf_usm_h(1)%ol, &
1718	surf_def_v(0)%ol, surf_lsm_v(0)%ol, surf_usm_v(0)%ol, &
1719	surf_def_v(1)%ol, surf_lsm_v(1)%ol, surf_usm_v(1)%ol, &
1720	surf_def_v(2)%ol, surf_lsm_v(2)%ol, surf_usm_v(2)%ol, &
1721	surf_def_v(3)%ol, surf_lsm_v(3)%ol, surf_usm_v(3)%ol )
1722	ELSE
1723	!
1724	!-- Output of averaged data
1725	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1726	surfaces%var_av(:,n_out) = 0.0_wp
1727
1728	ENDIF
1729
1730	CASE ( 'z0' )
1731	!
1732	!-- Output of instantaneous data
1733	IF ( av == 0 ) THEN
1734	CALL surface_data_output_collect( &
1735	surf_def_h(0)%z0, surf_lsm_h(0)%z0, surf_usm_h(0)%z0, &
1736	surf_def_h(1)%z0, surf_lsm_h(1)%z0, surf_usm_h(1)%z0, &
1737	surf_def_v(0)%z0, surf_lsm_v(0)%z0, surf_usm_v(0)%z0, &
1738	surf_def_v(1)%z0, surf_lsm_v(1)%z0, surf_usm_v(1)%z0, &
1739	surf_def_v(2)%z0, surf_lsm_v(2)%z0, surf_usm_v(2)%z0, &
1740	surf_def_v(3)%z0, surf_lsm_v(3)%z0, surf_usm_v(3)%z0 )
1741	ELSE
1742	!
1743	!-- Output of averaged data
1744	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1745	surfaces%var_av(:,n_out) = 0.0_wp
1746
1747	ENDIF
1748
1749	CASE ( 'z0h' )
1750	!
1751	!-- Output of instantaneous data
1752	IF ( av == 0 ) THEN
1753	CALL surface_data_output_collect( &
1754	surf_def_h(0)%z0h, surf_lsm_h(0)%z0h, surf_usm_h(0)%z0h, &
1755	surf_def_h(1)%z0h, surf_lsm_h(1)%z0h, surf_usm_h(1)%z0h, &
1756	surf_def_v(0)%z0h, surf_lsm_v(0)%z0h, surf_usm_v(0)%z0h, &
1757	surf_def_v(1)%z0h, surf_lsm_v(1)%z0h, surf_usm_v(1)%z0h, &
1758	surf_def_v(2)%z0h, surf_lsm_v(2)%z0h, surf_usm_v(2)%z0h, &
1759	surf_def_v(3)%z0h, surf_lsm_v(3)%z0h, surf_usm_v(3)%z0h )
1760	ELSE
1761	!
1762	!-- Output of averaged data
1763	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1764	surfaces%var_av(:,n_out) = 0.0_wp
1765
1766	ENDIF
1767
1768	CASE ( 'z0q' )
1769	!
1770	!-- Output of instantaneous data
1771	IF ( av == 0 ) THEN
1772	CALL surface_data_output_collect( &
1773	surf_def_h(0)%z0q, surf_lsm_h(0)%z0q, surf_usm_h(0)%z0q, &
1774	surf_def_h(1)%z0q, surf_lsm_h(1)%z0q, surf_usm_h(1)%z0q, &
1775	surf_def_v(0)%z0q, surf_lsm_v(0)%z0q, surf_usm_v(0)%z0q, &
1776	surf_def_v(1)%z0q, surf_lsm_v(1)%z0q, surf_usm_v(1)%z0q, &
1777	surf_def_v(2)%z0q, surf_lsm_v(2)%z0q, surf_usm_v(2)%z0q, &
1778	surf_def_v(3)%z0q, surf_lsm_v(3)%z0q, surf_usm_v(3)%z0q )
1779	ELSE
1780	!
1781	!-- Output of averaged data
1782	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1783	surfaces%var_av(:,n_out) = 0.0_wp
1784
1785	ENDIF
1786
1787	CASE ( 'theta1' )
1788	!
1789	!-- Output of instantaneous data
1790	IF ( av == 0 ) THEN
1791	CALL surface_data_output_collect( &
1792	surf_def_h(0)%pt1, surf_lsm_h(0)%pt1, surf_usm_h(0)%pt1, &
1793	surf_def_h(1)%pt1, surf_lsm_h(1)%pt1, surf_usm_h(1)%pt1, &
1794	surf_def_v(0)%pt1, surf_lsm_v(0)%pt1, surf_usm_v(0)%pt1, &
1795	surf_def_v(1)%pt1, surf_lsm_v(1)%pt1, surf_usm_v(1)%pt1, &
1796	surf_def_v(2)%pt1, surf_lsm_v(2)%pt1, surf_usm_v(2)%pt1, &
1797	surf_def_v(3)%pt1, surf_lsm_v(3)%pt1, surf_usm_v(3)%pt1 )
1798	ELSE
1799	!
1800	!-- Output of averaged data
1801	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1802	surfaces%var_av(:,n_out) = 0.0_wp
1803
1804	ENDIF
1805
1806	CASE ( 'qv1' )
1807	!
1808	!-- Output of instantaneous data
1809	IF ( av == 0 ) THEN
1810	CALL surface_data_output_collect( &
1811	surf_def_h(0)%qv1, surf_lsm_h(0)%qv1, surf_usm_h(0)%qv1, &
1812	surf_def_h(1)%qv1, surf_lsm_h(1)%qv1, surf_usm_h(1)%qv1, &
1813	surf_def_v(0)%qv1, surf_lsm_v(0)%qv1, surf_usm_v(0)%qv1, &
1814	surf_def_v(1)%qv1, surf_lsm_v(1)%qv1, surf_usm_v(1)%qv1, &
1815	surf_def_v(2)%qv1, surf_lsm_v(2)%qv1, surf_usm_v(2)%qv1, &
1816	surf_def_v(3)%qv1, surf_lsm_v(3)%qv1, surf_usm_v(3)%qv1 )
1817	ELSE
1818	!
1819	!-- Output of averaged data
1820	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1821	surfaces%var_av(:,n_out) = 0.0_wp
1822
1823	ENDIF
1824
1825	CASE ( 'thetav1' )
1826	!
1827	!-- Output of instantaneous data
1828	IF ( av == 0 ) THEN
1829	CALL surface_data_output_collect( &
1830	surf_def_h(0)%vpt1, surf_lsm_h(0)%vpt1, surf_usm_h(0)%vpt1, &
1831	surf_def_h(1)%vpt1, surf_lsm_h(1)%vpt1, surf_usm_h(1)%vpt1, &
1832	surf_def_v(0)%vpt1, surf_lsm_v(0)%vpt1, surf_usm_v(0)%vpt1, &
1833	surf_def_v(1)%vpt1, surf_lsm_v(1)%vpt1, surf_usm_v(1)%vpt1, &
1834	surf_def_v(2)%vpt1, surf_lsm_v(2)%vpt1, surf_usm_v(2)%vpt1, &
1835	surf_def_v(3)%vpt1, surf_lsm_v(3)%vpt1, surf_usm_v(3)%vpt1 )
1836	ELSE
1837	!
1838	!-- Output of averaged data
1839	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1840	surfaces%var_av(:,n_out) = 0.0_wp
1841
1842	ENDIF
1843
1844	CASE ( 'usws' )
1845	!
1846	!-- Output of instantaneous data
1847	IF ( av == 0 ) THEN
1848	CALL surface_data_output_collect( &
1849	surf_def_h(0)%usws, surf_lsm_h(0)%usws, surf_usm_h(0)%usws, &
1850	surf_def_h(1)%usws, surf_lsm_h(1)%usws, surf_usm_h(1)%usws, &
1851	surf_def_v(0)%usws, surf_lsm_v(0)%usws, surf_usm_v(0)%usws, &
1852	surf_def_v(1)%usws, surf_lsm_v(1)%usws, surf_usm_v(1)%usws, &
1853	surf_def_v(2)%usws, surf_lsm_v(2)%usws, surf_usm_v(2)%usws, &
1854	surf_def_v(3)%usws, surf_lsm_v(3)%usws, surf_usm_v(3)%usws, &
1855	momentumflux_output_conversion )
1856	ELSE
1857	!
1858	!-- Output of averaged data
1859	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1860	surfaces%var_av(:,n_out) = 0.0_wp
1861
1862	ENDIF
1863
1864	CASE ( 'vsws' )
1865	!
1866	!-- Output of instantaneous data
1867	IF ( av == 0 ) THEN
1868	CALL surface_data_output_collect( &
1869	surf_def_h(0)%vsws, surf_lsm_h(0)%vsws, surf_usm_h(0)%vsws, &
1870	surf_def_h(1)%vsws, surf_lsm_h(1)%vsws, surf_usm_h(1)%vsws, &
1871	surf_def_v(0)%vsws, surf_lsm_v(0)%vsws, surf_usm_v(0)%vsws, &
1872	surf_def_v(1)%vsws, surf_lsm_v(1)%vsws, surf_usm_v(1)%vsws, &
1873	surf_def_v(2)%vsws, surf_lsm_v(2)%vsws, surf_usm_v(2)%vsws, &
1874	surf_def_v(3)%vsws, surf_lsm_v(3)%vsws, surf_usm_v(3)%vsws, &
1875	momentumflux_output_conversion )
1876	ELSE
1877	!
1878	!-- Output of averaged data
1879	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1880	surfaces%var_av(:,n_out) = 0.0_wp
1881
1882	ENDIF
1883
1884	CASE ( 'shf' )
1885	!
1886	!-- Output of instantaneous data
1887	IF ( av == 0 ) THEN
1888	CALL surface_data_output_collect( &
1889	surf_def_h(0)%shf, surf_lsm_h(0)%shf, surf_usm_h(0)%shf, &
1890	surf_def_h(1)%shf, surf_lsm_h(1)%shf, surf_usm_h(1)%shf, &
1891	surf_def_v(0)%shf, surf_lsm_v(0)%shf, surf_usm_v(0)%shf, &
1892	surf_def_v(1)%shf, surf_lsm_v(1)%shf, surf_usm_v(1)%shf, &
1893	surf_def_v(2)%shf, surf_lsm_v(2)%shf, surf_usm_v(2)%shf, &
1894	surf_def_v(3)%shf, surf_lsm_v(3)%shf, surf_usm_v(3)%shf, &
1895	heatflux_output_conversion )
1896	ELSE
1897	!
1898	!-- Output of averaged data
1899	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1900	surfaces%var_av(:,n_out) = 0.0_wp
1901	ENDIF
1902
1903	CASE ( 'qsws' )
1904	!
1905	!-- Output of instantaneous data
1906	IF ( av == 0 ) THEN
1907	CALL surface_data_output_collect( &
1908	surf_def_h(0)%qsws, surf_lsm_h(0)%qsws, surf_usm_h(0)%qsws, &
1909	surf_def_h(1)%qsws, surf_lsm_h(1)%qsws, surf_usm_h(1)%qsws, &
1910	surf_def_v(0)%qsws, surf_lsm_v(0)%qsws, surf_usm_v(0)%qsws, &
1911	surf_def_v(1)%qsws, surf_lsm_v(1)%qsws, surf_usm_v(1)%qsws, &
1912	surf_def_v(2)%qsws, surf_lsm_v(2)%qsws, surf_usm_v(2)%qsws, &
1913	surf_def_v(3)%qsws, surf_lsm_v(3)%qsws, surf_usm_v(3)%qsws, &
1914	waterflux_output_conversion )
1915	ELSE
1916	!
1917	!-- Output of averaged data
1918	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1919	surfaces%var_av(:,n_out) = 0.0_wp
1920
1921	ENDIF
1922
1923	CASE ( 'ssws' )
1924	!
1925	!-- Output of instantaneous data
1926	IF ( av == 0 ) THEN
1927	CALL surface_data_output_collect( &
1928	surf_def_h(0)%ssws, surf_lsm_h(0)%ssws, surf_usm_h(0)%ssws, &
1929	surf_def_h(1)%ssws, surf_lsm_h(1)%ssws, surf_usm_h(1)%ssws, &
1930	surf_def_v(0)%ssws, surf_lsm_v(0)%ssws, surf_usm_v(0)%ssws, &
1931	surf_def_v(1)%ssws, surf_lsm_v(1)%ssws, surf_usm_v(1)%ssws, &
1932	surf_def_v(2)%ssws, surf_lsm_v(2)%ssws, surf_usm_v(2)%ssws, &
1933	surf_def_v(3)%ssws, surf_lsm_v(3)%ssws, surf_usm_v(3)%ssws )
1934	ELSE
1935	!
1936	!-- Output of averaged data
1937	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1938	surfaces%var_av(:,n_out) = 0.0_wp
1939
1940	ENDIF
1941
1942	CASE ( 'qcsws' )
1943	!
1944	!-- Output of instantaneous data
1945	IF ( av == 0 ) THEN
1946	CALL surface_data_output_collect( &
1947	surf_def_h(0)%qcsws, surf_lsm_h(0)%qcsws, surf_usm_h(0)%qcsws, &
1948	surf_def_h(1)%qcsws, surf_lsm_h(1)%qcsws, surf_usm_h(1)%qcsws, &
1949	surf_def_v(0)%qcsws, surf_lsm_v(0)%qcsws, surf_usm_v(0)%qcsws, &
1950	surf_def_v(1)%qcsws, surf_lsm_v(1)%qcsws, surf_usm_v(1)%qcsws, &
1951	surf_def_v(2)%qcsws, surf_lsm_v(2)%qcsws, surf_usm_v(2)%qcsws, &
1952	surf_def_v(3)%qcsws, surf_lsm_v(3)%qcsws, surf_usm_v(3)%qcsws )
1953	ELSE
1954	!
1955	!-- Output of averaged data
1956	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1957	surfaces%var_av(:,n_out) = 0.0_wp
1958
1959	ENDIF
1960
1961	CASE ( 'ncsws' )
1962	!
1963	!-- Output of instantaneous data
1964	IF ( av == 0 ) THEN
1965	CALL surface_data_output_collect( &
1966	surf_def_h(0)%ncsws, surf_lsm_h(0)%ncsws, surf_usm_h(0)%ncsws, &
1967	surf_def_h(1)%ncsws, surf_lsm_h(1)%ncsws, surf_usm_h(1)%ncsws, &
1968	surf_def_v(0)%ncsws, surf_lsm_v(0)%ncsws, surf_usm_v(0)%ncsws, &
1969	surf_def_v(1)%ncsws, surf_lsm_v(1)%ncsws, surf_usm_v(1)%ncsws, &
1970	surf_def_v(2)%ncsws, surf_lsm_v(2)%ncsws, surf_usm_v(2)%ncsws, &
1971	surf_def_v(3)%ncsws, surf_lsm_v(3)%ncsws, surf_usm_v(3)%ncsws )
1972	ELSE
1973	!
1974	!-- Output of averaged data
1975	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1976	surfaces%var_av(:,n_out) = 0.0_wp
1977
1978	ENDIF
1979
1980
1981	CASE ( 'qisws' )
1982	!
1983	!-- Output of instantaneous data
1984	IF ( av == 0 ) THEN
1985	CALL surface_data_output_collect( &
1986	surf_def_h(0)%qisws, surf_lsm_h(0)%qisws, surf_usm_h(0)%qisws, &
1987	surf_def_h(1)%qisws, surf_lsm_h(1)%qisws, surf_usm_h(1)%qisws, &
1988	surf_def_v(0)%qisws, surf_lsm_v(0)%qisws, surf_usm_v(0)%qisws, &
1989	surf_def_v(1)%qisws, surf_lsm_v(1)%qisws, surf_usm_v(1)%qisws, &
1990	surf_def_v(2)%qisws, surf_lsm_v(2)%qisws, surf_usm_v(2)%qisws, &
1991	surf_def_v(3)%qisws, surf_lsm_v(3)%qisws, surf_usm_v(3)%qisws )
1992	ELSE
1993	!
1994	!-- Output of averaged data
1995	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
1996	surfaces%var_av(:,n_out) = 0.0_wp
1997
1998	ENDIF
1999
2000	CASE ( 'nisws' )
2001	!
2002	!-- Output of instantaneous data
2003	IF ( av == 0 ) THEN
2004	CALL surface_data_output_collect( &
2005	surf_def_h(0)%nisws, surf_lsm_h(0)%nisws, surf_usm_h(0)%nisws, &
2006	surf_def_h(1)%nisws, surf_lsm_h(1)%nisws, surf_usm_h(1)%nisws, &
2007	surf_def_v(0)%nisws, surf_lsm_v(0)%nisws, surf_usm_v(0)%nisws, &
2008	surf_def_v(1)%nisws, surf_lsm_v(1)%nisws, surf_usm_v(1)%nisws, &
2009	surf_def_v(2)%nisws, surf_lsm_v(2)%nisws, surf_usm_v(2)%nisws, &
2010	surf_def_v(3)%nisws, surf_lsm_v(3)%nisws, surf_usm_v(3)%nisws )
2011	ELSE
2012	!
2013	!-- Output of averaged data
2014	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2015	surfaces%var_av(:,n_out) = 0.0_wp
2016
2017	ENDIF
2018
2019	CASE ( 'qrsws' )
2020	!
2021	!-- Output of instantaneous data
2022	IF ( av == 0 ) THEN
2023	CALL surface_data_output_collect( &
2024	surf_def_h(0)%qrsws, surf_lsm_h(0)%qrsws, surf_usm_h(0)%qrsws, &
2025	surf_def_h(1)%qrsws, surf_lsm_h(1)%qrsws, surf_usm_h(1)%qrsws, &
2026	surf_def_v(0)%qrsws, surf_lsm_v(0)%qrsws, surf_usm_v(0)%qrsws, &
2027	surf_def_v(1)%qrsws, surf_lsm_v(1)%qrsws, surf_usm_v(1)%qrsws, &
2028	surf_def_v(2)%qrsws, surf_lsm_v(2)%qrsws, surf_usm_v(2)%qrsws, &
2029	surf_def_v(3)%qrsws, surf_lsm_v(3)%qrsws, surf_usm_v(3)%qrsws )
2030	ELSE
2031	!
2032	!-- Output of averaged data
2033	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2034	surfaces%var_av(:,n_out) = 0.0_wp
2035
2036	ENDIF
2037
2038	CASE ( 'nrsws' )
2039	!
2040	!-- Output of instantaneous data
2041	IF ( av == 0 ) THEN
2042	CALL surface_data_output_collect( &
2043	surf_def_h(0)%nrsws, surf_lsm_h(0)%nrsws, surf_usm_h(0)%nrsws, &
2044	surf_def_h(1)%nrsws, surf_lsm_h(1)%nrsws, surf_usm_h(1)%nrsws, &
2045	surf_def_v(0)%nrsws, surf_lsm_v(0)%nrsws, surf_usm_v(0)%nrsws, &
2046	surf_def_v(1)%nrsws, surf_lsm_v(1)%nrsws, surf_usm_v(1)%nrsws, &
2047	surf_def_v(2)%nrsws, surf_lsm_v(2)%nrsws, surf_usm_v(2)%nrsws, &
2048	surf_def_v(3)%nrsws, surf_lsm_v(3)%nrsws, surf_usm_v(3)%nrsws )
2049	ELSE
2050	!
2051	!-- Output of averaged data
2052	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2053	surfaces%var_av(:,n_out) = 0.0_wp
2054
2055	ENDIF
2056
2057	CASE ( 'sasws' )
2058	!
2059	!-- Output of instantaneous data
2060	IF ( av == 0 ) THEN
2061	CALL surface_data_output_collect( &
2062	surf_def_h(0)%sasws, surf_lsm_h(0)%sasws, surf_usm_h(0)%sasws, &
2063	surf_def_h(1)%sasws, surf_lsm_h(1)%sasws, surf_usm_h(1)%sasws, &
2064	surf_def_v(0)%sasws, surf_lsm_v(0)%sasws, surf_usm_v(0)%sasws, &
2065	surf_def_v(1)%sasws, surf_lsm_v(1)%sasws, surf_usm_v(1)%sasws, &
2066	surf_def_v(2)%sasws, surf_lsm_v(2)%sasws, surf_usm_v(2)%sasws, &
2067	surf_def_v(3)%sasws, surf_lsm_v(3)%sasws, surf_usm_v(3)%sasws )
2068	ELSE
2069	!
2070	!-- Output of averaged data
2071	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2072	surfaces%var_av(:,n_out) = 0.0_wp
2073
2074	ENDIF
2075
2076	CASE ( 'q_surface' )
2077	!
2078	!-- Output of instantaneous data
2079	IF ( av == 0 ) THEN
2080	CALL surface_data_output_collect( &
2081	surf_def_h(0)%q_surface, surf_lsm_h(0)%q_surface, surf_usm_h(0)%q_surface, &
2082	surf_def_h(1)%q_surface, surf_lsm_h(1)%q_surface, surf_usm_h(1)%q_surface, &
2083	surf_def_v(0)%q_surface, surf_lsm_v(0)%q_surface, surf_usm_v(0)%q_surface, &
2084	surf_def_v(1)%q_surface, surf_lsm_v(1)%q_surface, surf_usm_v(1)%q_surface, &
2085	surf_def_v(2)%q_surface, surf_lsm_v(2)%q_surface, surf_usm_v(2)%q_surface, &
2086	surf_def_v(3)%q_surface, surf_lsm_v(3)%q_surface, surf_usm_v(3)%q_surface )
2087	ELSE
2088	!
2089	!-- Output of averaged data
2090	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2091	surfaces%var_av(:,n_out) = 0.0_wp
2092
2093	ENDIF
2094
2095	CASE ( 'theta_surface' )
2096	!
2097	!-- Output of instantaneous data
2098	IF ( av == 0 ) THEN
2099	CALL surface_data_output_collect( &
2100	surf_def_h(0)%pt_surface, surf_lsm_h(0)%pt_surface, surf_usm_h(0)%pt_surface, &
2101	surf_def_h(1)%pt_surface, surf_lsm_h(1)%pt_surface, surf_usm_h(1)%pt_surface, &
2102	surf_def_v(0)%pt_surface, surf_lsm_v(0)%pt_surface, surf_usm_v(0)%pt_surface, &
2103	surf_def_v(1)%pt_surface, surf_lsm_v(1)%pt_surface, surf_usm_v(1)%pt_surface, &
2104	surf_def_v(2)%pt_surface, surf_lsm_v(2)%pt_surface, surf_usm_v(2)%pt_surface, &
2105	surf_def_v(3)%pt_surface, surf_lsm_v(3)%pt_surface, surf_usm_v(3)%pt_surface )
2106	ELSE
2107	!
2108	!-- Output of averaged data
2109	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2110	surfaces%var_av(:,n_out) = 0.0_wp
2111
2112	ENDIF
2113
2114	CASE ( 'thetav_surface' )
2115	!
2116	!-- Output of instantaneous data
2117	IF ( av == 0 ) THEN
2118	CALL surface_data_output_collect( &
2119	surf_def_h(0)%vpt_surface, surf_lsm_h(0)%vpt_surface, surf_usm_h(0)%vpt_surface, &
2120	surf_def_h(1)%vpt_surface, surf_lsm_h(1)%vpt_surface, surf_usm_h(1)%vpt_surface, &
2121	surf_def_v(0)%vpt_surface, surf_lsm_v(0)%vpt_surface, surf_usm_v(0)%vpt_surface, &
2122	surf_def_v(1)%vpt_surface, surf_lsm_v(1)%vpt_surface, surf_usm_v(1)%vpt_surface, &
2123	surf_def_v(2)%vpt_surface, surf_lsm_v(2)%vpt_surface, surf_usm_v(2)%vpt_surface, &
2124	surf_def_v(3)%vpt_surface, surf_lsm_v(3)%vpt_surface, surf_usm_v(3)%vpt_surface )
2125	ELSE
2126	!
2127	!-- Output of averaged data
2128	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2129	surfaces%var_av(:,n_out) = 0.0_wp
2130
2131	ENDIF
2132
2133	CASE ( 'rad_net' )
2134	!
2135	!-- Output of instantaneous data
2136	IF ( av == 0 ) THEN
2137	CALL surface_data_output_collect( &
2138	surf_def_h(0)%rad_net, surf_lsm_h(0)%rad_net, surf_usm_h(0)%rad_net, &
2139	surf_def_h(1)%rad_net, surf_lsm_h(1)%rad_net, surf_usm_h(1)%rad_net, &
2140	surf_def_v(0)%rad_net, surf_lsm_v(0)%rad_net, surf_usm_v(0)%rad_net, &
2141	surf_def_v(1)%rad_net, surf_lsm_v(1)%rad_net, surf_usm_v(1)%rad_net, &
2142	surf_def_v(2)%rad_net, surf_lsm_v(2)%rad_net, surf_usm_v(2)%rad_net, &
2143	surf_def_v(3)%rad_net, surf_lsm_v(3)%rad_net, surf_usm_v(3)%rad_net )
2144	ELSE
2145	!
2146	!-- Output of averaged data
2147	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2148	surfaces%var_av(:,n_out) = 0.0_wp
2149
2150	ENDIF
2151
2152	CASE ( 'rad_lw_in' )
2153	!
2154	!-- Output of instantaneous data
2155	IF ( av == 0 ) THEN
2156	CALL surface_data_output_collect( &
2157	surf_def_h(0)%rad_lw_in, surf_lsm_h(0)%rad_lw_in, surf_usm_h(0)%rad_lw_in, &
2158	surf_def_h(1)%rad_lw_in, surf_lsm_h(1)%rad_lw_in, surf_usm_h(1)%rad_lw_in, &
2159	surf_def_v(0)%rad_lw_in, surf_lsm_v(0)%rad_lw_in, surf_usm_v(0)%rad_lw_in, &
2160	surf_def_v(1)%rad_lw_in, surf_lsm_v(1)%rad_lw_in, surf_usm_v(1)%rad_lw_in, &
2161	surf_def_v(2)%rad_lw_in, surf_lsm_v(2)%rad_lw_in, surf_usm_v(2)%rad_lw_in, &
2162	surf_def_v(3)%rad_lw_in, surf_lsm_v(3)%rad_lw_in, surf_usm_v(3)%rad_lw_in )
2163	ELSE
2164	!
2165	!-- Output of averaged data
2166	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2167	surfaces%var_av(:,n_out) = 0.0_wp
2168
2169	ENDIF
2170
2171	CASE ( 'rad_lw_out' )
2172	!
2173	!-- Output of instantaneous data
2174	IF ( av == 0 ) THEN
2175	CALL surface_data_output_collect( &
2176	surf_def_h(0)%rad_lw_out, surf_lsm_h(0)%rad_lw_out, surf_usm_h(0)%rad_lw_out, &
2177	surf_def_h(1)%rad_lw_out, surf_lsm_h(1)%rad_lw_out, surf_usm_h(1)%rad_lw_out, &
2178	surf_def_v(0)%rad_lw_out, surf_lsm_v(0)%rad_lw_out, surf_usm_v(0)%rad_lw_out, &
2179	surf_def_v(1)%rad_lw_out, surf_lsm_v(1)%rad_lw_out, surf_usm_v(1)%rad_lw_out, &
2180	surf_def_v(2)%rad_lw_out, surf_lsm_v(2)%rad_lw_out, surf_usm_v(2)%rad_lw_out, &
2181	surf_def_v(3)%rad_lw_out, surf_lsm_v(3)%rad_lw_out, surf_usm_v(3)%rad_lw_out )
2182	ELSE
2183	!
2184	!-- Output of averaged data
2185	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2186	surfaces%var_av(:,n_out) = 0.0_wp
2187
2188	ENDIF
2189
2190	CASE ( 'rad_sw_in' )
2191	!
2192	!-- Output of instantaneous data
2193	IF ( av == 0 ) THEN
2194	CALL surface_data_output_collect( &
2195	surf_def_h(0)%rad_sw_in, surf_lsm_h(0)%rad_sw_in, surf_usm_h(0)%rad_sw_in, &
2196	surf_def_h(1)%rad_sw_in, surf_lsm_h(1)%rad_sw_in, surf_usm_h(1)%rad_sw_in, &
2197	surf_def_v(0)%rad_sw_in, surf_lsm_v(0)%rad_sw_in, surf_usm_v(0)%rad_sw_in, &
2198	surf_def_v(1)%rad_sw_in, surf_lsm_v(1)%rad_sw_in, surf_usm_v(1)%rad_sw_in, &
2199	surf_def_v(2)%rad_sw_in, surf_lsm_v(2)%rad_sw_in, surf_usm_v(2)%rad_sw_in, &
2200	surf_def_v(3)%rad_sw_in, surf_lsm_v(3)%rad_sw_in, surf_usm_v(3)%rad_sw_in )
2201	ELSE
2202	!
2203	!-- Output of averaged data
2204	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2205	surfaces%var_av(:,n_out) = 0.0_wp
2206
2207	ENDIF
2208
2209	CASE ( 'rad_sw_out' )
2210	!
2211	!-- Output of instantaneous data
2212	IF ( av == 0 ) THEN
2213	CALL surface_data_output_collect( &
2214	surf_def_h(0)%rad_sw_out, surf_lsm_h(0)%rad_sw_out, surf_usm_h(0)%rad_sw_out, &
2215	surf_def_h(1)%rad_sw_out, surf_lsm_h(1)%rad_sw_out, surf_usm_h(1)%rad_sw_out, &
2216	surf_def_v(0)%rad_sw_out, surf_lsm_v(0)%rad_sw_out, surf_usm_v(0)%rad_sw_out, &
2217	surf_def_v(1)%rad_sw_out, surf_lsm_v(1)%rad_sw_out, surf_usm_v(1)%rad_sw_out, &
2218	surf_def_v(2)%rad_sw_out, surf_lsm_v(2)%rad_sw_out, surf_usm_v(2)%rad_sw_out, &
2219	surf_def_v(3)%rad_sw_out, surf_lsm_v(3)%rad_sw_out, surf_usm_v(3)%rad_sw_out )
2220	ELSE
2221	!
2222	!-- Output of averaged data
2223	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2224	surfaces%var_av(:,n_out) = 0.0_wp
2225
2226	ENDIF
2227
2228	CASE ( 'ghf' )
2229	!
2230	!-- Output of instantaneous data
2231	IF ( av == 0 ) THEN
2232	!
2233	!-- Sum up ground / wall heat flux. Note, for urban surfaces the wall heat flux is
2234	!-- aggregated from the different green, window and wall tiles.
2235	DO l = 0, 1
2236	DO m = 1, surf_usm_h(l)%ns
2237	surf_usm_h(l)%ghf(m) = surf_usm_h(l)%frac(m,ind_veg_wall) * &
2238	surf_usm_h(l)%wghf_eb(m) + &
2239	surf_usm_h(l)%frac(m,ind_pav_green) * &
2240	surf_usm_h(l)%wghf_eb_green(m) + &
2241	surf_usm_h(l)%frac(m,ind_wat_win) * &
2242	surf_usm_h(l)%wghf_eb_window(m)
2243	ENDDO
2244	ENDDO
2245	DO l = 0, 3
2246	DO m = 1, surf_usm_v(l)%ns
2247	surf_usm_v(l)%ghf(m) = surf_usm_v(l)%frac(m,ind_veg_wall) * &
2248	surf_usm_v(l)%wghf_eb(m) + &
2249	surf_usm_v(l)%frac(m,ind_pav_green) * &
2250	surf_usm_v(l)%wghf_eb_green(m) + &
2251	surf_usm_v(l)%frac(m,ind_wat_win) * &
2252	surf_usm_v(l)%wghf_eb_window(m)
2253	ENDDO
2254	ENDDO
2255
2256	CALL surface_data_output_collect( &
2257	surf_def_h(0)%ghf, surf_lsm_h(0)%ghf, surf_usm_h(0)%ghf, &
2258	surf_def_h(1)%ghf, surf_lsm_h(1)%ghf, surf_usm_h(1)%ghf, &
2259	surf_def_v(0)%ghf, surf_lsm_v(0)%ghf, surf_usm_v(0)%ghf, &
2260	surf_def_v(1)%ghf, surf_lsm_v(1)%ghf, surf_usm_v(1)%ghf, &
2261	surf_def_v(2)%ghf, surf_lsm_v(2)%ghf, surf_usm_v(2)%ghf, &
2262	surf_def_v(3)%ghf, surf_lsm_v(3)%ghf, surf_usm_v(3)%ghf )
2263	ELSE
2264	!
2265	!-- Output of averaged data
2266	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2267	surfaces%var_av(:,n_out) = 0.0_wp
2268
2269	ENDIF
2270
2271	CASE ( 'r_a' )
2272	!
2273	!-- Output of instantaneous data
2274	IF ( av == 0 ) THEN
2275	CALL surface_data_output_collect( &
2276	surf_def_h(0)%r_a, surf_lsm_h(0)%r_a, surf_usm_h(0)%r_a, &
2277	surf_def_h(1)%r_a, surf_lsm_h(1)%r_a, surf_usm_h(1)%r_a, &
2278	surf_def_v(0)%r_a, surf_lsm_v(0)%r_a, surf_usm_v(0)%r_a, &
2279	surf_def_v(1)%r_a, surf_lsm_v(1)%r_a, surf_usm_v(1)%r_a, &
2280	surf_def_v(2)%r_a, surf_lsm_v(2)%r_a, surf_usm_v(2)%r_a, &
2281	surf_def_v(3)%r_a, surf_lsm_v(3)%r_a, surf_usm_v(3)%r_a )
2282	ELSE
2283	!
2284	!-- Output of averaged data
2285	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2286	surfaces%var_av(:,n_out) = 0.0_wp
2287
2288	ENDIF
2289
2290	CASE ( 'r_soil' )
2291	!
2292	!-- Output of instantaneous data
2293	IF ( av == 0 ) THEN
2294	CALL surface_data_output_collect( &
2295	surf_def_h(0)%r_soil, surf_lsm_h(0)%r_soil, surf_usm_h(0)%r_soil, &
2296	surf_def_h(1)%r_soil, surf_lsm_h(1)%r_soil, surf_usm_h(1)%r_soil, &
2297	surf_def_v(0)%r_soil, surf_lsm_v(0)%r_soil, surf_usm_v(0)%r_soil, &
2298	surf_def_v(1)%r_soil, surf_lsm_v(1)%r_soil, surf_usm_v(1)%r_soil, &
2299	surf_def_v(2)%r_soil, surf_lsm_v(2)%r_soil, surf_usm_v(2)%r_soil, &
2300	surf_def_v(3)%r_soil, surf_lsm_v(3)%r_soil, surf_usm_v(3)%r_soil )
2301	ELSE
2302	!
2303	!-- Output of averaged data
2304	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2305	surfaces%var_av(:,n_out) = 0.0_wp
2306
2307	ENDIF
2308
2309	CASE ( 'r_canopy' )
2310	!
2311	!-- Output of instantaneous data
2312	IF ( av == 0 ) THEN
2313	CALL surface_data_output_collect( &
2314	surf_def_h(0)%r_canopy, surf_lsm_h(0)%r_canopy, surf_usm_h(0)%r_canopy, &
2315	surf_def_h(1)%r_canopy, surf_lsm_h(1)%r_canopy, surf_usm_h(1)%r_canopy, &
2316	surf_def_v(0)%r_canopy, surf_lsm_v(0)%r_canopy, surf_usm_v(0)%r_canopy, &
2317	surf_def_v(1)%r_canopy, surf_lsm_v(1)%r_canopy, surf_usm_v(1)%r_canopy, &
2318	surf_def_v(2)%r_canopy, surf_lsm_v(2)%r_canopy, surf_usm_v(2)%r_canopy, &
2319	surf_def_v(3)%r_canopy, surf_lsm_v(3)%r_canopy, surf_usm_v(3)%r_canopy )
2320	ELSE
2321	!
2322	!-- Output of averaged data
2323	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2324	surfaces%var_av(:,n_out) = 0.0_wp
2325
2326	ENDIF
2327
2328	CASE ( 'r_s' )
2329	!
2330	!-- Output of instantaneous data
2331	IF ( av == 0 ) THEN
2332	CALL surface_data_output_collect( &
2333	surf_def_h(0)%r_s, surf_lsm_h(0)%r_s, surf_usm_h(0)%r_s, &
2334	surf_def_h(1)%r_s, surf_lsm_h(1)%r_s, surf_usm_h(1)%r_s, &
2335	surf_def_v(0)%r_s, surf_lsm_v(0)%r_s, surf_usm_v(0)%r_s, &
2336	surf_def_v(1)%r_s, surf_lsm_v(1)%r_s, surf_usm_v(1)%r_s, &
2337	surf_def_v(2)%r_s, surf_lsm_v(2)%r_s, surf_usm_v(2)%r_s, &
2338	surf_def_v(3)%r_s, surf_lsm_v(3)%r_s, surf_usm_v(3)%r_s )
2339	ELSE
2340	!
2341	!-- Output of averaged data
2342	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2343	surfaces%var_av(:,n_out) = 0.0_wp
2344
2345	ENDIF
2346
2347	CASE ( 'rad_sw_dir' )
2348	!
2349	!-- Output of instantaneous data
2350	IF ( av == 0 ) THEN
2351	CALL surface_data_output_collect( &
2352	surf_def_h(0)%rad_sw_dir, surf_lsm_h(0)%rad_sw_dir, surf_usm_h(0)%rad_sw_dir, &
2353	surf_def_h(1)%rad_sw_dir, surf_lsm_h(1)%rad_sw_dir, surf_usm_h(1)%rad_sw_dir, &
2354	surf_def_v(0)%rad_sw_dir, surf_lsm_v(0)%rad_sw_dir, surf_usm_v(0)%rad_sw_dir, &
2355	surf_def_v(1)%rad_sw_dir, surf_lsm_v(1)%rad_sw_dir, surf_usm_v(1)%rad_sw_dir, &
2356	surf_def_v(2)%rad_sw_dir, surf_lsm_v(2)%rad_sw_dir, surf_usm_v(2)%rad_sw_dir, &
2357	surf_def_v(3)%rad_sw_dir, surf_lsm_v(3)%rad_sw_dir, surf_usm_v(3)%rad_sw_dir )
2358	ELSE
2359	!
2360	!-- Output of averaged data
2361	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2362	surfaces%var_av(:,n_out) = 0.0_wp
2363
2364	ENDIF
2365
2366	CASE ( 'rad_sw_dif' )
2367	!
2368	!-- Output of instantaneous data
2369	IF ( av == 0 ) THEN
2370	CALL surface_data_output_collect( &
2371	surf_def_h(0)%rad_sw_dif, surf_lsm_h(0)%rad_sw_dif, surf_usm_h(0)%rad_sw_dif, &
2372	surf_def_h(1)%rad_sw_dif, surf_lsm_h(1)%rad_sw_dif, surf_usm_h(1)%rad_sw_dif, &
2373	surf_def_v(0)%rad_sw_dif, surf_lsm_v(0)%rad_sw_dif, surf_usm_v(0)%rad_sw_dif, &
2374	surf_def_v(1)%rad_sw_dif, surf_lsm_v(1)%rad_sw_dif, surf_usm_v(1)%rad_sw_dif, &
2375	surf_def_v(2)%rad_sw_dif, surf_lsm_v(2)%rad_sw_dif, surf_usm_v(2)%rad_sw_dif, &
2376	surf_def_v(3)%rad_sw_dif, surf_lsm_v(3)%rad_sw_dif, surf_usm_v(3)%rad_sw_dif )
2377	ELSE
2378	!
2379	!-- Output of averaged data
2380	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2381	surfaces%var_av(:,n_out) = 0.0_wp
2382
2383	ENDIF
2384
2385	CASE ( 'rad_sw_ref' )
2386	!
2387	!-- Output of instantaneous data
2388	IF ( av == 0 ) THEN
2389	CALL surface_data_output_collect( &
2390	surf_def_h(0)%rad_sw_ref, surf_lsm_h(0)%rad_sw_ref, surf_usm_h(0)%rad_sw_ref, &
2391	surf_def_h(1)%rad_sw_ref, surf_lsm_h(1)%rad_sw_ref, surf_usm_h(1)%rad_sw_ref, &
2392	surf_def_v(0)%rad_sw_ref, surf_lsm_v(0)%rad_sw_ref, surf_usm_v(0)%rad_sw_ref, &
2393	surf_def_v(1)%rad_sw_ref, surf_lsm_v(1)%rad_sw_ref, surf_usm_v(1)%rad_sw_ref, &
2394	surf_def_v(2)%rad_sw_ref, surf_lsm_v(2)%rad_sw_ref, surf_usm_v(2)%rad_sw_ref, &
2395	surf_def_v(3)%rad_sw_ref, surf_lsm_v(3)%rad_sw_ref, surf_usm_v(3)%rad_sw_ref )
2396	ELSE
2397	!
2398	!-- Output of averaged data
2399	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2400	surfaces%var_av(:,n_out) = 0.0_wp
2401
2402	ENDIF
2403
2404	CASE ( 'rad_sw_res' )
2405	!
2406	!-- Output of instantaneous data
2407	IF ( av == 0 ) THEN
2408	CALL surface_data_output_collect( &
2409	surf_def_h(0)%rad_sw_res, surf_lsm_h(0)%rad_sw_res, surf_usm_h(0)%rad_sw_res, &
2410	surf_def_h(1)%rad_sw_res, surf_lsm_h(1)%rad_sw_res, surf_usm_h(1)%rad_sw_res, &
2411	surf_def_v(0)%rad_sw_res, surf_lsm_v(0)%rad_sw_res, surf_usm_v(0)%rad_sw_res, &
2412	surf_def_v(1)%rad_sw_res, surf_lsm_v(1)%rad_sw_res, surf_usm_v(1)%rad_sw_res, &
2413	surf_def_v(2)%rad_sw_res, surf_lsm_v(2)%rad_sw_res, surf_usm_v(2)%rad_sw_res, &
2414	surf_def_v(3)%rad_sw_res, surf_lsm_v(3)%rad_sw_res, surf_usm_v(3)%rad_sw_res )
2415	ELSE
2416	!
2417	!-- Output of averaged data
2418	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2419	surfaces%var_av(:,n_out) = 0.0_wp
2420
2421	ENDIF
2422
2423	CASE ( 'rad_lw_dif' )
2424	!
2425	!-- Output of instantaneous data
2426	IF ( av == 0 ) THEN
2427	CALL surface_data_output_collect( &
2428	surf_def_h(0)%rad_lw_dif, surf_lsm_h(0)%rad_lw_dif, surf_usm_h(0)%rad_lw_dif, &
2429	surf_def_h(1)%rad_lw_dif, surf_lsm_h(1)%rad_lw_dif, surf_usm_h(1)%rad_lw_dif, &
2430	surf_def_v(0)%rad_lw_dif, surf_lsm_v(0)%rad_lw_dif, surf_usm_v(0)%rad_lw_dif, &
2431	surf_def_v(1)%rad_lw_dif, surf_lsm_v(1)%rad_lw_dif, surf_usm_v(1)%rad_lw_dif, &
2432	surf_def_v(2)%rad_lw_dif, surf_lsm_v(2)%rad_lw_dif, surf_usm_v(2)%rad_lw_dif, &
2433	surf_def_v(3)%rad_lw_dif, surf_lsm_v(3)%rad_lw_dif, surf_usm_v(3)%rad_lw_dif )
2434	ELSE
2435	!
2436	!-- Output of averaged data
2437	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2438	surfaces%var_av(:,n_out) = 0.0_wp
2439
2440	ENDIF
2441
2442	CASE ( 'rad_lw_ref' )
2443	!
2444	!-- Output of instantaneous data
2445	IF ( av == 0 ) THEN
2446	CALL surface_data_output_collect( &
2447	surf_def_h(0)%rad_lw_ref, surf_lsm_h(0)%rad_lw_ref, surf_usm_h(0)%rad_lw_ref, &
2448	surf_def_h(1)%rad_lw_ref, surf_lsm_h(1)%rad_lw_ref, surf_usm_h(1)%rad_lw_ref, &
2449	surf_def_v(0)%rad_lw_ref, surf_lsm_v(0)%rad_lw_ref, surf_usm_v(0)%rad_lw_ref, &
2450	surf_def_v(1)%rad_lw_ref, surf_lsm_v(1)%rad_lw_ref, surf_usm_v(1)%rad_lw_ref, &
2451	surf_def_v(2)%rad_lw_ref, surf_lsm_v(2)%rad_lw_ref, surf_usm_v(2)%rad_lw_ref, &
2452	surf_def_v(3)%rad_lw_ref, surf_lsm_v(3)%rad_lw_ref, surf_usm_v(3)%rad_lw_ref )
2453	ELSE
2454	!
2455	!-- Output of averaged data
2456	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2457	surfaces%var_av(:,n_out) = 0.0_wp
2458
2459	ENDIF
2460
2461	CASE ( 'rad_lw_res' )
2462	!
2463	!-- Output of instantaneous data
2464	IF ( av == 0 ) THEN
2465	CALL surface_data_output_collect( &
2466	surf_def_h(0)%rad_lw_res, surf_lsm_h(0)%rad_lw_res, surf_usm_h(0)%rad_lw_res, &
2467	surf_def_h(1)%rad_lw_res, surf_lsm_h(1)%rad_lw_res, surf_usm_h(1)%rad_lw_res, &
2468	surf_def_v(0)%rad_lw_res, surf_lsm_v(0)%rad_lw_res, surf_usm_v(0)%rad_lw_res, &
2469	surf_def_v(1)%rad_lw_res, surf_lsm_v(1)%rad_lw_res, surf_usm_v(1)%rad_lw_res, &
2470	surf_def_v(2)%rad_lw_res, surf_lsm_v(2)%rad_lw_res, surf_usm_v(2)%rad_lw_res, &
2471	surf_def_v(3)%rad_lw_res, surf_lsm_v(3)%rad_lw_res, surf_usm_v(3)%rad_lw_res )
2472	ELSE
2473	!
2474	!-- Output of averaged data
2475	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2476	surfaces%var_av(:,n_out) = 0.0_wp
2477
2478	ENDIF
2479
2480	CASE ( 'uvw1' )
2481	!
2482	!-- Output of instantaneous data
2483	IF ( av == 0 ) THEN
2484	CALL surface_data_output_collect( &
2485	surf_def_h(0)%uvw_abs, surf_lsm_h(0)%uvw_abs, surf_usm_h(0)%uvw_abs, &
2486	surf_def_h(1)%uvw_abs, surf_lsm_h(1)%uvw_abs, surf_usm_h(1)%uvw_abs, &
2487	surf_def_v(0)%uvw_abs, surf_lsm_v(0)%uvw_abs, surf_usm_v(0)%uvw_abs, &
2488	surf_def_v(1)%uvw_abs, surf_lsm_v(1)%uvw_abs, surf_usm_v(1)%uvw_abs, &
2489	surf_def_v(2)%uvw_abs, surf_lsm_v(2)%uvw_abs, surf_usm_v(2)%uvw_abs, &
2490	surf_def_v(3)%uvw_abs, surf_lsm_v(3)%uvw_abs, surf_usm_v(3)%uvw_abs )
2491	ELSE
2492	!
2493	!-- Output of averaged data
2494	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2495	surfaces%var_av(:,n_out) = 0.0_wp
2496
2497	ENDIF
2498	!
2499	!-- Waste heat from indoor model
2500	CASE ( 'waste_heat' )
2501	!
2502	!-- Output of instantaneous data
2503	IF ( av == 0 ) THEN
2504	CALL surface_data_output_collect( &
2505	surf_def_h(0)%waste_heat, surf_lsm_h(0)%waste_heat, surf_usm_h(0)%waste_heat, &
2506	surf_def_h(1)%waste_heat, surf_lsm_h(1)%waste_heat, surf_usm_h(1)%waste_heat, &
2507	surf_def_v(0)%waste_heat, surf_lsm_v(0)%waste_heat, surf_usm_v(0)%waste_heat, &
2508	surf_def_v(1)%waste_heat, surf_lsm_v(1)%waste_heat, surf_usm_v(1)%waste_heat, &
2509	surf_def_v(2)%waste_heat, surf_lsm_v(2)%waste_heat, surf_usm_v(2)%waste_heat, &
2510	surf_def_v(3)%waste_heat, surf_lsm_v(3)%waste_heat, surf_usm_v(3)%waste_heat )
2511	ELSE
2512	!
2513	!-- Output of averaged data
2514	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2515	surfaces%var_av(:,n_out) = 0.0_wp
2516
2517	ENDIF
2518	!
2519	!-- Innermost building wall flux from indoor model
2520	CASE ( 'im_hf' )
2521	!
2522	!-- Output of instantaneous data
2523	IF ( av == 0 ) THEN
2524	CALL surface_data_output_collect( &
2525	surf_def_h(0)%iwghf_eb, surf_lsm_h(0)%iwghf_eb, surf_usm_h(0)%iwghf_eb, &
2526	surf_def_h(1)%iwghf_eb, surf_lsm_h(1)%iwghf_eb, surf_usm_h(1)%iwghf_eb, &
2527	surf_def_v(0)%iwghf_eb, surf_lsm_v(0)%iwghf_eb, surf_usm_v(0)%iwghf_eb, &
2528	surf_def_v(1)%iwghf_eb, surf_lsm_v(1)%iwghf_eb, surf_usm_v(1)%iwghf_eb, &
2529	surf_def_v(2)%iwghf_eb, surf_lsm_v(2)%iwghf_eb, surf_usm_v(2)%iwghf_eb, &
2530	surf_def_v(3)%iwghf_eb, surf_lsm_v(3)%iwghf_eb, surf_usm_v(3)%iwghf_eb )
2531	ELSE
2532	!
2533	!-- Output of averaged data
2534	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2535	surfaces%var_av(:,n_out) = 0.0_wp
2536
2537	ENDIF
2538	!
2539	!-- Surface albedo (tile approach)
2540	CASE ( 'albedo' )
2541	!
2542	!-- Output of instantaneous data
2543	IF ( av == 0 ) THEN
2544	CALL surface_data_output_collect( &
2545	surf_def_h(0)%albedo, surf_lsm_h(0)%albedo, surf_usm_h(0)%albedo, &
2546	surf_def_h(1)%albedo, surf_lsm_h(1)%albedo, surf_usm_h(1)%albedo, &
2547	surf_def_v(0)%albedo, surf_lsm_v(0)%albedo, surf_usm_v(0)%albedo, &
2548	surf_def_v(1)%albedo, surf_lsm_v(1)%albedo, surf_usm_v(1)%albedo, &
2549	surf_def_v(2)%albedo, surf_lsm_v(2)%albedo, surf_usm_v(2)%albedo, &
2550	surf_def_v(3)%albedo, surf_lsm_v(3)%albedo, surf_usm_v(3)%albedo )
2551	ELSE
2552	!
2553	!-- Output of averaged data
2554	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2555	surfaces%var_av(:,n_out) = 0.0_wp
2556
2557	ENDIF
2558	!
2559	!-- Surface emissivity (tile approach)
2560	CASE ( 'emissivity' )
2561	!
2562	!-- Output of instantaneous data
2563	IF ( av == 0 ) THEN
2564	CALL surface_data_output_collect( &
2565	surf_def_h(0)%emissivity, surf_lsm_h(0)%emissivity, surf_usm_h(0)%emissivity, &
2566	surf_def_h(1)%emissivity, surf_lsm_h(1)%emissivity, surf_usm_h(1)%emissivity, &
2567	surf_def_v(0)%emissivity, surf_lsm_v(0)%emissivity, surf_usm_v(0)%emissivity, &
2568	surf_def_v(1)%emissivity, surf_lsm_v(1)%emissivity, surf_usm_v(1)%emissivity, &
2569	surf_def_v(2)%emissivity, surf_lsm_v(2)%emissivity, surf_usm_v(2)%emissivity, &
2570	surf_def_v(3)%emissivity, surf_lsm_v(3)%emissivity, surf_usm_v(3)%emissivity )
2571	ELSE
2572	!
2573	!-- Output of averaged data
2574	surfaces%var_out(:) = surfaces%var_av(:,n_out) / REAL( average_count_surf, KIND=wp )
2575	surfaces%var_av(:,n_out) = 0.0_wp
2576
2577	ENDIF
2578	!
2579	!-- Add further variables:
2580	!-- 'css', 'cssws', 'qsws_liq', 'qsws_soil', 'qsws_veg'
2581
2582	END SELECT
2583	!
2584	!-- Write to binary file:
2585	!-- - surfaces%points ( 3, 1-npoints )
2586	!-- - surfaces%polygons ( 5, 1-ns )
2587	!-- - surfaces%var_out ( 1-ns, time )
2588	!-- - Dimension: 1-nsurfaces, 1-npoints - can be ordered consecutively
2589	!-- - Distinguish between average and non-average data
2590	IF ( to_vtk ) THEN
2591	DO i = 0, io_blocks - 1
2592	IF ( i == io_group ) THEN
2593	WRITE ( 25 + av ) LEN_TRIM( 'time' )
2594	WRITE ( 25 + av ) 'time'
2595	WRITE ( 25 + av ) time_since_reference_point
2596	WRITE ( 25 + av ) LEN_TRIM( trimvar )
2597	WRITE ( 25 + av ) TRIM( trimvar )
2598	WRITE ( 25 + av ) surfaces%var_out
2599	ENDIF
2600	#if defined( __parallel )
2601	CALL MPI_BARRIER( comm2d, ierr )
2602	#endif
2603	ENDDO
2604	ENDIF
2605
2606	IF ( to_netcdf ) THEN
2607	#if defined( __netcdf4_parallel )
2608	!
2609	!-- Write output array to file
2610	nc_stat = NF90_PUT_VAR( id_set_surf(av), id_var_dosurf(av,n_out), surfaces%var_out, &
2611	start = (/ surfaces%s(1), dosurf_time_count(av) /), &
2612	count = (/ surfaces%ns, 1 /) )
2613	CALL netcdf_handle_error( 'surface_data_output', 6667 )
2614	#endif
2615	ENDIF
2616
2617	ENDDO
2618
2619	!
2620	!-- If averaged output was written to NetCDF file, set the counter to zero
2621	IF ( av == 1 ) average_count_surf = 0
2622
2623	END SUBROUTINE surface_data_output
2624
2625	!--------------------------------------------------------------------------------------------------!
2626	! Description:
2627	! ------------
2628	!> Routine for controlling the data averaging.
2629	!--------------------------------------------------------------------------------------------------!
2630	SUBROUTINE surface_data_output_averaging
2631
2632	IMPLICIT NONE
2633
2634	CHARACTER(LEN=100) :: trimvar !< dummy variable for current output variable
2635
2636	INTEGER(iwp) :: l !< running index for surface orientation
2637	INTEGER(iwp) :: m !< running index for surface elements
2638	INTEGER(iwp) :: n_out !< counter variables for surface output
2639
2640	n_out = 0
2641	DO WHILE ( dosurf(1,n_out+1)(1:1) /= ' ' )
2642
2643	n_out = n_out + 1
2644	trimvar = TRIM( dosurf(1,n_out) )
2645
2646	SELECT CASE ( trimvar )
2647
2648	CASE ( 'us' )
2649	CALL surface_data_output_sum_up( &
2650	surf_def_h(0)%us, surf_lsm_h(0)%us, surf_usm_h(0)%us, &
2651	surf_def_h(1)%us, surf_lsm_h(1)%us, surf_usm_h(1)%us, &
2652	surf_def_v(0)%us, surf_lsm_v(0)%us, surf_usm_v(0)%us, &
2653	surf_def_v(1)%us, surf_lsm_v(1)%us, surf_usm_v(1)%us, &
2654	surf_def_v(2)%us, surf_lsm_v(2)%us, surf_usm_v(2)%us, &
2655	surf_def_v(3)%us, surf_lsm_v(3)%us, surf_usm_v(3)%us, n_out )
2656
2657	CASE ( 'ts' )
2658	CALL surface_data_output_sum_up( &
2659	surf_def_h(0)%ts, surf_lsm_h(0)%ts, surf_usm_h(0)%ts, &
2660	surf_def_h(1)%ts, surf_lsm_h(1)%ts, surf_usm_h(1)%ts, &
2661	surf_def_v(0)%ts, surf_lsm_v(0)%ts, surf_usm_v(0)%ts, &
2662	surf_def_v(1)%ts, surf_lsm_v(1)%ts, surf_usm_v(1)%ts, &
2663	surf_def_v(2)%ts, surf_lsm_v(2)%ts, surf_usm_v(2)%ts, &
2664	surf_def_v(3)%ts, surf_lsm_v(3)%ts, surf_usm_v(3)%ts, n_out )
2665
2666	CASE ( 'qs' )
2667	CALL surface_data_output_sum_up( &
2668	surf_def_h(0)%qs, surf_lsm_h(0)%qs, surf_usm_h(0)%qs, &
2669	surf_def_h(1)%qs, surf_lsm_h(1)%qs, surf_usm_h(1)%qs, &
2670	surf_def_v(0)%qs, surf_lsm_v(0)%qs, surf_usm_v(0)%qs, &
2671	surf_def_v(1)%qs, surf_lsm_v(1)%qs, surf_usm_v(1)%qs, &
2672	surf_def_v(2)%qs, surf_lsm_v(2)%qs, surf_usm_v(2)%qs, &
2673	surf_def_v(3)%qs, surf_lsm_v(3)%qs, surf_usm_v(3)%qs, n_out )
2674
2675	CASE ( 'ss' )
2676	CALL surface_data_output_sum_up( &
2677	surf_def_h(0)%ss, surf_lsm_h(0)%ss, surf_usm_h(0)%ss, &
2678	surf_def_h(1)%ss, surf_lsm_h(1)%ss, surf_usm_h(1)%ss, &
2679	surf_def_v(0)%ss, surf_lsm_v(0)%ss, surf_usm_v(0)%ss, &
2680	surf_def_v(1)%ss, surf_lsm_v(1)%ss, surf_usm_v(1)%ss, &
2681	surf_def_v(2)%ss, surf_lsm_v(2)%ss, surf_usm_v(2)%ss, &
2682	surf_def_v(3)%ss, surf_lsm_v(3)%ss, surf_usm_v(3)%ss, n_out )
2683
2684	CASE ( 'qcs' )
2685	CALL surface_data_output_sum_up( &
2686	surf_def_h(0)%qcs, surf_lsm_h(0)%qcs, surf_usm_h(0)%qcs, &
2687	surf_def_h(1)%qcs, surf_lsm_h(1)%qcs, surf_usm_h(1)%qcs, &
2688	surf_def_v(0)%qcs, surf_lsm_v(0)%qcs, surf_usm_v(0)%qcs, &
2689	surf_def_v(1)%qcs, surf_lsm_v(1)%qcs, surf_usm_v(1)%qcs, &
2690	surf_def_v(2)%qcs, surf_lsm_v(2)%qcs, surf_usm_v(2)%qcs, &
2691	surf_def_v(3)%qcs, surf_lsm_v(3)%qcs, surf_usm_v(3)%qcs, n_out )
2692
2693	CASE ( 'ncs' )
2694	CALL surface_data_output_sum_up( &
2695	surf_def_h(0)%ncs, surf_lsm_h(0)%ncs, surf_usm_h(0)%ncs, &
2696	surf_def_h(1)%ncs, surf_lsm_h(1)%ncs, surf_usm_h(1)%ncs, &
2697	surf_def_v(0)%ncs, surf_lsm_v(0)%ncs, surf_usm_v(0)%ncs, &
2698	surf_def_v(1)%ncs, surf_lsm_v(1)%ncs, surf_usm_v(1)%ncs, &
2699	surf_def_v(2)%ncs, surf_lsm_v(2)%ncs, surf_usm_v(2)%ncs, &
2700	surf_def_v(3)%ncs, surf_lsm_v(3)%ncs, surf_usm_v(3)%ncs, n_out )
2701
2702	CASE ( 'qis' )
2703	CALL surface_data_output_sum_up( &
2704	surf_def_h(0)%qis, surf_lsm_h(0)%qis, surf_usm_h(0)%qis, &
2705	surf_def_h(1)%qis, surf_lsm_h(1)%qis, surf_usm_h(1)%qis, &
2706	surf_def_v(0)%qis, surf_lsm_v(0)%qis, surf_usm_v(0)%qis, &
2707	surf_def_v(1)%qis, surf_lsm_v(1)%qis, surf_usm_v(1)%qis, &
2708	surf_def_v(2)%qis, surf_lsm_v(2)%qis, surf_usm_v(2)%qis, &
2709	surf_def_v(3)%qis, surf_lsm_v(3)%qis, surf_usm_v(3)%qis, n_out )
2710
2711	CASE ( 'nis' )
2712	CALL surface_data_output_sum_up( &
2713	surf_def_h(0)%nis, surf_lsm_h(0)%nis, surf_usm_h(0)%nis, &
2714	surf_def_h(1)%nis, surf_lsm_h(1)%nis, surf_usm_h(1)%nis, &
2715	surf_def_v(0)%nis, surf_lsm_v(0)%nis, surf_usm_v(0)%nis, &
2716	surf_def_v(1)%nis, surf_lsm_v(1)%nis, surf_usm_v(1)%nis, &
2717	surf_def_v(2)%nis, surf_lsm_v(2)%nis, surf_usm_v(2)%nis, &
2718	surf_def_v(3)%nis, surf_lsm_v(3)%nis, surf_usm_v(3)%nis, n_out )
2719
2720	CASE ( 'qrs' )
2721	CALL surface_data_output_sum_up( &
2722	surf_def_h(0)%qrs, surf_lsm_h(0)%qrs, surf_usm_h(0)%qrs, &
2723	surf_def_h(1)%qrs, surf_lsm_h(1)%qrs, surf_usm_h(1)%qrs, &
2724	surf_def_v(0)%qrs, surf_lsm_v(0)%qrs, surf_usm_v(0)%qrs, &
2725	surf_def_v(1)%qrs, surf_lsm_v(1)%qrs, surf_usm_v(1)%qrs, &
2726	surf_def_v(2)%qrs, surf_lsm_v(2)%qrs, surf_usm_v(2)%qrs, &
2727	surf_def_v(3)%qrs, surf_lsm_v(3)%qrs, surf_usm_v(3)%qrs, n_out )
2728
2729	CASE ( 'nrs' )
2730	CALL surface_data_output_sum_up( &
2731	surf_def_h(0)%nrs, surf_lsm_h(0)%nrs, surf_usm_h(0)%nrs, &
2732	surf_def_h(1)%nrs, surf_lsm_h(1)%nrs, surf_usm_h(1)%nrs, &
2733	surf_def_v(0)%nrs, surf_lsm_v(0)%nrs, surf_usm_v(0)%nrs, &
2734	surf_def_v(1)%nrs, surf_lsm_v(1)%nrs, surf_usm_v(1)%nrs, &
2735	surf_def_v(2)%nrs, surf_lsm_v(2)%nrs, surf_usm_v(2)%nrs, &
2736	surf_def_v(3)%nrs, surf_lsm_v(3)%nrs, surf_usm_v(3)%nrs, n_out )
2737
2738	CASE ( 'ol' )
2739	CALL surface_data_output_sum_up( &
2740	surf_def_h(0)%ol, surf_lsm_h(0)%ol, surf_usm_h(0)%ol, &
2741	surf_def_h(1)%ol, surf_lsm_h(1)%ol, surf_usm_h(1)%ol, &
2742	surf_def_v(0)%ol, surf_lsm_v(0)%ol, surf_usm_v(0)%ol, &
2743	surf_def_v(1)%ol, surf_lsm_v(1)%ol, surf_usm_v(1)%ol, &
2744	surf_def_v(2)%ol, surf_lsm_v(2)%ol, surf_usm_v(2)%ol, &
2745	surf_def_v(3)%ol, surf_lsm_v(3)%ol, surf_usm_v(3)%ol, n_out )
2746
2747	CASE ( 'z0' )
2748	CALL surface_data_output_sum_up( &
2749	surf_def_h(0)%z0, surf_lsm_h(0)%z0, surf_usm_h(0)%z0, &
2750	surf_def_h(1)%z0, surf_lsm_h(1)%z0, surf_usm_h(1)%z0, &
2751	surf_def_v(0)%z0, surf_lsm_v(0)%z0, surf_usm_v(0)%z0, &
2752	surf_def_v(1)%z0, surf_lsm_v(1)%z0, surf_usm_v(1)%z0, &
2753	surf_def_v(2)%z0, surf_lsm_v(2)%z0, surf_usm_v(2)%z0, &
2754	surf_def_v(3)%z0, surf_lsm_v(3)%z0, surf_usm_v(3)%z0, n_out )
2755
2756	CASE ( 'z0h' )
2757	CALL surface_data_output_sum_up( &
2758	surf_def_h(0)%z0h, surf_lsm_h(0)%z0h, surf_usm_h(0)%z0h, &
2759	surf_def_h(1)%z0h, surf_lsm_h(1)%z0h, surf_usm_h(1)%z0h, &
2760	surf_def_v(0)%z0h, surf_lsm_v(0)%z0h, surf_usm_v(0)%z0h, &
2761	surf_def_v(1)%z0h, surf_lsm_v(1)%z0h, surf_usm_v(1)%z0h, &
2762	surf_def_v(2)%z0h, surf_lsm_v(2)%z0h, surf_usm_v(2)%z0h, &
2763	surf_def_v(3)%z0h, surf_lsm_v(3)%z0h, surf_usm_v(3)%z0h, n_out )
2764
2765	CASE ( 'z0q' )
2766	CALL surface_data_output_sum_up( &
2767	surf_def_h(0)%z0q, surf_lsm_h(0)%z0q, surf_usm_h(0)%z0q, &
2768	surf_def_h(1)%z0q, surf_lsm_h(1)%z0q, surf_usm_h(1)%z0q, &
2769	surf_def_v(0)%z0q, surf_lsm_v(0)%z0q, surf_usm_v(0)%z0q, &
2770	surf_def_v(1)%z0q, surf_lsm_v(1)%z0q, surf_usm_v(1)%z0q, &
2771	surf_def_v(2)%z0q, surf_lsm_v(2)%z0q, surf_usm_v(2)%z0q, &
2772	surf_def_v(3)%z0q, surf_lsm_v(3)%z0q, surf_usm_v(3)%z0q, n_out )
2773
2774	CASE ( 'theta1' )
2775	CALL surface_data_output_sum_up( &
2776	surf_def_h(0)%pt1, surf_lsm_h(0)%pt1, surf_usm_h(0)%pt1, &
2777	surf_def_h(1)%pt1, surf_lsm_h(1)%pt1, surf_usm_h(1)%pt1, &
2778	surf_def_v(0)%pt1, surf_lsm_v(0)%pt1, surf_usm_v(0)%pt1, &
2779	surf_def_v(1)%pt1, surf_lsm_v(1)%pt1, surf_usm_v(1)%pt1, &
2780	surf_def_v(2)%pt1, surf_lsm_v(2)%pt1, surf_usm_v(2)%pt1, &
2781	surf_def_v(3)%pt1, surf_lsm_v(3)%pt1, surf_usm_v(3)%pt1, n_out )
2782
2783	CASE ( 'qv1' )
2784	CALL surface_data_output_sum_up( &
2785	surf_def_h(0)%qv1, surf_lsm_h(0)%qv1, surf_usm_h(0)%qv1, &
2786	surf_def_h(1)%qv1, surf_lsm_h(1)%qv1, surf_usm_h(1)%qv1, &
2787	surf_def_v(0)%qv1, surf_lsm_v(0)%qv1, surf_usm_v(0)%qv1, &
2788	surf_def_v(1)%qv1, surf_lsm_v(1)%qv1, surf_usm_v(1)%qv1, &
2789	surf_def_v(2)%qv1, surf_lsm_v(2)%qv1, surf_usm_v(2)%qv1, &
2790	surf_def_v(3)%qv1, surf_lsm_v(3)%qv1, surf_usm_v(3)%qv1, n_out )
2791
2792	CASE ( 'thetav1' )
2793	CALL surface_data_output_sum_up( &
2794	surf_def_h(0)%vpt1, surf_lsm_h(0)%vpt1, surf_usm_h(0)%vpt1, &
2795	surf_def_h(1)%vpt1, surf_lsm_h(1)%vpt1, surf_usm_h(1)%vpt1, &
2796	surf_def_v(0)%vpt1, surf_lsm_v(0)%vpt1, surf_usm_v(0)%vpt1, &
2797	surf_def_v(1)%vpt1, surf_lsm_v(1)%vpt1, surf_usm_v(1)%vpt1, &
2798	surf_def_v(2)%vpt1, surf_lsm_v(2)%vpt1, surf_usm_v(2)%vpt1, &
2799	surf_def_v(3)%vpt1, surf_lsm_v(3)%vpt1, surf_usm_v(3)%vpt1, n_out )
2800
2801	CASE ( 'usws' )
2802	CALL surface_data_output_sum_up( &
2803	surf_def_h(0)%usws, surf_lsm_h(0)%usws, surf_usm_h(0)%usws, &
2804	surf_def_h(1)%usws, surf_lsm_h(1)%usws, surf_usm_h(1)%usws, &
2805	surf_def_v(0)%usws, surf_lsm_v(0)%usws, surf_usm_v(0)%usws, &
2806	surf_def_v(1)%usws, surf_lsm_v(1)%usws, surf_usm_v(1)%usws, &
2807	surf_def_v(2)%usws, surf_lsm_v(2)%usws, surf_usm_v(2)%usws, &
2808	surf_def_v(3)%usws, surf_lsm_v(3)%usws, surf_usm_v(3)%usws, n_out, &
2809	momentumflux_output_conversion )
2810
2811	CASE ( 'vsws' )
2812	CALL surface_data_output_sum_up( &
2813	surf_def_h(0)%vsws, surf_lsm_h(0)%vsws, surf_usm_h(0)%vsws, &
2814	surf_def_h(1)%vsws, surf_lsm_h(1)%vsws, surf_usm_h(1)%vsws, &
2815	surf_def_v(0)%vsws, surf_lsm_v(0)%vsws, surf_usm_v(0)%vsws, &
2816	surf_def_v(1)%vsws, surf_lsm_v(1)%vsws, surf_usm_v(1)%vsws, &
2817	surf_def_v(2)%vsws, surf_lsm_v(2)%vsws, surf_usm_v(2)%vsws, &
2818	surf_def_v(3)%vsws, surf_lsm_v(3)%vsws, surf_usm_v(3)%vsws, n_out, &
2819	momentumflux_output_conversion )
2820
2821	CASE ( 'shf' )
2822	CALL surface_data_output_sum_up( &
2823	surf_def_h(0)%shf, surf_lsm_h(0)%shf, surf_usm_h(0)%shf, &
2824	surf_def_h(1)%shf, surf_lsm_h(1)%shf, surf_usm_h(1)%shf, &
2825	surf_def_v(0)%shf, surf_lsm_v(0)%shf, surf_usm_v(0)%shf, &
2826	surf_def_v(1)%shf, surf_lsm_v(1)%shf, surf_usm_v(1)%shf, &
2827	surf_def_v(2)%shf, surf_lsm_v(2)%shf, surf_usm_v(2)%shf, &
2828	surf_def_v(3)%shf, surf_lsm_v(3)%shf, surf_usm_v(3)%shf, n_out, &
2829	heatflux_output_conversion )
2830
2831	CASE ( 'qsws' )
2832	CALL surface_data_output_sum_up( &
2833	surf_def_h(0)%qsws, surf_lsm_h(0)%qsws, surf_usm_h(0)%qsws, &
2834	surf_def_h(1)%qsws, surf_lsm_h(1)%qsws, surf_usm_h(1)%qsws, &
2835	surf_def_v(0)%qsws, surf_lsm_v(0)%qsws, surf_usm_v(0)%qsws, &
2836	surf_def_v(1)%qsws, surf_lsm_v(1)%qsws, surf_usm_v(1)%qsws, &
2837	surf_def_v(2)%qsws, surf_lsm_v(2)%qsws, surf_usm_v(2)%qsws, &
2838	surf_def_v(3)%qsws, surf_lsm_v(3)%qsws, surf_usm_v(3)%qsws, n_out, &
2839	waterflux_output_conversion )
2840
2841	CASE ( 'ssws' )
2842	CALL surface_data_output_sum_up( &
2843	surf_def_h(0)%ssws, surf_lsm_h(0)%ssws, surf_usm_h(0)%ssws, &
2844	surf_def_h(1)%ssws, surf_lsm_h(1)%ssws, surf_usm_h(1)%ssws, &
2845	surf_def_v(0)%ssws, surf_lsm_v(0)%ssws, surf_usm_v(0)%ssws, &
2846	surf_def_v(1)%ssws, surf_lsm_v(1)%ssws, surf_usm_v(1)%ssws, &
2847	surf_def_v(2)%ssws, surf_lsm_v(2)%ssws, surf_usm_v(2)%ssws, &
2848	surf_def_v(3)%ssws, surf_lsm_v(3)%ssws, surf_usm_v(3)%ssws, n_out )
2849
2850	CASE ( 'qcsws' )
2851	CALL surface_data_output_sum_up( &
2852	surf_def_h(0)%qcsws, surf_lsm_h(0)%qcsws, surf_usm_h(0)%qcsws, &
2853	surf_def_h(1)%qcsws, surf_lsm_h(1)%qcsws, surf_usm_h(1)%qcsws, &
2854	surf_def_v(0)%qcsws, surf_lsm_v(0)%qcsws, surf_usm_v(0)%qcsws, &
2855	surf_def_v(1)%qcsws, surf_lsm_v(1)%qcsws, surf_usm_v(1)%qcsws, &
2856	surf_def_v(2)%qcsws, surf_lsm_v(2)%qcsws, surf_usm_v(2)%qcsws, &
2857	surf_def_v(3)%qcsws, surf_lsm_v(3)%qcsws, surf_usm_v(3)%qcsws, n_out )
2858
2859	CASE ( 'ncsws' )
2860	CALL surface_data_output_sum_up( &
2861	surf_def_h(0)%ncsws, surf_lsm_h(0)%ncsws, surf_usm_h(0)%ncsws, &
2862	surf_def_h(1)%ncsws, surf_lsm_h(1)%ncsws, surf_usm_h(1)%ncsws, &
2863	surf_def_v(0)%ncsws, surf_lsm_v(0)%ncsws, surf_usm_v(0)%ncsws, &
2864	surf_def_v(1)%ncsws, surf_lsm_v(1)%ncsws, surf_usm_v(1)%ncsws, &
2865	surf_def_v(2)%ncsws, surf_lsm_v(2)%ncsws, surf_usm_v(2)%ncsws, &
2866	surf_def_v(3)%ncsws, surf_lsm_v(3)%ncsws, surf_usm_v(3)%ncsws, n_out )
2867
2868	CASE ( 'qisws' )
2869	CALL surface_data_output_sum_up( &
2870	surf_def_h(0)%qisws, surf_lsm_h(0)%qisws, surf_usm_h(0)%qisws, &
2871	surf_def_h(1)%qisws, surf_lsm_h(1)%qisws, surf_usm_h(1)%qisws, &
2872	surf_def_v(0)%qisws, surf_lsm_v(0)%qisws, surf_usm_v(0)%qisws, &
2873	surf_def_v(1)%qisws, surf_lsm_v(1)%qisws, surf_usm_v(1)%qisws, &
2874	surf_def_v(2)%qisws, surf_lsm_v(2)%qisws, surf_usm_v(2)%qisws, &
2875	surf_def_v(3)%qisws, surf_lsm_v(3)%qisws, surf_usm_v(3)%qisws, n_out )
2876
2877	CASE ( 'nisws' )
2878	CALL surface_data_output_sum_up( &
2879	surf_def_h(0)%nisws, surf_lsm_h(0)%nisws, surf_usm_h(0)%nisws, &
2880	surf_def_h(1)%nisws, surf_lsm_h(1)%nisws, surf_usm_h(1)%nisws, &
2881	surf_def_v(0)%nisws, surf_lsm_v(0)%nisws, surf_usm_v(0)%nisws, &
2882	surf_def_v(1)%nisws, surf_lsm_v(1)%nisws, surf_usm_v(1)%nisws, &
2883	surf_def_v(2)%nisws, surf_lsm_v(2)%nisws, surf_usm_v(2)%nisws, &
2884	surf_def_v(3)%nisws, surf_lsm_v(3)%nisws, surf_usm_v(3)%nisws, n_out )
2885
2886	CASE ( 'qrsws' )
2887	CALL surface_data_output_sum_up( &
2888	surf_def_h(0)%qrsws, surf_lsm_h(0)%qrsws, surf_usm_h(0)%qrsws, &
2889	surf_def_h(1)%qrsws, surf_lsm_h(1)%qrsws, surf_usm_h(1)%qrsws, &
2890	surf_def_v(0)%qrsws, surf_lsm_v(0)%qrsws, surf_usm_v(0)%qrsws, &
2891	surf_def_v(1)%qrsws, surf_lsm_v(1)%qrsws, surf_usm_v(1)%qrsws, &
2892	surf_def_v(2)%qrsws, surf_lsm_v(2)%qrsws, surf_usm_v(2)%qrsws, &
2893	surf_def_v(3)%qrsws, surf_lsm_v(3)%qrsws, surf_usm_v(3)%qrsws, n_out )
2894
2895	CASE ( 'nrsws' )
2896	CALL surface_data_output_sum_up( &
2897	surf_def_h(0)%nrsws, surf_lsm_h(0)%nrsws, surf_usm_h(0)%nrsws, &
2898	surf_def_h(1)%nrsws, surf_lsm_h(1)%nrsws, surf_usm_h(1)%nrsws, &
2899	surf_def_v(0)%nrsws, surf_lsm_v(0)%nrsws, surf_usm_v(0)%nrsws, &
2900	surf_def_v(1)%nrsws, surf_lsm_v(1)%nrsws, surf_usm_v(1)%nrsws, &
2901	surf_def_v(2)%nrsws, surf_lsm_v(2)%nrsws, surf_usm_v(2)%nrsws, &
2902	surf_def_v(3)%nrsws, surf_lsm_v(3)%nrsws, surf_usm_v(3)%nrsws, n_out )
2903
2904	CASE ( 'sasws' )
2905	CALL surface_data_output_sum_up( &
2906	surf_def_h(0)%sasws, surf_lsm_h(0)%sasws, surf_usm_h(0)%sasws, &
2907	surf_def_h(1)%sasws, surf_lsm_h(1)%sasws, surf_usm_h(1)%sasws, &
2908	surf_def_v(0)%sasws, surf_lsm_v(0)%sasws, surf_usm_v(0)%sasws, &
2909	surf_def_v(1)%sasws, surf_lsm_v(1)%sasws, surf_usm_v(1)%sasws, &
2910	surf_def_v(2)%sasws, surf_lsm_v(2)%sasws, surf_usm_v(2)%sasws, &
2911	surf_def_v(3)%sasws, surf_lsm_v(3)%sasws, surf_usm_v(3)%sasws, n_out )
2912
2913	CASE ( 'q_surface' )
2914	CALL surface_data_output_sum_up( &
2915	surf_def_h(0)%q_surface, surf_lsm_h(0)%q_surface, surf_usm_h(0)%q_surface, &
2916	surf_def_h(1)%q_surface, surf_lsm_h(1)%q_surface, surf_usm_h(1)%q_surface, &
2917	surf_def_v(0)%q_surface, surf_lsm_v(0)%q_surface, surf_usm_v(0)%q_surface, &
2918	surf_def_v(1)%q_surface, surf_lsm_v(1)%q_surface, surf_usm_v(1)%q_surface, &
2919	surf_def_v(2)%q_surface, surf_lsm_v(2)%q_surface, surf_usm_v(2)%q_surface, &
2920	surf_def_v(3)%q_surface, surf_lsm_v(3)%q_surface, surf_usm_v(3)%q_surface, &
2921	n_out )
2922
2923	CASE ( 'theta_surface' )
2924	CALL surface_data_output_sum_up( &
2925	surf_def_h(0)%pt_surface, surf_lsm_h(0)%pt_surface, surf_usm_h(0)%pt_surface, &
2926	surf_def_h(1)%pt_surface, surf_lsm_h(1)%pt_surface, surf_usm_h(1)%pt_surface, &
2927	surf_def_v(0)%pt_surface, surf_lsm_v(0)%pt_surface, surf_usm_v(0)%pt_surface, &
2928	surf_def_v(1)%pt_surface, surf_lsm_v(1)%pt_surface, surf_usm_v(1)%pt_surface, &
2929	surf_def_v(2)%pt_surface, surf_lsm_v(2)%pt_surface, surf_usm_v(2)%pt_surface, &
2930	surf_def_v(3)%pt_surface, surf_lsm_v(3)%pt_surface, surf_usm_v(3)%pt_surface, &
2931	n_out )
2932
2933	CASE ( 'thetav_surface' )
2934	CALL surface_data_output_sum_up( &
2935	surf_def_h(0)%vpt_surface, surf_lsm_h(0)%vpt_surface, surf_usm_h(0)%vpt_surface, &
2936	surf_def_h(1)%vpt_surface, surf_lsm_h(1)%vpt_surface, surf_usm_h(1)%vpt_surface, &
2937	surf_def_v(0)%vpt_surface, surf_lsm_v(0)%vpt_surface, surf_usm_v(0)%vpt_surface, &
2938	surf_def_v(1)%vpt_surface, surf_lsm_v(1)%vpt_surface, surf_usm_v(1)%vpt_surface, &
2939	surf_def_v(2)%vpt_surface, surf_lsm_v(2)%vpt_surface, surf_usm_v(2)%vpt_surface, &
2940	surf_def_v(3)%vpt_surface, surf_lsm_v(3)%vpt_surface, surf_usm_v(3)%pt_surface, &
2941	n_out )
2942
2943	CASE ( 'rad_net' )
2944	CALL surface_data_output_sum_up( &
2945	surf_def_h(0)%rad_net, surf_lsm_h(0)%rad_net, surf_usm_h(0)%rad_net, &
2946	surf_def_h(1)%rad_net, surf_lsm_h(1)%rad_net, surf_usm_h(1)%rad_net, &
2947	surf_def_v(0)%rad_net, surf_lsm_v(0)%rad_net, surf_usm_v(0)%rad_net, &
2948	surf_def_v(1)%rad_net, surf_lsm_v(1)%rad_net, surf_usm_v(1)%rad_net, &
2949	surf_def_v(2)%rad_net, surf_lsm_v(2)%rad_net, surf_usm_v(2)%rad_net, &
2950	surf_def_v(3)%rad_net, surf_lsm_v(3)%rad_net, surf_usm_v(3)%rad_net, &
2951	n_out )
2952
2953	CASE ( 'rad_lw_in' )
2954	CALL surface_data_output_sum_up( &
2955	surf_def_h(0)%rad_lw_in, surf_lsm_h(0)%rad_lw_in, surf_usm_h(0)%rad_lw_in, &
2956	surf_def_h(1)%rad_lw_in, surf_lsm_h(1)%rad_lw_in, surf_usm_h(1)%rad_lw_in, &
2957	surf_def_v(0)%rad_lw_in, surf_lsm_v(0)%rad_lw_in, surf_usm_v(0)%rad_lw_in, &
2958	surf_def_v(1)%rad_lw_in, surf_lsm_v(1)%rad_lw_in, surf_usm_v(1)%rad_lw_in, &
2959	surf_def_v(2)%rad_lw_in, surf_lsm_v(2)%rad_lw_in, surf_usm_v(2)%rad_lw_in, &
2960	surf_def_v(3)%rad_lw_in, surf_lsm_v(3)%rad_lw_in, surf_usm_v(3)%rad_lw_in, &
2961	n_out )
2962
2963	CASE ( 'rad_lw_out' )
2964	CALL surface_data_output_sum_up( &
2965	surf_def_h(0)%rad_lw_out, surf_lsm_h(0)%rad_lw_out, surf_usm_h(0)%rad_lw_out, &
2966	surf_def_h(1)%rad_lw_out, surf_lsm_h(1)%rad_lw_out, surf_usm_h(1)%rad_lw_out, &
2967	surf_def_v(0)%rad_lw_out, surf_lsm_v(0)%rad_lw_out, surf_usm_v(0)%rad_lw_out, &
2968	surf_def_v(1)%rad_lw_out, surf_lsm_v(1)%rad_lw_out, surf_usm_v(1)%rad_lw_out, &
2969	surf_def_v(2)%rad_lw_out, surf_lsm_v(2)%rad_lw_out, surf_usm_v(2)%rad_lw_out, &
2970	surf_def_v(3)%rad_lw_out, surf_lsm_v(3)%rad_lw_out, surf_usm_v(3)%rad_lw_out, &
2971	n_out )
2972
2973	CASE ( 'rad_sw_in' )
2974	CALL surface_data_output_sum_up( &
2975	surf_def_h(0)%rad_sw_in, surf_lsm_h(0)%rad_sw_in, surf_usm_h(0)%rad_sw_in, &
2976	surf_def_h(1)%rad_sw_in, surf_lsm_h(1)%rad_sw_in, surf_usm_h(1)%rad_sw_in, &
2977	surf_def_v(0)%rad_sw_in, surf_lsm_v(0)%rad_sw_in, surf_usm_v(0)%rad_sw_in, &
2978	surf_def_v(1)%rad_sw_in, surf_lsm_v(1)%rad_sw_in, surf_usm_v(1)%rad_sw_in, &
2979	surf_def_v(2)%rad_sw_in, surf_lsm_v(2)%rad_sw_in, surf_usm_v(2)%rad_sw_in, &
2980	surf_def_v(3)%rad_sw_in, surf_lsm_v(3)%rad_sw_in, surf_usm_v(3)%rad_sw_in, &
2981	n_out )
2982
2983	CASE ( 'rad_sw_out' )
2984	CALL surface_data_output_sum_up( &
2985	surf_def_h(0)%rad_sw_out, surf_lsm_h(0)%rad_sw_out, surf_usm_h(0)%rad_sw_out, &
2986	surf_def_h(1)%rad_sw_out, surf_lsm_h(1)%rad_sw_out, surf_usm_h(1)%rad_sw_out, &
2987	surf_def_v(0)%rad_sw_out, surf_lsm_v(0)%rad_sw_out, surf_usm_v(0)%rad_sw_out, &
2988	surf_def_v(1)%rad_sw_out, surf_lsm_v(1)%rad_sw_out, surf_usm_v(1)%rad_sw_out, &
2989	surf_def_v(2)%rad_sw_out, surf_lsm_v(2)%rad_sw_out, surf_usm_v(2)%rad_sw_out, &
2990	surf_def_v(3)%rad_sw_out, surf_lsm_v(3)%rad_sw_out, surf_usm_v(3)%rad_sw_out, &
2991	n_out )
2992
2993	CASE ( 'ghf' )
2994	!
2995	!-- Sum up ground / wall heat flux. Note, for urban surfaces the wall heat flux is
2996	!-- aggregated from the different green, window and wall tiles.
2997	DO l = 0, 1
2998	DO m = 1, surf_usm_h(l)%ns
2999	surf_usm_h(l)%ghf(m) = surf_usm_h(l)%frac(m,ind_veg_wall) * &
3000	surf_usm_h(l)%wghf_eb(m) + &
3001	surf_usm_h(l)%frac(m,ind_pav_green) * &
3002	surf_usm_h(l)%wghf_eb_green(m) + &
3003	surf_usm_h(l)%frac(m,ind_wat_win) * &
3004	surf_usm_h(l)%wghf_eb_window(m)
3005	ENDDO
3006	ENDDO
3007	DO l = 0, 3
3008	DO m = 1, surf_usm_v(l)%ns
3009	surf_usm_v(l)%ghf(m) = surf_usm_v(l)%frac(m,ind_veg_wall) * &
3010	surf_usm_v(l)%wghf_eb(m) + &
3011	surf_usm_v(l)%frac(m,ind_pav_green) * &
3012	surf_usm_v(l)%wghf_eb_green(m) + &
3013	surf_usm_v(l)%frac(m,ind_wat_win) * &
3014	surf_usm_v(l)%wghf_eb_window(m)
3015	ENDDO
3016	ENDDO
3017
3018	CALL surface_data_output_sum_up( &
3019	surf_def_h(0)%ghf, surf_lsm_h(0)%ghf, surf_usm_h(0)%ghf, &
3020	surf_def_h(1)%ghf, surf_lsm_h(1)%ghf, surf_usm_h(1)%ghf, &
3021	surf_def_v(0)%ghf, surf_lsm_v(0)%ghf, surf_usm_v(0)%ghf, &
3022	surf_def_v(1)%ghf, surf_lsm_v(1)%ghf, surf_usm_v(1)%ghf, &
3023	surf_def_v(2)%ghf, surf_lsm_v(2)%ghf, surf_usm_v(2)%ghf, &
3024	surf_def_v(3)%ghf, surf_lsm_v(3)%ghf, surf_usm_v(3)%ghf, n_out )
3025
3026	CASE ( 'r_a' )
3027	CALL surface_data_output_sum_up( &
3028	surf_def_h(0)%r_a, surf_lsm_h(0)%r_a, surf_usm_h(0)%r_a, &
3029	surf_def_h(1)%r_a, surf_lsm_h(1)%r_a, surf_usm_h(1)%r_a, &
3030	surf_def_v(0)%r_a, surf_lsm_v(0)%r_a, surf_usm_v(0)%r_a, &
3031	surf_def_v(1)%r_a, surf_lsm_v(1)%r_a, surf_usm_v(1)%r_a, &
3032	surf_def_v(2)%r_a, surf_lsm_v(2)%r_a, surf_usm_v(2)%r_a, &
3033	surf_def_v(3)%r_a, surf_lsm_v(3)%r_a, surf_usm_v(3)%r_a, n_out )
3034
3035	CASE ( 'r_soil' )
3036	CALL surface_data_output_sum_up( &
3037	surf_def_h(0)%r_soil, surf_lsm_h(0)%r_soil, surf_usm_h(0)%r_soil, &
3038	surf_def_h(1)%r_soil, surf_lsm_h(1)%r_soil, surf_usm_h(1)%r_soil, &
3039	surf_def_v(0)%r_soil, surf_lsm_v(0)%r_soil, surf_usm_v(0)%r_soil, &
3040	surf_def_v(1)%r_soil, surf_lsm_v(1)%r_soil, surf_usm_v(1)%r_soil, &
3041	surf_def_v(2)%r_soil, surf_lsm_v(2)%r_soil, surf_usm_v(2)%r_soil, &
3042	surf_def_v(3)%r_soil, surf_lsm_v(3)%r_soil, surf_usm_v(3)%r_soil, n_out )
3043
3044	CASE ( 'r_canopy' )
3045	CALL surface_data_output_sum_up( &
3046	surf_def_h(0)%r_canopy, surf_lsm_h(0)%r_canopy, surf_usm_h(0)%r_canopy, &
3047	surf_def_h(1)%r_canopy, surf_lsm_h(1)%r_canopy, surf_usm_h(1)%r_canopy, &
3048	surf_def_v(0)%r_canopy, surf_lsm_v(0)%r_canopy, surf_usm_v(0)%r_canopy, &
3049	surf_def_v(1)%r_canopy, surf_lsm_v(1)%r_canopy, surf_usm_v(1)%r_canopy, &
3050	surf_def_v(2)%r_canopy, surf_lsm_v(2)%r_canopy, surf_usm_v(2)%r_canopy, &
3051	surf_def_v(3)%r_canopy, surf_lsm_v(3)%r_canopy, surf_usm_v(3)%r_canopy, &
3052	n_out )
3053
3054	CASE ( 'r_s' )
3055	CALL surface_data_output_sum_up( &
3056	surf_def_h(0)%r_s, surf_lsm_h(0)%r_s, surf_usm_h(0)%r_s, &
3057	surf_def_h(1)%r_s, surf_lsm_h(1)%r_s, surf_usm_h(1)%r_s, &
3058	surf_def_v(0)%r_s, surf_lsm_v(0)%r_s, surf_usm_v(0)%r_s, &
3059	surf_def_v(1)%r_s, surf_lsm_v(1)%r_s, surf_usm_v(1)%r_s, &
3060	surf_def_v(2)%r_s, surf_lsm_v(2)%r_s, surf_usm_v(2)%r_s, &
3061	surf_def_v(3)%r_s, surf_lsm_v(3)%r_s, surf_usm_v(3)%r_s, n_out )
3062
3063
3064	CASE ( 'rad_sw_dir' )
3065	CALL surface_data_output_sum_up( &
3066	surf_def_h(0)%rad_sw_dir, surf_lsm_h(0)%rad_sw_dir, surf_usm_h(0)%rad_sw_dir, &
3067	surf_def_h(1)%rad_sw_dir, surf_lsm_h(1)%rad_sw_dir, surf_usm_h(1)%rad_sw_dir, &
3068	surf_def_v(0)%rad_sw_dir, surf_lsm_v(0)%rad_sw_dir, surf_usm_v(0)%rad_sw_dir, &
3069	surf_def_v(1)%rad_sw_dir, surf_lsm_v(1)%rad_sw_dir, surf_usm_v(1)%rad_sw_dir, &
3070	surf_def_v(2)%rad_sw_dir, surf_lsm_v(2)%rad_sw_dir, surf_usm_v(2)%rad_sw_dir, &
3071	surf_def_v(3)%rad_sw_dir, surf_lsm_v(3)%rad_sw_dir, surf_usm_v(3)%rad_sw_dir, &
3072	n_out )
3073	CASE ( 'rad_sw_dif' )
3074	CALL surface_data_output_sum_up( &
3075	surf_def_h(0)%rad_sw_dif, surf_lsm_h(0)%rad_sw_dif, surf_usm_h(0)%rad_sw_dif, &
3076	surf_def_h(1)%rad_sw_dif, surf_lsm_h(1)%rad_sw_dif, surf_usm_h(1)%rad_sw_dif, &
3077	surf_def_v(0)%rad_sw_dif, surf_lsm_v(0)%rad_sw_dif, surf_usm_v(0)%rad_sw_dif, &
3078	surf_def_v(1)%rad_sw_dif, surf_lsm_v(1)%rad_sw_dif, surf_usm_v(1)%rad_sw_dif, &
3079	surf_def_v(2)%rad_sw_dif, surf_lsm_v(2)%rad_sw_dif, surf_usm_v(2)%rad_sw_dif, &
3080	surf_def_v(3)%rad_sw_dif, surf_lsm_v(3)%rad_sw_dif, surf_usm_v(3)%rad_sw_dif, &
3081	n_out )
3082
3083	CASE ( 'rad_sw_ref' )
3084	CALL surface_data_output_sum_up( &
3085	surf_def_h(0)%rad_sw_ref, surf_lsm_h(0)%rad_sw_ref, surf_usm_h(0)%rad_sw_ref, &
3086	surf_def_h(1)%rad_sw_ref, surf_lsm_h(1)%rad_sw_ref, surf_usm_h(1)%rad_sw_ref, &
3087	surf_def_v(0)%rad_sw_ref, surf_lsm_v(0)%rad_sw_ref, surf_usm_v(0)%rad_sw_ref, &
3088	surf_def_v(1)%rad_sw_ref, surf_lsm_v(1)%rad_sw_ref, surf_usm_v(1)%rad_sw_ref, &
3089	surf_def_v(2)%rad_sw_ref, surf_lsm_v(2)%rad_sw_ref, surf_usm_v(2)%rad_sw_ref, &
3090	surf_def_v(3)%rad_sw_ref, surf_lsm_v(3)%rad_sw_ref, surf_usm_v(3)%rad_sw_ref, &
3091	n_out )
3092
3093	CASE ( 'rad_sw_res' )
3094	CALL surface_data_output_sum_up( &
3095	surf_def_h(0)%rad_sw_res, surf_lsm_h(0)%rad_sw_res, surf_usm_h(0)%rad_sw_res, &
3096	surf_def_h(1)%rad_sw_res, surf_lsm_h(1)%rad_sw_res, surf_usm_h(1)%rad_sw_res, &
3097	surf_def_v(0)%rad_sw_res, surf_lsm_v(0)%rad_sw_res, surf_usm_v(0)%rad_sw_res, &
3098	surf_def_v(1)%rad_sw_res, surf_lsm_v(1)%rad_sw_res, surf_usm_v(1)%rad_sw_res, &
3099	surf_def_v(2)%rad_sw_res, surf_lsm_v(2)%rad_sw_res, surf_usm_v(2)%rad_sw_res, &
3100	surf_def_v(3)%rad_sw_res, surf_lsm_v(3)%rad_sw_res, surf_usm_v(3)%rad_sw_res, &
3101	n_out )
3102
3103	CASE ( 'rad_lw_dif' )
3104	CALL surface_data_output_sum_up( &
3105	surf_def_h(0)%rad_lw_dif, surf_lsm_h(0)%rad_lw_dif, surf_usm_h(0)%rad_lw_dif, &
3106	surf_def_h(1)%rad_lw_dif, surf_lsm_h(1)%rad_lw_dif, surf_usm_h(1)%rad_lw_dif, &
3107	surf_def_v(0)%rad_lw_dif, surf_lsm_v(0)%rad_lw_dif, surf_usm_v(0)%rad_lw_dif, &
3108	surf_def_v(1)%rad_lw_dif, surf_lsm_v(1)%rad_lw_dif, surf_usm_v(1)%rad_lw_dif, &
3109	surf_def_v(2)%rad_lw_dif, surf_lsm_v(2)%rad_lw_dif, surf_usm_v(2)%rad_lw_dif, &
3110	surf_def_v(3)%rad_lw_dif, surf_lsm_v(3)%rad_lw_dif, surf_usm_v(3)%rad_lw_dif, &
3111	n_out )
3112
3113	CASE ( 'rad_lw_ref' )
3114	CALL surface_data_output_sum_up( &
3115	surf_def_h(0)%rad_lw_ref, surf_lsm_h(0)%rad_lw_ref, surf_usm_h(0)%rad_lw_ref, &
3116	surf_def_h(1)%rad_lw_ref, surf_lsm_h(1)%rad_lw_ref, surf_usm_h(1)%rad_lw_ref, &
3117	surf_def_v(0)%rad_lw_ref, surf_lsm_v(0)%rad_lw_ref, surf_usm_v(0)%rad_lw_ref, &
3118	surf_def_v(1)%rad_lw_ref, surf_lsm_v(1)%rad_lw_ref, surf_usm_v(1)%rad_lw_ref, &
3119	surf_def_v(2)%rad_lw_ref, surf_lsm_v(2)%rad_lw_ref, surf_usm_v(2)%rad_lw_ref, &
3120	surf_def_v(3)%rad_lw_ref, surf_lsm_v(3)%rad_lw_ref, surf_usm_v(3)%rad_lw_ref, &
3121	n_out )
3122
3123	CASE ( 'rad_lw_res' )
3124	CALL surface_data_output_sum_up( &
3125	surf_def_h(0)%rad_lw_res, surf_lsm_h(0)%rad_lw_res, surf_usm_h(0)%rad_lw_res, &
3126	surf_def_h(1)%rad_lw_res, surf_lsm_h(1)%rad_lw_res, surf_usm_h(1)%rad_lw_res, &
3127	surf_def_v(0)%rad_lw_res, surf_lsm_v(0)%rad_lw_res, surf_usm_v(0)%rad_lw_res, &
3128	surf_def_v(1)%rad_lw_res, surf_lsm_v(1)%rad_lw_res, surf_usm_v(1)%rad_lw_res, &
3129	surf_def_v(2)%rad_lw_res, surf_lsm_v(2)%rad_lw_res, surf_usm_v(2)%rad_lw_res, &
3130	surf_def_v(3)%rad_lw_res, surf_lsm_v(3)%rad_lw_res, surf_usm_v(3)%rad_lw_res, &
3131	n_out )
3132
3133	CASE ( 'uvw1' )
3134	CALL surface_data_output_sum_up( &
3135	surf_def_h(0)%uvw_abs, surf_lsm_h(0)%uvw_abs, surf_usm_h(0)%uvw_abs, &
3136	surf_def_h(1)%uvw_abs, surf_lsm_h(1)%uvw_abs, surf_usm_h(1)%uvw_abs, &
3137	surf_def_v(0)%uvw_abs, surf_lsm_v(0)%uvw_abs, surf_usm_v(0)%uvw_abs, &
3138	surf_def_v(1)%uvw_abs, surf_lsm_v(1)%uvw_abs, surf_usm_v(1)%uvw_abs, &
3139	surf_def_v(2)%uvw_abs, surf_lsm_v(2)%uvw_abs, surf_usm_v(2)%uvw_abs, &
3140	surf_def_v(3)%uvw_abs, surf_lsm_v(3)%uvw_abs, surf_usm_v(3)%uvw_abs, n_out )
3141
3142	CASE ( 'waste_heat' )
3143	CALL surface_data_output_sum_up( &
3144	surf_def_h(0)%waste_heat, surf_lsm_h(0)%waste_heat, surf_usm_h(0)%waste_heat,&
3145	surf_def_h(1)%waste_heat, surf_lsm_h(1)%waste_heat, surf_usm_h(1)%waste_heat,&
3146	surf_def_v(0)%waste_heat, surf_lsm_v(0)%waste_heat, surf_usm_v(0)%waste_heat,&
3147	surf_def_v(1)%waste_heat, surf_lsm_v(1)%waste_heat, surf_usm_v(1)%waste_heat,&
3148	surf_def_v(2)%waste_heat, surf_lsm_v(2)%waste_heat, surf_usm_v(2)%waste_heat,&
3149	surf_def_v(3)%waste_heat, surf_lsm_v(3)%waste_heat, surf_usm_v(3)%waste_heat,&
3150	n_out )
3151
3152	CASE ( 'im_hf' )
3153	CALL surface_data_output_sum_up( &
3154	surf_def_h(0)%iwghf_eb, surf_lsm_h(0)%iwghf_eb, surf_usm_h(0)%iwghf_eb, &
3155	surf_def_h(1)%iwghf_eb, surf_lsm_h(1)%iwghf_eb, surf_usm_h(1)%iwghf_eb, &
3156	surf_def_v(0)%iwghf_eb, surf_lsm_v(0)%iwghf_eb, surf_usm_v(0)%iwghf_eb, &
3157	surf_def_v(1)%iwghf_eb, surf_lsm_v(1)%iwghf_eb, surf_usm_v(1)%iwghf_eb, &
3158	surf_def_v(2)%iwghf_eb, surf_lsm_v(2)%iwghf_eb, surf_usm_v(2)%iwghf_eb, &
3159	surf_def_v(3)%iwghf_eb, surf_lsm_v(3)%iwghf_eb, surf_usm_v(3)%iwghf_eb, &
3160	n_out )
3161
3162	CASE ( 'albedo' )
3163	CALL surface_data_output_sum_up( &
3164	surf_def_h(0)%albedo, surf_lsm_h(0)%albedo, surf_usm_h(0)%albedo, &
3165	surf_def_h(1)%albedo, surf_lsm_h(1)%albedo, surf_usm_h(1)%albedo, &
3166	surf_def_v(0)%albedo, surf_lsm_v(0)%albedo, surf_usm_v(0)%albedo, &
3167	surf_def_v(1)%albedo, surf_lsm_v(1)%albedo, surf_usm_v(1)%albedo, &
3168	surf_def_v(2)%albedo, surf_lsm_v(2)%albedo, surf_usm_v(2)%albedo, &
3169	surf_def_v(3)%albedo, surf_lsm_v(3)%albedo, surf_usm_v(3)%albedo, n_out )
3170
3171
3172	CASE ( 'emissivity' )
3173	CALL surface_data_output_sum_up( &
3174	surf_def_h(0)%emissivity, surf_lsm_h(0)%emissivity, surf_usm_h(0)%emissivity,&
3175	surf_def_h(1)%emissivity, surf_lsm_h(1)%emissivity, surf_usm_h(1)%emissivity,&
3176	surf_def_v(0)%emissivity, surf_lsm_v(0)%emissivity, surf_usm_v(0)%emissivity,&
3177	surf_def_v(1)%emissivity, surf_lsm_v(1)%emissivity, surf_usm_v(1)%emissivity,&
3178	surf_def_v(2)%emissivity, surf_lsm_v(2)%emissivity, surf_usm_v(2)%emissivity,&
3179	surf_def_v(3)%emissivity, surf_lsm_v(3)%emissivity, surf_usm_v(3)%emissivity,&
3180	n_out )
3181
3182	END SELECT
3183	ENDDO
3184
3185
3186	END SUBROUTINE surface_data_output_averaging
3187
3188	!--------------------------------------------------------------------------------------------------!
3189	! Description:
3190	! ------------
3191	!> Sum-up the surface data for average output variables.
3192	!--------------------------------------------------------------------------------------------------!
3193	SUBROUTINE surface_data_output_sum_up_1d( var_def_h0, var_lsm_h0, var_usm_h0, &
3194	var_def_h1, var_lsm_h1, var_usm_h1, &
3195	var_def_v0, var_lsm_v0, var_usm_v0, &
3196	var_def_v1, var_lsm_v1, var_usm_v1, &
3197	var_def_v2, var_lsm_v2, var_usm_v2, &
3198	var_def_v3, var_lsm_v3, var_usm_v3, &
3199	n_out, fac )
3200
3201	IMPLICIT NONE
3202
3203	INTEGER(iwp) :: k !< height index of surface element
3204	INTEGER(iwp) :: m !< running index for surface elements
3205	INTEGER(iwp) :: n_out !< index for output variable
3206	INTEGER(iwp) :: n_surf !< running index for surface elements
3207
3208	REAL(wp), DIMENSION(:), OPTIONAL :: fac !< passed output conversion factor for heatflux output
3209	REAL(wp), DIMENSION(nzb:nzt+1) :: conversion_factor !< effective array for output conversion factor
3210
3211	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h0 !< output variable at upward-facing default-type surfaces
3212	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_h0 !< output variable at upward-facing natural-type surfaces
3213	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_h0 !< output variable at upward-facing urban-type surfaces
3214	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h1 !< output variable at downward-facing default-type surfaces
3215	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_h1 !< output variable at downward-facing natural-type surfaces
3216	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_h1 !< output variable at downward-facing urban-type surfaces
3217	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v0 !< output variable at northward-facing default-type surfaces
3218	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v1 !< output variable at southward-facing default-type surfaces
3219	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v2 !< output variable at eastward-facing default-type surfaces
3220	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v3 !< output variable at westward-facing default-type surfaces
3221	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v0 !< output variable at northward-facing natural-type surfaces
3222	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v1 !< output variable at southward-facing natural-type surfaces
3223	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v2 !< output variable at eastward-facing natural-type surfaces
3224	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v3 !< output variable at westward-facing natural-type surfaces
3225	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v0 !< output variable at northward-facing urban-type surfaces
3226	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v1 !< output variable at southward-facing urban-type surfaces
3227	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v2 !< output variable at eastward-facing urban-type surfaces
3228	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v3 !< output variable at westward-facing urban-type surfaces
3229
3230	!
3231	!-- Set conversion factor to one if not present
3232	IF ( .NOT. PRESENT( fac ) ) THEN
3233	conversion_factor = 1.0_wp
3234	ELSE
3235	conversion_factor = fac
3236	ENDIF
3237	!
3238	!-- Set counter variable to zero before the variable is written to the output array.
3239	n_surf = 0
3240
3241	!
3242	!-- Write the horizontal surfaces.
3243	!-- Before each variable is written to the output data structure, first check if the variable
3244	!-- for the respective surface type is defined. If a variable is not defined, skip the block and
3245	!-- increment the counter variable by the number of surface elements of this type. Usually this is
3246	!-- zero, however, there might be the situation that e.g. urban surfaces are defined but the
3247	!-- respective variable is not allocated for this surface type. To write the data on the exact
3248	!-- position, increment the counter.
3249	IF ( ALLOCATED( var_def_h0 ) ) THEN
3250	DO m = 1, surf_def_h(0)%ns
3251	n_surf = n_surf + 1
3252	k = surf_def_h(0)%k(m)
3253	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_def_h0(m) * &
3254	conversion_factor(k)
3255	ENDDO
3256	ELSE
3257	n_surf = n_surf + surf_def_h(0)%ns
3258	ENDIF
3259	IF ( ALLOCATED( var_lsm_h0 ) ) THEN
3260	DO m = 1, surf_lsm_h(0)%ns
3261	n_surf = n_surf + 1
3262	k = surf_lsm_h(0)%k(m)
3263	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_lsm_h0(m) * &
3264	conversion_factor(k)
3265	ENDDO
3266	ELSE
3267	n_surf = n_surf + surf_lsm_h(0)%ns
3268	ENDIF
3269	IF ( ALLOCATED( var_usm_h0 ) ) THEN
3270	DO m = 1, surf_usm_h(0)%ns
3271	n_surf = n_surf + 1
3272	k = surf_usm_h(0)%k(m)
3273	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_usm_h0(m) * &
3274	conversion_factor(k)
3275	ENDDO
3276	ELSE
3277	n_surf = n_surf + surf_usm_h(0)%ns
3278	ENDIF
3279	IF ( ALLOCATED( var_def_h1 ) ) THEN
3280	DO m = 1, surf_def_h(1)%ns
3281	n_surf = n_surf + 1
3282	k = surf_def_h(1)%k(m)
3283	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_def_h1(m) * &
3284	conversion_factor(k)
3285	ENDDO
3286	ELSE
3287	n_surf = n_surf + surf_def_h(1)%ns
3288	ENDIF
3289	IF ( ALLOCATED( var_lsm_h1 ) ) THEN
3290	DO m = 1, surf_lsm_h(1)%ns
3291	n_surf = n_surf + 1
3292	k = surf_lsm_h(1)%k(m)
3293	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_lsm_h1(m) * &
3294	conversion_factor(k)
3295	ENDDO
3296	ELSE
3297	n_surf = n_surf + surf_lsm_h(1)%ns
3298	ENDIF
3299	IF ( ALLOCATED( var_usm_h1 ) ) THEN
3300	DO m = 1, surf_usm_h(1)%ns
3301	n_surf = n_surf + 1
3302	k = surf_usm_h(1)%k(m)
3303	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_usm_h1(m) * &
3304	conversion_factor(k)
3305	ENDDO
3306	ELSE
3307	n_surf = n_surf + surf_usm_h(1)%ns
3308	ENDIF
3309	!
3310	!-- Write northward-facing
3311	IF ( ALLOCATED( var_def_v0 ) ) THEN
3312	DO m = 1, surf_def_v(0)%ns
3313	n_surf = n_surf + 1
3314	k = surf_def_v(0)%k(m)
3315	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_def_v0(m) * &
3316	conversion_factor(k)
3317	ENDDO
3318	ELSE
3319	n_surf = n_surf + surf_def_v(0)%ns
3320	ENDIF
3321	IF ( ALLOCATED( var_lsm_v0 ) ) THEN
3322	DO m = 1, surf_lsm_v(0)%ns
3323	n_surf = n_surf + 1
3324	k = surf_lsm_v(0)%k(m)
3325	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_lsm_v0(m) * &
3326	conversion_factor(k)
3327	ENDDO
3328	ELSE
3329	n_surf = n_surf + surf_lsm_v(0)%ns
3330	ENDIF
3331	IF ( ALLOCATED( var_usm_v0 ) ) THEN
3332	DO m = 1, surf_usm_v(0)%ns
3333	n_surf = n_surf + 1
3334	k = surf_usm_v(0)%k(m)
3335	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_usm_v0(m) * &
3336	conversion_factor(k)
3337	ENDDO
3338	ELSE
3339	n_surf = n_surf + surf_usm_v(0)%ns
3340	ENDIF
3341	!
3342	!-- Write southward-facing
3343	IF ( ALLOCATED( var_def_v1 ) ) THEN
3344	DO m = 1, surf_def_v(1)%ns
3345	n_surf = n_surf + 1
3346	k = surf_def_v(1)%k(m)
3347	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_def_v1(m) * &
3348	conversion_factor(k)
3349	ENDDO
3350	ELSE
3351	n_surf = n_surf + surf_def_v(1)%ns
3352	ENDIF
3353	IF ( ALLOCATED( var_lsm_v1 ) ) THEN
3354	DO m = 1, surf_lsm_v(1)%ns
3355	n_surf = n_surf + 1
3356	k = surf_lsm_v(1)%k(m)
3357	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_lsm_v1(m) * &
3358	conversion_factor(k)
3359	ENDDO
3360	ELSE
3361	n_surf = n_surf + surf_lsm_v(1)%ns
3362	ENDIF
3363	IF ( ALLOCATED( var_usm_v1 ) ) THEN
3364	DO m = 1, surf_usm_v(1)%ns
3365	n_surf = n_surf + 1
3366	k = surf_usm_v(1)%k(m)
3367	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) + var_usm_v1(m) * &
3368	conversion_factor(k)
3369	ENDDO
3370	ELSE
3371	n_surf = n_surf + surf_usm_v(1)%ns
3372	ENDIF
3373	!
3374	!-- Write eastward-facing
3375	IF ( ALLOCATED( var_def_v2 ) ) THEN
3376	DO m = 1, surf_def_v(2)%ns
3377	n_surf = n_surf + 1
3378	k = surf_def_v(2)%k(m)
3379	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3380	+ var_def_v2(m) * conversion_factor(k)
3381	ENDDO
3382	ELSE
3383	n_surf = n_surf + surf_def_v(2)%ns
3384	ENDIF
3385	IF ( ALLOCATED( var_lsm_v2 ) ) THEN
3386	DO m = 1, surf_lsm_v(2)%ns
3387	n_surf = n_surf + 1
3388	k = surf_lsm_v(2)%k(m)
3389	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3390	+ var_lsm_v2(m) * conversion_factor(k)
3391	ENDDO
3392	ELSE
3393	n_surf = n_surf + surf_lsm_v(2)%ns
3394	ENDIF
3395	IF ( ALLOCATED( var_usm_v2 ) ) THEN
3396	DO m = 1, surf_usm_v(2)%ns
3397	n_surf = n_surf + 1
3398	k = surf_usm_v(2)%k(m)
3399	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3400	+ var_usm_v2(m) * conversion_factor(k)
3401	ENDDO
3402	ELSE
3403	n_surf = n_surf + surf_usm_v(2)%ns
3404	ENDIF
3405	!
3406	!-- Write westward-facing
3407	IF ( ALLOCATED( var_def_v3 ) ) THEN
3408	DO m = 1, surf_def_v(3)%ns
3409	n_surf = n_surf + 1
3410	k = surf_def_v(3)%k(m)
3411	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3412	+ var_def_v3(m) * conversion_factor(k)
3413	ENDDO
3414	ELSE
3415	n_surf = n_surf + surf_def_v(3)%ns
3416	ENDIF
3417	IF ( ALLOCATED( var_lsm_v3 ) ) THEN
3418	DO m = 1, surf_lsm_v(3)%ns
3419	n_surf = n_surf + 1
3420	k = surf_lsm_v(3)%k(m)
3421	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3422	+ var_lsm_v3(m) * conversion_factor(k)
3423	ENDDO
3424	ELSE
3425	n_surf = n_surf + surf_lsm_v(3)%ns
3426	ENDIF
3427	IF ( ALLOCATED( var_usm_v3 ) ) THEN
3428	DO m = 1, surf_usm_v(3)%ns
3429	n_surf = n_surf + 1
3430	k = surf_usm_v(3)%k(m)
3431	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3432	+ var_usm_v3(m) * conversion_factor(k)
3433	ENDDO
3434	ELSE
3435	n_surf = n_surf + surf_usm_v(3)%ns
3436	ENDIF
3437
3438	END SUBROUTINE surface_data_output_sum_up_1d
3439
3440	!--------------------------------------------------------------------------------------------------!
3441	! Description:
3442	! ------------
3443	!> Sum-up the surface data for average output variables for properties which are defined using tile
3444	!> approach.
3445	!--------------------------------------------------------------------------------------------------!
3446	SUBROUTINE surface_data_output_sum_up_2d( var_def_h0, var_lsm_h0, var_usm_h0, &
3447	var_def_h1, var_lsm_h1, var_usm_h1, &
3448	var_def_v0, var_lsm_v0, var_usm_v0, &
3449	var_def_v1, var_lsm_v1, var_usm_v1, &
3450	var_def_v2, var_lsm_v2, var_usm_v2, &
3451	var_def_v3, var_lsm_v3, var_usm_v3, &
3452	n_out, fac )
3453
3454	IMPLICIT NONE
3455
3456	INTEGER(iwp) :: k !< height index of surface element
3457	INTEGER(iwp) :: m !< running index for surface elements
3458	INTEGER(iwp) :: n_out !< index for output variable
3459	INTEGER(iwp) :: n_surf !< running index for surface elements
3460
3461	REAL(wp), DIMENSION(:), OPTIONAL :: fac !< passed output conversion factor for heatflux output
3462	REAL(wp), DIMENSION(nzb:nzt+1) :: conversion_factor !< effective array for output conversion factor
3463
3464	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_h0 !< output variable at upward-facing default-type surfaces
3465	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_h0 !< output variable at upward-facing natural-type surfaces
3466	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_h0 !< output variable at upward-facing urban-type surfaces
3467	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_h1 !< output variable at downward-facing default-type surfaces
3468	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_h1 !< output variable at downward-facing natural-type surfaces
3469	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_h1 !< output variable at downward-facing urban-type surfaces
3470	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v0 !< output variable at northward-facing default-type surfaces
3471	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v1 !< output variable at southward-facing default-type surfaces
3472	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v2 !< output variable at eastward-facing default-type surfaces
3473	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v3 !< output variable at westward-facing default-type surfaces
3474	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v0 !< output variable at northward-facing natural-type surfaces
3475	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v1 !< output variable at southward-facing natural-type surfaces
3476	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v2 !< output variable at eastward-facing natural-type surfaces
3477	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v3 !< output variable at westward-facing natural-type surfaces
3478	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v0 !< output variable at northward-facing urban-type surfaces
3479	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v1 !< output variable at southward-facing urban-type surfaces
3480	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v2 !< output variable at eastward-facing urban-type surfaces
3481	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v3 !< output variable at westward-facing urban-type surfaces
3482
3483	!
3484	!-- Set conversion factor to one if not present
3485	IF ( .NOT. PRESENT( fac ) ) THEN
3486	conversion_factor = 1.0_wp
3487	ELSE
3488	conversion_factor = fac
3489	ENDIF
3490	!
3491	!-- Set counter variable to zero before the variable is written to the output array.
3492	n_surf = 0
3493
3494	!
3495	!-- Write the horizontal surfaces.
3496	!-- Before each variable is written to the output data structure, first check if the variable
3497	!-- for the respective surface type is defined. If a variable is not defined, skip the block and
3498	!-- increment the counter variable by the number of surface elements of this type. Usually this is
3499	!-- zero, however, there might be the situation that e.g. urban surfaces are defined but the
3500	!-- respective variable is not allocated for this surface type. To write the data on the exact
3501	!-- position, increment the counter.
3502	IF ( ALLOCATED( var_def_h0 ) ) THEN
3503	DO m = 1, surf_def_h(0)%ns
3504	n_surf = n_surf + 1
3505	k = surf_def_h(0)%k(m)
3506	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3507	+ SUM ( surf_def_h(0)%frac(m,:) * &
3508	var_def_h0(m,:) ) * conversion_factor(k)
3509	ENDDO
3510	ELSE
3511	n_surf = n_surf + surf_def_h(0)%ns
3512	ENDIF
3513	IF ( ALLOCATED( var_def_h1 ) ) THEN
3514	DO m = 1, surf_def_h(1)%ns
3515	n_surf = n_surf + 1
3516	k = surf_def_h(1)%k(m)
3517	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3518	+ SUM ( surf_def_h(1)%frac(m,:) * &
3519	var_def_h1(m,:) ) * conversion_factor(k)
3520	ENDDO
3521	ELSE
3522	n_surf = n_surf + surf_def_h(1)%ns
3523	ENDIF
3524	IF ( ALLOCATED( var_lsm_h0 ) ) THEN
3525	DO m = 1, surf_lsm_h(0)%ns
3526	n_surf = n_surf + 1
3527	k = surf_lsm_h(0)%k(m)
3528	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3529	+ SUM ( surf_lsm_h(0)%frac(m,:) * &
3530	var_lsm_h0(m,:) ) * conversion_factor(k)
3531	ENDDO
3532	ELSE
3533	n_surf = n_surf + surf_lsm_h(0)%ns
3534	ENDIF
3535	IF ( ALLOCATED( var_lsm_h1 ) ) THEN
3536	DO m = 1, surf_lsm_h(1)%ns
3537	n_surf = n_surf + 1
3538	k = surf_lsm_h(1)%k(m)
3539	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3540	+ SUM ( surf_lsm_h(1)%frac(m,:) * &
3541	var_lsm_h1(m,:) ) * conversion_factor(k)
3542	ENDDO
3543	ELSE
3544	n_surf = n_surf + surf_lsm_h(1)%ns
3545	ENDIF
3546	IF ( ALLOCATED( var_usm_h0 ) ) THEN
3547	DO m = 1, surf_usm_h(0)%ns
3548	n_surf = n_surf + 1
3549	k = surf_usm_h(0)%k(m)
3550	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3551	+ SUM ( surf_usm_h(0)%frac(m,:) * &
3552	var_usm_h0(m,:) ) * conversion_factor(k)
3553	ENDDO
3554	ELSE
3555	n_surf = n_surf + surf_usm_h(0)%ns
3556	ENDIF
3557	IF ( ALLOCATED( var_usm_h1 ) ) THEN
3558	DO m = 1, surf_usm_h(1)%ns
3559	n_surf = n_surf + 1
3560	k = surf_usm_h(1)%k(m)
3561	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3562	+ SUM ( surf_usm_h(1)%frac(m,:) * &
3563	var_usm_h1(m,:) ) * conversion_factor(k)
3564	ENDDO
3565	ELSE
3566	n_surf = n_surf + surf_usm_h(1)%ns
3567	ENDIF
3568	!
3569	!-- Write northward-facing
3570	IF ( ALLOCATED( var_def_v0 ) ) THEN
3571	DO m = 1, surf_def_v(0)%ns
3572	n_surf = n_surf + 1
3573	k = surf_def_v(0)%k(m)
3574	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3575	+ SUM ( surf_def_v(0)%frac(m,:) * &
3576	var_def_v0(m,:) ) * conversion_factor(k)
3577	ENDDO
3578	ELSE
3579	n_surf = n_surf + surf_def_v(0)%ns
3580	ENDIF
3581	IF ( ALLOCATED( var_lsm_v0 ) ) THEN
3582	DO m = 1, surf_lsm_v(0)%ns
3583	n_surf = n_surf + 1
3584	k = surf_lsm_v(0)%k(m)
3585	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3586	+ SUM ( surf_lsm_v(0)%frac(m,:) * &
3587	var_lsm_v0(m,:) ) * conversion_factor(k)
3588	ENDDO
3589	ELSE
3590	n_surf = n_surf + surf_lsm_v(0)%ns
3591	ENDIF
3592	IF ( ALLOCATED( var_usm_v0 ) ) THEN
3593	DO m = 1, surf_usm_v(0)%ns
3594	n_surf = n_surf + 1
3595	k = surf_usm_v(0)%k(m)
3596	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3597	+ SUM ( surf_usm_v(0)%frac(m,:) * &
3598	var_usm_v0(m,:) ) * conversion_factor(k)
3599	ENDDO
3600	ELSE
3601	n_surf = n_surf + surf_usm_v(0)%ns
3602	ENDIF
3603	!
3604	!-- Write southward-facing
3605	IF ( ALLOCATED( var_def_v1 ) ) THEN
3606	DO m = 1, surf_def_v(1)%ns
3607	n_surf = n_surf + 1
3608	k = surf_def_v(1)%k(m)
3609	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3610	+ SUM ( surf_def_v(1)%frac(m,:) * &
3611	var_def_v1(m,:) ) * conversion_factor(k)
3612	ENDDO
3613	ELSE
3614	n_surf = n_surf + surf_def_v(1)%ns
3615	ENDIF
3616	IF ( ALLOCATED( var_lsm_v1 ) ) THEN
3617	DO m = 1, surf_lsm_v(1)%ns
3618	n_surf = n_surf + 1
3619	k = surf_lsm_v(1)%k(m)
3620	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3621	+ SUM ( surf_lsm_v(1)%frac(m,:) * &
3622	var_lsm_v1(m,:) ) * conversion_factor(k)
3623	ENDDO
3624	ELSE
3625	n_surf = n_surf + surf_lsm_v(1)%ns
3626	ENDIF
3627	IF ( ALLOCATED( var_usm_v1 ) ) THEN
3628	DO m = 1, surf_usm_v(1)%ns
3629	n_surf = n_surf + 1
3630	k = surf_usm_v(1)%k(m)
3631	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3632	+ SUM ( surf_usm_v(1)%frac(m,:) * &
3633	var_usm_v1(m,:) ) * conversion_factor(k)
3634	ENDDO
3635	ELSE
3636	n_surf = n_surf + surf_usm_v(1)%ns
3637	ENDIF
3638	!
3639	!-- Write eastward-facing
3640	IF ( ALLOCATED( var_def_v2 ) ) THEN
3641	DO m = 1, surf_def_v(2)%ns
3642	n_surf = n_surf + 1
3643	k = surf_def_v(2)%k(m)
3644	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3645	+ SUM ( surf_def_v(2)%frac(m,:) * &
3646	var_def_v2(m,:) ) * conversion_factor(k)
3647	ENDDO
3648	ELSE
3649	n_surf = n_surf + surf_def_v(2)%ns
3650	ENDIF
3651	IF ( ALLOCATED( var_lsm_v2 ) ) THEN
3652	DO m = 1, surf_lsm_v(2)%ns
3653	n_surf = n_surf + 1
3654	k = surf_lsm_v(2)%k(m)
3655	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3656	+ SUM ( surf_lsm_v(2)%frac(m,:) * &
3657	var_lsm_v2(m,:) ) * conversion_factor(k)
3658	ENDDO
3659	ELSE
3660	n_surf = n_surf + surf_lsm_v(2)%ns
3661	ENDIF
3662	IF ( ALLOCATED( var_usm_v2 ) ) THEN
3663	DO m = 1, surf_usm_v(2)%ns
3664	n_surf = n_surf + 1
3665	k = surf_usm_v(2)%k(m)
3666	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3667	+ SUM ( surf_usm_v(2)%frac(m,:) * &
3668	var_usm_v2(m,:) ) * conversion_factor(k)
3669	ENDDO
3670	ELSE
3671	n_surf = n_surf + surf_usm_v(2)%ns
3672	ENDIF
3673	!
3674	!-- Write westward-facing
3675	IF ( ALLOCATED( var_def_v3 ) ) THEN
3676	DO m = 1, surf_def_v(3)%ns
3677	n_surf = n_surf + 1
3678	k = surf_def_v(3)%k(m)
3679	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3680	+ SUM ( surf_def_v(3)%frac(m,:) * &
3681	var_def_v3(m,:) ) * conversion_factor(k)
3682	ENDDO
3683	ELSE
3684	n_surf = n_surf + surf_def_v(3)%ns
3685	ENDIF
3686	IF ( ALLOCATED( var_lsm_v3 ) ) THEN
3687	DO m = 1, surf_lsm_v(3)%ns
3688	n_surf = n_surf + 1
3689	k = surf_lsm_v(3)%k(m)
3690	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3691	+ SUM ( surf_lsm_v(3)%frac(m,:) * &
3692	var_lsm_v3(m,:) ) * conversion_factor(k)
3693	ENDDO
3694	ELSE
3695	n_surf = n_surf + surf_lsm_v(3)%ns
3696	ENDIF
3697	IF ( ALLOCATED( var_usm_v3 ) ) THEN
3698	DO m = 1, surf_usm_v(3)%ns
3699	n_surf = n_surf + 1
3700	k = surf_usm_v(3)%k(m)
3701	surfaces%var_av(n_surf,n_out) = surfaces%var_av(n_surf,n_out) &
3702	+ SUM ( surf_usm_v(3)%frac(m,:) * &
3703	var_usm_v3(m,:) ) * conversion_factor(k)
3704	ENDDO
3705	ELSE
3706	n_surf = n_surf + surf_usm_v(3)%ns
3707	ENDIF
3708
3709	END SUBROUTINE surface_data_output_sum_up_2d
3710
3711	!--------------------------------------------------------------------------------------------------!
3712	! Description:
3713	! ------------
3714	!> Collect the surface data from different types and different orientation.
3715	!--------------------------------------------------------------------------------------------------!
3716	SUBROUTINE surface_data_output_collect_1d( var_def_h0, var_lsm_h0, var_usm_h0, &
3717	var_def_h1, var_lsm_h1, var_usm_h1, &
3718	var_def_v0, var_lsm_v0, var_usm_v0, &
3719	var_def_v1, var_lsm_v1, var_usm_v1, &
3720	var_def_v2, var_lsm_v2, var_usm_v2, &
3721	var_def_v3, var_lsm_v3, var_usm_v3, fac )
3722
3723	IMPLICIT NONE
3724
3725	INTEGER(iwp) :: k !< height index of surface element
3726	INTEGER(iwp) :: m !< running index for surface elements
3727	INTEGER(iwp) :: n_surf !< running index for surface elements
3728
3729	REAL(wp), DIMENSION(:), OPTIONAL :: fac !< passed output conversion factor for heatflux output
3730	REAL(wp), DIMENSION(nzb:nzt+1) :: conversion_factor !< effective array for output conversion factor
3731
3732	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h0 !< output variable at upward-facing default-type surfaces
3733	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_h1 !< output variable at downward-facing default-type surfaces
3734	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_h0 !< output variable at upward-facing natural-type surfaces
3735	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_h1 !< output variable at downward-facing natural-type surfaces
3736	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_h0 !< output variable at upward-facing urban-type surfaces
3737	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_h1 !< output variable at downward-facing urban-type surfaces
3738	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v0 !< output variable at northward-facing default-type surfaces
3739	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v1 !< output variable at southward-facing default-type surfaces
3740	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v2 !< output variable at eastward-facing default-type surfaces
3741	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_def_v3 !< output variable at westward-facing default-type surfaces
3742	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v0 !< output variable at northward-facing natural-type surfaces
3743	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v1 !< output variable at southward-facing natural-type surfaces
3744	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v2 !< output variable at eastward-facing natural-type surfaces
3745	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_lsm_v3 !< output variable at westward-facing natural-type surfaces
3746	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v0 !< output variable at northward-facing urban-type surfaces
3747	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v1 !< output variable at southward-facing urban-type surfaces
3748	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v2 !< output variable at eastward-facing urban-type surfaces
3749	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: var_usm_v3 !< output variable at westward-facing urban-type surfaces
3750
3751	!
3752	!-- Set conversion factor to one if not present
3753	IF ( .NOT. PRESENT( fac ) ) THEN
3754	conversion_factor = 1.0_wp
3755	ELSE
3756	conversion_factor = fac
3757	ENDIF
3758	!
3759	!-- Set counter variable to zero before the variable is written to the output array.
3760	n_surf = 0
3761	!
3762	!-- Write the horizontal surfaces.
3763	!-- Before each variable is written to the output data structure, first check if the variable
3764	!-- for the respective surface type is defined. If a variable is not defined, skip the block and
3765	!-- increment the counter variable by the number of surface elements of this type. Usually this is
3766	!-- zero, however, there might be the situation that e.g. urban surfaces are defined but the
3767	!-- respective variable is not allocated for this surface type. To write the data on the exact
3768	!-- position, increment the counter.
3769	IF ( ALLOCATED( var_def_h0 ) ) THEN
3770	DO m = 1, surf_def_h(0)%ns
3771	n_surf = n_surf + 1
3772	k = surf_def_h(0)%k(m)
3773	surfaces%var_out(n_surf) = var_def_h0(m) * conversion_factor(k)
3774	ENDDO
3775	ELSE
3776	n_surf = n_surf + surf_def_h(0)%ns
3777	ENDIF
3778	IF ( ALLOCATED( var_def_h1 ) ) THEN
3779	DO m = 1, surf_def_h(1)%ns
3780	n_surf = n_surf + 1
3781	k = surf_def_h(1)%k(m)
3782	surfaces%var_out(n_surf) = var_def_h1(m) * conversion_factor(k)
3783	ENDDO
3784	ELSE
3785	n_surf = n_surf + surf_def_h(1)%ns
3786	ENDIF
3787	IF ( ALLOCATED( var_lsm_h0 ) ) THEN
3788	DO m = 1, surf_lsm_h(0)%ns
3789	n_surf = n_surf + 1
3790	k = surf_lsm_h(0)%k(m)
3791	surfaces%var_out(n_surf) = var_lsm_h0(m) * conversion_factor(k)
3792	ENDDO
3793	ELSE
3794	n_surf = n_surf + surf_lsm_h(0)%ns
3795	ENDIF
3796	IF ( ALLOCATED( var_lsm_h1 ) ) THEN
3797	DO m = 1, surf_lsm_h(1)%ns
3798	n_surf = n_surf + 1
3799	k = surf_lsm_h(1)%k(m)
3800	surfaces%var_out(n_surf) = var_lsm_h1(m) * conversion_factor(k)
3801	ENDDO
3802	ELSE
3803	n_surf = n_surf + surf_lsm_h(1)%ns
3804	ENDIF
3805	IF ( ALLOCATED( var_usm_h0 ) ) THEN
3806	DO m = 1, surf_usm_h(0)%ns
3807	n_surf = n_surf + 1
3808	k = surf_usm_h(0)%k(m)
3809	surfaces%var_out(n_surf) = var_usm_h0(m) * conversion_factor(k)
3810	ENDDO
3811	ELSE
3812	n_surf = n_surf + surf_usm_h(0)%ns
3813	ENDIF
3814	IF ( ALLOCATED( var_usm_h1 ) ) THEN
3815	DO m = 1, surf_usm_h(1)%ns
3816	n_surf = n_surf + 1
3817	k = surf_usm_h(1)%k(m)
3818	surfaces%var_out(n_surf) = var_usm_h1(m) * conversion_factor(k)
3819	ENDDO
3820	ELSE
3821	n_surf = n_surf + surf_usm_h(1)%ns
3822	ENDIF
3823	!
3824	!-- Write northward-facing
3825	IF ( ALLOCATED( var_def_v0 ) ) THEN
3826	DO m = 1, surf_def_v(0)%ns
3827	n_surf = n_surf + 1
3828	k = surf_def_v(0)%k(m)
3829	surfaces%var_out(n_surf) = var_def_v0(m) * conversion_factor(k)
3830	ENDDO
3831	ELSE
3832	n_surf = n_surf + surf_def_v(0)%ns
3833	ENDIF
3834	IF ( ALLOCATED( var_lsm_v0 ) ) THEN
3835	DO m = 1, surf_lsm_v(0)%ns
3836	n_surf = n_surf + 1
3837	k = surf_lsm_v(0)%k(m)
3838	surfaces%var_out(n_surf) = var_lsm_v0(m) * conversion_factor(k)
3839	ENDDO
3840	ELSE
3841	n_surf = n_surf + surf_lsm_v(0)%ns
3842	ENDIF
3843	IF ( ALLOCATED( var_usm_v0 ) ) THEN
3844	DO m = 1, surf_usm_v(0)%ns
3845	n_surf = n_surf + 1
3846	k = surf_usm_v(0)%k(m)
3847	surfaces%var_out(n_surf) = var_usm_v0(m) * conversion_factor(k)
3848	ENDDO
3849	ELSE
3850	n_surf = n_surf + surf_usm_v(0)%ns
3851	ENDIF
3852	!
3853	!-- Write southward-facing
3854	IF ( ALLOCATED( var_def_v1 ) ) THEN
3855	DO m = 1, surf_def_v(1)%ns
3856	n_surf = n_surf + 1
3857	k = surf_def_v(1)%k(m)
3858	surfaces%var_out(n_surf) = var_def_v1(m) * conversion_factor(k)
3859	ENDDO
3860	ELSE
3861	n_surf = n_surf + surf_def_v(1)%ns
3862	ENDIF
3863	IF ( ALLOCATED( var_lsm_v1 ) ) THEN
3864	DO m = 1, surf_lsm_v(1)%ns
3865	n_surf = n_surf + 1
3866	k = surf_lsm_v(1)%k(m)
3867	surfaces%var_out(n_surf) = var_lsm_v1(m) * conversion_factor(k)
3868	ENDDO
3869	ELSE
3870	n_surf = n_surf + surf_lsm_v(1)%ns
3871	ENDIF
3872	IF ( ALLOCATED( var_usm_v1 ) ) THEN
3873	DO m = 1, surf_usm_v(1)%ns
3874	n_surf = n_surf + 1
3875	k = surf_usm_v(1)%k(m)
3876	surfaces%var_out(n_surf) = var_usm_v1(m) * conversion_factor(k)
3877	ENDDO
3878	ELSE
3879	n_surf = n_surf + surf_usm_v(1)%ns
3880	ENDIF
3881	!
3882	!-- Write eastward-facing
3883	IF ( ALLOCATED( var_def_v2 ) ) THEN
3884	DO m = 1, surf_def_v(2)%ns
3885	n_surf = n_surf + 1
3886	k = surf_def_v(2)%k(m)
3887	surfaces%var_out(n_surf) = var_def_v2(m) * conversion_factor(k)
3888	ENDDO
3889	ELSE
3890	n_surf = n_surf + surf_def_v(2)%ns
3891	ENDIF
3892	IF ( ALLOCATED( var_lsm_v2 ) ) THEN
3893	DO m = 1, surf_lsm_v(2)%ns
3894	n_surf = n_surf + 1
3895	k = surf_lsm_v(2)%k(m)
3896	surfaces%var_out(n_surf) = var_lsm_v2(m) * conversion_factor(k)
3897	ENDDO
3898	ELSE
3899	n_surf = n_surf + surf_lsm_v(2)%ns
3900	ENDIF
3901	IF ( ALLOCATED( var_usm_v2 ) ) THEN
3902	DO m = 1, surf_usm_v(2)%ns
3903	n_surf = n_surf + 1
3904	k = surf_usm_v(2)%k(m)
3905	surfaces%var_out(n_surf) = var_usm_v2(m) * conversion_factor(k)
3906	ENDDO
3907	ELSE
3908	n_surf = n_surf + surf_usm_v(2)%ns
3909	ENDIF
3910	!
3911	!-- Write westward-facing
3912	IF ( ALLOCATED( var_def_v3 ) ) THEN
3913	DO m = 1, surf_def_v(3)%ns
3914	n_surf = n_surf + 1
3915	k = surf_def_v(3)%k(m)
3916	surfaces%var_out(n_surf) = var_def_v3(m) * conversion_factor(k)
3917	ENDDO
3918	ELSE
3919	n_surf = n_surf + surf_def_v(3)%ns
3920	ENDIF
3921	IF ( ALLOCATED( var_lsm_v3 ) ) THEN
3922	DO m = 1, surf_lsm_v(3)%ns
3923	n_surf = n_surf + 1
3924	k = surf_lsm_v(3)%k(m)
3925	surfaces%var_out(n_surf) = var_lsm_v3(m) * conversion_factor(k)
3926	ENDDO
3927	ELSE
3928	n_surf = n_surf + surf_lsm_v(3)%ns
3929	ENDIF
3930	IF ( ALLOCATED( var_usm_v3 ) ) THEN
3931	DO m = 1, surf_usm_v(3)%ns
3932	n_surf = n_surf + 1
3933	k = surf_usm_v(3)%k(m)
3934	surfaces%var_out(n_surf) = var_usm_v3(m) * conversion_factor(k)
3935	ENDDO
3936	ELSE
3937	n_surf = n_surf + surf_usm_v(3)%ns
3938	ENDIF
3939
3940	END SUBROUTINE surface_data_output_collect_1d
3941
3942	!--------------------------------------------------------------------------------------------------!
3943	! Description:
3944	! ------------
3945	!> Collect the surface data from different types and different orientation for properties which are
3946	!> defined using tile approach.
3947	!--------------------------------------------------------------------------------------------------!
3948	SUBROUTINE surface_data_output_collect_2d( var_def_h0, var_lsm_h0, var_usm_h0, &
3949	var_def_h1, var_lsm_h1, var_usm_h1, &
3950	var_def_v0, var_lsm_v0, var_usm_v0, &
3951	var_def_v1, var_lsm_v1, var_usm_v1, &
3952	var_def_v2, var_lsm_v2, var_usm_v2, &
3953	var_def_v3, var_lsm_v3, var_usm_v3, fac )
3954
3955	IMPLICIT NONE
3956
3957	INTEGER(iwp) :: k !< height index of surface element
3958	INTEGER(iwp) :: m !< running index for surface elements
3959	INTEGER(iwp) :: n_surf !< running index for surface elements
3960
3961	REAL(wp), DIMENSION(:), OPTIONAL :: fac !< passed output conversion factor for heatflux output
3962	REAL(wp), DIMENSION(nzb:nzt+1) :: conversion_factor !< effective array for output conversion factor
3963
3964	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_h0 !< output variable at upward-facing default-type surfaces
3965	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_h1 !< output variable at downward-facing default-type surfaces
3966	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_h0 !< output variable at upward-facing natural-type surfaces
3967	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_h1 !< output variable at downward-facing natural-type surfaces
3968	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_h0 !< output variable at upward-facing urban-type surfaces
3969	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_h1 !< output variable at downward-facing urban-type surfaces
3970	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v0 !< output variable at northward-facing default-type surfaces
3971	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v1 !< output variable at southward-facing default-type surfaces
3972	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v2 !< output variable at eastward-facing default-type surfaces
3973	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_def_v3 !< output variable at westward-facing default-type surfaces
3974	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v0 !< output variable at northward-facing natural-type surfaces
3975	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v1 !< output variable at southward-facing natural-type surfaces
3976	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v2 !< output variable at eastward-facing natural-type surfaces
3977	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_lsm_v3 !< output variable at westward-facing natural-type surfaces
3978	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v0 !< output variable at northward-facing urban-type surfaces
3979	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v1 !< output variable at southward-facing urban-type surfaces
3980	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v2 !< output variable at eastward-facing urban-type surfaces
3981	REAL(wp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN) :: var_usm_v3 !< output variable at westward-facing urban-type surfaces
3982
3983	!
3984	!-- Set conversion factor to one if not present
3985	IF ( .NOT. PRESENT( fac ) ) THEN
3986	conversion_factor = 1.0_wp
3987	ELSE
3988	conversion_factor = fac
3989	ENDIF
3990	!
3991	!-- Set counter variable to zero before the variable is written to the output array.
3992	n_surf = 0
3993	!
3994	!-- Write the horizontal surfaces.
3995	!-- Before each variable is written to the output data structure, first check if the variable
3996	!-- for the respective surface type is defined. If a variable is not defined, skip the block and
3997	!-- increment the counter variable by the number of surface elements of this type. Usually this is
3998	!-- zero, however, there might be the situation that e.g. urban surfaces are defined but the
3999	!-- respective variable is not allocated for this surface type. To write the data on the exact
4000	!-- position, increment the counter.
4001	IF ( ALLOCATED( var_def_h0 ) ) THEN
4002	DO m = 1, surf_def_h(0)%ns
4003	n_surf = n_surf + 1
4004	k = surf_def_h(0)%k(m)
4005	surfaces%var_out(n_surf) = SUM ( surf_def_h(0)%frac(m,:) * var_def_h0(m,:) ) * &
4006	conversion_factor(k)
4007	ENDDO
4008	ELSE
4009	n_surf = n_surf + surf_def_h(0)%ns
4010	ENDIF
4011	IF ( ALLOCATED( var_def_h1 ) ) THEN
4012	DO m = 1, surf_def_h(1)%ns
4013	n_surf = n_surf + 1
4014	k = surf_def_h(1)%k(m)
4015	surfaces%var_out(n_surf) = SUM ( surf_def_h(1)%frac(m,:) * var_def_h1(m,:) ) * &
4016	conversion_factor(k)
4017	ENDDO
4018	ELSE
4019	n_surf = n_surf + surf_def_h(1)%ns
4020	ENDIF
4021	IF ( ALLOCATED( var_lsm_h0 ) ) THEN
4022	DO m = 1, surf_lsm_h(0)%ns
4023	n_surf = n_surf + 1
4024	k = surf_lsm_h(0)%k(m)
4025	surfaces%var_out(n_surf) = SUM ( surf_lsm_h(0)%frac(m,:) * var_lsm_h0(m,:) ) * &
4026	conversion_factor(k)
4027	ENDDO
4028	ELSE
4029	n_surf = n_surf + surf_lsm_h(0)%ns
4030	ENDIF
4031	IF ( ALLOCATED( var_lsm_h1 ) ) THEN
4032	DO m = 1, surf_lsm_h(1)%ns
4033	n_surf = n_surf + 1
4034	k = surf_lsm_h(1)%k(m)
4035	surfaces%var_out(n_surf) = SUM ( surf_lsm_h(1)%frac(m,:) * var_lsm_h1(m,:) ) * &
4036	conversion_factor(k)
4037	ENDDO
4038	ELSE
4039	n_surf = n_surf + surf_lsm_h(0)%ns
4040	ENDIF
4041	IF ( ALLOCATED( var_usm_h0 ) ) THEN
4042	DO m = 1, surf_usm_h(0)%ns
4043	n_surf = n_surf + 1
4044	k = surf_usm_h(0)%k(m)
4045	surfaces%var_out(n_surf) = SUM ( surf_usm_h(0)%frac(m,:) * var_usm_h0(m,:) ) &
4046	* conversion_factor(k)
4047	ENDDO
4048	ELSE
4049	n_surf = n_surf + surf_usm_h(0)%ns
4050	ENDIF
4051	IF ( ALLOCATED( var_usm_h1 ) ) THEN
4052	DO m = 1, surf_usm_h(1)%ns
4053	n_surf = n_surf + 1
4054	k = surf_usm_h(1)%k(m)
4055	surfaces%var_out(n_surf) = SUM ( surf_usm_h(1)%frac(m,:) * var_usm_h1(m,:) ) &
4056	* conversion_factor(k)
4057	ENDDO
4058	ELSE
4059	n_surf = n_surf + surf_usm_h(1)%ns
4060	ENDIF
4061	!
4062	!-- Write northward-facing
4063	IF ( ALLOCATED( var_def_v0 ) ) THEN
4064	DO m = 1, surf_def_v(0)%ns
4065	n_surf = n_surf + 1
4066	k = surf_def_v(0)%k(m)
4067	surfaces%var_out(n_surf) = SUM ( surf_def_v(0)%frac(m,:) * var_def_v0(m,:) ) * &
4068	conversion_factor(k)
4069	ENDDO
4070	ELSE
4071	n_surf = n_surf + surf_def_v(0)%ns
4072	ENDIF
4073	IF ( ALLOCATED( var_lsm_v0 ) ) THEN
4074	DO m = 1, surf_lsm_v(0)%ns
4075	n_surf = n_surf + 1
4076	k = surf_lsm_v(0)%k(m)
4077	surfaces%var_out(n_surf) = SUM ( surf_lsm_v(0)%frac(m,:) * var_lsm_v0(m,:) ) * &
4078	conversion_factor(k)
4079	ENDDO
4080	ELSE
4081	n_surf = n_surf + surf_lsm_v(0)%ns
4082	ENDIF
4083	IF ( ALLOCATED( var_usm_v0 ) ) THEN
4084	DO m = 1, surf_usm_v(0)%ns
4085	n_surf = n_surf + 1
4086	k = surf_usm_v(0)%k(m)
4087	surfaces%var_out(n_surf) = SUM ( surf_usm_v(0)%frac(m,:) * var_usm_v0(m,:) ) * &
4088	conversion_factor(k)
4089	ENDDO
4090	ELSE
4091	n_surf = n_surf + surf_usm_v(0)%ns
4092	ENDIF
4093	!
4094	!-- Write southward-facing
4095	IF ( ALLOCATED( var_def_v1 ) ) THEN
4096	DO m = 1, surf_def_v(1)%ns
4097	n_surf = n_surf + 1
4098	k = surf_def_v(1)%k(m)
4099	surfaces%var_out(n_surf) = SUM ( surf_def_v(1)%frac(m,:) * var_def_v1(m,:) ) * &
4100	conversion_factor(k)
4101	ENDDO
4102	ELSE
4103	n_surf = n_surf + surf_def_v(1)%ns
4104	ENDIF
4105	IF ( ALLOCATED( var_lsm_v1 ) ) THEN
4106	DO m = 1, surf_lsm_v(1)%ns
4107	n_surf = n_surf + 1
4108	k = surf_lsm_v(1)%k(m)
4109	surfaces%var_out(n_surf) = SUM ( surf_lsm_v(1)%frac(m,:) * var_lsm_v1(m,:) ) * &
4110	conversion_factor(k)
4111	ENDDO
4112	ELSE
4113	n_surf = n_surf + surf_lsm_v(1)%ns
4114	ENDIF
4115	IF ( ALLOCATED( var_usm_v1 ) ) THEN
4116	DO m = 1, surf_usm_v(1)%ns
4117	n_surf = n_surf + 1
4118	k = surf_usm_v(1)%k(m)
4119	surfaces%var_out(n_surf) = SUM ( surf_usm_v(1)%frac(m,:) * var_usm_v1(m,:) ) * &
4120	conversion_factor(k)
4121	ENDDO
4122	ELSE
4123	n_surf = n_surf + surf_usm_v(1)%ns
4124	ENDIF
4125	!
4126	!-- Write eastward-facing
4127	IF ( ALLOCATED( var_def_v2 ) ) THEN
4128	DO m = 1, surf_def_v(2)%ns
4129	n_surf = n_surf + 1
4130	k = surf_def_v(2)%k(m)
4131	surfaces%var_out(n_surf) = SUM ( surf_def_v(2)%frac(m,:) * var_def_v2(m,:) ) * &
4132	conversion_factor(k)
4133	ENDDO
4134	ELSE
4135	n_surf = n_surf + surf_def_v(2)%ns
4136	ENDIF
4137	IF ( ALLOCATED( var_lsm_v2 ) ) THEN
4138	DO m = 1, surf_lsm_v(2)%ns
4139	n_surf = n_surf + 1
4140	k = surf_lsm_v(2)%k(m)
4141	surfaces%var_out(n_surf) = SUM ( surf_lsm_v(2)%frac(m,:) * var_lsm_v2(m,:) ) * &
4142	conversion_factor(k)
4143	ENDDO
4144	ELSE
4145	n_surf = n_surf + surf_lsm_v(2)%ns
4146	ENDIF
4147	IF ( ALLOCATED( var_usm_v2 ) ) THEN
4148	DO m = 1, surf_usm_v(2)%ns
4149	n_surf = n_surf + 1
4150	k = surf_usm_v(2)%k(m)
4151	surfaces%var_out(n_surf) = SUM ( surf_usm_v(2)%frac(m,:) * var_usm_v2(m,:) ) * &
4152	conversion_factor(k)
4153	ENDDO
4154	ELSE
4155	n_surf = n_surf + surf_usm_v(2)%ns
4156	ENDIF
4157	!
4158	!-- Write westward-facing
4159	IF ( ALLOCATED( var_def_v3 ) ) THEN
4160	DO m = 1, surf_def_v(3)%ns
4161	n_surf = n_surf + 1
4162	k = surf_def_v(3)%k(m)
4163	surfaces%var_out(n_surf) = SUM ( surf_def_v(3)%frac(m,:) * var_def_v3(m,:) ) * &
4164	conversion_factor(k)
4165	ENDDO
4166	ELSE
4167	n_surf = n_surf + surf_def_v(3)%ns
4168	ENDIF
4169	IF ( ALLOCATED( var_lsm_v3 ) ) THEN
4170	DO m = 1, surf_lsm_v(3)%ns
4171	n_surf = n_surf + 1
4172	k = surf_lsm_v(3)%k(m)
4173	surfaces%var_out(n_surf) = SUM ( surf_lsm_v(3)%frac(m,:) * var_lsm_v3(m,:) ) * &
4174	conversion_factor(k)
4175	ENDDO
4176	ELSE
4177	n_surf = n_surf + surf_lsm_v(3)%ns
4178	ENDIF
4179	IF ( ALLOCATED( var_usm_v3 ) ) THEN
4180	DO m = 1, surf_usm_v(3)%ns
4181	n_surf = n_surf + 1
4182	k = surf_usm_v(3)%k(m)
4183	surfaces%var_out(n_surf) = SUM ( surf_usm_v(3)%frac(m,:) * var_usm_v3(m,:) ) * &
4184	conversion_factor(k)
4185	ENDDO
4186	ELSE
4187	n_surf = n_surf + surf_usm_v(3)%ns
4188	ENDIF
4189
4190	END SUBROUTINE surface_data_output_collect_2d
4191
4192	!--------------------------------------------------------------------------------------------------!
4193	! Description:
4194	! ------------
4195	!> Parin for output of surface parameters
4196	!--------------------------------------------------------------------------------------------------!
4197	SUBROUTINE surface_data_output_parin
4198
4199	IMPLICIT NONE
4200
4201	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
4202
4203
4204	NAMELIST /surface_data_output_parameters/ averaging_interval_surf, data_output_surf, &
4205	dt_dosurf, dt_dosurf_av, skip_time_dosurf, &
4206	skip_time_dosurf_av, to_netcdf, to_vtk
4207
4208	line = ' '
4209
4210	!
4211	!-- Try to find the namelist
4212	REWIND ( 11 )
4213	line = ' '
4214	DO WHILE ( INDEX( line, '&surface_data_output_parameters' ) == 0 )
4215	READ ( 11, '(A)', END=14 ) line
4216	ENDDO
4217	BACKSPACE ( 11 )
4218
4219	!
4220	!-- Read namelist
4221	READ ( 11, surface_data_output_parameters, ERR = 10 )
4222	!
4223	!-- Set flag that indicates that surface data output is switched on
4224	surface_output = .TRUE.
4225	GOTO 14
4226
4227	10 BACKSPACE( 11 )
4228	READ( 11 , '(A)') line
4229	CALL parin_fail_message( 'surface_data_output_parameters', line )
4230
4231	14 CONTINUE
4232
4233
4234	END SUBROUTINE surface_data_output_parin
4235
4236
4237	!--------------------------------------------------------------------------------------------------!
4238	! Description:
4239	! ------------
4240	!> Check the input parameters for consistency. Further pre-process the given output variables, i.e.
4241	!> separate them into average and non-average output variables and map them onto internal output
4242	!> array.
4243	!--------------------------------------------------------------------------------------------------!
4244	SUBROUTINE surface_data_output_check_parameters
4245
4246	USE control_parameters, &
4247	ONLY: averaging_interval, &
4248	dt_data_output, &
4249	indoor_model, &
4250	initializing_actions, &
4251	message_string
4252
4253	USE pegrid, &
4254	ONLY: numprocs_previous_run
4255
4256	IMPLICIT NONE
4257
4258	CHARACTER(LEN=100) :: trimvar !< dummy for single output variable
4259	CHARACTER(LEN=100) :: unit !< dummy for unit of output variable
4260
4261	INTEGER(iwp) :: av !< id indicating average or non-average data output
4262	INTEGER(iwp) :: ilen !< string length
4263	INTEGER(iwp) :: n_out !< running index for number of output variables
4264	!
4265	!-- Check if any output file type is selected
4266	IF ( .NOT. to_vtk .AND. .NOT. to_netcdf ) THEN
4267	WRITE( message_string, * ) 'no output file type selected for surface-data output!&' // &
4268	'Set at least either "to_vtk" or "to_netcdf" to .TRUE.'
4269	CALL message( 'surface_data_output_check_parameters', 'PA0662', 1, 2, 0, 6, 0 )
4270	ENDIF
4271	!
4272	!-- Check the average interval
4273	IF ( averaging_interval_surf == 9999999.9_wp ) THEN
4274	averaging_interval_surf = averaging_interval
4275	ENDIF
4276	!
4277	!-- Set the default data-output interval dt_data_output if necessary
4278	IF ( dt_dosurf == 9999999.9_wp ) dt_dosurf = dt_data_output
4279	IF ( dt_dosurf_av == 9999999.9_wp ) dt_dosurf_av = dt_data_output
4280
4281	IF ( averaging_interval_surf > dt_dosurf_av ) THEN
4282	WRITE( message_string, * ) 'averaging_interval_surf = ', averaging_interval_surf, &
4283	' must be <= dt_dosurf_av = ', dt_dosurf_av
4284	CALL message( 'surface_data_output_check_parameters', 'PA0536', 1, 2, 0, 6, 0 )
4285	ENDIF
4286
4287	#if ! defined( __netcdf4_parallel )
4288	!
4289	!-- Surface output via NetCDF requires parallel NetCDF
4290	IF ( to_netcdf ) THEN
4291	message_string = 'to_netcdf = .True. requires parallel NetCDF'
4292	CALL message( 'surface_data_output_check_parameters', 'PA0116', 1, 2, 0, 6, 0 )
4293	ENDIF
4294	#endif
4295	!
4296	!-- In case of parallel NetCDF output the output timestep must not be zero. This is because the
4297	!-- number of requiered output timesteps is pre-calculated, which is not possible with zero output
4298	!-- timestep.
4299	IF ( netcdf_data_format > 4 ) THEN
4300	IF ( dt_dosurf == 0.0_wp ) THEN
4301	message_string = 'dt_dosurf = 0.0 while using a variable timestep and parallel ' // &
4302	'netCDF4 is not allowed.'
4303	CALL message( 'surface_data_output_check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
4304	ENDIF
4305
4306	IF ( dt_dosurf_av == 0.0_wp ) THEN
4307	message_string = 'dt_dosurf_av = 0.0 while using a variable timestep and parallel ' // &
4308	'netCDF4 is not allowed.'
4309	CALL message( 'surface_data_output_check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
4310	ENDIF
4311	ENDIF
4312
4313	!
4314	!-- In case of restart runs, check it the number of cores has been changed.
4315	!-- With surface output this is not allowed.
4316	IF ( TRIM( initializing_actions ) == 'read_restart_data' .AND. &
4317	numprocs_previous_run /= numprocs ) THEN
4318	message_string = 'The number of cores has been changed between restart runs. ' // &
4319	'This is not allowed when surface data output is used.'
4320	CALL message( 'surface_data_output_check_parameters', 'PA0585', 1, 2, 0, 6, 0 )
4321	ENDIF
4322	!
4323	!-- Count number of output variables and separate output strings for average and non-average output
4324	!-- variables.
4325	n_out = 0
4326	DO WHILE ( data_output_surf(n_out+1)(1:1) /= ' ' )
4327
4328	n_out = n_out + 1
4329	ilen = LEN_TRIM( data_output_surf(n_out) )
4330	trimvar = TRIM( data_output_surf(n_out) )
4331
4332	!
4333	!-- Check for data averaging
4334	av = 0
4335	IF ( ilen > 3 ) THEN
4336	IF ( data_output_surf(n_out)(ilen-2:ilen) == '_av' ) THEN
4337	trimvar = data_output_surf(n_out)(1:ilen-3)
4338	av = 1
4339	ENDIF
4340	ENDIF
4341
4342	dosurf_no(av) = dosurf_no(av) + 1
4343	dosurf(av,dosurf_no(av)) = TRIM( trimvar )
4344
4345	!
4346	!-- Check if all output variables are known and assign a unit
4347	unit = 'not set'
4348	SELECT CASE ( TRIM( trimvar ) )
4349
4350	CASE ( 'css', 'cssws', 'qsws_liq', 'qsws_soil', 'qsws_veg' )
4351	message_string = TRIM( trimvar ) // ' is not yet implemented in the surface output'
4352	CALL message( 'surface_data_output_check_parameters', 'PA0537', 1, 2, 0, 6, 0 )
4353
4354	CASE ( 'us', 'uvw1' )
4355	unit = 'm/s'
4356
4357	CASE ( 'ss', 'qcs', 'ncs', 'qis', 'nis', 'qrs', 'nrs' )
4358	unit = '1'
4359
4360	CASE ( 'z0', 'z0h', 'z0q', 'ol' )
4361	unit = 'm'
4362
4363	CASE ( 'ts', 'theta1', 'thetav1', 'theta_surface', 'thetav_surface' )
4364	unit = 'K'
4365
4366	CASE ( 'usws', 'vsws' )
4367	unit = 'm2/s2'
4368
4369	CASE ( 'qcsws', 'ncsws', 'qisws', 'nisws', 'qrsws', 'nrsws', 'sasws' )
4370
4371	CASE ( 'shf' )
4372	unit = 'K m/s'
4373
4374	CASE ( 'qsws' )
4375	unit = 'kg/kg m/s'
4376
4377	CASE ( 'ssws' )
4378	unit = 'kg/m2/s'
4379
4380	CASE ( 'qs', 'q_surface', 'qv1' )
4381	unit = 'kg/kg'
4382
4383	CASE ( 'rad_net' )
4384	unit = 'W/m2'
4385
4386	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_dif', 'rad_lw_ref', 'rad_lw_res' )
4387	unit = 'W/m2'
4388
4389	CASE ( 'rad_sw_in', 'rad_sw_out', 'rad_sw_dif', 'rad_sw_ref', 'rad_sw_res', 'rad_sw_dir' )
4390	unit = 'W/m2'
4391
4392	CASE ( 'ghf' )
4393	unit = 'W/m2'
4394
4395	CASE ( 'r_a', 'r_canopy', 'r_soil', 'r_s' )
4396	unit = 's/m'
4397
4398	CASE ( 'waste_heat', 'im_hf' )
4399	IF ( .NOT. indoor_model ) THEN
4400	message_string = TRIM( trimvar ) // ' requires the indoor model'
4401	CALL message( 'surface_data_output_check_parameters', 'PA0588', 1, 2, 0, 6, 0 )
4402	ENDIF
4403
4404	unit = 'W/m2'
4405
4406	CASE ( 'albedo', 'emissivity' )
4407	unit = '1'
4408
4409	CASE DEFAULT
4410	message_string = TRIM( trimvar ) // ' is not part of the surface output'
4411	CALL message( 'surface_data_output_check_parameters', 'PA0538', 1, 2, 0, 6, 0 )
4412	END SELECT
4413
4414	dosurf_unit(av,dosurf_no(av)) = unit
4415
4416	ENDDO
4417
4418	END SUBROUTINE surface_data_output_check_parameters
4419
4420
4421	!--------------------------------------------------------------------------------------------------!
4422	! Description:
4423	! ------------
4424	!> Last action.
4425	!--------------------------------------------------------------------------------------------------!
4426	SUBROUTINE surface_data_output_last_action( av )
4427
4428	USE control_parameters, &
4429	ONLY: io_blocks, &
4430	io_group
4431
4432	#if defined( __parallel )
4433	USE pegrid, &
4434	ONLY: comm2d, &
4435	ierr
4436	#endif
4437
4438	IMPLICIT NONE
4439
4440	INTEGER(iwp) :: av !< id indicating average or non-average data output
4441	INTEGER(iwp) :: i !< loop index
4442
4443	!
4444	!--Return, if nothing to output
4445	IF ( dosurf_no(av) == 0 ) RETURN
4446	!
4447	!--If output to VTK files is enabled, check if files are open and write an end-of-file statement.
4448	IF ( to_vtk ) THEN
4449	CALL check_open( 25 + av )
4450	!
4451	!-- Write time coordinate
4452	DO i = 0, io_blocks - 1
4453	IF ( i == io_group ) THEN
4454	WRITE ( 25 + av ) LEN_TRIM( 'END' )
4455	WRITE ( 25 + av ) 'END'
4456	ENDIF
4457	#if defined( __parallel )
4458	CALL MPI_BARRIER( comm2d, ierr )
4459	#endif
4460	ENDDO
4461	ENDIF
4462
4463	END SUBROUTINE surface_data_output_last_action
4464
4465
4466	!--------------------------------------------------------------------------------------------------!
4467	! Description:
4468	! ------------
4469	!> Read module-specific global restart data (Fortran binary format).
4470	!---------------------------------------------------------------------------------------------------!
4471	SUBROUTINE surface_data_output_rrd_global_ftn( found )
4472
4473
4474	USE control_parameters, &
4475	ONLY: length, &
4476	restart_string
4477
4478	IMPLICIT NONE
4479
4480	LOGICAL, INTENT(OUT) :: found !< flag indicating if variable was found
4481
4482	found = .TRUE.
4483
4484	SELECT CASE ( restart_string(1:length) )
4485
4486	CASE ( 'average_count_surf' )
4487	READ ( 13 ) average_count_surf
4488	CASE ( 'time_dosurf_av' )
4489	READ ( 13 ) time_dosurf_av
4490
4491	CASE DEFAULT
4492
4493	found = .FALSE.
4494
4495	END SELECT
4496
4497
4498	END SUBROUTINE surface_data_output_rrd_global_ftn
4499
4500
4501	!--------------------------------------------------------------------------------------------------!
4502	! Description:
4503	! ------------
4504	!> Read module-specific global restart data (MPI-IO).
4505	!--------------------------------------------------------------------------------------------------!
4506	SUBROUTINE surface_data_output_rrd_global_mpi
4507
4508	CALL rrd_mpi_io( 'average_count_surf', average_count_surf )
4509	CALL rrd_mpi_io( 'time_dosurf_av', time_dosurf_av )
4510
4511	END SUBROUTINE surface_data_output_rrd_global_mpi
4512
4513
4514	!--------------------------------------------------------------------------------------------------!
4515	! Description:
4516	! ------------
4517	!> Read module-specific local restart data arrays (Fortran binary format).
4518	!--------------------------------------------------------------------------------------------------!
4519	SUBROUTINE surface_data_output_rrd_local_ftn( found )
4520
4521
4522	USE control_parameters, &
4523	ONLY: length, &
4524	restart_string
4525
4526	IMPLICIT NONE
4527
4528	LOGICAL, INTENT(OUT) :: found
4529
4530
4531	found = .TRUE.
4532
4533	SELECT CASE ( restart_string(1:length) )
4534
4535	CASE ( 'surfaces%var_av' )
4536	READ ( 13 ) surfaces%var_av
4537
4538	CASE DEFAULT
4539
4540	found = .FALSE.
4541
4542	END SELECT
4543
4544
4545	END SUBROUTINE surface_data_output_rrd_local_ftn
4546
4547
4548	!--------------------------------------------------------------------------------------------------!
4549	! Description:
4550	! ------------
4551	!> Read module-specific local restart data arrays (MPI-IO).
4552	!--------------------------------------------------------------------------------------------------!
4553	SUBROUTINE surface_data_output_rrd_local_mpi
4554
4555	IMPLICIT NONE
4556
4557	CHARACTER(LEN=3) :: dum !< dummy string to create output-variable name
4558
4559	INTEGER(iwp) :: i !< grid index in x-direction
4560	INTEGER(iwp) :: j !< grid index in y-direction
4561	INTEGER(iwp) :: l !< running index surface orientation
4562	INTEGER(iwp) :: m !< running index surface elements
4563	INTEGER(iwp) :: n !< counting variable
4564	INTEGER(iwp) :: nv !< running index over number of variables
4565
4566	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< end index of surface data at (j,i)
4567	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: global_start_index !< index array for surface data (MPI-IO)
4568	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: num_surf !< number of surface data at (j,i)
4569	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index of surface data at (j,i)
4570
4571	LOGICAL :: array_found !< flag indicating whether variable is in restart data or not
4572	LOGICAL :: ldum !< dummy variable
4573
4574	REAL(wp), DIMENSION(:), ALLOCATABLE :: surf_in !< input array in expected restart format
4575
4576	!
4577	!-- Note, surface data which is written to file is organized in a different way than
4578	!-- the output surface data. The output surface data is a concatenated array of the
4579	!-- different surface types and orientations, while the mpi-io expects surface data that
4580	!-- is consecutive in terms of start- and end-index, i.e. organized along the (j,i)
4581	!-- grid index. Hence, data need to be tranformed back to the output surface data.
4582	ALLOCATE( end_index(nys:nyn,nxl:nxr) )
4583	ALLOCATE( num_surf(nys:nyn,nxl:nxr) )
4584	ALLOCATE( start_index(nys:nyn,nxl:nxr) )
4585	ALLOCATE( global_start_index(nys:nyn,nxl:nxr) )
4586
4587	ALLOCATE( surf_in(1:surfaces%ns) )
4588
4589	CALL rd_mpi_io_check_array( 'surfaces%start_index', found = array_found )
4590	IF ( array_found ) CALL rrd_mpi_io( 'surfaces%start_index', start_index )
4591
4592	CALL rd_mpi_io_check_array( 'surfaces%end_index', found = array_found )
4593	IF ( array_found ) CALL rrd_mpi_io( 'surfaces%end_index', end_index )
4594
4595	CALL rd_mpi_io_check_array( 'surfaces%global_start_index', found = array_found )
4596	IF ( array_found ) CALL rrd_mpi_io( 'surfaces%global_start_index', global_start_index )
4597
4598	CALL rd_mpi_io_surface_filetypes( start_index, end_index, ldum, global_start_index )
4599
4600	DO nv = 1, dosurf_no(1)
4601	WRITE( dum, '(I3.3)' ) nv
4602
4603	CALL rd_mpi_io_check_array( 'surfaces%var_av' // TRIM( dum ), found = array_found )
4604
4605	IF ( array_found ) THEN
4606
4607	CALL rrd_mpi_io_surface( 'surfaces%var_av' // TRIM(dum), surf_in )
4608	!
4609	!-- Write temporary input variable back to surface-output data array.
4610	n = 0
4611	num_surf = 0
4612	DO l = 0, 1
4613	DO m = 1, surf_def_h(l)%ns
4614	i = surf_def_h(l)%i(m)
4615	j = surf_def_h(l)%j(m)
4616	n = n + 1
4617	surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
4618	num_surf(j,i) = num_surf(j,i) + 1
4619	ENDDO
4620	DO m = 1, surf_lsm_h(l)%ns
4621	i = surf_lsm_h(l)%i(m)
4622	j = surf_lsm_h(l)%j(m)
4623	n = n + 1
4624	surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
4625	num_surf(j,i) = num_surf(j,i) + 1
4626	ENDDO
4627	DO m = 1, surf_usm_h(l)%ns
4628	i = surf_usm_h(l)%i(m)
4629	j = surf_usm_h(l)%j(m)
4630	n = n + 1
4631	surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
4632	num_surf(j,i) = num_surf(j,i) + 1
4633	ENDDO
4634	ENDDO
4635
4636	DO l = 0, 3
4637	DO m = 1, surf_def_v(l)%ns
4638	i = surf_def_v(l)%i(m)
4639	j = surf_def_v(l)%j(m)
4640	n = n + 1
4641	surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
4642	num_surf(j,i) = num_surf(j,i) + 1
4643	ENDDO
4644	DO m = 1, surf_lsm_v(l)%ns
4645	i = surf_lsm_v(l)%i(m)
4646	j = surf_lsm_v(l)%j(m)
4647	n = n + 1
4648	surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
4649	num_surf(j,i) = num_surf(j,i) + 1
4650	ENDDO
4651	DO m = 1, surf_usm_v(l)%ns
4652	i = surf_usm_v(l)%i(m)
4653	j = surf_usm_v(l)%j(m)
4654	n = n + 1
4655	surfaces%var_av(n,nv) = surf_in(start_index(j,i)+num_surf(j,i))
4656	num_surf(j,i) = num_surf(j,i) + 1
4657	ENDDO
4658	ENDDO
4659	ENDIF
4660	ENDDO
4661
4662
4663	END SUBROUTINE surface_data_output_rrd_local_mpi
4664
4665
4666	!--------------------------------------------------------------------------------------------------!
4667	! Description:
4668	! ------------
4669	!> This routine writes the respective restart data.
4670	!--------------------------------------------------------------------------------------------------!
4671	SUBROUTINE surface_data_output_wrd_global
4672
4673	IMPLICIT NONE
4674
4675	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
4676
4677	CALL wrd_write_string( 'average_count_surf' )
4678	WRITE ( 14 ) average_count_surf
4679
4680	CALL wrd_write_string( 'time_dosurf_av' )
4681	WRITE ( 14 ) time_dosurf_av
4682
4683	ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN
4684
4685	CALL wrd_mpi_io( 'average_count_surf', average_count_surf )
4686	CALL wrd_mpi_io( 'time_dosurf_av', time_dosurf_av )
4687
4688	ENDIF
4689
4690	END SUBROUTINE surface_data_output_wrd_global
4691
4692
4693	!--------------------------------------------------------------------------------------------------!
4694	! Description:
4695	! ------------
4696	!> This routine writes restart data which individual on each PE
4697	!--------------------------------------------------------------------------------------------------!
4698	SUBROUTINE surface_data_output_wrd_local
4699
4700	IMPLICIT NONE
4701
4702	CHARACTER(LEN=3) :: dum !< dummy string to create output-variable name
4703
4704	INTEGER(iwp) :: i !< grid index in x-direction
4705	INTEGER(iwp) :: j !< grid index in y-direction
4706	INTEGER(iwp) :: l !< running index surface orientation
4707	INTEGER(iwp) :: m !< running index surface elements
4708	INTEGER(iwp) :: n !< counting variable
4709	INTEGER(iwp) :: nv !< running index over number of variables
4710	INTEGER(iwp) :: start_index_aggregated !< sum of start-index at (j,i) over all surface types
4711
4712	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: end_index !< end index of surface data at (j,i)
4713	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: global_start_index !< index array for surface data (MPI-IO)
4714	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: num_surf !< number of surface data at (j,i)
4715	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index of surface data at (j,i)
4716
4717	LOGICAL :: surface_data_to_write !< switch for MPI-I/O if PE has surface data to write
4718
4719	REAL(wp), DIMENSION(:), ALLOCATABLE :: surf_out !< surface data in expected restart format
4720
4721
4722	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
4723
4724	IF ( ALLOCATED( surfaces%var_av ) ) THEN
4725	CALL wrd_write_string( 'surfaces%var_av' )
4726	WRITE ( 14 ) surfaces%var_av
4727	ENDIF
4728
4729	ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN
4730	!
4731	!-- Note, surface data which is written to file is organized in a different way than
4732	!-- the output surface data. The output surface data is a concatenated array of the
4733	!-- different surface types and orientations, while the mpi-io expects surface data that
4734	!-- is consecutive in terms of start- and end-index, i.e. organized along the (j,i)
4735	!-- grid index. Hence, data need to be tranformed before it can be written to file.
4736	IF ( ALLOCATED( surfaces%var_av ) ) THEN
4737	ALLOCATE( end_index(nys:nyn,nxl:nxr) )
4738	ALLOCATE( num_surf(nys:nyn,nxl:nxr) )
4739	ALLOCATE( start_index(nys:nyn,nxl:nxr) )
4740	ALLOCATE( global_start_index(nys:nyn,nxl:nxr) )
4741	ALLOCATE( surf_out(1:surfaces%ns) )
4742	!
4743	!-- Determine the start and end index at each (j,i)-pair and resort the surface data
4744	start_index = 1
4745	end_index = 0
4746	start_index_aggregated = 1
4747	num_surf = 0
4748	DO l = 0, 1
4749	DO m = 1, surf_def_h(l)%ns
4750	i = surf_def_h(l)%i(m)
4751	j = surf_def_h(l)%j(m)
4752	num_surf(j,i) = num_surf(j,i) + 1
4753	ENDDO
4754	DO m = 1, surf_lsm_h(l)%ns
4755	i = surf_lsm_h(l)%i(m)
4756	j = surf_lsm_h(l)%j(m)
4757	num_surf(j,i) = num_surf(j,i) + 1
4758	ENDDO
4759	DO m = 1, surf_usm_h(l)%ns
4760	i = surf_usm_h(l)%i(m)
4761	j = surf_usm_h(l)%j(m)
4762	num_surf(j,i) = num_surf(j,i) + 1
4763	ENDDO
4764	ENDDO
4765
4766	DO l = 0, 3
4767	DO m = 1, surf_def_v(l)%ns
4768	i = surf_def_v(l)%i(m)
4769	j = surf_def_v(l)%j(m)
4770	num_surf(j,i) = num_surf(j,i) + 1
4771	ENDDO
4772	DO m = 1, surf_lsm_v(l)%ns
4773	i = surf_lsm_v(l)%i(m)
4774	j = surf_lsm_v(l)%j(m)
4775	num_surf(j,i) = num_surf(j,i) + 1
4776	ENDDO
4777	DO m = 1, surf_usm_v(l)%ns
4778	i = surf_usm_v(l)%i(m)
4779	j = surf_usm_v(l)%j(m)
4780	num_surf(j,i) = num_surf(j,i) + 1
4781	ENDDO
4782	ENDDO
4783
4784	start_index = 0
4785	end_index = 0
4786	start_index_aggregated = 1
4787	DO i = nxl, nxr
4788	DO j = nys, nyn
4789	start_index(j,i) = start_index_aggregated
4790	end_index(j,i) = start_index(j,i) + num_surf(j,i) - 1
4791	start_index_aggregated = start_index_aggregated + num_surf(j,i)
4792	ENDDO
4793	ENDDO
4794
4795	CALL rd_mpi_io_surface_filetypes( start_index, end_index, surface_data_to_write, &
4796	global_start_index )
4797	CALL wrd_mpi_io( 'surfaces%start_index', start_index )
4798	CALL wrd_mpi_io( 'surfaces%end_index', end_index )
4799	CALL wrd_mpi_io( 'surfaces%global_start_index', global_start_index )
4800
4801	DO nv = 1, dosurf_no(1)
4802	n = 0
4803	num_surf = 0
4804	DO l = 0, 1
4805	DO m = 1, surf_def_h(l)%ns
4806	i = surf_def_h(l)%i(m)
4807	j = surf_def_h(l)%j(m)
4808	n = n + 1
4809	surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
4810	num_surf(j,i) = num_surf(j,i) + 1
4811	ENDDO
4812	DO m = 1, surf_lsm_h(l)%ns
4813	i = surf_lsm_h(l)%i(m)
4814	j = surf_lsm_h(l)%j(m)
4815	n = n + 1
4816	surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
4817	num_surf(j,i) = num_surf(j,i) + 1
4818	ENDDO
4819	DO m = 1, surf_usm_h(l)%ns
4820	i = surf_usm_h(l)%i(m)
4821	j = surf_usm_h(l)%j(m)
4822	n = n + 1
4823	surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
4824	num_surf(j,i) = num_surf(j,i) + 1
4825	ENDDO
4826	ENDDO
4827
4828	DO l = 0, 3
4829	DO m = 1, surf_def_v(l)%ns
4830	i = surf_def_v(l)%i(m)
4831	j = surf_def_v(l)%j(m)
4832	n = n + 1
4833	surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
4834	num_surf(j,i) = num_surf(j,i) + 1
4835	ENDDO
4836	DO m = 1, surf_lsm_v(l)%ns
4837	i = surf_lsm_v(l)%i(m)
4838	j = surf_lsm_v(l)%j(m)
4839	n = n + 1
4840	surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
4841	num_surf(j,i) = num_surf(j,i) + 1
4842	ENDDO
4843	DO m = 1, surf_usm_v(l)%ns
4844	i = surf_usm_v(l)%i(m)
4845	j = surf_usm_v(l)%j(m)
4846	n = n + 1
4847	surf_out(start_index(j,i)+num_surf(j,i)) = surfaces%var_av(n,nv)
4848	num_surf(j,i) = num_surf(j,i) + 1
4849	ENDDO
4850	ENDDO
4851
4852	WRITE( dum, '(I3.3)' ) nv
4853
4854	CALL wrd_mpi_io_surface( 'surfaces%var_av' // TRIM( dum ), surf_out )
4855	ENDDO
4856
4857	ENDIF
4858
4859	ENDIF
4860
4861	END SUBROUTINE surface_data_output_wrd_local
4862
4863
4864	END MODULE surface_data_output_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |