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

Last change on this file since 3817 was 3817, checked in by suehring, 2 years ago

Correct output coordinates for vertical surface-element data

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