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

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

unused_variables removed, bugfix in im_define_netcdf_grid argument list, trim added to avoid truncation compiler warnings, save attribute added to local targets to avoid outlive pointer target warning, first argument removed from module_interface_rrd_*, file module_interface reformatted with respect to coding standards, bugfix in surface_data_output_rrd_local (variable k removed)

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