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

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

remove unused variables in surface_data_output_mod and use compiler diretives to declare variables only if netcdf4-output is used

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