Home

Context Navigation

surface_data_output_mod.f90 @ 4495

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago
restart data handling with MPI-IO added, first part
Property svn:keywords set to `Id`
File size: 246.2 KB

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