source: palm/trunk/SOURCE/surface_data_output_mod.f90

Last change on this file was 4893, checked in by raasch, 7 months ago

revised output of surface data via MPI-IO for better performance

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