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

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

plant canopy: Missing working precision + bugfix in calculation of wind speed; surface data output: Correct x,y-coordinates of vertical surfaces in netcdf output; Change definition of azimuth angle, reference is north 0 degree

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