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

Last change on this file since 4517 was 4517, checked in by raasch, 14 months ago

added restart with MPI-IO for reading local arrays

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