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

Last change on this file since 4313 was 4286, checked in by resler, 20 months ago

Merge branch resler into trunk

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