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

Last change on this file since 4598 was 4577, checked in by raasch, 10 months ago

further re-formatting to follow the PALM coding standard

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