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

Last change on this file since 4313 was 4286, checked in by resler, 21 months ago

Merge branch resler into trunk

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