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

Last change on this file since 3735 was 3735, checked in by dom_dwd_user, 3 years ago

biometeorology_mod.f90:
(N) Fixed auto-setting of thermal index calculation flags by output as
originally proposed by resler.
(C) removed bio_pet and outher configuration variables.
(C) Updated namelist.
(B) Forcing initialization of tmrt_av_grid to avoid mysterious mrt
values at i==0, j==0

module_interface_mod.f90:
(C) Receiving parameter j (averaging 0==.F./1==.T.) in
module_interface_check_data_output from check_parameters.f90.
(C) Passing j to bio_check_parameters.

check_parameters.f90:
(C) Passing j to module_interface_check_data_output

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