Home

Context Navigation

surface_data_output_mod.f90 @ 3881

Last change on this file since 3881 was 3881, checked in by suehring, 5 years ago
Make level 3 initialization of urban-surfaces consistent to input data standard; revise flagging of ground-floor level surfaces at buidlings; bugfixes in level 3 initialization of albedo; bugfix in OpenMP directive; check for zero output timestep in surface output; assure that Obukhov lenght does not become zero
Property svn:keywords set to `Id`
File size: 244.5 KB

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