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

Last change on this file since 4444 was 4444, checked in by raasch, 5 years ago

bugfix: cpp-directives for serial mode added

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