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

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

document last commit

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