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

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

Consider restart data in time-averaged surface output; revise error message; split initialization of surface-output module

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