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

Last change on this file since 4129 was 4129, checked in by gronemeier, 2 years ago

changes in surface_data_output_mod:

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