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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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