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

Last change on this file since 3881 was 3881, checked in by suehring, 2 years ago

Make level 3 initialization of urban-surfaces consistent to input data standard; revise flagging of ground-floor level surfaces at buidlings; bugfixes in level 3 initialization of albedo; bugfix in OpenMP directive; check for zero output timestep in surface output; assure that Obukhov lenght does not become zero

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