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

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