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

Last change on this file since 4029 was 4029, checked in by raasch, 2 years ago

bugfix: decycling of chemistry species after nesting data transfer, exchange of ghost points and boundary conditions separated for chemical species and SALSA module, nest_chemistry option removed, netcdf variable NF90_NOFILL is used as argument instead of 1 in calls to NF90_DEF_VAR_FILL

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