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

Last change on this file since 3744 was 3744, checked in by suehring, 3 years ago

Coupling of indoor model to atmosphere; output of indoor temperatures and waste heat; enable restarts with indoor model; bugfix plant transpiration; bugfix - missing calculation of 10cm temperature

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