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

Last change on this file since 3727 was 3727, checked in by gronemeier, 3 years ago

New: surface output available in NetCDF format (Makefile, netcdf_interface_mod, surface_data_output_mod)

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